A Formal Methods Approach to Pattern Synthesis
in Reaction Diffusion Systems

arXiv:1409.5671v1 [cs.AI] 12 Sep 2014

Ebru Aydin Gol, Ezio Bartocci and Calin Belta
Abstract— We propose a technique to detect and generate
patterns in a network of locally interacting dynamical systems.
Central to our approach is a novel spatial superposition logic,
whose semantics is defined over the quad-tree of a partitioned
image. We show that formulas in this logic can be efficiently
learned from positive and negative examples of several types of
patterns. We also demonstrate that pattern detection, which is
implemented as a model checking algorithm, performs very
well for test data sets different from the learning sets. We
define a quantitative semantics for the logic and integrate the
model checking algorithm with particle swarm optimization in
a computational framework for synthesis of parameters leading
to desired patterns in reaction-diffusion systems.

I. INTRODUCTION
From the stripes of a zebra and the spots on a leopard
to the filaments (Anabaena) [1], spirals, squares (Thiopedia
rosea), and vortex (Paenibacillus) [2] formed by singlecell organisms, patterns can be found everywhere in nature.
Pattern formation is at the very origin of morphogenesis and
developmental biology, and it is at the core of technologies
such as self-assembly, tissue engineering, and amorphous
computing. Even though it received a lot of attention from
diverse communities such as biology, computer science, and
physics, the problem of pattern formation is still not well
understood.
Pattern recognition is usually formulated as a machine
learning problem [3], in which patterns are characterized
either statistically [4] or through a structural relationship
among their features [5]. Despite its success in several
application areas [6], pattern recognition still lacks a formal
foundation. Can patterns be specified in a formal language
with well-defined syntax and semantics? Can we develop
algorithms for pattern detection from specification given
in such a language? Given a large collection of locally
interacting agents, can we design parameter synthesis rules,
control and interaction strategies guaranteeing the emergence
of global patterns? In this paper, by drawing inspiration from
model checking [7], [8], we provide partial answers to these
questions.
We address the following problem: Given a network of
locally interacting dynamical systems, and given sets of positive and negative examples of a desired pattern, find parameter values that guarantee the occurrence of the pattern in
the network at steady state. Our approach is based on a novel
spatial superposition logic, called Tree Spatial Superposition
Logic (TSSL), whose semantics is defined over quad-trees of
Ebru Aydin Gol (ebru@bu.edu) and Calin Belta (cbelta@bu.edu) are
with Boston University. Ezio Bartocci (ezio.bartocci@tuwien.ac.at) is with
Vienna University of Technology.

partitioned images. The decision of whether a pattern exists
in an image becomes a model checking problem. A pattern
descriptor is a TSSL formula, and we employ machinelearning techniques to infer such a formula from the given
positive and negative examples of the pattern. To synthesize
parameters of the original networked system leading to a
desired pattern, we use a particle swarm optimization (PSO)
algorithm. The optimization fitness function is given by a
measure of satisfaction induced by the quantitative semantics
that we introduce for the logic. We present examples showing
that formulas in the proposed logic are good classifiers for
some commonly encountered patterns. While the overall
algorithm can, in principle, be applied to any network of
locally interacting systems, in this paper we focus on the
Turing reaction-diffusion system [9], and show that patternproducing parameters can be automatically generated with
our method.
The rest of the paper is organized as follows. In Section II
we discuss the work. In Section III we formulate the problem
and outline our approach. We define the syntax and semantics
of TSSL in Section IV. A machine learning technique to
learn TSSL formulas from positive and negative examples of
desired patterns is developed in Section V. The solution to
the pattern generation problem is presented in Section VI as
a supervised, iterative procedure that integrates quantitative
model checking and optimization. We conclude with final
remarks and directions for future work in Section VII.
II. RELATED WORK
Pattern recognition is a well-established technique in
machine learning. Given a data set and a set of classes, the
goal is to assign each data to one class, or to provide a
“most likely” matching of the data to the classes. The two
main steps in pattern recognition are: (a) to extract distinctive
features [10], [11], [12], [13] with relevant information from
a set of input data representing the pattern of interest and (b)
to build, using one of the several available machine learning
techniques (see [14] for a detailed survey), an accurate
classifier trained with the extracted features. The descriptor
chosen in feature extraction phase depends on the application
domain and the specific problem.
This work is related to pattern recognition in computer
vision, where these descriptors may assume different forms.
Feature descriptors such as Textons [10] and Histograms of
Oriented Gradients (HoG) [11] are concerned with statistical
information of color distribution or of intensity gradients
and edge directions. The scale-invariant feature transform
(SIFT), proposed by Lowe in [13], is based on the appearance

of an object at particular interest points, and is invariant to
image scale and rotation. The shape context [12] is another
feature descriptor intended to describe the shape of an object
by the points of its contours and the surrounding context.
In this paper we establish an interesting connection between verification and pattern recognition. Both classical
verification [15], [16], [17], [18], [19] and pattern recognition
techniques aim to verify (and possibly quantify) the emergence of a behavioral pattern. We propose logic formulas
as pattern descriptors and verification techniques as pattern
classifiers. The logical nature of such pattern descriptors
allows to reason about the patterns and to infer interesting
properties. For example, in [20], the spatial modalities are
used to characterize self-similar (fractal) texture. Furthermore, combining different pattern descriptors using both
modal and logical operators is quite intuitive.
This paper is inspired by the original work on morphogenesis by Alan Turing [9], and is closely related to [20]. In the
latter, the authors introduced a Linear Spatial Superposition
Logic (LSSL), whose formulas were interpreted, as in this
paper, over quad-tree partitions. The existence of a pattern
in an image corresponded to the existence of a path in the
corresponding tree from the root to the leaf corresponding
to a representative point in the image. As a consequence,
the method was shown to work for spirals, for which the
center was chosen as the representative point. The tree logic
proposed here is more general as it does not depend on the
choice of such a point and captures the pattern “globally”.
For example, the patterns considered in this paper cannot be
expressed in LSSL, because they rely on a tree representation
rather than a path representation.
As opposed to [20], we also define a quantitative semantics
for the logic, and use the distance to satisfaction as a fitness
function while searching for pattern-producing parameters.
This quantitative semantics and the discounted model checking on a computational tree are inspired from [21], with the
notable difference that we do not need a metric distance, but
rather a measure of satisfiability. Such measures have also
been used in [15], [16], [17], [18], [19]. The main novelty
of this paper, compared to the other pattern recognition
approaches, is that we can quantify “how far” a system is
from producing a desired pattern. This, which is possible due
to the quantitative semantics of our logic, enables the use
of optimization algorithms to search for pattern-producing
parameters.
III. PROBLEM FORMULATION
Notation. We use R, R+ , N and N+ to denote the set of
real numbers, non-negative reals, integer numbers, and nonnegative integers, respectively. For any c ∈ R and set S ⊆ R,
S>c := {x ∈ S | x > c}, and for any a, b ∈ R, S[a,b] := {x ∈ S |
a ≤ x ≤ b}.
A reaction-diffusion system S is modeled as a spatially
distributed and locally interacting K × K rectangular grid of
identical systems, where each location (i, j) ∈ N[1,K] × N[1,K]

corresponds to a system:
(n)

dxi, j

(n)

(n)

(n)

1
(n)
∑ xv ,
|νi, j | v∈ν
i, j

= Dn (ui, j − xi, j ) + fn (xi, j , R),

n = 1, . . . , N,
(III.1)
(1)
(N)
where xi, j = [xi, j , . . . , xi, j ] is the state vector of system Si, j ,
which captures the concentrations of all species of interest.
D and R are the parameters of system S. D = [D1 , . . . , DN ] ∈
RN+ is the vector of diffusion coefficients. R ∈ RP−N is the
vector of parameters that defines the local dynamics fn : RN+ ×
RP−N → R for each of the species n = 1, . . . , N. Note that
the parameters and dynamics are the same for all systems
Si, j , (i, j) ∈ N[1,K] × N[1,K] . The diffusion coefficient is strictly
positive for diffusible species and it is 0 for non-diffusible
(1)
(N)
species. Finally, ui, j = [ui, j , . . . , ui, j ] is the input of system
Si, j from the neighboring systems:
Si, j :

dt

ui, j =

where νi, j denotes the set of indices of systems adjacent to
Si, j .
Given a parameter vector p = [D, R] ∈ RP , we use S(p)
to denote an instantiation of a reaction-diffusion system. We
use x(t) ∈ RK×K×N
to denote the state of system S(p) at time
+
(p)
t, and xi, j (t) ∈ RN+ to denote the state of system Si, j at time t.
While the model captures the dynamics of concentrations of
all species of interest, we assume that a subset {n1 , . . . , no } ⊆
{1, . . . , N} of the species is observable through:
H : RK×K×N
→ RK×K×o
:
+
[0,b]

y = H(x),

for some b ∈ R+ . For example, a subset of the genes in
a gene network are tagged with fluorescent reporters. The
relative concentrations of the corresponding proteins can be
inferred by using fluorescence microscopy.
We are interested in analyzing the observations generated
by system (III.1) in steady state. Therefore, we focus on
parameters that generate steady state behavior, which can be
easily checked through a running average:
K

K

N

(n)

(n)

∑ ∑ ∑ | xi, j (t) − xi, j

|< ε,

(III.2)

i=1 j=1 n=1
(n)

(n)

T
where xi, j = t−T
xi, j (τ)dτ/T for some T ≤ t. The system
is said to be in steady state at time t¯, if (III.2) holds for
all t ≥ t¯. In the rest of the paper, we will simply call the
observation of a trajectory at steady state as the observation
of the trajectory, and denote it as H(x(t¯)).
Example 3.1: We consider a 32 × 32 reaction-diffusion
system with two species (i.e. K = 32, N = 2):

R

(1)

dxi, j

dt
(2)
dxi, j



(1)
(1)
(1) (2)
(1)
= D1 ui, j − xi, j + R1 xi, j xi, j − xi, j + R2 ,


(2)
(2)
(1) (2)
= D2 ui, j − xi, j + R3 xi, j xi, j + R4 .

(III.3)
dt
The system is inspired from Turing’s reaction-diffusion system and is presented in [22] as a model of the skin pigments

of an animal. At a cell (location (i, j)), the concentration
(1)
of species 1, xi, j , depends on the concentration of species
1 in this cell and in its neighbors (if D1 > 0), and the
(2)
concentration of species 2 in this cell only, i.e. xi, j . Similarly,
(2)

xi, j depends on the concentration of species 2 in this cell
(1)

and in its neighbors (if D2 > 0), and xi, j (if R3 6= 0).
We assume that species 1 is observable through mapping
H : R32×32×2
→ R32×32
+
[0,1] given by:
(1)

y = H(x), where yi, j =

xi, j

(1)

.

maxm,n xm,n

We simulate the system from random initial conditions with
parameters R = [1, −12, −1, 16], and different diffusion parameters D1 = [5.6, 24.5], D2 = [0.2, 20], and D3 = [1.4, 5.3].
The observed concentrations of species 1 at different time
points are shown in Figure 1. At time t = 50, all trajectories
are in steady state. Note that, in all three cases, the spatial
distribution of the steady state concentrations of species 1 has
some regularity, i.e. it forms a “pattern”. We will use large
spots (LS), fine patches (FP), and small spots (SS) to refer to
the patterns corresponding to D1 , D2 , and D3 , respectively.
(a)

(b)

(c)
t=0

t=5

t=10

t=20

t=30

t=40

t=50

t=60

Fig. 1. Observations generated by system (III.3) with parameters R and
(a) D1 , (b) D2 , and (c) D3 from Example 3.1 (the concentration of species
1 is represented with shades of red). The steady state observations produce
(a) large spots (LS), (b) fine patches (FP), and (c) small spots (SS).

Problem 3.1: Given a reaction-diffusion system S as defined in (III.1), a finite set of initial conditions X0 ⊂ RK×K×N ,
ranges of the design parameters P = P1 ×. . .× PP , Pi ⊂ R, i =
1, . . . , P, a set of observations Y+ = {yi }i=1,...,N+ that contain
a desired pattern, a set of observations Y− = {yi }i=1,...,N−
that do not contain the pattern, find parameters p∗ ∈ P such
∗
that the trajectories of system S(p ) originating from X0 are
guaranteed to produce observations similar to the ones from
the set Y+ .
To solve Problem 3.1, we need to perform two steps:
•
•

Design a mechanism that decides whether an observation contains a pattern.
Develop a search algorithm over the state space of the
design parameters to find p∗ .

The first step requires to define a pattern descriptor. To
this goal, we develop a new spatial logic over spatialsuperposition trees obtained from the observations, and treat
the decision problem as a model checking problem. The new
logic and the superposition trees are explained in Section
IV. Then, finding a pattern descriptor reduces to finding a

formula of the new logic that specifies the desired pattern. We
employ machine-learning techniques to learn such a formula
from the given sets of observations Y+ and Y− .
The second step is the synthesis of parameters p∗ such that
the observations produced by the corresponding reaction∗
diffusion system S(p ) satisfy the formula learned in the first
step. To this end, we introduce quantitative semantics for
the new logic, which assigns a positive valuation only to the
superposition-trees that satisfy the formula. This quantitative
valuation is treated as a measure of satisfaction, and is
used as the fitness function in a particle swarm optimization
(PSO) algorithm. The choice of PSO is motivated by its
inherent distributed nature, and its ability to operate on
irregular search spaces, i.e. it does not require a differentiable
fitness function. Finally, we propose a supervised, iterative
procedure to find p∗ that solves Problem 3.1. The procedure
involves iterative applications of steps one and two, and
an update of the set Y− until a parameter set that solves
Problem 3.1 is found, which is decided by the user.
IV. TREE SPATIAL SUPERPOSITION LOGIC
A. Quad-tree spatial representation
We represent the observations of a reaction-diffusion system as a matrix Ak,k of 2k × 2k elements ai, j with k ∈ N>0 .
Each element corresponds to a small region in the space
(1)
(o)
and is defined as a tuple ai, j = hai, j , · · · , ai, j i of values
representing the concentration of the observable species
(c)
within an interval ai, j ∈ [0, b], with b ∈ R+ . Given a matrix
Ak,k , we use Ak,k [is , ie ; js , je ] to denote the sub-matrix formed
by selecting the rows with indices from is to ie and the
columns with indices from js to je .
Definition 4.1: A quad-tree Q = (V, R) is a quaternary
tree [23] representation of Ak,k where each vertex v ∈ V
represents a sub-matrix of Ak,k and the relation R ⊂ V × V
defines the four children of each node v that is not a leaf.
A vertex v is a leaf when all the elements of the sub-matrix
that it represents have the same values.
Figure 2 shows an example of a quadtree, where node
v0 represents the entire matrix; child v1 represents the submatrix {1, · · · , 2k−1 } × {1, · · · , 2k−1 }; child v7 represents the
sub-matrix {2k−2 + 1, · · · , 2k−1 } × {2k−2 + 1, · · · , 2k−1 }; etc.
In Figure 2, we also label each edge in the quad-tree with the
direction of the sub-matrix represented by the child: north
west (NW), north east (NE), south west (SW), south east
(SE).
v0
v1

NW$
NW$

NE$

NE$

v5

v8

v6

v7

SW$

NW$ SW$

v4
NW$

SW$

v5

v3

v2
(a)$

Fig. 2.

SE$

SW$

v6

SE$

v3

v2

v1

SE$

v0

SE$

NE$

v8

v7
(b)$

Quad-tree representation (b) of a matrix (a).

NE$

v4

Definition 4.2: We define the mean function µc : V →
[0, b] for sub-matrix Ak,k [is , ie ; js , je ] represented by the vertex
v ∈ V of the quad-tree Q = (V, R) as follows:
1
(c)
ai, j
µc (v) =
(ie − is + 1)( je − js + 1) i, j∈{i ,··· ,i∑
s
e }×{ js ,··· , je }
The function µc provides the expected value for an observable variable with index c, 1 ≤ c ≤ o in a particular region
of the space represented by the vertex v.
Definition 4.3: Two vertices va , vb ∈ V are said to be
equivalent when the mean function applied to the elements of
the sub-matrices that they represent produce the same values:
va ≡ vb ⇐⇒ µc (va ) = µc (vb ), ∀c, 1 ≤ c ≤ o
We use the mean of the concentration of the observable
species as a spatial abstraction (superposition) of the observations in a particular region of the system, avoiding in
this way to enumerate the observations of all locations. This
approach is inspired by previous papers [20], [24], where
the authors aim to combat the state-explosion problem that
would stem otherwise.
Proposition 4.1: Given a vertex v ∈ V of a quad-tree Q =
(V, R) and its four children vNE , vNW , vSE , vSW the following
property holds:
µc (vNE ) + µc (vNW ) + µc (vSE ) + µc (vSW )
µc (v) =
4
Proof: The proof can be easily derived by expanding
the terms of Definition 4.2.
Proposition 4.2: The number of vertices needed for the
quad-tree representation Q = (V, R) of a matrix Ak,k is upper
bounded by ∑ki=0 22i .
Proof: The proof follows from the fact that the worst
case scenario is when all the elements have different values.
In this case the cardinality of the set V is equal to the cardinality of a full and complete quaternary tree. For example,
to represent the matrix A3,3 , it would require a max number
of vertices |V | ≤ 1 + 4 + 16 + 64 = 85.
NW,NE,SE,SW

v0
B =1/2
W =1/2$

(a)$
v2

v1

v7

v8

v9

v10

v11

v12

v13

v14

v15

(b)$

v16

v17

v18 v19

$

NE,NW,SE,SW

SE

v6

B =1
W = 0$

$ sι

v4

SE

SE

SE

v5

v3

B =1/2
W =1/2$

v20

$ s1
(c)$

NE,SW

s2

B =0
W = 1$

$

s3

NW,NE,SE,SW

Fig. 3. A checkerboard pattern as a matrix of pixels (a), the quad-tree
representation (b) and the derived quad transition system (c), where B and
W denote black and white, respectively.

Algorithm B UILDING Q UAD T RANSITION S YSTEM
(1)

(o)

Input:

Matrix Ak,k of 2k × 2k of elements ai, j = hai, j , · · · , ai, j i,

Output:

its quad-tree Q = (V, R), the root v0 ∈ V , and a labeling
function LQ : R → D = {NW, NE, SE, SW }
Quad Transition System QT S = (S, sι , τ, Σ, [.], L)

1:
2:
3:
4:

Σ := {m1 , · · · , mo } . Initialize the set of variables Σ of QT S .
τ = 0/ . Initialize the set τ of the transition relation τ of QT S .

5:
6:

LF := {v ∈ V | 6 ∃t ∈ V : (v,t) ∈ R} . LF is the set of leaves of Q
PLF := {Pi ⊆ LF, 1 ≤ i ≤ n | Pi 6= 0/ ∧ ∀va , vb ∈ Pi ,

S := {sι } . Initialize the set of states S of QT S .
T S := {hsι , {v0 }i}

.

Each tuple in TS contains a state in S and a set of vertices in V.

∀vc ∈ Pj6=i , va ≡ vb ∧ va 6≡ vc }

.
7:

for each P̂ ∈ PLF do

.
.
8:
9:
10:
11:

14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:

For each partition element, create a state s0 with a self-loop and
a transition to the state sι if P̂ contains a child of v0 .

add new state s0 to S and a tuple hs0 , P̂i to T S
τ := τ ∪ {(s0 , s0 )} ∪ {(s, s0 ) : hs,V Si ∈ T S,
∃v ∈ V S, ∃v0 ∈ P̂ : (v, v0 ) ∈ R}
end for
FS := {v ∈ V |(v0 , v) ∈ R}\LF

.
12:
13:

PLF is a partition of LF with equivalent leaves.

explore the children of v0 that are not leaves.

while FS 6= 0/ do . FS contains the frontier vertices to be explored.
LFS := {v ∈ FS | ∀v0 ∈ V : (v, v0 ) ∈ R :
∃hs,V Si ∈ T S ∧ v0 ∈ V S}
/ Pi 6= 0,
/ ∀va , vb ∈ Pi ,
PLFS := {Pi∈I ⊆ LFS | I 6= 0,
∀vc ∈ Pj6=i , va ≡ vb ∧ va 6≡ vc }
for each P̂ ∈ PLFS do
add new state s0 to S and a tuple hs0 , P̂i to T S
S
τ := ( s:hs,V Si∈T S:∃v∈P̂,∃v0 ∈V S,(v,v0 )∈R (s0 , s)) ∪ τ
if ∃v ∈ P̂ ∧ ∃hs,V Si : ∃v0 ∈ V S ∧ (v0 , v) ∈ R then
τ := τ ∪ {(s, s0 )}
end if
end for
for each v̂ ∈ FS\LFS do
add new state s0 to S and a tuple hs0 , {v̂}i to TS
S
τ := ( s:hs,V Si∈T S:∃v0 ∈V S,(v̂,v0 )∈R (s0 , s)) ∪ τ
if ∃hs,V Si : ∃v0 ∈ V S ∧ (v0 , v̂) ∈ R then
τ := τ ∪ {(s, s0 )}
end if
end for
FS := {v ∈ V | ∃v̄ ∈ FS, (v̄, v) ∈ R}\LF
end while
define func [.] as [c̄](s̄) := µc̄ (vs̄ ), c̄ ∈ {1, · · · , o},
vs̄ ∈ V S : hs̄,V Si ∈ T S
define func L as L(s,t) := (t = s)?D :
S

ṽ∈V˜S,v̄∈V¯T :hs,V˜Si,ht,V¯T i∈T S,(ṽ,v̄)∈R LQ(ṽ, v̄)

B. Quad Transition System
We now introduce the notion of quad transition system
that extends the classical quad-tree structure, allowing for a
more compact exploration for model checking.
Definition 4.4: A Quad Transition System (QTS) is a
tuple QT S = (S, sι , τ, Σ, [.], L), where:
1) S is a finite set of states with sι ∈ S the initial state;

33:

return S, sι , τ, Σ, [.], L

2) τ ⊆ S × S is the transition relation. We require τ to be
non-blocking and bounded-branching: ∀s ∈ S, ∃t ∈ S :
(s,t) ∈ τ and ∀s ∈ S, if T (s) = {t : (s,t) ∈ τ} is the set
of all successors of s, the cardinality of |T (s)| ≤ 4;

3) Σ is a finite set of variables;
4) [.] is a function [.] : S → (Σ → [0, b]) that assigns to
each state s ∈ S and a variable m ∈ Σ a rational value
[s](m) in [0, b] with b ∈ R+ ;
5) L is a labeling function for the transition L : τ → 2D
with D = {NW, NE, SE, SW } and with the property
that ∀(s,t), (s,t 0 ) ∈ τ, with t 6= t 0 it holds that L(s,t) ∩
S
/ ∀t∈S:(s,t)∈τ L(s,t) = D .
L(s,t 0 ) = 0,
The B UILDING Q UAD T RANSITION S YSTEM algorithm
shows how to generate a QTS starting from a quad-tree
representation Q = (V, R) of a a matrix Ak,k and a labeling
function LQ : R → D .
Proposition 4.3: A quad transition system (QTS) QT S =
(S, sι , τ, Σ, [.], L) generated by the B UILDING Q UAD T RANSI TION S YSTEM algorithm has always a least fixed point, that
is ∃s ∈ S : (s, s) ∈ τ.
Proof: This property holds because the algorithm
generates a state with a self-loop transition for each partition
of equivalent leaves in the quad-tree.
Definition 4.5 (Labeled paths): Given a set B of labels
representing the spatial directions, a labeled path (lpath) of a
QTS Q is an infinite sequence πB = s0 s1 s2 · · · of states such
/ ∀i ∈ N. Given a state
that (si , si+1 ) ∈ τ ∧ L(si , si+1 ) ∩ B 6= 0,
s, we denote LPathsB (s) the set of all labeled paths starting
in s, and with πBi the i-th element of a path πB ∈ LPathsB (s).
For example, in Figure 3, LPathsB (sι ) = {sι s1 s2 s2 · · · } if B =
{NW, SE}.
C. TSSL Syntax and Semantics
Definition 4.6 (TSSL syntax): The syntax of TSSL is
defined as follows:
ϕ ::= > | ⊥ | m ∼ d | ¬ϕ | ϕ1 ∧ ϕ2 | ∃B

ϕ |∀B

ϕ |∃B ϕ1 Uk ϕ2 |∀B ϕ1 Uk ϕ2

/
with ∼∈ { ≤ , ≥}, d ∈ [0, b], b ∈ R+ , k ∈ N>0 , B ⊆ D : B 6= 0,
and m ∈ Σ, with Σ the set of variables.
From this basic syntax one can derive other two temporal
operators: the exist eventually operator ∃B Fk , the forall
eventually operator ∀B Fk , the exist globally operator ∃B Gk ,
and the forall globally operator ∀B Gk defined such that:
∃B Fk ϕ := ∃B > Uk ϕ
∀B Fk ϕ := ∀B > Uk ϕ

∃B Gk ϕ := ¬∀B Fk ¬ϕ.
∀B Gk ϕ := ¬∃B Fk ¬ϕ.

The TSSL logic resembles the classic CTL logic [25], with
the main difference that the next and until are not temporal,
but spatial operators meaning a change of resolution (or
zoom in). The set B selects the spatial directions in which
the operator is allowed to work and the parameter k limits
the until to operate on a finite sequence of states. In the
following we provide the TSSL qualitative semantics that,
given a spatial model and a formula representing the pattern
to detect, provides a yes/no answer.
Definition 4.7 (TSSL Qualitative Semantics): Let Q =
(S, sι , τ, Σ, [.], L) be a QTS, Then, Q satisfies a TSSL formula
ϕ, written Q |= ϕ, if and only if Q , sι |= ϕ, where:

Q , s |= >
Q , s |= m ∼ d
Q , s |= ¬ϕ
Q , s |= ϕ1 ∧ ϕ2
Q , s |= ∃B ϕ
Q , s |= ∀B ϕ
Q , s |= ∃B ϕ1 Uk ϕ2

and
⇔
⇔
⇔
⇔
⇔
⇔

Q , s |= ∀B ϕ1 Uk ϕ2

⇔

Q, s 6|= ⊥
[s](m) ∼ d
Q , s 6|= ϕ
Q , s |= ϕ1 ∧ Q , s |= ϕ2
∃s0 : (s, s0 ) ∈ τ ∧ L(s, s0 ) ∩ B 6= 0/ ∧ Q , s0 |= ϕ
∀s0 : (s, s0 ) ∈ τ ∧ L(s, s0 ) ∩ B 6= 0/ ∧ Q , s0 |= ϕ
∃πB ∈ LPathsB (s) : ∃i, 0 < i ≤ k :
(Q, πBi |= ϕ2 ) ∧ (∀ j < i, (Q, π j |= ϕ1 ))
∀πB ∈ LPathsB (s) : ∃i, 0 < i ≤ k :
(Q, πBi |= ϕ2 ) ∧ (∀ j < i, (Q, π j |= ϕ1 ))

Example 4.1: Checkerboard pattern. The checkerboard
pattern from Fig 3 a) can be characterized with the following
TSSL formula (B∗ = {SW, NE, NW, SE}):
∀B∗

(∀B∗

((∀{SW,NE}

(m ≥ 1)) ∧ (∀{NW,SE}

(m ≤ 0)))).

The “eventually” operator can be used to define all the
possible checkerboards of different sizes less or equal than
42 as follows:
∀B∗ F2 ((∀{SW,NE}

(m ≥ 1)) ∧ (∀{NW,SE}

(m ≤ 0)))

The qualitative semantics is useful to check if a given spatial
model violates or satisfies a pattern expressed in TSSL. However, it does not provide any information about how much
the property is violated or satisfied. This information may
be useful to guide a simulation-based parameter exploration
for pattern generation. For this reason we equip our logic
also with a quantitative valuation that provides a measure of
satisfiability in the same spirit of [17]. Since the valuation
of a TSSL formula with spatial operators requires to traverse
and to compare regions of space at different resolution,
we apply a discount factor of 14 on the result each time a
transition is taken in QTS.
Definition 4.8 ( TSSL Quantitative Semantics): Let
Q = (S, sι , τ, Σ, [.], L) be a QTS. The quantitative valuation
JϕK : S → [−b, b] of a TSSL formula ϕ is defined as follows:
J>K(s) = b

J⊥K(s) = −b

Jm ∼ dK(s) = (∼ is ≥) ? ([m](s) − d) : (d − [m](s))
J¬ϕK(s) = −JϕK(s)

Jϕ1 ∧ ϕ2 K(s) = min(Jϕ1 K(s), Jϕ2 K(s))
1
J∃B ϕK(s) =
max J ϕK(πB1 )
4 πB ∈LPathsB (s)
J∀B

ϕK(s) =

1
min
J ϕK(πB1 )
4 πB ∈LPathsB (s)

J∃B ϕ1 Uk ϕ2 K(s) =

πB ∈LPathsB (s)

J∀B ϕ1 Uk ϕ2 K(s) =

πB ∈LPathsB (s)

sup
inf

{min(
{min(

1
1
Jϕ2 K(πBi ), inf{ j Jϕ1 K(πBj ) | j < i}) | 0 < i ≤ k}}
4i
4
1
1
Jϕ2 K(πBi ), inf{ j Jϕ1 K(πBj ) | j < i}) | 0 < i ≤ k}
4i
4

Theorem 4.1 (Soundness): Let Q = (S, sι , τ, Σ, [.], L) be
a QTS, s ∈ S a state of Q , and ϕ a TSSL formula. Then, the
following properties hold for the two semantics:
JϕK(s) > 0 =⇒ Q , s |= ϕ
JϕK(s) < 0 =⇒ Q , s 6|= ϕ

Proof: The proof can be derived by structural induction
on the operational semantics.
Remark 4.1: Theorem 4.1 provides the basis of the techniques for pattern generation discussed in the following
sections. It is worth to note that, in the case JϕK(s) = 0,
it is not possible to infer whether Q violates or satisfies a
TSSL formula ϕ and only in this particular case we need to
resort to the qualitative semantics for determining it.

V. TSSL PATTERN CLASSIFIERS
A QTS can be seen in the context of multi-resolution
representation, since the nodes that appear at deeper levels
provide information for higher resolutions. Therefore, a
TSSL formula can effectively capture properties of an image.
However, it is difficult to write a formula that describes a
desired property, such as a pattern. Here, we propose to
use machine-learning techniques to find such a formula from
given sets of positive (Y+ ) and negative (Y− ) examples.
We first define a labeled data set from the given data sets
Y+ and Y− as

L = {(Qy , +) | y ∈ Y+ } ∪ {(Qy , −) | y ∈ Y− },
where Qy is the QTS generated from y. Then, we separate
the data set L into disjoint training and testing sets LL , LT . In
machine-learning, the training set is used to learn a classifier
for a target class, e.g. +, and the testing set is used to measure
the accuracy of the classifier. We employ RIPPER [26], a rule
based learner, to learn a classifier from LL , and then translate
the classifier into a TSSL formula characterizing +. Each
rule obtained from the learning algorithm is described as
ri : Ci ⇒∼i ,
where Ci is a boolean formula over linear predicates over
the variables of the states of a QTS, e.g. [s](m) > d, and ∼i
takes values from the label set {+, −}. A linear predicate for
a state s ∈ S can be written as a TSSL formula via the QTS
path from the root sι to s. Therefore, each Ci can be translated
into an equivalent TSSL formula Φi . The classification rules
are interpreted as nested if-else statements. Hence, a logically
equivalent TSSL formula for the desired property is defined
as follows:
!
Φ+ :=

_
j∈R+

Φj ∧

^

¬Φi ,

(V.4)

i=1,..., j−1

where R+ is the set of indices of rules ri with ∼i = +, and
Φi is the TSSL formula obtained from Ci .
(a)

(b)

(1)

(1)

Fig. 4. Sample sets of images from the sets (a) Y+ and (b) Y− for the
LS pattern.

Example 5.1: LS pattern. For the LS pattern from Ex(1)
ample 3.1, we generate a data set Y+ containing 8000
positive examples by simulating the reaction-diffusion system (III.3) from random initial conditions with parameters
(1)
R and D1 . Similarly, to generate the data set Y− containing
8000 negative examples, we simulate system (III.3) from
random initial conditions. However, in this case we use R and
randomly choose the diffusion coefficients from R2[0,30] . As
stated before, we only consider the observation of a system
in steady-state, for this reason, simulated trajectories that

do not reach steady state-in 60 time units are discarded. A
(1)
(1)
sample set of images from the sets Y+ and Y− is shown
(1)
in Figure 4. We generate a labeled set L
of QTS from
(1)
(1)
(1)
these sets, and separate L into LL , LT . We use RIPPER
algorithm implemented in Weka [27] to learn a classifier
(1)
from LL . The learning step took 228.5sec on an iMac with
a Intel Core i5 processor at 2.8GHz with 8GB of memory.
The classifier consists of 24 rules. The first rule is
r1 :(R ≥ 0.59) ∧ (R ≤ 0.70) ∧ (R.NW.NW.NW.SE ≤ 0.75)∧
(R.NW.NW.NW.NW ≥ 0.45) ⇒ +,
where R denotes the root of a QTS, and the labels of the
children are explained in Figure 2. Rule r1 translates to the
following TSSL formula:
Φ1 :(m ≥ 0.59) ∧ (m ≤ 0.70) ∧ (∃NW ∃NW ∃NW
(∃NW ∃NW ∃NW ∃NW m ≥ 0.45).

∃SE

m ≥ 0.75)∧

(1)

We define the TSSL formula Φ+ characterizing the pat(1)
(1)
tern as in (V.4), and model check QTSs from LT (|LT | =
(1)
8000) against Φ+ , which yields a high prediction accuracy
(96.11%) with 311 miss-classified QTSs.
FP and SS patterns. We follow the above explained
(i)
(i)
steps to generate data sets Y+ , Y− , generate labeled data
(i)
(i)
(i)
sets LL , LT , and finally learn formulas Φ+ for the FP and
SS patterns corresponding to diffusion coefficient vectors Di ,
i = 2, 3 from Example 3.1. Due to the space limitations, we
only present the results on the test sets. The model checking
of the QTSs from the corresponding test sets yields high
(2)
prediction accuracies 98.01%, and 93.13% for Φ+ , and
(3)
Φ+ , respectively.
VI. PARAMETER SYNTHESIS FOR
PATTERN GENERATION
In this section we present the solution to Problem 3.1, i.e.
a framework to synthesize parameters p ∈ P of a reactiondiffusion system S (III.1) such that the observations of
system S(p) satisfy a given TSSL formula Φ. First, we
show that the parameters of a reaction-diffusion system that
produce trajectories satisfying the TSSL formula can be
found by optimizing quantitative model checking results.
Second, we include the optimization in a supervised iterative
procedure for parameter synthesis.
We slightly abuse the terminology and say that a trajectory
x(t),t ≥ 0 of system S(p) satisfies Φ if the QTS Q =
(S, sι , τ, Σ, [.], L) of the corresponding observation, H(x(t¯)),
satisfies Φ, i.e Q |= Φ, or JΦK(sι ) > 0.
We first define an induced quantitative valuation of a
system S(p) and a set of initial conditions X0 from a TSSL
formula Φ as:
JΦK(S(p) ) = min {JΦK(sι ) | Q = (S, sι , τ, Σ, [.], L) is QTS of H(x(t¯)), x(0) = x0 }
x0 ∈X0

(VI.5)

The definition of the induced valuation of a system S(p)
implies that all trajectories of S(p) originating from X0 satisfy
Φ if JΦK(S(p) ) > 0. Therefore, it is sufficient to find p
that maximizes (VI.5). It is assumed that the ranges P =
P1 × . . . × PP of the design parameters are known. Therefore,

the parameters maximizing (VI.5) can be found with a greedy
search on a quantization of P . However, the computation of
JΦK(S(p) ) for a given p ∈ P is expensive, since it requires to
perform the following steps for each x0 ∈ X0 : simulating the
system S(p) from x0 , generating QTS Q of the corresponding
observation, and quantitative model checking of Q against
Φ. Here, we use the particle swarm optimization (PSO)
algorithm [28] over P with (VI.5) as the fitness function.
The choice of PSO is motivated by its inherent distributed
nature, and its ability to operate on irregular search spaces.
In particular, PSO does not require a differentiable fitness
function.
Example 6.1: LS pattern. We consider the reactiondiffusion system from Example 3.1 and the TSSL formula
(1)
Φ+ corresponding to the LS pattern from Example 5.1.
We assume that the parameters of the local dynamics are
known, R = [1, −12, −1, 16], and the diffusion coefficients
D1 and D2 are set as the design parameters with P = R2[0,30] .
We implement PSO to find p ∈ P maximizing the induced
valuation (VI.5). The PSO computation was distributed on 16
processors at 2.1GHz on a cluster, and the running time was
around 18 minutes. The optimized parameters are D1 = 2.25
and D2 = 29.42, and the valuation of the system is 0.0023.
A set of observations obtained by simulating S([2.25,29.42]) is
shown in Figure 6-(a). Note that, while all the observations
have some spatial periodicity indicating the presence of a
pattern, they are still different from the desired LS pattern.
FP and SS patterns. We also apply the PSO algorithm on
the same setting explained above to maximize the induced
(2)
valuation (VI.5) for the TSSL formulas Φ+ (FP pattern)
(3)
and Φ+ (SS pattern) from Example 5.1. The optimized
(2)
parameters are [0.083, 11.58] and [1.75, 7.75] for Φ+ and
(3)
Φ+ , respectively. Sets of observations obtained by simulating systems S([0.083,11.58]) and S([1.75,7.75]) are shown in
Figure 5. In contrast with the LS pattern, the observations
are similar to the ones from the corresponding data sets i.e.
(2)
(3)
Y+ and Y+ .
(a)

(b)

Fig. 5. Sample set of observations obtained by simulating (a) S([0.083,11.58])
and (b) S([1.75,7.75]) .

(a)

(b)

(c)

Fig. 6. Sample set of observations obtained by simulating (a) S([2.25,29.42]) ,
(b) S([3.75,28.75]) , and (c) S([6.25,29.42]) .

Remark 6.1: In this paper, we consider the observations
generated from a given set of initial conditions X0 . However,

the initial condition can be set as a design parameter and
optimized in PSO over a given domain RK×K×N
.
[a,b]
As seen in Example 6.1, it is possible that simulations
of the system corresponding to optimized parameters do
not necessarily lead to desired patterns. This should not
be unexpected, as the formula reflects the original training
set of positive and negative examples, and was not “aware”
that these new simulations are not good patterns. A natural
extension of our method should allow to add the newly
obtained simulations to the negative training set, and to
reiterate the whole procedure. This approach is summarized
in the I NTERACTIVE D ESIGN algorithm.
Algorithm I NTERACTIVE D ESIGN
Input:

Parametric reaction-diffusion system S, ranges of parameters P ,
a set of initial states X0 , sets of observations Y+ and Y−

Output: Optimized parameters p, the corresponding valuation γ
(no solution if γ < 0)

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:

while True do
Φ = Learning(Y+ , Y− )
{p, γ} = Optimization(S, X0 , Φ)
. γ is the induced valuation of S(p)
if γ < 0 then return p, γ
end if
UserQuery: Show observations of trajectories
of S(p) originating from X0 .
if User approves then return p, γ
else
Y− = Y− ∪ {H(x(t¯)) | x(t),t ≥ 0,
is generated by S(p) , x(0) ∈ X0 }.
end if
end while

We start with the user defined sets of observations Y+
and Y− , and learn a TSSL formula Φ from the QTS
representations of the observations (Section V). Then, in
the optimization step, we find a set of parameters p that
maximizes γ = JΦK(S(p) ). If γ < 0, then we terminate the
algorithm as parameters producing observations similar to
the ones from the set Y+ with respect to the TSSL formula
Φ could not be found. If γ ≥ 0, then the observations
of system S(p) satisfy Φ. Finally, the user inspects the
observations generated from the reaction-diffusion system
with the optimized set of parameters S(p) . If the observations
are similar to the ones from the set Y+ , then we find a
solution. If, however, the user decides that the observations
do not contain the pattern, then we add observations obtained
from system S(p) to Y− , and repeat the process, i.e learn a
new formula, run the optimization until the user terminates
the process or the optimization step fails (γ < 0).
Example 6.2: LS pattern. We apply I NTERAC TIVE D ESIGN algorithm to the system from Example 6.1. A
sample set of observations obtained in the first iteration is
shown in Figure 6-(a). We decide that these observations are
(1)
not similar to the ones from the set Y+ shown in Figure 4-

(a), and add these 250 observations generated with the
(1)
optimized parameters to Y− (line 9). In the second iteration,
the optimized parameters are D1 = 3.75 and D2 = 28.75,
and the observations obtained by simulating S([3.75,28.75])
are shown in Figure 6-(b). We continue by adding these to
(1)
Y− . The parameters computed in the third iteration are
D1 = 6.25 and D2 = 29.42. The observations obtained by
simulating S([6.25,29.42]) are shown in Figure 6-(c). Although
the optimized parameters are different from D1 , which was
(1)
used to generate Y+ , the observations of S([6.25,29.42]) are
(1)
similar to the ones from the set Y+ and we terminate the
algorithm.
VII. CONCLUSION AND FUTURE WORK
We defined a tree spatial superposition logic (TSSL)
whose semantics is naturally interpreted over quad trees
of partitioned images. We showed that formulas in this
logic can be efficiently learned from positive and negative
examples. We defined a quantitative semantics for TSSL
and combined with an optimization algorithm to develop
a supervised, iterative procedure for synthesis of patternproducing parameters.
While the experiments show that the current version of the
logic works quite well and can accommodate translational
and rotational symmetries commonly found in biology patterns, there are several directions of future work. First, we
expect that even better results could be obtained if more
statistical moments were used, rather than just the mean
as in the current version of this work. Second, we do not
exploit the full semantics of the logic in this paper. In
future work, we plan to investigate reasoning about multiple
branches and using the “until” operator. Third, we plan
to apply this method to more realistic networks, such as
populations of locally interacting engineered cells. We expect
that experimental techniques from synthetic biology can be
used to “tune” existing synthetic gene circuits to produce
global desired patterns.
R EFERENCES
[1] J. Golden and H. Yoon, “Heterocyst formation in anabaena,” Curr
Opin Microbiol., vol. 1, no. 6, pp. 623–629, 1998.
[2] R. Scherrer and V. Shull, “Structure, partial elemental composition,
and size of thiopedia rosea cells and platelets,” Can J Microbiol.,
vol. 32, no. 7, pp. 607–610, 1986.
[3] C. M. Bishop, Pattern Recognition and Machine Learning. Springer,
2006.
[4] A. Jain, R. Duin, and J. Mao, “Statistical pattern recognition: A
review,” IEEE Transactions on pattern analysis and machine learning,
vol. 22, pp. 4–37, 2000.
[5] T. Pavlidis, Structural Pattern Recognition. Springer-Verlag, 1980.
[6] R. C. Veltkamp and M. Hagedoorn, “State-of-the-art in shape matching,” Principles of Visual Information Retrieval, Tech. Rep., 1999.
[7] E. A. Emerson, “Temporal and modal logic,” in Handbook of Theoretical Computer Science: Formal Models and Semantics, J. van Leeuwen,
Ed. North-Holland Pub. Co./MIT Press, 1990, vol. B, pp. 995–1072.
[8] E. M. M. Clarke, D. Peled, and O. Grumberg, Model checking. MIT
Press, 1999.
[9] A. M. Turing, “The chemical basis of morphogenesis,” Philosophical
Transactions of the Royal Society of London, vol. 327, pp. 37–72,
1952.
[10] B. Julesz, “Textons, the elements of texture perception, and their
interactions,” Nature, vol. 290, pp. 91–97, 1981.

[11] N. Dalal and B. Triggs, “Histograms of oriented gradients for human
detection,” in Proc. of CVPR 2005: the IEEE Computer Society
Conference on Computer Vision and Pattern Recognition, vol. 1, June
2005, pp. 886–893.
[12] S. Belongie, J. Malik, and J. Puzicha, “Shape Matching and Object
Recognition Using Shape Contexts,” IEEE Transactions on Pattern
Analysis and Machine Intelligence, vol. 24, pp. 509–521, 2002.
[13] D. G. Lowe, “Object recognition from local scale-invariant features,”
in Proc. of the International Conference on Computer Vision, vol. 2,
1999, pp. 1150–1157.
[14] S. J. Russell and P. Norvig, Artificial Intelligence: A Modern Approach.
Prentice Hall, 2002.
[15] A. Rizk and F. Fages, “From Model-Checking to Temporal Logic
Constraint Solving,” in Proc. of CP 2009: the 15th International
Conference on Principles and Practice of Constraint Programming,
Lisbon, Portugal, 20-24 September, ser. Lecture Notes in Computer
Science, vol. 5732. Springer, 2009, pp. 319–334.
[16] A. Donzé and O. Maler, “Robust satisfaction of temporal logic over
real-valued signals,” in Proc. of FORMATS 2010, the 8th International
Conference on Formal Modeling and Analysis of Timed Systems,
Klosterneuburg, Austria, September 8–10, vol. 6246, 2010, pp. 92–
106.
[17] A. Donzé, E. Fanchon, L. M. Gattepaille, O. Maler, and P. Tracqui,
“Robustness analysis and behavior discrimination in enzymatic reaction networks,” PLoS One, vol. 6, no. 9, p. e24246, 2011.
[18] G. Fainekos and G. Pappas, “Robust Sampling for MITL specifications,” in Proc. of FORMATS 2007, the 5th International Conference
on Formal Modeling and Analysis of Timed Systems, ser. Lecture Notes
in Computer Science, 2007, vol. 8044, pp. 264–279.
[19] G. E. Fainekos and G. J. Pappas, “Robustness of temporal logic
specifications for continuous-time signals,” Theor. Comput. Sci., vol.
410, no. 42, pp. 4262–4291, 2009.
[20] R. Grosu, S. Smolka, F. Corradini, A. Wasilewska, E. Entcheva, and
E. Bartocci, “Learning and detecting emergent behavior in networks
of cardiac myocytes,” Communications of the ACM, vol. 52, no. 3, pp.
97–105, 2009.
[21] L. de Alfaro, M. Faella, and M. Stoeling, “Linear and branching
system metrics,” IEEE Transactions on Software Engineering, vol. 35,
no. 2, pp. 258–273, 2009.
[22] R. Collantes, “Algorithm alley. Dr. Dobb’s journal,” December 1996.
[23] R. Finkel and J. Bentley, “Quad trees a data structure for retrieval on
composite keys,” Acta Informatica, vol. 4, no. 1, pp. 1–9, 1974.
[24] Y. Kwon and G. Agha, “Scalable modeling and performance evaluation
of wireless sensor networks,” in Real-Time and Embedded Technology
and Applications Symposium, 2006. Proceedings of the 12th IEEE,
2006, pp. 49–58.
[25] E. M. Clarke and E. A. Emerson, “Design and synthesis of synchronization skeletons using branching-time temporal logic,” in Proc. of
Logic of Programs Workshop, ser. Lecture Notes in Computer Science,
vol. 131, 1982, pp. 52–71.
[26] W. W. Cohen, “Fast effective rule induction,” in Proceedings of the
Twelfth International Conference on Machine Learning.
Morgan
Kaufmann, 1995, pp. 115–123.
[27] M. Hall, E. Frank, G. Holmes, B. Pfahringer, P. Reutemann, and I. H.
Witten, “The weka data mining software: An update,” SIGKDD Explor.
Newsl., vol. 11, no. 1, pp. 10–18, Nov. 2009.
[28] J. Kennedy and R. Eberhart, “Particle swarm optimization,” in Proceedings of the IEEE International Conference on Neural Networks,
vol. 4, 1995, pp. 1942–1948.

