arXiv:1609.02063v1 [] 6 Sep 2016

Zuse Institute Berlin

I SABEL B ECKENBACH L EON E IFLER
KONSTANTIN FACKELDEY A MBROS G LEIXNER
A NDREAS G REVER M ARCUS W EBER JAKOB W ITZIG

Mixed-Integer Programming for Cycle
Detection in Non-reversible Markov
Processes

A version of this paper is submitted to Multiscale Modeling and Simulation: A SIAM Interdisciplinary Journal. Date of submission: 25.08.2016.

ZIB Report 16-39 (August 2016)

Takustr. 7
14195 Berlin
Germany

Zuse Institute Berlin
Takustr. 7
D-14195 Berlin
Telefon: +49 30-84185-0
Telefax: +49 30-84185-125
e-mail: bibliothek@zib.de
URL: http://www.zib.de
ZIB-Report (Print) ISSN 1438-0064
ZIB-Report (Internet) ISSN 2192-7782

Mixed-Integer Programming for Cycle Detection
in Non-reversible Markov Processes∗
Isabel Beckenbach1, Leon Eifler1, Konstantin Fackeldey2, Ambros Gleixner1,
Andreas Grever2, Marcus Weber2, and Jakob Witzig1
1

Zuse Institute Berlin, Department Optimization, Takustr. 7, 14195 Berlin,
Germany, {beckenbach,eifler,gleixner,witzig}@zib.de
2

Zuse Institute Berlin, Department Numerical Mathematics, Takustr. 7,
14195 Berlin, Germany, {fackeldey,grever,weber}@zib.de

September 8, 2016

Abstract
In this paper, we present a new, optimization-based method to exhibit
cyclic behavior in non-reversible stochastic processes. While our method
is general, it is strongly motivated by discrete simulations of ordinary
differential equations representing non-reversible biological processes, in
particular molecular simulations. Here, the discrete time steps of the
simulation are often very small compared to the time scale of interest,
i.e., of the whole process. In this setting, the detection of a global cyclic
behavior of the process becomes difficult because transitions between individual states may appear almost reversible on the small time scale of
the simulation.
We address this difficulty using a mixed-integer programming model
that allows us to compute a cycle of clusters with maximum net flow, i.e.,
large forward and small backward probability. For a synthetic genetic
regulatory network consisting of a ring-oscillator with three genes, we
show that this approach can detect the most productive overall cycle,
outperforming classical spectral analysis methods. Our method applies to
general non-equilibrium steady state systems such as catalytic reactions,
for which the objective value computes the effectiveness of the catalyst.

1

Introduction

Simulation data stemming from chemical or biological processes typically lead
to a huge amount of data points (time series) in some high dimensional space.
Often a direct interpretation of these data for prediction or understanding of
the underlying physical process is almost impossible due to the range of spatial
and temporal scales. This led to the development of coarse graining methods
that provide relevant information of the system on a level with less complexity.
∗ A version of this paper is submitted to Multiscale Modeling and Simulation: A SIAM
Interdisciplinary Journal. Date of submission: 25.08.2016.

1

One example widely used in the context of biological and chemical processes
are Markov State Models (MSM) [3, 7, 8, 22]. In a MSM the underlying long
time series is described by a Markov chain on some low dimensional space, i.e.,
there exists a stochastic transition matrix P , whose entries pij can be interpreted
as the portion of the system that will transit from state i to state j in one time
step. If the vector v(t) represents a distribution at time step t, the matrix vector
multiplication v(t)T P = v(t + 1)T is a propagation of that distribution for one
time step. The stationary distribution meets the condition π T = π T P , which
means that π is a steady state.
A Markov chain with n states is called reversible if and only if the detailed
balance condition
πi pij = πj pji
(1)
holds for all i, j = 1, . . . , n. MSMs are well understood if they are applied to
reversible processes, e.g., simulation of a molecule in water, for which spectral
clustering is a commonly used coarse graining method.
In practice, however, many chemical and biological processes are not reversible. If the detailed balance condition is not met for states i and j, we can
define a net flow between these two states by taking the difference πi pij − πj pji .
It follows immediately that for n > 2 we can find cycles of positive flow as long
as the process is non-reversible.
Non-reversible processes with a stationary distribution are sometimes called
non-equilibrium steady state (NESS) processes (see [10]). A catalytic process
is one example of a NESS. A catalytic process is a chemical process in which
the rate of chemical reactions is increased due to the presence of a catalyst. If
one considers the (ensemble distribution of) states of the catalyst, they usually
undergo a cycle in the conformational space, returning to its initial state at
the end of the process. Thus, the catalyst is in a NESS. Each cycle transforms
educts into products as illustrated in Figure 1. The faster the productive cycle,
the more effective the catalyst. For processes with such a behavior we want to
find a clustering that maximizes the net flow of the cycle.
Spectral clustering is applied in two cases. If the process consists of m
metastabilities, i.e., there is a set of states between which jumps only rarely
occur, it is assumed that the spectrum of P has m leading eigenvalues (see [7]).
If the process consists of a dominant cycle of m states, i.e., there is an ordered
set of states and the process jumps to the next state with a high probability, it
is assumed that the spectrum of P has m complex eigenvalues close to the unit
circle (see [10]). Though processes involving different time scales like simulation
of molecular dynamics are not likely to produce a dominant cycle. To analyze
general cyclic behavior a more flexible approach is needed.
Contribution. Our main contribution is two-fold. First, we develop a new
technique to detect cycles in Markov processes called cycle clustering. It is
general in the sense that the inherent non-reversibility may not be dominating,
i.e., the average “speed” of the cycle may be small relative to the overall time
scale of interest. In contrast, previous methods such as [10] were designed for
dominant cycles only. The computational results show that we are able to find
cycles that are far from being dominant, which makes our method applicable
to biological processes like catalysis. As such, our model could, for instance, be
used to compute the effectiveness of the catalyst.

2

Figure 1: In a catalytic process the catalyst always returns to its initial state

Second, as a more abstract methodological contribution, we model and solve
the underlying clustering problem using mixed-integer optimization. Compared
to some classical methods, this approach requires few assumptions on the input data: transition matrices do not need to be reversible and there are no
assumptions on its spectrum. The separation of clustering model and the underlying solution algorithm helps to focus on the semantics of the clustering
independently of how they can be computed. The optimality guarantees given
by mixed-integer programming solvers enable us to compute “best” clusterings
and cycles with respect to the specified objective function.
The article is organized as follows. In Section 2, we give brief introductions
to spectral clustering and the current state-of-the-art methods as well as an
overview of mixed-integer programming techniques. In Section 3, we develop
our new, optimization-based method to exhibit cyclic behavior in non-reversible
Markov processes. Moreover, we describe the underlying mixed-integer programming model and prove its complexity. In Section 4, we investigate the
viability of our approach by computational experiments both on small-scale, artificially created instances and the repressilator system from [11]. In Section 5,
we give concluding remarks.

2

Background

The clustering method presented in this paper is based on mixed-integer optimization. As such it distinguishes itself from the spectral approach to cluster
analysis most prevalent in molecular dynamics. In the following, we will briefly
describe the main ideas of spectral clustering and mixed-integer programming.

3

2.1

Spectral Clustering

Spectral clustering is a common method to identify metastable sets in Markov
State Models. Given an undirected graph G := (V, E) defined by a set of n
vertices V , a set of edges E ⊆ V2 , and a weighted adjacency matrix W ∈
[0, 1]n×n , where wij = wji > 0 if and only if there exists P
an edge eij ∈ E.
Moreover, let D be a diagonal degree matrix, where dii = j wij > 0 for all
i ∈ V . The transition matrix P that describes a random walk on G is given
w
by P := D−1 W , i.e., the probability to jump from node i to j is given by diiij .
If the graph is connected and non-bipartite, there exists a unique stationary
distribution vector π ∈ [0, 1]n , such that π T P = π T (cf. [26]).
For each set of nodes C ⊆ V , the cut between C and its complement V \ C
is defined by the set of edges eij ∈ E with i ∈ C andPj ∈ V \ C. The weight of
the cut between C and V \ C is defined by ω(C) = i∈C,j∈V \C wij . Now, the
aim is to find a partition C1 , . . . , Ck of V such that the edges between clusters
have small weight. A min-cut problem can be formulated as minimizing
k

cut(C1 , . . . , Ck ) :=

1X
ω(Ci , V \ Ci ).
2 i=1

(2)

In practice, similar objective functions like RatioCut [16] or Ncut [25] are used.
More general information on spectral clustering can be found in [26].
A relaxation of the optimization problem leads to an eigenvalue problem
for the first k eigenvectors of a graph Laplacian L = I − D−1 W = I − P .
Eigenvectors corresponding to eigenvalues that are close to the Perron root
λ1 = 1 are the basis vectors of a special invariant subspace of P . There exists
a transformed basis of this invariant subspace with the following property. The
corresponding transformed basis vectors can be interpreted as the membership
vectors of the metastable subsets of the state space. The method which finds
this linear basis transformation is called PCCA+ [8].
The analysis using a spectral approach is well understood when it comes
to reversible Markov chains, i.e., when the evolution of the process is invariant
under time reversal. If the process is non-reversible then the weight matrix W is
not symmetric and the eigenvectors are not orthogonal. The canonical approach
cannot be applied anymore. In [28] a variation on the PCCA+ algorithm has
been developed that can be applied to non-reversible matrices, called G-PCCA.
Instead of the spectral decomposition, G-PCCA uses the Schur decomposition
such that the former follows as a special case.
Along a different line of research, some recent articles have tried to identify
dominant structures such as cycles in non-reversible transition matrices, e.g., [9,
10]. In these methods, the cycles are assumed to be “dominant”, which means
that there is a high probability inside the Markov chain to follow these cycles.
In this setting, complex pairs of eigenvalues of the transition matrices are clearly
identifiable.
In constrast, the method presented in this paper aims at finding cycles in
non-reversible transition matrices that are not necessarily dominant, but exhibit
only rare circular jumps between metastable sets. This type of transitions is typical for catalytic and other biological processes. On the timescale of simulation
the process is nearly reversible, but it has a small (in terms of probabilities) tendency towards the direction of a cyclic behavior. In these cases, non-reversibility
4

is inherent to the process and does not stem from, e.g.,, truncation errors due
to a finite sampling of a reversible Markov chain. Hence, methods that try to
make such input amenable to analytical tools for reversible Markov chains, e.g.,
by computing the nearest reversible matrix as in [21], are not applicable because
they destroy the characteristics of the process.

2.2

Mixed-Integer Programming

A mixed-integer program (MIP) is an optimization problem that can be written
in the form
(P )

zM IP = min{cT x | Ax ≥ b, ` ≤ x ≤ u, x ∈ Zl × Rn−l },

with objective function c ∈ Rn , constraint matrix A ∈ Rm×n , constraint righthand side b ∈ Rm , and lower and upper bound vectors `, u ∈ (R ∪ {±∞})n on
the variables. When omitting the integrality conditions, we obtain the linear
program (LP)
zLP = min{cT x | Ax ≥ b, ` ≤ x ≤ u, x ∈ Rn }.
It constitutes a relaxation of the corresponding MIP and provides a lower bound
on its optimum, i.e., zLP ≤ zM IP . This fact plays an important role in the
LP-based branch-and-bound algorithm [6, 18], the most widely used general
algorithm to solve MIPs to global optimality.
LP-based branch-and-bound is a divide-and-conquer method which starts
by solving the LP relaxation of the problem to compute a lower bound and a
solution candidate x? . If x? fulfills the integrality restrictions, the problem is
solved to optimality; if not, it is split into (typically two) disjoint subproblems,
thereby removing x? from the feasible region of both LPs. Typically, an integer
variable xi with fractional solution value x?i is selected and the restrictions
xi ≥ dx?i e and xi ≤ bx?i c are added to the two subproblems, respectively. This
step is called branching. As this process is iterated, we store and update the
best solution x̃ found so far whenever one of the subproblems has an integral
LP solution.
The key observation is that a subproblem can be disregarded when its lower
bound is greater or equal than the objective value of x̃. This is called bounding.
The branch-and-bound process is typically illustrated as a tree, cf. Figure 2. The
root node represents the original problem and the two subproblems created by
the branching step correspond to two child nodes being created for the current
node.
In modern MIP solvers the general branch-and-bound scheme is extended
by various algorithms to enhance the performance, see, e.g., [1, 2, 23] and many
more. Nevertheless, mixed-integer programming is complex both in theory
(N P-hard, see, e.g., [15]) and in practice: As explained, state-of-the-art solvers
eventually rely on enumerative search over an exponentially large solution space
and may converge slowly. However, even when terminated early for hard problem instances, they typically provide good solutions and give proven guarantees
on the quality of the solutions returned.

5

P
x i1 ≤

bx?,P
i1 c

xi1 ≥ dx?,P
i1 e

P1
xi2 ≤

1
bx?,P
i2 c

P2
x i2 ≥

P3
5
xi3 ≤ bx?,P
i3 c

P5

1
dx?,P
i2 e

P4

xi4 ≤

2
bx?,P
i4 c

P7

2
xi4 ≥ dx?,P
i4 e

P8

5
xi3 ≥ dx?,P
i3 e

P6

Figure 2: Illustration of a branch-and-bound tree. Solving the LP relaxation of P ,
P1 , and P3 led to branching steps on variables i1 , i2 , and i3 , respectively. First feasible
solution x̃ obtained by solving the LP relaxation of P5 (dashed node). Subproblem
P6 has an infeasible LP relaxation after branching on i3 (dotted node). Afterwards,
subproblems P4 , P7 , and P8 can be disregarded due to bounding (gray shaped), i.e.,
the lower bound of each subproblem is not smaller than the objective value of x̃. Since
all subproblems are processed, x̃ is an optimal solution of P .

3

Cycle Clustering

In this section we present a new method for detecting global cyclic behavior
of a Markov process by partitioning its state space into ordered clusters. The
technique, which we call cycle clustering, can be applied to any discrete-time
Markov process over a discrete (or discretized) state space with a stationary
distribution. We show how an “optimal” cycle clustering can be computed using
algorithms from mixed-integer programming.

3.1

Setting

We consider a Markov process with a finite set of states B = {1, . . . , n}. We will
call states also bins in order to indicate that they might stem from a discretization of a continuous space as it is common in molecular simulations. However,
this is not an underlying assumption of our method.
Let P ∈ [0, 1]n×n be the matrix of conditional transition probabilities, where
pij is equal to the probability of moving from bin i to bin j in one time step
when initially in bin i. Then the only requirements for applying our method are
1. that P is right stochastic, i.e., all row sums are one, and
2. that P has a stationary distribution π, i.e.,Pπ ∈ [0, 1]n is a left-hand
eigenvector for eigenvalue 1, π T P = π T , and i=1,...,n πi = 1.
Using the stationary distribution, we can compute the matrix W = diag(π)P
of unconditional transition probabilities, which is the main input data for our
method. The entry qij = πi pij equals the probability or intense of transitions
from i to j in the whole ensemble of transitions.

3.2

Net Flow and Coherence

While in many biological processes such as catalysis it may be intuitively clear
how to define cyclic behavior using application-specific interpretations of the
6

state space, it is much less obvious in the abstract formalism of a Markov state
model. Because of the probabilistic aspect it is not meaningful to look for
a sequence of states that are visited in one fixed cyclical order. Moreover, a
trajectory may not even visit all states as it completes one “iteration” of the
cycle.
We address these difficulties by combining cycle detection with clustering.
Our goal is to partition the set of states into a fixed number of clusters and
order them in form of a cycle such that with high probability we will encounter
• transitions from one cluster to the next cluster in cycle direction, or
• transitions within one cluster, but
• no transitions between clusters in backward direction.
To quantify this, we introduce the following measure of non-reversibility between
two sets of states.
Definition 1 (net flow). Given two disjoint sets of states A, B ⊆ B, A ∩ B = ∅,
we call
X
X
f (A, B) :=
(πi pij − πj pji ) =
(qij − qji )
i∈A,j∈B

i∈A,j∈B

the net flow from set A to set B.
The net flow f (A, B) corresponds to the portion of particles transiting
from A to B minus the portion transiting backwards from B to A in one time
step. Its value equals the sum of the deviations from the detailed balance conditions (1). The net flow is signed and by definition, f (A, B) = −f (B, A) for
all A, B.
To ensure that, at the same time, clusters contain related groups of states,
we use the following definition.
Definition 2 (coherence). Given a set of states A ⊆ B, we call
X
X
g(A) :=
πi pij =
qij
i,j∈A

i,j∈A

the coherence of set A.
The coherence g(A) is hence equal to the unconditional probability of residing and remaining within A given the stationary distribution π. As such, it can
be interpreted as a proxy for measuring closeness in the original state space.

3.3

Clustering Model

By an m-cycle clustering we denote a partitioning of the set of states into
m pairwise disjoint clusters,
B=

m
[

Ck , Ck ∩ C` = ∅ for all k 6= `,

k=1

7

endowed with the cyclic order C1 → C2 → . . . → Cm → C1 . We are interested
in cycle clusterings with large total net flow between consecutive clusters,
m−1
X

f (Ck , Ck+1 ) + f (Cm , C1 ),

k=1

and large total coherence of the individual clusters, i.e.,
m
X

g(Ck ).

k=1

To use matrix notation, we can encode a clustering in an assignment matrix X ∈ {0, 1}n×m ,
(
1 for i ∈ Ck ,
Xik =
0 otherwise.
Then our objective can be expressed in terms of the projected matrix of unconditional transition probabilities
W := X T W X ∈ [0, 1]m×m .
Its diagonal entries wkk carry the coherences g(Ck ). The net flow values f (Ck , C` )
equal wk` − w`k , i.e., they can be read from the off diagonal entries of
∆ := W − W

T

= X T (W − W T )X.

By construction, the diagonal entries of ∆ are zero. Furthermore, from π T P =
π T it follows that the row sum vector of W equals its column sum vector, i.e.,
the row sums and column sums of ∆ are all zero.
As a first consequence, ∆ is the zero matrix if m ≤ 2, i.e., a clustering
into a 2-cycle cannot exhibit any non-reversibility. In this sense, the smallest
interesting case is a 3-cycle clustering. As a second consequence, for any 3-cycle
clustering, ∆ must have the special structure


0
ε −ε
ε ,
∆ = −ε 0
ε −ε 0
where we may assume ε ≥ 0 after reordering. Hence, in a 3-cycle clustering,
the net flow between each two clusters, ε, is identical. Maximizing the total
net flow 3ε is equivalent to maximizing pairwise non-reversibility. In our experiments later, we will also focus on this prototypical case of a cycle clustering
with m = 3.
Finally, note that the two objectives, non-reversibility in terms of net flow
and coherence, are not necessarily aligned: In general, we cannot assume that
there will be one clustering that maximizes both criteria at the same time.
Hence, we combine both and use a scaling parameter α > 0 to control the
emphasis on coherence to obtain the weighted objective
m−1
X

f (Ck , Ck+1 ) + f (Cm , C1 ) + α

k=1

m
X
k=1

8

g(Ck ).

(3)

4
0.45

0.55

1

0.1

9
0.55

0.45

8

3

5
0.45
0.1

0.1

0.45

2
0.45

0.45

7

6
0.55

Figure 3: A cyclic process with non-dominant cycle. The edges weights are the
unconditional transition probabilities πi pij . For better readability they have been
scaled by a factor of nine.

In our experiments, we will use a small value of α = 0.001 as a default.
The following gives an example where both criteria are important to detect
a meaningful 3-cycle.
Example 3.1. Consider the process shown in Figure 3, which has nine bins
grouped into three “natural” clusters. The unconditional transition probabilities
are attached to the edges, scaled by a factor of nine, the number of bins, for
better readability. The example is constructed symmetrically such that there is
a small cyclic flow between bins 1, 2, and 3. Each of these bins is connected to
two more bins in a reversible fashion.
To obtain the highest possible total net flow of 0.3/9, a 3-cycle clustering
must assign the states 1, 2, and 3 to different clusters. The assignment of
the remaining states has no influence on the net flow and could therefore be
performed arbitrarily—unless we take into account coherence. Coherence is
maximized by clustering 4 and 5 with 1, 6 and 7 with 2, and 8 and 9 with 3,
thus detecting the natural cyclic structure of the process.
Remark 3. In earlier models, we had also experimented with the natural idea
of maximizing some measure of reversibility within each cluster. However, this
incurs difficulties when the backward and forward probability between two states
is zero or very small. First, this introduces a connection between two states
that might not be at all or are only very weakly related. Second, it is inherently
instable to classify such edges as irreversible or reversible because of the small
differences. Coherence, as used here, avoids all these disadvantages.

3.4

Mixed-Integer Programming Formulation

In order to actually compute optimal cycle clusterings w.r.t. (3), we use the
following MIP model. For each bin i ∈ B = {1, . . . , n} and cluster k ∈ C =
{1, . . . , m}, we introduce a binary decision variables xik with
xik = 1 ⇐⇒ i ∈ Ck ⇐⇒ bin i is assigned to cluster k,

9

which correspond to the entries of the assignment matrix X used before. Furthermore, with slight abuse of notation, we introduce continuous variables fk ∈
R≥0 for the net flow from cluster k to φ(k), where
(
k + 1 if k < m,
φ : C 7→ C, φ(k) =
1
otherwise,
and continuous variables gk ∈ R≥0 for the coherence of cluster k ∈ C. Then
computing an optimal cycle clustering for a fixed number of m clusters can be
expressed as the following MIP model:
X
X
max
fk + α ·
gk
(4)
k∈C

s.t.

X

k∈C

xik = 1

for all i ∈ B

(5)

xik ≥ 1

for all k ∈ C

(6)

qij (xik xjφ(k) − xiφ(k) xjk )

for all k ∈ C

(7)

qij xik xjk

for all k ∈ C

(8)

for all i ∈ B, k ∈ C

(9)

k∈C

X
i∈B

fk =

X
i,j∈B

gk =

X
i,j∈B

xik ∈ {0, 1}
fk , gk ∈ R≥0

for all k ∈ C

(10)

Constraints of type (5) ensure that each bin i is assigned to exactly one
cluster k. Constraints (6) assert that there are no empty clusters. The net
flow between two consecutive clusters k and φ(k) is described by constraints
of type (7). The coherence within each cluster k is modeled by constraints of
type (8).
The products of binary variables appearing in constraints (7) and (8) are
nonlinear. We have applied a standard reformulation technique [13] to obtain a
mixed-integer linear programming formulation that can be solved by standard
state-of-the-art MIP solvers. This requires the introduction of additional auxiliary variables, but yields significantly lower solution times than using global
mixed-integer nonlinear programming solvers.
Remark 4. We want to point out that the MIP model is even more general
than our initial development of the cycle clustering approach. It only requires a
non-negative matrix W as input and does not rely on the form W = diag(π)P
with P being a stochastic matrix and π its stationary distribution vector.

3.5

Complexity of Cycle Clustering

While it is known that mixed-integer programming is N P-hard in general [24],
special subclasses still may be easier. In this section we discuss the complexity of
cycle clustering and show that it is N P-hard by a reduction from the multiway
cut problem [5].
Definition 5 (multiway cut). Let G = (V, E) be a graph with non-negative
edge weights c(e) ≥ 0 and a set of specified vertices S = {s1 , . . . , sm } ⊆ V
10

called terminals. A multiway cut is a subset of edges E 0 ⊆ E that separates the
terminals s1 , . . . , sm in the sense that there exists no path from any terminal to
any other terminal in (V, E \ E 0 ). The multiway cut problem is that of finding
a weight-minimal multiway cut.
The multiway cut problem is N P-hard for any fixed m ≥ 3 [5], which allows
us to prove the following.
Theorem 6. The cycle clustering problem as defined in Section 3.3 is N P-hard
for any m ∈ N, m ≥ 3, and any α > 0.
Proof. Suppose we are given a multiway cut instance over an undirected graph
G = (V, E) with the set of nodes V = {1, . . . , n}, the set of edges E ⊆ V2 , and
edge weights c(e) ≥ 0. Let S = {s1 , . . . , sm } be the set of terminals for m ≥ 3.
Because edges between terminals contribute a constant offset to
 the objective
value of any multiway cut, we may assume w.l.o.g. that E ∩ S2 = ∅. Because
disconnected P
non-terminal nodes can be assigned arbitrarily, we may assume
w.l.o.g. that
v∈V : {uv}∈E c({uv}) > 0 for all u ∈ V \ S.
In order to show that this instance can be solved as a cycle clustering problem, we construct a directed graph D = (V, A) with arc set
A = {(uv) ∈ V × V | {uv} ∈ E, u < v}

(forward arcs)

∪ {(uv) ∈ V × V | {uv} ∈ E, u > v}

(backward arcs)

∪ {(si sφ(i) ) ∈ S × S | {si sφ(i) } ∈
/ E}.

(auxiliary arcs)

The first two parts split edges into forward and backward arcs and the third
set, which may be empty, ensures that the cycle s1 → s2 → . . . → sm → s1 is
present in D.
To define arc weights, we partition the edges of G into
Ě = {{uv} ∈ E | |{u, v} ∩ S| ≤ 1} and Ê = {{uv} ∈ E | u, v ∈ S}
and set


for all a = (si sφ(i) ) ∈ A,
M
d(a) = 0
for all a = (si sj ) ∈ A with j 6= φ(i),


c(e)/2 for all a = (uv) ∈ A, {uv} ∈ Ě.
Here M is chosen sufficiently large such that each terminal will later be forced
into a different cluster, i.e.,
X
X
αc(e) ≥
αc(e).
M>
e∈E

e∈Ě

From this construction we derive a weighted adjacency matrix Q ∈ Rn×n via
(
d(a) if a = (uv) ∈ A,
quv =
0
otherwise.
P
0
Normalization of each row by its row sum kQu· k1 =
u0 ∈V quu , which by
assumption is non-zero, gives a stochastic matrix P ∈ [0, 1]n×n with entries
puv =

quv
.
kQu· k1
11

This transition matrix P has a unique stationary distribution π given by
kQu· k1
0
u0 ∈V kQu · k1

πu = P

for u = 1, . . . , n. The corresponding matrix W of unconditional transition
probabilities has entries
kQu· k1
quv
quv
·
=P
0 0
kQu0 · k1 kQu· k1
0
u ,v 0 ∈V qu v

wuv = πu puv = P

u0 ∈V

for u, v = 1, . . . , n, and thus

wuv − wvu



M
= −M


0

for u = si , v = sφ(i) ,
for v = sφ(i) , u = si ,
otherwise.

Hence, only edges between consecutive terminals violate the detailed balance
condition and contribute to the net flow.
Now let C1 , . . . , Cm be an optimal solution of the cycle clustering problem
w.r.t. the constructed matrix P and the stationary distribution π. W.l.o.g.,
assume s1 ∈ C1 . Due to the choice of M , in any optimal clustering, terminal sk
must be in cluster Ck .
Finally, since the assignment of non-terminal nodes does not affect the net
flow, they must be assigned such as to maximize coherence. The following calculation shows that maximizing this remaining part of the objective function is
equivalent to minimizing the weight of the edges in the corresponding multiway
cut,
α

m
X

g(Ck ) = α

m
X
X

wuv = α

k=1 u,v∈Ck

k=1

=P
=P

α

u0 ,v 0 ∈V

|

{z

qu0 v0

quv
P

k=1 u,v∈Ck
m
X
X

α

u0 ,v 0 ∈V

m
X
X

quv = P

k=1 u,v∈Ck

X

u0 ,v 0 ∈V

u0 ,v 0 ∈V

qu0 v0

X

c(e)

k=1 e∈E∩(Ck )
2

X

constant

m
X

α

c(e) −
qu0 v0
Sm
Ck
} e∈E
k=1 ( 2 )
| {z } e∈E\
|
{z

constant >0

qu0 v0


c(e)
}

multiway cut weight

To summarize, we gave a polynomial reduction of the multiway cut problem
to cycle clustering, proving that cycle clustering is N P-hard.

4

Computational Experiments

To evaluate our new clustering approach we used both synthetic instances and
a well-known system of differential equations that model the interaction of
genes [11]. In the first part of this section we describe the set of instances
we have used for the computational experiments in more detail. In the second
part we discuss the solving environment and the software we have used. Finally,
we present our computational results.
12

Algorithm 1 HMC with drift
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:

Input: start vector x0 , inverse temperature β, N , drift d, random vectors r1 , . . . , rN , uniformly distributed numbers u1 , . . . , uN ∈ [0, 1]
Output: trajectory x0 , . . . , xN −1
for i = 1 to N − 1 do
xnew ← xi−1 + ri + d
if exp(−β(Ω(xnew ) − Ω(xi−1 ))) < ui then
xi ← xnew
d ← update(d)
else
xi ← xi−1
end if
end for

4.1

Testset

Catalytic Cycle. To create data sets that resemble molecular dynamical simulations, a hybrid Monte-Carlo method (HMC) was applied to a synthetic, twodimensional energy landscape Ω as described in [4, 12]. In this variant of an
HMC, the system is propagated with a drift to one of the minima defined by
the potential Ω and an additional random value, followed by a Metropolis-like
acceptance step that assures the convergence of the distribution defined by the
function Ω. If the system enters a predefined set the drift is updated for the
next state of the cycle. This creates the dynamics of a metastable system with
rare asymmetric jumps, that one would expect from a catalytic cycle.
Algorithm 1 yields a sampling of length N , which is used to compute the
transition matrix. We extracted n vectors c1 , . . . , cn from the sampling such
that the fill distance
h := max

min kxj − ci k2

j=1,...,N i=1,...,n

was minimized. The vectors ci are the centers of regions (bins) and following [27], radial basis functions
exp(−kx − ci k22 )
Φi (x) := Pn
2
k=1 exp(−kx − ck k2 )
with values in (0, 1) were used as membership function, i.e., instead of assigning
every x to one bin, it is assigned to bin i with the fraction Φi (x). The transition
matrix P was then defined as
PN
Φi (xk )Φj (x̃k )
pij := k=0
PN
k=0 Φi (xk )
where the notation x̃k refers to a propagation of the system by a time step t,
i.e., x̃k is the state of xk after t steps. In the examples we set t = 1.
The method was applied to examples with three, four, and six minima to
show the influence of the applied drift on the clustering and the value of maximal
flow. The potential functions have their minima arranged on a circle around a
single maximum and are of the form
13



Ω3 (x, y) = 6 exp −3(x2 + y 2 )




−8 exp −(x − x? )2 − (y + y ? )2 − 8 exp −(x + x? )2 − (y + y ? )2


−8 exp −x2 − (y − 1)2


Ω4 (x, y) = 4 exp −3(x2 + y 2 )




−8 exp −x2 − (y − 1.5)2 − 8 exp −(x − 1)2 − y 2




−8 exp −(x + 1)2 − y 2 − 8 exp −x2 − (y + 1.5)2


Ω6 (x, y) = 4 exp −3(x2 + y 2 )




−8 exp −(x − 2x? )2 − (y + 2y ? )2 − 8 exp −(x + 2x? )2 − (y + 2y ? )2




−8 exp −(x + 2x? )2 − (y − 2y ? )2 − 8 exp −(x − 2x? )2 − (y − 2y ? )2




−8 exp −(x + 2)2 − (y − 1)2 − 8 exp −(x − 2)2 − y 2
√
with x? = 0.5 and y ? = 0.5 3. The minima resemble the metastable macroscopic states of the system. Trajectories with N = 10000 steps were simulated
to construct matrices with n = 20 bins.
Repressilator. The repressilator is a system of differential equations that
describes a synthetic genetic regulatory network [11] consisting of three genes,
TetR, γcI, and LacI. Each of the genes produces a protein p that represses the
production of mRNA m of another gene. The symmetric system was described
in [11] by the equations
dmA
v
+ v0
= −mA +
dt
1 + phC
dmB
v
= −mB +
+ v0
dt
1 + phA
dmC
v
+ v0
= −mC +
dt
1 + phB

dpA
= −β(pA − mA )
dt
dpB
= −β(pB − mB )
dt
dpC
= −β(pC − mC )
dt

with v = 298.2 transcriptions per second, β = 1/5 the ratio of protein decay rate
to mRNA decay rate, a growth constant v0 = 0.03, and a Hill coefficient h = 2.
Trajectories were started at 200 points in the six-dimensional cube of the
interval [0, 20] and simulated for 1.5 seconds using the ode45 function by Matlab.
The starting points were generated by a Niederreiter sequence (see [19, 20])
with values from [0,1] and scaled afterwards. The transition matrix was defined
similar as above by
exp(−0.2kxstart
− xend
i
j k)
pij = P200
,
start
− xend
k=1 exp(−0.2kxi
k k)

4.2

i, j = 1, . . . , 200.

Testing Environment

For computing optimal cycle clusterings we used the mixed-integer programming
solver SCIP, which is free for academic purposes [14].
Although in theory SCIP could solve general MIP formulations out-of-thebox, our cycle clustering MIPs proved hard in practice. Hence, to speed up the
14

instances
Ω3
Ω3
Ω4
Ω6
Rep (MIP)
Rep (G-PCCA)

drift best obj. obj. bound net flow coherence time [s]
0.1
0.2
0.1
0.1

0.0018
0.0057
0.0037
0.0056

0.0018
0.0057
0.0037
0.0056

0.0013
0.0052
0.0033
0.0051

0.5374
0.5360
0.3704
0.4595

0.1
0.1
1.4
6.9

–
–

0.1395
0.0163

0.2345
–

0.1391
0.0159

0.3636
0.3754

3600
<1

Table 1: Summary of results for all test instances. The repressilator MIP could not
be solved to optimality within the time limit of 3600 seconds. Note that the G-PCCA
result in the last line are added as a reference and computed by a different clustering
algorithm without guaranteed bound on the objective, see [28].

solving process of SCIP, we implemented three problem-specific heuristics that
try to find good primal solutions. They are called at the root node and during
the branch-and-bound search:
• First, we implemented a greedy heuristic to construct a feasible clustering by iteratively assigning the bins to clusters. The heuristic starts with
the assignment of bin 1 to cluster 1 and assigns all remaining bins iteratively. Therefore, the best possible assignment w.r.t. non-reversibility and
coherence is computed in each assignment step.
• Second, inspired by the approach with Schur vectors [28], we use the solution of the LP relaxation at each node within the branch-and-bound
procedure as a starting point for a rounding heuristic that knows the specific problem structure.
• Third, we implemented an improvement heuristic similar to [17] that iteratively tries to identify a bin that can be moved to a different cluster such
that the objective function is increased.
Note that none of these procedures is guaranteed to find an optimal clustering
or even to succeed at all at finding a feasible solution. However, applied regularly as part of the global solution process of SCIP they help to accelerate the
convergence of the primal and dual bound significantly.
All tests were run sequentially on identical machines with an Intel Xeon
Quad-core with 3.2 GHz and 48 GB of RAM. To balance net flow and coherence
in the objective function, we used a value of α = 0.001.

4.3

Computational Results

We analyzed results for the two different experimental settings of Section 4.1.
In the Catalytic Cycle example, there exist strong metastabilities with a weak
non-reversible net flow between them. Our results show that our method is
able to identify non-dominant cycles between metastable clusters. The results
for the second example Repressilator show that we are also able to find dominating cyclic structures. Moreover, our novel MIP formulation provides biologically meaningful clusterings and outperforms state-of-the-art approaches in
that sense.
15

2

2

1.5

1.5

1

1

0.5

0.5

3

2

1

0

0

−0.5

−0.5

0

−1
−1

−1

−1.5

−1.5

−2

−2
−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

−2
−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

−3
−3

−2

−1

0

1

2

3

Figure 4: Visualization of the solutions obtained with a 3-cycle, 4-cycle, and 6-cycle
clustering for the synthetic potentials Ω3 , Ω4 , and Ω6 , respectively.

Catalytic Cycle. In the first four rows of Table 1 and in Figure 4 one can
see the results of the HMC simulations of potential energy surfaces different
number of metastabilities (3, 4, and 6) and a small or larger drift (0.1 and 0.2).
Column “best obj.” states the objective function value of the best
P cycle clustering that was computed and is composed
of
column
“net
flow”
(
k f (Ck , Cφ(k) ))
P
plus 0.001 times column “coherence” ( k g(Ck )). Because SCIP could compute
proven optimal cycle clusterings for these four instances, the “best obj.” values
equal the “obj. bound” values that state the proven upper bound.
In all four cases our approach was able to correctly identify the metastabilities and the direction of the drift.
Furthermore, we can observe that the coherence within the clusters is large
although α was chosen small. Note that it equals 1 minus the probability of seeing transitions between clusters. By construction, the net flow is comparatively
small, but could still be detected consistently. The next example shows that
our method is also able to analyze systems with a fast and productive cycle.
Repressilator. For the repressilator model, the last two rows of Table 1 compare the cycle clustering solution (MIP) with a solution of the spectral clustering
algorithm (G-PCCA) explained in Section 2 and [10, 28]. Note that here the
“drift” column does not apply and G-PCCA was not designed with this objective function in mind and by its nature does not compute any proven objective
bound. In this sense we want to emphasize that this comparison does not provide any kind of benchmark.
The resulting MIP problem was significantly larger and computationally
more challenging than for the synthetic instances, and SCIP could not solve
it to optimality within one hour. Nevertheless, the best clustering solution
returned by SCIP exhibits both large net flow and coherence. Furthermore, the
objective bound guarantees that it is at most a factor of 1.8 from the objective
of an optimal cycle clustering.
To analyze the clusterings in more detail, consider the projected transition
matrices in the notation of Chapter 3.3,


0.126 0.134 0.075
W MIP = 0.088 0.150 0.123
0.121 0.077 0.107

16

f1

f2

15
10
5
0 p m p m p m
A
A B
B C
C

pA mA pB mB pC mC

pA mA pB mB pC mC

f3

(a) Solution obtained by G-PCCA.

f1

f2

15
10
5
pA mA pB mB pC mC 0 pA mA pB mB pC mC

pA mA pB mB pC mC

f3

(b) Solution obtained by cycle clustering.

Figure 5: Average values of the concentrations in each cluster.

with


W MIP −

T
W MIP

0
= 0.046 · −1
1

1
0
−1


−1
1 .
0

Here, the cycle C1 → C2 → C3 can be clearly identified. In comparison to this,
the matrix


0.082 0.090 0.103
W G-PCCA = 0.084 0.113 0.123
0.108 0.117 0.180
does not have such a clear cyclic order and the net flow is one order of magnitude
smaller,


0
1 −1
T
1 .
W G-PCCA − W G-PCCA = 0.005 · −1 0
1 −1 0
In Figure 5a, we have plotted the arithmetic averages of the protein and
mRNA concentrations in each cluster both for the G-PCCA clustering and the
cycle clustering. In the cycle clustering solution, one can clearly identify a peak
17

of pA , mA in cluster C1 , a peak of pB , mB in cluster C2 , and a peak of pC , mC in
cluster C3 . In contrast, the protein and mRNA concentrations in the G-PCCA
clusters seem more uniform and—besides a slightly decreased concentration of
mC , mA , and mB in cluster C1 , C2 , and C3 , respectively—make it difficult to
detect a particular structure.
In this sense, cycle clustering succeeds in separating states in a biologically meaningful way, using only simulation data, without additional knowledge
about the defining dynamical system. Each identified cluster of the system corresponds to a biological entity of a protein together with its mRNA. In constrast,
the G-PCCA solution does not seem to be interpretable in terms of (separate)
biological entities. The high productivity of the system is not reflected by the
objective value of the G-PCCA solution, while the MIP solution accounts well
for the cyclic nature of the repressilator.
MIP Performance. As can be seen from the running time in the last column
of Table 1, the improved clustering solutions come at the price of increased
computing times. While G-PCCA took less than a second, SCIP could not
solve the MIP problem to optimality within the time limit of one hour. This
is not surprising given Theorem 6, which proves that cycle clustering is a hard
combinatorial optimization problem.
However, good cycle clustering solutions are usually computed very early
during the solution process and for the repressilator solution the dual bound
proves that the best clustering w.r.t. our objective function can be at most
2.6 times as good as the solution stated in Table 1. In future research, we
will focus on improving the performance of MIP solvers for cycle clustering by
further dedicated techniques.

5

Conclusion

Many non-reversible biological processes seem to be mainly reversible on a small
time-scale. In order to identify the global cyclic behavior of the system, we
formulated an optimization problem for partitioning the state space into certain
cluster which are “visited” in a non-reversible manner. Standard approaches
which use spectral information of the transition matrix P are constructed such
that they find dominant cycles or the strongest metastabilities, but they do not
account for non-reversible cycles if these cycles are hidden, i.e., not dominant.
We prove that the new clustering method amounts to solving an N P-hard
combinatorial optimization problem. However, computational experiments show
that our solution strategy, which uses a mixed-integer programming formulation,
effectively finds optimal or near-optimal clusterings. The results for a genetic
regulatory network demonstrate that the identified clusters are meaningful in
the biological context. This clustering was not found by standard or spectral
approaches. One reason for that could be, that the problem of finding a global
cycle with maximal net flow is a more complex problem than determining a
dominant cycle with the help of spectral analysis, i.e. Schur-decomposition.

18

Acknowledgements
The work for this article has been partly conducted within the Research Campus
Modal funded by the German Federal Ministry of Education and Research (fund
number 05M14ZAM).

References
[1] T. Achterberg, T. Berthold, T. Koch, and K. Wolter. Constraint integer
programming: A new approach to integrate cp and mip. In International
Conference on Integration of Artificial Intelligence (AI) and Operations Research (OR) Techniques in Constraint Programming, pages 6–20. Springer,
2008.
[2] D. Applegate, R. Bixby, W. Cook, and V. Chvátal. On the solution of
travelling salesman problems. Universität Bonn. Institut für Ökonometrie
und Operations Research, 1998.
[3] G. Bowman, V. Pande, and F. Noé, editors. An Introduction to Markov
State Models and Their Application to Long Timescale Molecular Simulation, volume 797 of Advances in Experimental Medicine and Biology.
Springer Berlin Heidelberg, 2014.
[4] S. Brooks, A. Gelman, G. Jones, and X.-L. Meng. Handbook of Markov
Chain Monte Carlo. CRC press, 2011.
[5] E. Dahlhaus, D. S. Johnson, C. H. Papadimitriou, P. D. Seymour, and
M. Yannakakis. The complexity of multiway cuts. In Proceedings of the
twenty-fourth annual ACM symposium on Theory of computing, pages 241–
251. ACM, 1992.
[6] R. J. Dakin. A tree-search algorithm for mixed integer programming problems. The computer journal, 8(3):250–255, 1965.
[7] P. Deuflhard, W. Huisinga, A. Fischer, and C. Schütte. Identification of
almost invariant aggregates in reversible nearly uncoupled markov chains.
Linear Algebra and its Applications, 315(1-3):39 – 59, 2000.
[8] P. Deuflhard and M. Weber. Robust perron cluster analysis in conformation
dynamics. Linear Algebra and its Applications, 398:161 – 184, 2005. Special
Issue on Matrices and Mathematical Biology.
[9] N. Djurdjevac-Conrad, R. Banisch, and C. Schütte. Modularity of directed
networks: Cycle decomposition approach. Journal of Computational Dynamics 2 (2015) pp. 1-24, 2014.
[10] N. Djurdjevac-Conrad, M. Weber, and C. Schütte. Finding dominant structures of nonreversible markov processes. Multiscale Modeling and Simulation, 2016. accepted for publication.
[11] M. Elowitz and S. Leibler. A synthetic oscillatory network of transcriptional
regulators. Nature, 403(20):335–338, 2000.

19

[12] K. Fackeldey and M. Weber. Local refinements in classical molecular dynamics simulations. 2nd International Conference on Mathematical Modeling in Physical Sciences 2013, Journal of Physics: Conference Series
490:012016, 2014.
[13] R. Fortet. L’algebre de boole et ses applications en recherche opérationnelle.
Trabajos de Estadistica y de Investigación Operativa, 11(2):111–118, 1960.
[14] G. Gamrath, T. Fischer, T. Gally, A. M. Gleixner, G. Hendel, T. Koch, S. J.
Maher, M. Miltenberger, B. Müller, M. E. Pfetsch, C. Puchert, D. Rehfeldt,
S. Schenker, R. Schwarz, F. Serrano, Y. Shinano, S. Vigerske, D. Weninger,
M. Winkler, J. T. Witt, and J. Witzig. The SCIP Optimization Suite 3.2.
ZIB-Report 15-60, Zuse Institute Berlin, 2016.
[15] M. R. Garey and D. S. Johnson. Computers and Intractability: A Guide
to the Theory of NP-Completeness. W. H. Freeman & Co., New York, NY,
USA, 1979.
[16] L. Hagen and A. B. Kahng. New spectral methods for ratio cut partitioning
and clustering. IEEE transactions on computer-aided design of integrated
circuits and systems, 11(9):1074–1085, 1992.
[17] B. Kernighan and S. Lin. An Efficient Heuristic Procedure for Partitioning
Graphs. The Bell Systems Technical Journal, 49(2), 1970.
[18] A. H. Land and A. G. Doig. An automatic method of solving discrete programming problems. Econometrica: Journal of the Econometric Society,
pages 497–520, 1960.
[19] R. Lidl and H. Niederreiter. Finite Fields. Cambridge University Press,
Cambridge, 1984.
[20] H. Niederreiter. Low-discrepancy and low-dispersion sequences. Journal of
Number Theory, (30):51–78, 1988.
[21] A. Nielsen and M. Weber. Computing the nearest reversible markov chain.
Numerical Linear Algebra with Applications, 22(3):483 – 499, 2015.
[22] J.-H. Prinz, H. Wu, M. Sarich, B. Keller, M. Senne, M. Held, J. D. Chodera,
C. Schütte, and F. Noé. Markov models of molecular kinetics: Generation
and validation. The Journal of Chemical Physics, 134(17), 2011.
[23] M. W. Savelsbergh. Preprocessing and probing techniques for mixed integer
programming problems. ORSA Journal on Computing, 6(4):445–454, 1994.
[24] A. Schrijver. Combinatorial Optimization: Polyhedra and Efficiency, volume 24. Springer, 2003.
[25] J. Shi and J. Malik. Normalized cuts and image segmentation. IEEE
Transactions on pattern analysis and machine intelligence, 22(8):888–905,
2000.
[26] U. von Luxburg. Tutorial on spectral clustering. Max Planck Institute for
Biological Cybernetics, Technical Report No. 149, Tüubingen, 2006.

20

[27] M. Weber. A Subspace Approach to Molecular Markov State Models via a
New Infinitesimal Generator. Habilitation thesis, FU Berlin, 2011.
[28] M. Weber and K. Fackeldey. G-PCCA: Spectral clustering for nonreversible markov chains. ZIB-Report 15-35, Zuse Institute Berlin, 2015.

21

