Kinetic Path Summation, Multi–Sheeted Extension of Master Equation, and Evaluation of Ergodicity Coefficient

# Kinetic Path Summation, Multi–Sheeted Extension of Master Equation, and Evaluation of Ergodicity Coefficient

A. N. Gorban University of Leicester, UK
###### Abstract

We study the Master equation with time–dependent coefficients, a linear kinetic equation for the Markov chains or for the monomolecular chemical kinetics. For the solution of this equation a path summation formula is proved. This formula represents the solution as a sum of solutions for simple kinetic schemes (kinetic paths), which are available in explicit analytical form. The relaxation rate is studied and a family of estimates for the relaxation time and the ergodicity coefficient is developed. To calculate the estimates we introduce the multi–sheeted extensions of the initial kinetics. This approach allows us to exploit the internal (“micro”)structure of the extended kinetics without perturbation of the base kinetics.

###### keywords:
Path summation, Master Equation, ergodicity coefficient, transition graph, reaction network, kinetics, relaxation time, replica
###### Pacs:

Physica A 390 (2011) 1009–1025

Corresponding author: University of Leicester, LE1 7RH, UK

## 1 Introduction

### 1.1 The problem

First-order kinetics form the simplest and well-studied class of kinetic systems. It includes the continuous-time Markov chains [1, 2] (the Master Equation [3]), kinetics of monomolecular and pseudomonomolecular reactions [4], provides a natural language for description of fluxes in networks and has many other applications, from physics and chemistry to biology, engineering, sociology, and even political science.

At the same time, the first-order kinetics are very fundamental and provide the background for kinetic description of most of nonlinear systems: we almost always start from the Master Equation (it may be very high-dimensional) and then reduce the description to a lower level but with nonlinear kinetics.

For the description of the first order kinetics we select the species–concentration language of chemical kinetics, which is completely equivalent to the states–probabilities language of the Markov chains theory and is a bit more flexible in the normalization choice: the sum of concentration could be any positive number, while for the Markov chains we have to introduce special “incomplete states”.

The first-order kinetic system is weakly ergodic if it allows the only conservation law: the sum of concentration. Such a system forgets its initial condition: the distance between any two trajectories with the same value of the conservation law tends to zero when time goes to infinity. Among all possible distances, the distance () plays a special role: it decreases monotonically in time for any first order kinetic system. Further in this paper, we use the norm on the space of concentrations.

Straightforward analysis of the relaxation rate for a linear system includes computation of the spectrum of the operator of the shift in time. For an autonomous system, we have to find the “slowest” nonzero eigenvalue of the kinetic (generator) matrix. For a system with time–dependent coefficients, we have to solve the linear differential equations for the fundamental operator (the shift in time). After that, we have to analyze the spectrum of this operator. Beyond the simplest particular cases there exist no analytical formulas for such calculations.

Nevertheless, there exists the method for evaluation of the contraction rate for the first order kinetics, based on the analysis of transition graph. For this evaluation, we need to solve kinetic equations for some irreversible acyclic subsystems which we call the kinetic paths (10). These kinetic paths are combined from simple fragments of the initial kinetic systems. For such systems, it is trivial to solve the kinetic equations in quadratures even if the coefficients are time–dependent. The explicit recurrent formulas for these solutions are given (12).

We construct the explicit formula for the solution of the kinetic equation for an arbitrary system with time–dependent coefficients by the summation of solutions of an infinite number of kinetic paths (15).

On the basis of this summation formula we produce a representation of the contraction rate for weakly ergodic systems (23). Because of monotonicity, any partial sum of this formula gives an estimate for this contraction.

To calculate the estimates we introduce the multi–sheeted extensions of the initial kinetics. Such a multi–sheeted extension is a larger Markov chain together with a projection of its (the larger) state space on the initial state space and the following property: the projection of the extended random walk is a random walk for the initial chain (Section 4.2).

This approach allows us to exploit the internal (“micro”)structure of the extended kinetics without perturbation of the base kinetics.

It is difficult to find, who invented the kinetic path approach. We have used it in 1980s [5], but consider this idea as a scientific “folklore”.

In this paper we study the backgrounds of the kinetic path methods. This return to backgrounds is inspired, in particular, by the series of work [6, 7], where the kinetic path summation formula was introduced (independently, on another material and with different argumentation) and applied to analysis of large stochastic systems. The method was compared to the kinetic Gillespie algorithm [8] and on model systems it was demonstrated [7] that for ensembles of rare trajectories far from equilibrium, the path sampling method performs better.

For the linear chains of reversible semi-Markovian processes with nearest neighbors hopping, the path summation formula was developed with counting all possible trajectories in Laplace space [9]. Higher order propagators and the first passage time were also evaluated. This problem statement was inspired, in particular, by the evolving field of single molecules (for more detail see [10]).

The idea of kinetic path with selection of the dominant paths gives an effective generalization of the limiting step approximation in chemical kinetics [11, 12].

## 2 Basic Notions

Let us recall the basic facts about the first-order kinetics. We consider a general network of linear reactions. This network is represented as a directed graph (digraph) ([13, 14]): vertices correspond to components (, edges correspond to reactions (). For the set of vertices we use notation , and for the set of edges notation . For each vertex, , a positive real variable (concentration) is defined. Each reaction is represented by a pair of numbers , . For each reaction a nonnegative continuous bounded function, the reaction rate coefficient (the variable “rate constant”) is given. To follow the standard notation of the matrix multiplication, the order of indexes in is always inverse with respect to reaction: it is , where the arrow shows the direction of the reaction. The kinetic equations have the form

 dcidt=∑j, j≠i(kij(t)cj−kji(t)ci), (1)

or in the vector form: . The quantities are concentrations of and is a vector of concentrations. We don’t assume any special relation between constants, and consider them as independent quantities.

For each , the matrix of kinetic coefficients has the following properties:

• non-diagonal elements of are non-negative;

• diagonal elements of are non-positive;

• elements in each column of have zero sum.

This family of matrices coincides with the family of generators of finite Markov chains in continuous time ([1, 2]).

A linear conservation law is a linear function defined on the concentrations , whose value is preserved by the dynamics (1). Equation (1) always has a linear conservation law: .

Another important and simple property of this equation is the preservation of positivity for the solution of (1) : if for all then for .

For many technical reasons it is useful to discuss not only positive solutions to (1) and further we do not automatically assume that .

The time shift operator which transforms into is . This is a column-stochastic matrix:

 uij(t,t0)≥0 ,  ∑iuij(t,t0)=1  (t≥t0) .

This matrix satisfies the equation:

 dU(t,t0)dt=KU(t,t0)  or  duildt=∑j(kij(t)ujl−kji(t)uil) (2)

with initial conditions , where is the unit operator ().

Every stochastic in column operator is a contraction in the norm on the invariant hyperplanes . It is sufficient to study the restriction of on the invariant subspace :

 ∥Ux∥≤δ∥x∥if∑ixi=0

for some . The minimum of such is , the norm of the operator restricted to its invariant subspace . One of the definitions of weak ergodicity is [15]. The unit ball of the norm restricted to the subspace is a polyhedron with vertices

 gij=12(ei−ej),i≠j , (3)

where are the standard basis vectors in : , is the Kronecker delta. For a norm with the polyhedral unit ball, the norm of the operator is

 maxv∈V∥U(v)∥ ,

where is the set of vertices of the unit ball. Therefore, for a ball with vertices (3)

 δU=∥U∥=12maxi,j∑k|uki−ukj|≤1 . (4)

This is a half of the maximum of the distances between columns of . The ergodicity coefficient, , is zero for a matrix with unit norm and one if transforms any two vectors with the same sum of coordinates in one vector ().

The contraction coefficient (4) is a norm of operator and therefore has a “submultiplicative” property: for two stochastic in column operators the coefficient could be estimated through a product of the coefficients

 δUW≤δUδW . (5)

We will systematically use this property in such a way. In many estimates we find an upper border , . In such a case, exponentially with . Nevertheless, the estimate may originally have a positive limit when . In this situation we can use for bounded and for exploit the multiplicative estimate (5). The moment may be defined, for example, by maximization of the negative Lyapunov exponent:

 τ1=argmaxτ>0{−ln(δ(τ))τ} . (6)

For a system with external fluxes the kinetic equation has the form

 dcidt=∑j(kij(t)cj−kji(t)ci)+Πi(t) . (7)

The Duhamel integral gives for this system with initial condition :

 c(t)=U(t,t0)c(t0)+∫tt0U(t,τ)Π(τ) dτ ,

where is the vector of fluxes with components .

In particular, for stochastic in column operators this formula gives: an identity for the linear conservation law

 ∑ici(t)=∑ici(t0)+∫tt0∑iΠi(τ) dτ , (8)

and an inequality for the norm

 ∥c(t)∥≤∥U(t,t0)c(t0)∥+∫tt0∑i∥Π(τ)∥ dτ≤∥c(t0)∥+∫tt0∑i∥Π(τ)∥ dτ . (9)

We need the last formula for the estimation of contraction coefficients when the vector is not positive.

## 3 Kinetic Paths

Two vertices are called adjacent if they share a common edge. A directed path is a sequence of adjacent edges where each step goes in direction of an edge. A vertex is reachable from a vertex , if there exists a directed path from to .

Formally, a path in a reaction graph is any finite sequence of indexes (a multiindex) (, ) such that for all (i.e. there exists a reaction ). The number of the vertices in the path may be any natural number (including 1), and any vertex can be included in the path several times. If then we call the one-vertex path degenerated. There is a natural order on the set of paths: if is continuation of , i.e. and . In this order, the antecedent element (or the parent) for each is , the path which we produce from by deletion of the last step. With this definition of parents , the set of the paths with a given start point is a rooted tree.

###### Definition 1

For each path we define an auxiliary set of reaction, the kinetic path :

 BI1(i1)ki2i1−−−−→B22(i2)ki3i2−−−−→…kiqiq−1−−−−→BIq(iq)⏐⏐↓κi1¯¯¯¯i2⏐⏐↓κi2¯¯¯¯i3⏐⏐↓κiq (10)

The vertices of the kinetic path (10) are auxiliary components. Each of them is determined by the path multiindex and the position in the path . There is a correspondence between the auxiliary component and the component of the original network. The coefficient is a sum of the reaction rate coefficients for all outgoing reactions from the vertex of the original network, and the coefficient is this sum without the term which corresponds to the reaction :

 κi=∑l, l≠ikli,κi¯j=∑l, l≠i,jkli .

A quantity, the concentration , corresponds to any vertex of the kinetic path and a kinetic equation of the standard form can be written for this path. The end vertex, , plays a special role in the further consideration and we use the special notations: , , , is the reaction rate coefficient of the last outgoing reactions in (10) (the last vertical arrow) and is the reaction rate coefficient of the last incoming reaction in (10) (the last horizontal arrow).

We use for the incoming flux for the terminal vertex of the kinetic path (10) and for the outgoing flux for this vertex.

Let us consider the set of all paths with the same start point and the solutions of all the correspondent kinetic equations with initial conditions:

 bI1(i1)=1,bIl(il)=0forl>1 .

For the concentrations of the terminal vertices this self-consistent set of initial conditions gives the infinite chain (or, to be more precise, the tree) of simple kinetic equations for the set of variables , :

 ˙ς1=−κ1(t)ς1,˙ςI=−κI(t)ςI+kI(t)ςI− , (11)

where index 1 corresponds to the degenerated path which consists of one vertex (the start point only) and corresponds to .

This simple chain of equations with initial conditions, and for , has a recurrent representation of solution:

 ς1(t)=exp(−∫tt0κ1(τ)dτ),ςI(t)=∫tt0exp(−∫tθκI(τ)dτ)kI(θ)ςI−(θ)dθ . (12)

The analogues of the Kirchhoff rules from the theory of electric or hydraulic circuits are useful for outgoing flux of a path and for incoming fluxes of the paths which are the one-step continuations of this path (i.e. ):

 κJςJ=∑I, I−=JkIςI− . (13)

Let us rewrite this formula as a relation between the outgoing flux from the last vertex of and incoming fluxes for the last vertices of paths ():

 P−J=∑I, I−=JP+I . (14)

The Kirchhoff rule (14) together with the kinetic equation for given initial conditions immediately implies the following summation formula.

###### Theorem 1

Let us consider the solution to the initial kinetic equations (1) with the initial conditions . Then

 cj(t)=∑I∈I1, iI=jςI(t) (15)

Proof. To prove this formula let us prove that the sum from the right hand side (i) exists (ii) satisfies the initial kinetic equations (1) and (iii) satisfies the selected initial conditions.

Convergence of the series with positive terms follows from the boundedness of the set of the partial sums, which follows from the Kirchhoff rules. According to them,

 ∑I∈I1ςI(t)≡1

because consists of the paths with the selected initial point only.

The sum

 Cj=∑I∈I1, iI=jςI

satisfies the kinetic equation (1). Indeed, let be the set of all paths from to . Let us find the set of all paths of the form . This set (we call it ) consists of all paths to all points which are connected to by a reaction:

 I−1j=⋃(l,j)∈EI1l .

From this identity and the chain of the kinetic equations (11) we get immediately that

 dCidt=∑j, j≠i(kij(t)Cj−kji(t)Ci), (16)

The coincidence of the initial conditions for and is obvious. Hence, because of the uniqueness theorem for equations (1) we proved that .

It is convenient to reformulate Theorem 1 in the terms of the fundamental operator . The th column of is a solution of (1) with initial conditions . Therefore, we have proved the following theorem. Let be the set of all paths with the initial vertex and the end vertex and be the solutions of the chain (11) for with initial conditions: and for .

###### Theorem 2
 uji(t,t0)=∑I∈IijςI(t) .     □ (17)

Remark 1. It is important that all the terms in the sum (17) are non-negative, and any partial sum gives the approximation to from below.

Remark 2. If the kinetic coefficients are constant then the Laplace transform gives a very simple representation for solution to the chain (11) (see also computations in [9, 6]). The kinetic path (10) is a sequence of elementary links

 …kirir−1−−−−→Brr(ir)kir+1ir−−−−→…⏐⏐↓κir¯¯¯¯¯¯¯¯ir+1 (18)

The transfer function for the link (18) is the ratio of the output Laplace Transform to the input Laplace Transform for the equation. Let the input be a function and the output be , where is the solution to equation

 ˙bi1=−κi1bir+Xi1(t);˙bir=−κirbir+kirir−1Xir(t)(r>1)

with zero initial conditions. The Laplace transform gives

 Wi1=1p+κi1,Wir=kirir−1p+κir(r>1)

for a link (18) and for the whole path (10) we get

 WI=1p+κi1q∏r=2kirir−1p+κir. (19)

(compare, for example, to formula (9) in [6]). It is worth to mention commutativity of this product: it does not change after a permutation of internal links. For the infinite chain (11) with initial conditions, and for , the Laplace transformation of solutions is

 LςI=WI (20)

## 4 Evaluation of Ergodicity Coefficient

### 4.1 Preliminaries: Weak Ergodicity and Annihilation Formula

#### 4.1.1 Geometric Criterion of Weak Ergodicity

In this Subsection, let us consider a reaction kinetic system (1) with constant coefficients for .

A set is positively invariant with respect to the kinetic equations (1), if any solution that starts in at time () belongs to for ( if ). It is straightforward to check that the standard simplex is a positively invariant set for kinetic equation (1): just check that if for some , and all then . This simple fact immediately implies the following properties of :

• All eigenvalues of have non-positive real parts, , because solutions cannot leave in positive time;

• If then , because the intersection of with any plane is a polygon, and a polygon cannot be invariant with respect to rotations to sufficiently small angles;

• The Jordan cell of that corresponds to the zero eigenvalue is diagonal – because all solutions should be bounded in for positive time.

• The shift in time operator is a contraction in the norm for : there exists such a monotonically decreasing (non-increasing) function (, , that for any two solutions of (1)

 ∑i|ci(t)−c′i(t)|≤δ(t)∑i|ci(0)−c′i(0)|. (21)

Moreover, if for the values of all linear conservation laws coincide then monotonically when .

The first-order kinetic system is weakly ergodic if it allows only the conservation law: the sum of concentration. Such a system forgets its initial condition: distance between any two trajectories with the same value of the conservation law tends to zero when time goes to infinity.

The difference between weakly ergodic and ergodic systems is in obligatory existence of a strictly positive stationary distribution: for an ergodic system, in addition, a strictly positive steady state exists: and all . Examples of weakly ergodic but not ergodic systems: a chain of reactions and symmetric random walk on an infinite lattice.

The weak ergodicity of the network follows from its topological properties.

###### Theorem 3

The following properties are equivalent (and each one of them can be used as an alternative definition of weak ergodicity):

1. There exists a unique independent linear conservation law for kinetic equations (this is ).

2. For any normalized initial state () there exists a limit state

 c∗=limt→∞exp(Kt)c(0)

that is the same for all normalized initial conditions: For all ,

 limt→∞exp(Kt)c=b0(c)c∗.
3. For each two vertices we can find such a vertex that is reachable both from and from . This means that the following structure exists:

 Ai→…→Ak←…←Aj . (22)

One of the paths can be degenerated: it may be or .

4. For operator is a strong contraction in the invariant subspace in the norm: , the function is strictly monotonic and when     .

The proof of this theorem could be extracted from detailed books about Markov chains and networks ([1, 17]). In its present form it was published in [5] with explicit estimations of the ergodicity coefficients.

Let us demonstrate how to prove the geometric criterion of weak ergodicity, the equivalence .

Let us assume that there are several linearly independent conservation laws, linear functionals , . The linear transform maps the standard simplex in (, ) onto a polyhedron . Because of linear independence of the system , , this has nonempty interior. Hence, it has no less than vertices , .

The preimage of every point in is a positively invariant subset with respect to kinetic equations because the standard simplex is positively invariant and the functionals are the conservation laws. In particular, preimage of every vertex is a positively invariant face of , ; if .

Each vertex of the standard simplex corresponds to a component : at this vertex and other there. Let the vertices from correspond to the components which form a set ; if .

For any and every reaction the component also belongs to because is positively invariant and a solution to kinetic equations cannot leave this face. Therefore, if , and then there is no such vertex that is reachable both from and from . We proved the implication .

Now, let us assume that the statement 3 is wrong and there exist two such components and that no components are reachable both from and . Let and be the sets of components reachable from and (including themselves), respectively; .

For every concentration vector a limit exists (because all eigenvalues of have non-positive real part and the Jordan cell of that corresponds to the zero eigenvalue is diagonal). The operator is linear operator in . Let us define two linear conservation laws:

 bi(c)=∑Ar∈Sic∗r(c),  bj(c)=∑Ar∈Sjc∗r(c) .

These functionals are linearly independent because for a vector with coordinates we get , and for a vector with coordinates we get , . Hence, the system has at least two linearly independent linear conservation laws. Therefore, .

#### 4.1.2 Annihilation Formula

In this Section, we find an exact expression for the contraction coefficients for the time evolution operator in norm on the invariant subspace . The unit -ball in this subspace is a polyhedron with vertices , where are the standard basic vectors in (3). The contraction coefficient of an operator is its norm on that subspace (4), this is half of the maximum of the distances between columns of .

The kinetic path summation formula (17) estimates the matrix elements of from below, but this does not give the possibility to evaluate the difference between these elements. To use the summation formula efficiently, we need another expression for the contraction coefficient.

The th column of is a solution of the kinetic equations (1) with initial conditions . For each let us introduce the incoming flux for the vertex in this solution:

 Πij(t)=∑qkjq(t)cq(t)

(the upper index indicates the number of column in , the lower index corresponds to the number of vertex ).

Formula (4) for the contraction coefficient gives

 δ(t,t0)=12maxi,j∥U(t,t0)(ei−ej)∥ .

is a solution to the kinetic equation (1) with initial conditions , and for . This is the difference between two solutions, and . Let us use the notation

 Gij(t)=12U(t,t0)(ei−ej) .

For each we define

 Π+q=∑l,c+l>c−lkql(c+l−c−l),Π−q=∑l,c+l

The decrease in the norm of can be represented as an annihilation of a flux with an equal amount of concentration from the vertex by the following rules:

1. If then the flux annihilates with an equal amount of positive concentration stored at vertex (Fig. (a)a);

2. If then the flux annihilates with an equal amount of negative concentration stored at vertex (Fig. (b)b);

3. If then the flux annihilates with the equal amount from the opposite flux (Fig. (c)c).

Let us summarize these rules in one formula:

###### Proposition 1
 ddt∥Gij(t)∥l1=−∑q, c+q>c−qΠ−q(t)−∑q, c+q

Immediately from (23) we obtain the following integral formula

 1−∥Gij(t)∥l1=∫tt0⎛⎜⎝∑q, c+q>c−qΠ−q(τ)+∑q, c+q

The annihilation formula gives us a better understanding of the nature of contraction but is not fully constructive because it uses fluxes from solutions to the initial kinetic equation (1).

### 4.2 Multi–Sheeted Extensions of Kinetic System

Let us introduce a multi–sheeted extension of a kinetic system.

###### Definition 2

The vertices of a multi–sheeted extension of the system (1) are where is a finite or countable set. An individual vertex is (). The corresponding concentration is . The reaction rate constant for is . This system is a multi–sheeted extension of the initial system if an identity holds:

 ∑rk(j,r)(i,l)=kji  for all  l . (25)

This means that the flux from each vertex is distributed between sheets, but the sum through sheets is the same as for the initial system. We call the kinetic behavior of the sum the base kinetics.

A simple proposition is important for further consideration.

###### Proposition 2

If is a solution to the extended multi–sheeted system then

 ci(t)=∑lc(i,l)(t) (26)

is a solution to the initial system and

 ∑il|c(i,l)(t)|≥∑i|ci(t)| . (27)

(Here we do not assume positivity of all ).

Formula (25) allows us to redirect reactions from one sheet to another (Fig. 2) without any change of the base kinetics. In the next section we show how to use this possibility for effective calculations.

Formula (26) means that kinetics of the extended system in projection on the initial space is the base kinetics: the components are projected in the projected vector of concentrations is and the projected kinetics is given by the initial Master equation with the projected coefficients . “Recharging” is any change of the non-negative extended coefficients which does not change the projected coefficients.

The key role in the further estimates plays formula (27). We will apply this formula to the solutions with the zero sums of coordinates, they are differences between the normalized positive solutions.

### 4.3 Fluxes and Mixers

In this Subsection, we present the system of estimates for the contraction coefficient. The main idea is based on the following property which can be used as an alternative definition of weak ergodicity (Theorem 3): For each two vertices we can find a vertex that is reachable both from and from . This means that the following structure exists:

 Ai→…→Aq←…←Aj.

One of the paths can be degenerated: it may be or . The positive flux from meets the negative flux from at point and one of them annihilates with the equal amount of the concentration of opposite sign.

Let us generalize this construction. Let us fix three different vertices: (the “positive source”), (the “negative source”) and (the “mixing point”). The degenerated case or we discuss separately. Let be such a system of vertices that , and there exists an oriented path in from to . Analogously, let be such a system of vertices that , and there exists an oriented path in from to . We assume that .

With each subset of vertices we associate a kinetic system (subsystem): for

 ˙cr=∑l, Al∈S, r≠lkrlcl−n∑p=1kprcr . (28)

In this subsystem, we retain all the outgoing reaction for and delete the reactions which lead to vertices in from “abroad”.

The flux from to is

 Π+S=∑r, Ar∈S+kqrcr(t) ,

where is a component of the solution of (28) for with initial conditions . Analogously, we define the flux

 Π−S=∑r, Ar∈S−kqrcr(t) ,

where is a component of the solution of (28) for with initial conditions . Decrease of the norm is estimated by the following theorem.

The system we call a mixer, that is a device for mixing. An elementary mixer consists of two kinetic paths (22) with the corespondent outgoing reactions:

 Ai1ki2i1−−−−→…kirir−1−−−−→Airkirir+1←−−−−…kir+l−1ir+l←−−−−−−Air+l⏐⏐↓κi1¯¯¯¯i2⏐⏐↓κirκir+l¯¯¯¯¯¯¯¯¯¯¯¯ir+l−1⏐⏐↓ (29)

where , , .

The degenerated elementary mixer consists of one kinetic path:

 Ai1<