# Structural polyhedral stability of a biochemical network
is equivalent to finiteness of the associated generalised Petri net

Franco Blanchinia, Carlos Andrés Deviab and Giulia Giordanoc The work of
C.A.D. was supported by the Delft Technology Fellowship grant awarded to G.G.,
who also acknowledges the support of the Strategic Grant MOSES at the
University of Trento.a Dipartimento di Scienze Matematiche, Informatiche e
Fisiche, Università degli Studi di Udine, Via delle Scienze 206, 33100 Udine,
Italy<EMAIL_ADDRESS>Delft Center for Systems and Control, Delft
University of Technology, Mekelweg 2, 2628 CD Delft, The Netherlands.
<EMAIL_ADDRESS>Department of Industrial Engineering, University
of Trento, Via Sommarive 9, 38123 Trento, Italy<EMAIL_ADDRESS>
###### Abstract

We consider biochemical systems associated with a generalised class of Petri
nets with possibly negative token numbers. We show that the existence of a
structural polyhedral Lyapunov function for the biochemical system is
equivalent to the boundedness of the associated Petri net evolution or,
equivalently, to the finiteness of the number of states reachable from each
initial condition. For networks that do not admit a polyhedral Lyapunov
function, we investigate whether it is possible to enforce polyhedral
structural stability by applying a strong negative feedback on some _pinned_
nodes: in terms of the Petri net, this is equivalent to turning pinned nodes
into _black holes_ that clear any positive or negative incoming token. If such
nodes are chosen so that the transformed Petri net has bounded discrete
trajectories, then there exists a stabilising pinning control: the biochemical
network becomes Lyapunov stable if a sufficiently strong local negative
feedback is applied to the pinned nodes. These results allow us to
structurally identify the critical nodes to be locally controlled so as to
ensure the stability of the whole network.

##  I Introduction and Motivation

Structural analysis investigates how several systems often encountered in
nature enjoy important properties in view of their interconnection structure,
regardless of parameter values. Here, we consider structural stability and
_stabilisation_ of biochemical networks [12, 18, 23], adopting piecewise-
linear Lyapunov functions [8, 9, 10], which – along with the complementary
piecewise-linear-in-rate Lyapunov functions [1, 2, 3, 9] that can be seen as
their dual [11] – have proven effective in the stability analysis of chemical
reaction networks. A recent contribution [4] shows that this type of functions
can be very useful to detect, more in general, non-oscillatory behaviours.

Chemical reaction networks have been often analysed resorting to discrete-
event frameworks, employing for instance Petri nets [5, 37]: a chemical
reaction is seen as a process that assembles the needed number of reactant
molecules and releases the proper number of product molecules.

In this paper, we show how the existence of a piecewise-linear Lyapunov
function for the large class of _unitary_ (bio)chemical networks can be
interpreted as the boundedness of the evolution of a suitable _generalised_
Petri net, which can have both positive and negative token numbers.

This generalisation of Petri nets has been widely investigated in the
literature, e.g. under the name of _lending Petri net_ [6], and _negative
tokens_ [32] have been also called _anti-tokens_ [25, 36] or _debit tokens_
[7, 28]. However, to the best of the authors’ knowledge, this is the first
time that such a concept is associated with the stability of biochemical
networks.

We also consider a _structural_ pinning control problem. Pinning some nodes
means applying a strong feedback to these nodes with the goal of controlling
the whole network.

Pinning control has been extensively investigated in past years [29, 38], with
one of the main driving questions being how many nodes to pin and which ones
[33]. The approach has been used to address network control problems ranging
from asymptotic convergence [21] and noise rejection [16] to consensus [17]
and synchronisation [34]. Pinning techniques have been applied in different
areas such as circuits [21], power grids [33], protein networks and gene
regulatory networks [16]. Pinning control of a (bio)chemical reaction network
system can be seen as equivalent to the conversion of the pinned nodes into
_black holes_ , which swallow any incoming token (either positive or
negative), in the associated generalised Petri net.

The main contributions of this paper are summarised next.

  * •

We formulate the structural stability problem for biochemical reaction
networks (Section II) and we associate a biochemical network with a
Generalised Petri Net (GPN), with possibly negative tokens (Section III).

  * •

The GPN is fully determined by the network structure, and it does not depend
on the (monotonic) reaction rate functions.

  * •

The boundedness of all possible evolutions of the GPN is equivalent to the
existence of a polyhedral Lyapunov function (PLF) for the biochemical system;
such a PLF can be computed based on the efficient numerical procedure proposed
in [8] (Section III-A).

  * •

The equivalence with the GPN suggests more efficient stopping criteria for the
numerical procedure (Section III-C).

  * •

For networks that do not admit a PLF, we study how to convert some nodes into
_black holes_ that swallow any incoming token (Section III-D), so as to ensure
the stability of the network.

  * •

We show that converting some nodes into black holes is equivalent to _pinning_
them, i.e. virtually fixing their state variables to imposed values (Section
IV); after pinning appropriately chosen nodes, the network can admit a PLF.

  * •

We illustrate our results by assessing the structural stability, possibly
after pinning suitably chosen nodes, of some examples from the biochemical
literature, including transcription and translation models (Section V).

##  II Problem Formulation

Consider the general class of biochemical systems

$\dot{x}(t)=Sg(x(t))+g_{0},$ (1)

where the state vector $x(t)\in\mathbb{R}^{n}_{+}$ includes the concentrations
of the involved biochemical species, $S\in\mathbb{Z}^{n\times m}$ is the
stoichiometric matrix, the vector function
$g:\mathbb{R}^{n}_{+}\to\mathbb{R}^{m}$ represents reaction rates, and
$g_{0}\in\mathbb{R}^{n}_{+}$ is a vector of constant influxes. We make the
following standing assumptions.

###### Assumption 1

System (1) admits the equilibrium $\bar{x}\in\mathbb{R}^{n}_{+}$, such that

$0=Sg(\bar{x})+g_{0}.$ (2)

###### Assumption 2

The network is _unitary_ [8], namely, each of the entries of matrix $S$ is
either $1$, $0$ or $-1$.

###### Assumption 3

Function $g_{k}(x)$ is nonnegative and strictly monotonic (either increasing
or decreasing) in each of its arguments; it depends on variable $x_{i}$ if and
only if $S_{ik}=-1$, and it is zero if and only if one of its arguments is
zero.

###### Remark 1

Unitary networks, as in Assumption 2, cannot include multi-molecular
reactions, such as $2X_{1}+X_{2}\xrightharpoonup{\hbox{\scriptsize$$}}X_{3}$.
However, multi-molecular reactions are known to occur in fact as chains of bi-
molecular reactions (e.g.,
$X_{1}+X_{2}\xrightharpoonup{\hbox{\scriptsize$$}}X_{4}$ and
$X_{1}+X_{4}\xrightharpoonup{\hbox{\scriptsize$$}}X_{3}$), which do lead to a
unitary network. Note also that Assumption 3 rules out autocatalytic reactions
(e.g. of the form $X_{1}\xrightharpoonup{\hbox{\scriptsize$$}}2X_{1}$), for
which structural stability could never be guaranteed.

As shown in [8, 10, 12], the variable shift $z(t)=x(t)-\bar{x}$ allows us to
write the system in the equivalent form

$\dot{z}(t)=BD(z(t))Cz(t),$ (3)

where $D$ is a diagonal matrix with positive diagonal entries,

$D(z)=\mbox{diag}\\{D_{1}(z),D_{2}(z),\dots,D_{q}(z)\\},~{}~{}~{}~{}D_{k}(z)>0.$

Matrices $B$, $D$ and $C$ are derived from (1) as follows

  * •

The diagonal entries of $D$ are related to the absolute values of the $q$
nonzero partial derivatives $\partial g_{j}/\partial x_{i}$, with
$j\in\\{1,\dots,m\\}$, $i\in\\{1,\dots,n\\}$, arbitrarily ordered.

  * •

The $k$th column of $B$, $B_{k}$, corresponding to $D_{kk}$ related to
$|\partial g_{j}/\partial x_{i}|$, is equal to the column $S_{j}$ of $S$.

  * •

The $k$th row of $C$, $C_{k}^{\top}$, corresponding to $D_{kk}$ related to
$|\partial g_{j}/\partial x_{i}|$, has a single nonzero entry, the $i$th,
equal to the sign of $\partial g_{j}/\partial x_{i}$.

The proof [10, 12] relies on the fact that the Jacobian of $Sg(x)$ can be
written as $J(x)=B\Delta C$, where the diagonal matrix $\Delta$ includes the
absolute value of all partial derivatives, and on the integral formula

$\displaystyle
Sg(x)-Sg(\bar{x})=\left[\int_{0}^{1}~{}J(\bar{x}+\sigma(x-\bar{x}))d\sigma\right](x-\bar{x})$
$\displaystyle=$ $\displaystyle
B\left[\int_{0}^{1}~{}\Delta(\bar{x}+\sigma(x-\bar{x}))d\sigma\right]C(x-\bar{x})=BD(x)C(x-\bar{x}).$

_Structural_ stability (which needs to hold regardless of the numerical values
and functional expressions within matrix $D$, [8, 10, 12]) can be studied by
absorbing system (3) in a differential inclusion

$\dot{z}(t)=BD(t)Cz(t),$ (4)

where $D(t)$ is a positive definite diagonal matrix of size $q$.

(a)  (b)

Figure 1: Graph and GPN representation of the system in Example 1.

###### Example 1

Consider the biochemical reaction network
$\emptyset\xrightharpoonup{\hbox{\scriptsize$u_{1}$}}X_{1}$,
$X_{1}\xrightharpoonup{\hbox{\scriptsize$g_{1}$}}X_{2}+X_{3}$,
$X_{1}+X_{3}\xrightharpoonup{\hbox{\scriptsize$g_{13}$}}\emptyset$,
$X_{2}\xrightharpoonup{\hbox{\scriptsize$g_{2}$}}\emptyset$, admitting the
graph representation in Fig. 1, left. The system evolution is described by the
differential equations

$\begin{cases}\dot{x}_{1}=u_{1}-g_{1}(x_{1})-g_{13}(x_{1},x_{3})\\\
\dot{x}_{2}=g_{1}(x_{1})-g_{2}(x_{2})\\\
\dot{x}_{3}=g_{1}(x_{1})-g_{13}(x_{1},x_{3})\end{cases}$

which can be recast in the form (3) with

$B=\left[\begin{smallmatrix}-1&~{}~{}0&-1&-1\\\ ~{}~{}1&-1&~{}~{}0&~{}~{}0\\\
~{}~{}1&~{}~{}0&-1&-1\end{smallmatrix}\right]\quad\mbox{and}\quad
C=\left[\begin{smallmatrix}1&0&0\\\ 0&1&0\\\ 0&0&1\\\
1&0&0\end{smallmatrix}\right].$

###### Remark 2

For unitary networks, satisfying Assumption 2, $C_{i}^{\top}B_{i}=-1$ for all
$i=1,\dots,q$.

Henceforth we work under an additional, mild assumption.

###### Assumption 4

There exist unknown bounds $D^{-}_{k}>0$ and $D^{+}_{k}>0$ such that

$D^{-}_{k}\leq D_{k}\leq D^{+}_{k},\quad k=1,\dots,q.$ (5)

###### Definition 1

The system (3) under Assumptions 1-4 is _structurally stable_ if the
differential inclusion (4) is Lyapunov stable (possibly marginally), while it
is _structurally asymptotically stable_ if (4) is asymptotically stable.

###### Remark 3

The purely technical Assumption 4 does not change the structural nature of our
investigation. The upper bound $D^{+}_{k}$ ensures compactness and can be
arbitrarily large. The arbitrarily small $D^{-}_{k}$ ensures $D_{k}$ to be
bounded away from zero (cf. the $\epsilon$-perturbation in [4, Definition 2]).
For instance, the scalar system $\dot{z}(t)=-D(t)z(t)$ is not necessarily
asymptotically stable with the weaker bound $D>0$ (e.g., if $D(t)=e^{-t}$,
$\dot{z}=-D(t)z(t)$ does not converge to $0$), while stability is asymptotic
if $D(t)\geq D^{-}>0$.

##  III Polyhedral Lyapunov Functions and Generalised Petri Net Boundedness

###### Definition 2

Given an uncertain dynamical system

$\dot{x}=f(x,w),~{}~{}~{}w\in{\cal W},$

where ${\cal W}$ is a closed set, the positively homogeneous convex function
$V(x)$ is a _Lyapunov function (LF)_ for the system $\dot{x}=f(x,w)$ if, for
some $\beta\geq 0$, the (generalised) Lyapunov derivative

$D^{+}V(x,w)=\limsup_{h\rightarrow 0}\frac{V(x+hf(x,w))-V(x))}{h}\leq-\beta
V(x)$

for all $x$ and $w\in{\cal W}$. The LF is _weak_ if the inequality holds for
$\beta=0$, _strong_ if $\beta>0$. The LF $V(x)$ is _polyhedral (PLF)_ if it
can be written as

$V(x)=\|Fx\|_{\infty},$ (6)

where matrix $F$ has full column rank, or

$V(x)=\min\\{\|p\|_{1}:~{}x=Xp\\},$ (7)

where matrix $X$ has full row rank. _Polyhedral stability_ means that the
system admits a PLF.

The procedure in [8] to generate a PLF for the system associates the original
differential inclusion (4) with a discrete difference inclusion

$y(k+1)=\Phi(k)y(k),~{}~{}\Phi(k)\in\mathcal{F},$ (8)

where

$\mathcal{F}=\\{\Phi_{i}:\Phi_{i}=I+B_{i}C_{i}^{\top},\,\,i=1,\dots,q\\}.$ (9)

The procedure iterates over polyhedral sets, starting from the unit ball of
the $1$-norm: $\mathcal{Y}^{0}=\mbox{conv}\\{[-I~{}I]\\}$, where conv denotes
the convex hull.

###### Procedure 1

[8]

  1. 1.

$Y^{0}=I$, $\mathcal{Y}^{0}=\mbox{conv}\\{[-Y^{0}~{}Y^{0}]\\}$;

  2. 2.

$Y^{k+1}:=[\Phi_{1}~{}\Phi_{2}~{}\dots~{}\Phi_{q}]Y^{k}$;

  3. 3.

$\mathcal{Y}^{k+1}:=\mbox{conv}\\{[-Y^{k+1}~{}Y^{k+1}]\\}$,

  4. 4.

IF $\mathcal{Y}^{k+1}=\mathcal{Y}^{k}$, set
$\bar{\mathcal{Y}}=\mathcal{Y}^{k}$ and STOP; ELSE go to step $2$.

For a practical implementation, further stopping conditions should be added
before the ELSE statement at Step 4 of the procedure; otherwise, as currently
stated, the procedure might fail to stop. Convergence issues are one of the
aspects we will investigate: stopping criteria with a negative outcome (i.e.,
no structural PLF exists) will be discussed in Section III-C for numerical
purposes.

If Procedure 1 stops, the polytope $\bar{\mathcal{Y}}=\mathcal{Y}^{k}$, with
vertices $X=[-Y^{k}~{}Y^{k}]$, is the unit ball of a PLF as in (7). If we
apply the same procedure to the dual system $\dot{z}=C^{\top}DB^{\top}z$, by
considering $\Phi_{i}^{\top}$, under convergence assumptions, we obtain the
PLF in the dual form in (6) with $F=Y^{\top}$, where $Y=[-Y^{k}~{}Y^{k}]$. The
efficient implementation of the procedure requires removing the redundant
columns at each iteration $k$; see [8] for details. The sequence
$\mathcal{Y}^{k}$ being formed by integer vectors drastically improves
computability and provides efficient stopping criteria when the procedure
fails to converge.

As proven in [8], the stability of the differential inclusion (4) is
equivalent to the stability of (8); if (4) admits a (weak) PLF, then it is
marginally stable, which implies that system (3), with $D(z)$ continuous, is
stable, and is asymptotically stable if and only if its Jacobian $BDC$ is
structurally non-singular [10]. Moreover, if $C_{i}^{\top}B_{i}=-1$, the
stability of (4) is equivalent to the existence of a PLF for both (4) and (8)
and also equivalent to the fact that Procedure 1 successfully stops in finite
time. If $C_{i}^{\top}B_{i}=-1$, systems (3) and (4) admit a structural (weak)
Lyapunov function if and only if they admit a (weak) structural PLF. Hence,
the existence of a PLF guarantees stability, and even asymptotic stability
under structural non-singularity assumptions, as summarised in the following
result.

###### Theorem 1

[10] Assume that system (4), under Assumptions 1-4, admits a (weak) PLF. Then,
it is asymptotically stable if and only if matrix $BDC$ is non-singular for
all possible matrices $D$ satisfying (5).

Asymptotic stability is shown to be exponential in [22].

###### Remark 4

Structural non-singularity is easy to check, as shown in [26]: it is
equivalent to $\det[-B\hat{D}C]>0$ for all matrices $\hat{D}$ on the vertices
of the hyper-rectangle in (5).

###  III-A Generalised Petri net model

Procedure 1 can be interpreted as the evolution of a Generalised Petri Net
(GPN), a discrete-event system that, albeit similar to a Petri net, does not
work as a Petri net because the number of tokens at each node can be either
positive or negative [32]. The integer vector $\tau(k)$ denotes the presence
of $\tau_{\ell}(k)$ tokens in position $\ell$ at time $k$.

###### Example 2

The system in Example 1 can be associated with the difference inclusion
$y(k+1)=\Phi(k)y(k)$ where, at each step $k$, $\Phi(k)$ is one of the
following matrices:

$\Phi_{1}=\left[\begin{smallmatrix}0&0&0\\\ 1&1&0\\\
1&0&1\end{smallmatrix}\right],\,\Phi_{2}=\left[\begin{smallmatrix}1&0&0\\\
0&0&0\\\
0&0&1\end{smallmatrix}\right],\,\Phi_{3}=\left[\begin{smallmatrix}1&0&-1\\\
0&1&~{}~{}0\\\
0&0&~{}~{}0\end{smallmatrix}\right],\,\Phi_{4}=\left[\begin{smallmatrix}~{}~{}0&0&0\\\
~{}~{}0&1&0\\\ -1&0&1\end{smallmatrix}\right].$

We start from the unit ball of the $1$-norm,
$\mathbb{X}^{0}=[-I~{}~{}I]=[-v_{1}~{}-v_{2}~{}-v_{3}~{}~{}v_{1}~{}~{}v_{2}~{}~{}v_{3}]$,
and consider just the positive vertices (the evolution of the others can be
immediately obtained, being the opposite). Vertex $v_{1}=[1~{}0~{}0]^{\top}$
is transformed into $\Phi_{1}v_{1}=[0~{}1~{}1]^{\top}=v_{4}$,
$\Phi_{2}v_{1}=\Phi_{3}v_{1}=[1~{}0~{}0]^{\top}=v_{1}$,
$\Phi_{4}v_{1}=[0~{}0~{}-1]^{\top}=-v_{3}$; vertex $v_{2}=[0~{}1~{}0]^{\top}$
is transformed into
$\Phi_{1}v_{2}=\Phi_{3}v_{2}=\Phi_{4}v_{2}=[0~{}1~{}0]^{\top}=v_{2}$,
$\Phi_{2}v_{2}=[0~{}0~{}0]^{\top}$; vertex $v_{3}=[0~{}0~{}1]^{\top}$ is
transformed into
$\Phi_{1}v_{3}=\Phi_{2}v_{3}=\Phi_{4}v_{3}=[0~{}0~{}1]^{\top}=v_{3}$,
$\Phi_{3}v_{3}=[-1~{}0~{}0]^{\top}=-v_{1}$. The sole newly generated vertex is
$v_{4}=[0~{}1~{}1]^{\top}$ (and its opposite). The procedure applied to
$v_{4}$ gives $\Phi_{1}v_{4}=\Phi_{4}v_{4}=[0~{}1~{}1]^{\top}=v_{4}$,
$\Phi_{2}v_{4}=[0~{}0~{}1]^{\top}=v_{3}$,
$\Phi_{3}v_{1}=[-1~{}1~{}0]^{\top}=v_{5}$. Applying the procedure to the only
new vertex, $v_{5}=[-1~{}1~{}0]^{\top}$, gives
$\Phi_{1}v_{5}=[0~{}0~{}-1]^{\top}=-v_{3}$,
$\Phi_{2}v_{5}=[-1~{}0~{}0]^{\top}=-v_{1}$,
$\Phi_{3}v_{5}=[-1~{}1~{}0]^{\top}=v_{5}$,
$\Phi_{4}v_{5}=[0~{}1~{}1]^{\top}=v_{4}$. No new vertices are generated at
this step, hence the procedure stops successfully: the system admits a PLF
having unit ball $\mbox{conv}\\{[-X~{}X]\\}$, with
$X=[v_{1}~{}~{}v_{2}~{}~{}v_{3}~{}~{}v_{4}~{}~{}v_{5}]$.

(a) First step, initial marking $v_{1}=[1~{}0~{}0]^{\top}$.

(b) First step, initial marking $v_{2}=[0~{}1~{}0]^{\top}$.

(c) First step, initial marking $v_{3}=[0~{}0~{}1]^{\top}$.

(d) Second step, initial marking $v_{4}=[0~{}1~{}1]^{\top}$.

(e) Third step, initial marking $v_{5}=[-1~{}1~{}0]^{\top}$.

Figure 2: Discrete-event interpretation of the procedure evolution. On the
left side, enabled transitions are in green and non-enabled transitions in
red. On the right, the activated transition is in green.

The evolution of the discrete-time system in the numerical Procedure 1 can be
related to the evolution of a particular discrete-event system, in which the
initial conditions represent an _initial marking_ $\tau(0)$, assigning an
integer number of tokens to each node (associated with each of the chemical
species), and a transition (associated with each of the black rectangles in
Fig. 1, right, i.e. with each of the reactions occurring in the network) is
_enabled_ whenever at least one of the starting nodes of the transition
contains a non-zero number of tokens. When either the number of tokens in the
starting node is positive, or the starting nodes are two and the number of
tokens is non-negative in both (and positive in at least one of them), then
the transition takes tokens from the starting node(s) and moves them to the
arrival node(s), if explicitly present (otherwise, they simply disappear).
When either the starting node is one only and the number of tokens therein is
negative, or the starting nodes are two and the number of tokens is non-
positive in both (and negative in at least one of them), then the transition
takes negative tokens from the starting node(s) and moves them to the arrival
node(s), again, if explicitly present. When a transition is enabled, and
performed, all the tokens are moved until one of the starting nodes has zero
tokens. If one token number is negative and the other positive, their effects
are superposed.

###### Example 3

For the system in Examples 1 and 2, the GPN evolution is shown in Fig. 2: on
the left side, the initial marking is illustrated (enabled transitions are in
green, non-enabled transitions in red); on the right side, the new marking
generated by the action of each of the enabled transitions is illustrated.
Reaction I corresponds to the evolution matrix $\Phi_{1}$ (and generates the
same outcome for the same initial conditions), reaction II corresponds to
matrix $\Phi_{2}$ (and again generates the same outcome for the same initial
conditions), while reaction III corresponds to matrices $\Phi_{3}$ and
$\Phi_{4}$. In this latter case, the outcome of (at least) one of the two
evolution matrices is the same as that of the transition related to reaction
III in the discrete-event evolution, while the outcome of the other (if
different) is always the unchanged input vertex. Interestingly, the evolution
is the very same as that of the numerical procedure, and the same new vertices
are generated at each iteration.

###### Proposition 1

If $C_{i}^{\top}B_{i}=-1$ $\forall i=1,\dots,q$, the sequence generated by the
discrete-time switching system (8)-(9) uniquely corresponds to the evolution
of a GPN discrete-event system if the initial vector $y(0)=\tau(0)$ has
integer components.

###### Proof:

To keep the notation simple, we assume that all reactions are functions of at
most two variables. The integer operator $\Phi_{h}$ corresponds to the $h$th
non-zero derivative, say $\partial g_{i}/\partial x_{j}$. Consider the $i$th
reaction $X_{j}+X_{w}\xrightharpoonup{\hbox{\scriptsize$g_{i}$}}X_{\ell}$.
Then,

$\Phi_{h}=I+B_{h}C_{h}^{\top}=I+S_{i}e_{j}^{\top}$

where $S_{i}$ is the $i$th column of $S$ and $e_{j}^{\top}$ is the $j$th
canonical row vector. Hence, $\Phi_{h}$ has all columns equal to the identity
matrix excluding the $j$th. Column $j$ has the same entries as $B_{h}=S_{i}$,
excluding the diagonal entry $(j,j)$ which is $0$ (see Example 2). Operator
$\Phi_{h}$ applied to any integer vector $\tau$ as $\tau^{+}=\Phi_{h}\tau$, a
shorthand notation for $\tau(k+1)=\Phi_{h}\tau(k)$, corresponds to the
following transition $T_{h}$ performed on the Petri net:
1) $\tau^{+}_{j}=0$ (remove $\tau_{j}$, i.e. all, tokens from node $j$);
2) $\tau^{+}_{\ell}=\tau_{\ell}+B_{\ell h}\tau_{j}=\tau_{\ell}+\tau_{j}$ (add
$\tau_{j}$ tokens to node $\ell$);
3) $\tau^{+}_{w}=\tau_{w}+B_{wh}\tau_{j}=\tau_{w}-\tau_{j}$ (remove $\tau_{j}$
tokens from $w$).

This defines a one-to-one correspondence between matrix $\Phi_{h}$ applied on
integer vectors and the corresponding transition $T_{h}$ (such that
$\tau^{+}=T_{h}\tau$). ∎

Henceforth, we denote by $\Phi(h)$, $h=0,1,\dots$, the generic matrix sequence
$\Phi(h)\in\\{\Phi_{1},\Phi_{2},\dots,\Phi_{q}\\}=\mathcal{F}$. The family of
all their products,

$\Pi(\mathcal{F})\doteq\left\\{\prod_{h=0}^{K}\Phi(h),~{}K\geq
0,~{}~{}\Phi(h)\in\mathcal{F}\right\\},$ (10)

is an algebraic semigroup under the multiplication operation, i.e., if both
$P_{1}$ and $P_{2}$ are in $\Pi(\mathcal{F})$, then also
$P_{1}P_{2}\in\Pi(\mathcal{F})$.

###### Theorem 2

If $C_{i}^{\top}B_{i}=-1$ $\forall i=1,\dots,q$, the following statements are
equivalent:
(i) for any initial integer token distribution vector $\tau(0)=\tau_{0}$, the
set of possible evolutions of the GPN $\tau(k+1)=T_{k}\tau(k)$,

${\mathcal{R}}(\tau_{0})=\\{\tau=T_{\ell}\circ T_{\ell-1}\dots\circ
T_{0}~{}\tau_{0},\ell\geq 0,~{}T_{i}~{}\mbox{arbitrary}\\},$

i.e. the reachable set from $\tau_{0}$, is finite.
(ii) The semigroup $\Pi(\mathcal{F})$ in (10) is finite.
(iii) Procedure 1 stops in finite time.

###### Proof:

(i) $\Rightarrow$ (ii): If the GPN generates a finite number of configurations
given any initial token distribution $\tau(0)$, then from Proposition 1, for
any integer $y(0)$, the number of states reached by $y(k+1)=\Phi(k)y(k)$, is
finite, hence bounded. Since the generic trajectory given $y(0)$ has the form

$y(k)=\Phi(k)\Phi(k-1)\dots\Phi(1)\Phi(0)y(0),$

this implies that the products in the set $\Pi(\mathcal{F})$ are uniformly
bounded by some constant $\mu>0$: $\|\prod_{h=0}^{K}\Phi(h)\|\leq\mu$. Then
(ii) follows because all these products are integer matrices.
(ii) $\Rightarrow$ (iii): Procedure 1 does stop in finite time, because all
the columns of the matrices $Y^{k}$ are generated as
$\prod_{h=0}^{K}\Phi(h)Y^{0}$, which are in a finite number.
(iii) $\Rightarrow$ (i): If Procedure 1 stops in finite time, the matrix
family $\mathcal{F}$ admits the PLF induced by the final set
$\bar{\mathcal{Y}}$. Hence, for any integer initial condition $y(0)$ the
sequence $y(k)$ is bounded and, being integer, it is finite. Statement (i)
then follows from Proposition 1. ∎

###  III-B Interpretation of the results

Theorem 2 has an interesting biochemical interpretation, where tokens can be
seen as molecules. The transition operators $T_{k}$ defined in the proof of
Proposition 1 remove the tokens (if any) from some nodes (associated with
species) so as to generate tokens at other nodes. The existence of a PLF is
equivalent to the fact that no (infinite) sequence of these transitions can
drive the token count to infinity at some node.

Consider the important special case of mono-molecular reaction networks, where
all the internal reactions have the form
$X_{i}\xrightharpoonup{\hbox{\scriptsize$g_{k}$}}X_{j}$. Any transition just
moves all tokens in a node (possibly a negative number) to another node
leaving unchanged the total amount. For instance, if we initialise the network
with just a token at node $1$, the set of all possibly reached states
corresponds to a single token at some node. All mono-molecular reaction
networks are associated with bounded GPNs: indeed they are nonlinear
compartmental systems, well known to be structurally stable [30].

Besides reactions of the form
$X_{i}\xrightharpoonup{\hbox{\scriptsize$g_{k}$}}X_{j}$, let us consider
internal reactions of the form
$X_{h}+X_{r}\xrightharpoonup{\hbox{\scriptsize$g_{s}$}}\emptyset$ and
$X_{w}\xrightharpoonup{\hbox{\scriptsize$g_{u}$}}\emptyset$. The former new
reaction introduces operators such that, if $x_{h}$ tokens are present at node
$X_{h}$, they are removed ($x_{h}^{+}=0$) and the opposite amount appears at
node $X_{r}$, $x_{r}^{+}=x_{r}-x_{h}$. The latter new reaction removes all the
tokens present at node $X_{w}$ ($x_{w}^{+}=0$). Although the total amount of
tokens now is not conserved, it cannot increase: hence, these GPNs are also
bounded and the corresponding networks are structurally stable, in agreement
with [8].

As a simple unbounded case, consider the reactions
$\emptyset\xrightharpoonup{\hbox{\scriptsize$u_{1}$}}X_{1}$,
$X_{1}\xrightharpoonup{\hbox{\scriptsize$g_{1}$}}X_{2}+X_{3}$,
$X_{2}\xrightharpoonup{\hbox{\scriptsize$g_{2}$}}X_{1}$,
$X_{3}\xrightharpoonup{\hbox{\scriptsize$g_{3}$}}X_{1}$,
$X_{2}\xrightharpoonup{\hbox{\scriptsize$\tilde{g}_{2}$}}\emptyset$,
$X_{3}\xrightharpoonup{\hbox{\scriptsize$\tilde{g}_{3}$}}\emptyset$,
associated with

$\begin{cases}\dot{x}_{1}=u_{1}-g_{1}(x_{1})+g_{2}(x_{2})+g_{3}(x_{3})\\\
\dot{x}_{2}=g_{1}(x_{1})-g_{2}(x_{2})-\tilde{g}_{2}(x_{2})\\\
\dot{x}_{3}=g_{1}(x_{1})-g_{3}(x_{3})-\tilde{g}_{3}(x_{3})\end{cases}$

The corresponding Petri net is not bounded. Start with just one token at
$X_{1}$. Then $\partial g_{1}/\partial x_{1}$ can act producing two tokens,
one at $X_{2}$ and one at $X_{3}$. Then $g_{2}$ and $g_{3}$ can both act to
transfer the two tokens back at $X_{1}$. Repeating the argument, we see an
unbounded increase of tokens at node $1$. Indeed, the Jacobian $J=BDC$, where

$B=\left[\begin{smallmatrix}-1&~{}1&~{}1&~{}0~{}&~{}0\\\
~{}1&-1&~{}0&-1~{}&~{}0\\\
~{}1&~{}0&-1&~{}0~{}&-1\end{smallmatrix}\right]~{}\mbox{and}~{}C=\left[\begin{smallmatrix}1&0&0\\\
0&1&0\\\ 0&0&1\\\ 0&1&0\\\ 0&0&1\end{smallmatrix}\right],$

is not structurally Hurwitz: $\det[-J]=D_{1}(D_{4}D_{5}-D_{2}D_{3})$, the
constant term of the characteristic polynomial, can be negative.

###### Remark 5

(Stoichiometric compatibility class.) If the network evolves in a proper
stoichiometric compatibility class, namely $\sigma^{\top}S=0$ for some vector
$\sigma\neq 0$, so that $\sigma^{\top}B=0$, we have a conservation law:
$\sigma^{\top}z(t)$ is constant. This property is preserved by the discrete
operators $\Phi_{h}$:
$\sigma^{\top}y^{+}=\sigma^{\top}[I+B_{h}C^{\top}_{h}]y=\sigma^{\top}y$. For
our analysis, we can reduce the system by applying a state transformation that
turns $BDC$ into $T^{-1}BDCT\doteq\hat{B}D\hat{C}$ and then neglecting some of
the variables (see Example 7); the key condition
$C_{h}^{\top}B_{h}=\hat{C}_{h}^{\top}\hat{B}_{h}=-1$ is invariant. Therefore,
if $T$ and its inverse are integer matrices, as usually happens, the proposed
theory applies without changes.

Interestingly, the convergence of Procedure 1 (i.e., the fact that the
procedure stops in finite time) implies that the joint spectral radius [27] of
the matrix family $\mathcal{F}$ is equal to one:

$\sigma(\mathcal{F})\doteq\lim_{k\to
0}\max_{\Phi(\cdot)\in\mathcal{F}}\|\Phi(1)\Phi(2)\Phi(3)\dots\Phi(k)\|^{\frac{1}{k}}=1.$

Indeed, the convergence of Procedure 1 implies the boundedness of the
trajectories, hence $\sigma(\mathcal{F})\leq 1$. Conversely, since
$C_{k}^{\top}B_{k}=-1$, matrix $I+B_{k}C_{k}$ admits $1$ as an eigenvalue,
hence the spectral radius cannot be smaller than $1$.

Boundedness of the GPN evolution is equivalent to polyhedral stability, and
the GPN evolves according to an asynchronous mechanism: hence, the different
time scales of the system components play no role in defining its stability
properties. To formalise this concept, we modify system (1) as

$\Theta\dot{x}=Sg(x)+g_{0},$ (11)

where $\Theta=\mbox{diag}\\{\theta_{1},\theta_{2},\dots,\theta_{n}\\}\succ 0$
is a diagonal matrix of positive time constants.

###### Proposition 2

Structural polyhedral stability, guaranteed when Procedure 1 stops in finite
time, implies the structural stability of system (11) for any arbitrary
diagonal $\Theta\succ 0$.

###### Proof:

Scaling the state variable as $y=\Theta x$ turns equation (11) into
$\dot{y}=Sg(\Theta^{-1}y)+g_{0}$. The proof follows immediately by noticing
that matrices $B$ and $C$ are the same regardless of $\Theta$, while the
derivatives are scaled as $\partial g_{i}/\partial
y_{j}=\Theta_{j}^{-1}\partial g_{i}/\partial x_{j}$, which does not alter
their sign. ∎

###  III-C Stopping criteria for Procedure 1

The procedure may fail to converge; in this case, the system does not admit
any structural PLF. A possible stopping criterion, proposed in [8], is to
interrupt the procedure when either the size of the region $\mathcal{Y}^{k}$
reaches a bound $\mu$, or an assigned maximum number of steps, $n_{max}$, is
reached.

We discuss here other possible criteria. The procedure will never converge if,
for some $k$, $\mathcal{Y}^{k}$ includes the original region $\mathcal{Y}^{0}$
in its interior [13, 14]. The inclusion of the initial polytope, with vertex
matrix $[-I~{}I]$, should be checked at each step; this can be done as
follows.

###### Proposition 3

If the polytope $\mathcal{Y}=\mbox{conv}\\{[-Y,Y]\\}$ has a non-empty
interior, then it includes $\mathcal{Y}^{0}=\mbox{conv}\\{[-I,I]\\}$ in its
interior if and only if

$\nu_{i}=\min\\{\|p_{i}\|_{1}:Yp_{i}=e_{i}\\}<1,~{}~{}~{}i=1,2,\dots,n,$ (12)

where $\|p\|_{1}=\sum|p_{i}|$ is the $1$-norm and $e_{i}$ is the $i$th vector
of the canonical basis.

###### Proof:

From expression (7) we see that $\nu_{i}=V_{Y}(e_{i})$, where $V_{Y}$ is the
polyhedral norm with unit ball $\mathcal{Y}$. Vector $e_{i}$ (along with its
opposite $-e_{i}$) is in the interior of $\mathcal{Y}$ iff
$\nu_{i}=V_{Y}(e_{i})<1$. Moreover, $\mathcal{Y}^{0}=\mbox{conv}\\{[-I,I]\\}$
is in the interior of $\mathcal{Y}$ iff it vertices, i.e. $\pm e_{i}$, are in
the interior. ∎

Note that problem (12) can be solved via linear programming and allows to
efficiently stop the procedure at an early stage when the system does not
admit a structural PLF.

Another stopping criterion relies on eventually periodic matrices. A square
matrix $M$ is said _eventually periodic_ if there exist a non-negative integer
$m$ and a positive integer $p$ such that

$M^{m}=M^{m+kp}\qquad~{}~{}\mbox{for all integer}~{}k\geq 0.$ (13)

###### Proposition 4

If the square matrix $M$ is eventually periodic, then its eigenvalues are
either zero or roots of the unity.

###### Proof:

Take the eigenpair $(\lambda,v)$, $Mv=\lambda v$, $v\neq 0$. If $M$ is
eventually periodic,

$(M^{m}-M^{m+kp})v=\lambda^{m}(1-\lambda^{kp})v=0.$

Then $\lambda^{m}(1-\lambda^{kp})=0$, i.e., $\lambda$ must be either $0$ or a
root of the unity. ∎

###### Proposition 5

The set $\Pi(\mathcal{F})$ is a finite set only if each matrix in
$\Pi(\mathcal{F})$ is eventually periodic.

###### Proof:

By contradiction, if a matrix in $\Pi(\mathcal{F})$ is not eventually
periodic, then its powers form an infinite sequence of different matrices,
hence the set $\Pi(\mathcal{F})$ is infinite. ∎

The previous condition is not sufficient: even if all matrices are eventually
periodic, there can be an infinite sequence of products among them. Combining
Propositions 4 and 5 yields the following corollary.

###### Corollary 1

If $\Pi(\mathcal{F})$ includes a matrix whose eigenvalues are not either zero
or roots of the unity, then it has infinite cardinality.

###### Remark 6

To check the condition, there is no need to compute _all_ the products in
$\Pi(\mathcal{F})$. In fact,

  * •

if $(I+B_{i}C_{i}^{\top})\in\mathcal{F}$, then
$(I+B_{i}C_{i}^{\top})^{n}=(I+B_{i}C_{i}^{\top})$ for any positive integer
$n$. Indeed, since $C_{i}^{\top}B_{i}=-1$,

$(I+B_{i}C_{i}^{\top})^{2}=I+2B_{i}C_{i}^{\top}+B_{i}(C_{i}^{\top}B_{i})C_{i}^{\top}=I+B_{i}C_{i}^{\top};$

  * •

if $\Phi_{1}=I+B_{1}C_{1}^{\top}$ and $\Phi_{2}=I+B_{2}C_{2}^{\top}$ with
$C_{1}=C_{2}$, then $\Phi_{1}\Phi_{2}=\Phi_{2}$, as it can be seen in a
similar way.

Based on Corollary 1, an alternative stopping criterion is achieved by
computing the sequence of products $\Pi(\mathcal{F})$ of increasing order and
stopping whenever one of them has an eigenvalue that is neither zero nor a
root of the unity. This produces, in principle, an exponentially growing list
of matrices. Yet, extensive numerical experiments have shown that, in most
cases, the stopping condition is quickly reached.

###  III-D Turning nodes into black holes

We introduce a new type of node, called _black hole_ , in which any incoming
token (either positive or negative) is cleared, so that the black hole
contains zero tokens throughout the system evolution. If a node is replaced by
a black hole, then the GPN associated with the reaction network is transformed
and behaves differently.

###### Example 4

Consider the reaction network in Fig. 3, where
$\emptyset\xrightharpoonup{\hbox{\scriptsize$$}}X_{1}$,
$\emptyset\xrightharpoonup{\hbox{\scriptsize$$}}X_{2}$,
$X_{1}+X_{2}\xrightharpoonup{\hbox{\scriptsize$$}}X_{3}+X_{4}$,
$X_{4}\xrightharpoonup{\hbox{\scriptsize$$}}X_{2}$,
$X_{1}+X_{3}\xrightharpoonup{\hbox{\scriptsize$$}}\emptyset$. The GPN
associated with this network is unbounded, i.e., no structural PLF exists for
the system. If we turn $X_{2}$ into a black hole, then we virtually have
$X_{2}:=\emptyset$ and the transformed reaction network becomes:
$\emptyset\xrightharpoonup{\hbox{\scriptsize$$}}X_{1}$,
$X_{1}\xrightharpoonup{\hbox{\scriptsize$$}}X_{3}+X_{4}$,
$X_{4}\xrightharpoonup{\hbox{\scriptsize$$}}\emptyset$,
$X_{1}+X_{3}\xrightharpoonup{\hbox{\scriptsize$$}}\emptyset$. The GPN
associated with the transformed network is bounded. Boundedness is achieved
also if node $X_{4}$ is turned into a black hole, instead of $X_{2}$.
Conversely, turning either node $X_{1}$ or node $X_{3}$ into a black hole does
not yield boundedness.

Figure 3: Graph of the chemical reaction network corresponding to the system
in Example 4.

As we show in the next section, replacing a node with a black hole can be
regarded as applying a strong feedback to that node, thus enforcing a pinning
control.

##  IV Pinning Control

Pinning control strategies control just some of the state variables, each by
means of a (strong) local feedback, so as to stabilise the whole network.
Without loss of generality, we split the state vector as
$z=[z_{1}^{\top}\,\,z_{2}^{\top}]^{\top}$ and we assume that a feedback
control action with $\gamma>0$ is applied to $z_{1}$, the first $p$ entries of
$z$:

$\dot{z}_{1}=S_{1}(g(x)-g(\bar{x}))-\gamma z_{1},$

where $S_{1}$ contains the first $p$ rows of $S$. The system (3) can then be
split as

$\begin{bmatrix}\dot{z}_{1}\\\
\dot{z}_{2}\end{bmatrix}=\underbrace{\begin{bmatrix}B_{1}D(z)C_{1}-\gamma
I&B_{1}D(z)C_{2}\\\
B_{2}D(z)C_{1}&B_{2}D(z)C_{2}\end{bmatrix}}_{A(D(z),\gamma)}\begin{bmatrix}z_{1}\\\
z_{2}\end{bmatrix},$ (14)

where $B_{1}$ and $B_{2}$ contain the first $p$ and the last $n-p$ rows of
$B$, respectively, while $C_{1}$ and $C_{2}$ contain the first $p$ and the
last $n-p$ columns of $C$.

We have a first main result.

###### Theorem 3

Under Assumptions 1-4, the following conditions are equivalent:

  * (i)

The $(n-p)$-dimensional subsystem

$\dot{z}_{2}(t)=B_{2}D(t)C_{2}z_{2}(t)$ (15)

admits a structural PLF $U(z_{2})$, in the strong sense of Definition 2.

  * (ii)

There exists a polyhedral function $V(z)$ and $\bar{\gamma}>0$ such that
$V(z)$ is a structural strong LF for system (14) if $\gamma\geq\bar{\gamma}$ .

###### Remark 7

The structural result needs Assumption 4. Consider the linear differential
inclusion with matrix

$\begin{bmatrix}-\gamma+D_{1}(t)&1\\\ 1&-D_{2}(t)\end{bmatrix}.$

No matter how large $\gamma>0$ is, the system becomes unstable if $D_{1}(t)>0$
diverges and/or $D_{2}(t)>0$ converges to $0$; conversely, under Assumption 4,
a stabilising $\gamma$ always exists, which depends on the bounds in (5).

Proof of Theorem 3

(i) $\Rightarrow$ (ii). Due to the compactness Assumption 4, the matrix family
in system (14) is polytopic and we can write

$A(D,\gamma)=\begin{bmatrix}B_{1}DC_{1}&B_{1}DC_{2}\\\
B_{2}DC_{1}&B_{2}DC_{2}\end{bmatrix}=\sum_{k=1}^{M}\lambda_{k}\begin{bmatrix}A^{(k)}_{11}&A^{(k)}_{12}\\\
A^{(k)}_{21}&A^{(k)}_{22}\end{bmatrix},$

with $\sum_{k=1}^{M}\lambda_{k}=1$, $\lambda_{k}>0$, because each of the four
matrix blocks in $A(D,\gamma)$ can be expressed as a polytopic matrix
$\sum_{k=1}^{M}\lambda_{k}A^{(k)}_{ij}$ [14].

By assumption, system (15) admits a structural PLF, hence
$[A^{(k)}_{22}]X_{2}=X_{2}P_{22}^{(k)}$, where $P_{22}^{(k)}$ is strictly
column diagonally dominant and $X_{2}$ has full row rank [14, 15, 31]. Then we
can take

$\hat{X}=\begin{bmatrix}\rho I&0\\\ 0&X_{2}\end{bmatrix},$

where $\rho>0$ is a parameter to be selected, and write the
$\gamma$-parametrised equation

$\displaystyle\underbrace{\left[\begin{smallmatrix}-\gamma
I+A^{(k)}_{11}&A^{(k)}_{12}\\\
A^{(k)}_{21}&A^{(k)}_{22}\end{smallmatrix}\right]}_{\hat{A}^{(k)}}\hat{X}=\hat{X}\underbrace{\left[\begin{smallmatrix}-\gamma
I+A^{(k)}_{11}&A_{12}^{(k)}X_{2}/\rho\\\
M_{21}^{(k)}\rho&P_{22}^{(k)}\end{smallmatrix}\right]}_{\hat{P}^{(k)}}$ (16)

where $M_{21}^{(k)}$ are matrices that satisfy
$X_{2}M_{21}^{(k)}=A_{21}^{(k)}$, and always exist because $X_{2}$ has full
row rank.

Since $P_{22}^{(k)}$ is strictly diagonally dominant, we can take $\rho>0$
large enough, so that the diagonal dominance is preserved in the last $n-p$
columns of $\hat{P}^{(k)}$. Now, for any choice of $\rho$, there exists a
$\bar{\gamma}$ large enough such that, for $\gamma\geq\bar{\gamma}$, diagonal
dominance of the first $p$ columns of $\hat{P}^{(k)}$ is ensured. The
resulting equations $\hat{A}^{(k)}\hat{X}=\hat{X}\hat{P}^{(k)}$, with
$P^{(k)}$ diagonally dominant, ensure that $[-\hat{X},\hat{X}]$ are the
vertices of a PLF for system (14).

(ii) $\Rightarrow$ (i). A technical Lemma is required.

###### Lemma 1

Assume that the convex and compact set ${\cal S}$ including the origin as an
interior point (C-set) is positively invariant for the linear time-invariant
system

$\begin{bmatrix}\dot{z}_{1}\\\
\dot{z}_{2}\end{bmatrix}=\begin{bmatrix}F_{11}-\mu I&F_{12}\\\
F_{21}&F_{22}\end{bmatrix}\begin{bmatrix}z_{1}\\\ z_{2}\end{bmatrix}$ (17)

for all $\mu>0$. Then the intersection ${\cal S}_{2}=\\{z\in{\cal S}\colon
z_{1}=0\\}$ is a C-set in the subspace $[z_{1}^{\top}~{}z_{2}^{\top}]^{\top}$
with $z_{1}=0$ and is positively invariant for the subsystem
$\dot{z}_{2}=F_{22}z_{2}$.

The proof of Lemma 1 is in the appendix.

By assumption, ${\cal S}=\\{z\colon V(z)\leq 1\\}$, the unit ball of $V(z)$,
is an invariant set for system (14) for all $\gamma\geq\bar{\gamma}$. Let us
perturb system (14) and write it as

$\dot{z}=[A(D,\gamma)+\epsilon I]z\doteq F_{\epsilon}(D,\gamma)z,$ (18)

with $\epsilon>0$ small enough to ensure that exponential stability is
preserved: if $D^{+}V(z)\leq-\beta V(z)$, just take $0<\epsilon<\beta$.

For any fixed $\bar{D}$, if $\gamma\geq\bar{\gamma}$, ${\cal S}$ is invariant
for the linear time invariant system with state matrix
$F_{\epsilon}(\bar{D},\gamma)$. Take $\mu=\gamma-\bar{\gamma}$ and apply Lemma
1. The intersection ${\cal S}_{2}=\\{z\in{\cal S}\colon z_{1}=0\\}$, which is
a polyhedral C-set in the $z_{2}$-space, is positively invariant for the
system $\dot{z}_{2}=[A_{22}(\bar{D},\gamma)+\epsilon I]z_{2}$. Since this
claim is true for any choice of $D$, the C-set ${\cal S}_{2}$ is robustly
positively invariant for the differential inclusion
$\dot{z}_{2}=[A_{22}(D(t),\gamma)+\epsilon I]z_{2}$, which is thus at least
marginally stable. As a consequence, $\dot{z}_{2}=A_{22}(D,\gamma)z_{2}$ is
exponentially stable, hence it admits a PLF [14, 15, 31]. ∎

Although the same local feedback parameter $\gamma$ is considered in Theorem 3
for all pinned nodes, different parameters $\gamma_{i}$ could be adopted for
the $p$ nodes, provided that, for all $i$, $\gamma_{i}\geq\bar{\gamma}$ is
large enough to ensure diagonal dominance of the first $p$ columns in the last
matrix in (16). Furthermore, the result easily extends to nonlinear feedback
strategies $k_{i}(z_{i})z_{i}$, provided that $k_{i}(z_{i})\geq\bar{\gamma}$
for all $z_{i}$.

We now need to face a technical issue. Procedure 1 can be adopted to find a
structural PLF for system (15), and all the results in Section III, including
the stopping criterion, remain valid. Unfortunately, the procedure provides a
_weak_ structural PLF $V(z_{2})$, and not a _strong_ one as required by
Theorem 3. To fix the problem we consider three facts.

  * •

In view of Theorem 1, if we find a weak PLF for the differential inclusion
(15), we can claim its robust asymptotic stability if (and only if) matrix
$BDC$ is structurally non-singular.

  * •

A classical result [14, 15, 31] ensures that, if (15) with compact bounds (5)
(introduced exactly for this technical reason) is asymptotically stable, then
it is also exponentially stable and it admits a strong PLF $U(z_{2})$.

  * •

Hence, if we find $V(z_{2})$ and $BDC$ is structurally non-singular, we know
that a strong PLF $U(z_{2})$ exists. This allows us to apply Theorem 3;
fortunately, we do not need to compute $U(z_{2})$.

###### Corollary 2

Assume that $V(z_{2})$ is a weak structural PLF for (15) and that
$\det[B_{2}DC_{2}]\neq 0$ for all $D$ in (5). Then, there exists
$\bar{\gamma}$ such that, for all $\gamma>\bar{\gamma}$, system (14) is
structurally exponentially stable.

Structural non-singularity of $BDC$ can be checked as discussed in Remark 4.

###### Remark 8

All the results presented in this section so far hold also for non-unitary
networks: we do not need to assume $C_{k}^{\top}B_{k}=-1$.

###  IV-A Lyapunov function for the free variables $z_{2}$

If we find a PLF for the $z_{2}$-subsystem of dimension $n-p$, and we have
structural non-singularity of $BDC$, then the stability of the overall system
is ensured for large enough $\gamma$, according to Corollary 2. How can we
exploit the GPN to this aim?

###### Theorem 4

If $C_{i}^{\top}B_{i}=-1$ $\forall i=1,\dots,q$, the following statements are
equivalent:

  * i)

The $z_{2}$-subsystem (15) admits a weak structural PLF.

  * ii)

The evolution of the GPN where nodes $1,\dots,p$ have been turned into black
holes is bounded for any integer initial marking.

###### Proof:

We show that a (weak) PLF exists for (15) if and only if a (weak) PLF exists
for the differential inclusion

$\begin{bmatrix}\dot{z}_{1}\\\ \dot{z}_{2}\end{bmatrix}=\begin{bmatrix}0&0\\\
0&B_{2}D(t)C_{2}\end{bmatrix}\begin{bmatrix}z_{1}\\\ z_{2}\end{bmatrix}.$ (19)

Indeed, let $B_{2}DC_{2}=\sum_{k=1}^{q}~{}D_{k}B_{2,k}C_{2,k}^{\top}$. A weak
PLF for (19) exists if and only if the equation [14]

$\begin{bmatrix}0&0\\\ 0&A_{22}^{(k)}\end{bmatrix}\begin{bmatrix}X_{1}\\\
X_{2}\end{bmatrix}=\begin{bmatrix}X_{1}\\\ X_{2}\end{bmatrix}P^{(k)}$ (20)

holds with $X=[X_{1}^{\top}X_{2}^{\top}]^{\top}$ full row rank and $P^{(k)}$
weakly diagonally dominant. Hence, $X_{2}$ has full row rank and

$A_{22}^{(k)}X_{2}=X_{2}P_{22}^{(k)}$ (21)

holds, which is equivalent to the existence of a weak PLF for the subsystem
(15). Conversely, if (21) holds with $X_{2}$ full row rank and $P_{22}^{(k)}$
diagonally dominant, then (20) holds with $X$ and $P^{(k)}$ as follows

$\begin{bmatrix}0&0\\\ 0&A_{22}^{(k)}\end{bmatrix}\begin{bmatrix}I&0\\\
0&X_{2}\end{bmatrix}=\begin{bmatrix}I&0\\\
0&X_{2}\end{bmatrix}\begin{bmatrix}0&0\\\ 0&P_{22}^{(k)}\end{bmatrix}.$

To complete the proof, note that, after having turned nodes $1,\dots,p$ into
black holes, the evolution of the GPN is represented by the integer operators

$\Psi(k)=\begin{bmatrix}0&0\\\ 0&I+B_{2,h}C_{2,h}^{\top}\end{bmatrix}$

If we apply Procedure 1, initialised with $X_{0}=[I,-I]$, we generate the
matrices $\mathbb{X}_{k}$

$\mathbb{X}^{0}=\left[\begin{array}[]{rr|rr}I&0&-I&0\\\
0&I&0&-I\end{array}\right]~{}\mbox{and}~{}\mathbb{X}_{k}=\begin{bmatrix}0\\\
\mathbb{X}_{2}^{k}\end{bmatrix},~{}~{}k>0$

where the first $p$ rows are zero, while the remaining $n-p$ rows are exactly
those we get by applying the procedure to subsystem (15). On the other hand,
Procedure 1 applied to the original system with the pinned nodes converges if
and only if a (weak) structural PLF exists for (15), as shown in [8]. ∎

The conditions of Theorem 4 are equivalent to the fact that Procedure 1
converges after having zeroed the rows of $B$ and the columns of $C$
corresponding to the pinned nodes. From a computational standpoint, however,
applying the procedure on the subsystem (15) is more convenient, because the
system size is reduced and the stopping criteria discussed in Section III-C
remain valid for the restricted subspace, while they are no longer valid for
the original state space.

###### Remark 9

_(Strong convergence.)_ According to Theorem 4, the differential inclusion can
be reduced as in (19). If matrix $B_{2}DC_{2}$ is structurally non-singular,
then pinning the first $p$ nodes makes the differential inclusion _strongly
convergent_ , according to the definition in [22].

It is worth stressing that a _dual procedure_ can be adopted as well. We can
derive a Lyapunov function (if it exists) defined in terms of planes, as in
(6). As shown in [8, 9], this is equivalent to applying the Procedure 1 to the
dual system $\dot{w}(t)=C^{\top}D(t)B^{\top}w(t)$. In view of duality
properties [11, 14, 31], a PLF exists for the primal system if and only if it
exists for its dual.

###  IV-B Arc pinning: Regulating the reactions

Node pinning means imposing a strong feedback to some nodes. By arc pinning,
we mean that a strong feedback is imposed to some flows, hence

$g_{j}(x):=g_{j}(x)-\gamma(g_{j}(x)-\bar{r}_{j}),$

so that, roughly speaking, the flow is forced to have a prescribed nominal
value.

This _dual_ arc-pinning problem can be solved by writing the system in
reaction coordinates, as in [1, 2, 3], and then adopting the
$EDF$-decomposition [9]. If $\bar{x}$ is an equilibrium point of the system,
under suitable conditions, we can define a transformation from concentration
coordinates $(x_{i})_{i=1}^{n}$ to reaction coordinates $(r_{j})_{j=1}^{m}$ as
$r(x(t))=g(x(t))-g(\bar{x})$. After this transformation, the system becomes

$\dot{r}(t)=\bigg{[}\frac{\partial g}{\partial x}\bigg{]}Sr(t).$ (22)

A procedure similar to the $BDC$-decomposition [9] transforms system (22) into
the linear differential inclusion

$\dot{r}(t)=ED(t)Fr(t),$ (23)

where $D$ is a diagonal matrix with positive diagonal entries. The theory
remains completely unchanged. Pinning reactions $1,\dots,p$ is equivalent to
zeroing the first $p$ rows of $E$ and the first $p$ columns of $F$.

###  IV-C Periodic forcing input

Assume that $g_{0}(t)$ is a periodic input [19, 24, 35] and $x_{p}(t)$ is a
periodic target trajectory, having the same period as $g_{0}(t)$,
corresponding to $g_{0}(t)$. Let $x(t)$ be any other trajectory. Then, we can
write

$\dot{x}(t)=Sg(x(t))+g_{0}(t)\,\,\mbox{and}\,\,\dot{x}_{p}(t)=Sg(x_{p}(t))+g_{0}(t).$

Denoting $z(t)=x(t)-x_{p}(t)$, the $BDC$-decomposition leads to the dynamical
system

$\dot{z}(t)=BD(z(t),t)Cz(t).$

Since our analysis considers a differential inclusion with _arbitrary_ time-
varying $D(\cdot)$, the stability – or stabilisation via pinning control – of
such a differential inclusion implies $z(t)\rightarrow 0$, hence the stability
– or stabilisation – of the periodic trajectory.

##  V Pinning Control of Reaction Networks

Pinning control, as mentioned, consists in applying strong local feedback
actions to some nodes (or arcs) with the aim of regulating the whole network.
How can one select the node(s) to be pinned in order to achieve the control
objective?

Based on our results, we can re-formulate the question as: Which are the nodes
that, if converted into black holes, ensure global boundedness of the GPN
evolution, hence leading to Lyapunov stability of the overall system once they
are subject to a sufficiently strong local feedback?

We provide here some examples of chemical reaction networks: the readers are
invited to have a preliminary look at the network graphs (shown in Figures 4
and 5) and see if they can spot immediately which nodes are the _most
important_ ones to be governed so as to rule all the others; the authors of
this paper often failed to guess these nodes in advance.

Figure 4: Graph of the chemical reaction network corresponding to the system
in Example 5.

###### Example 5

For the network in Fig. 4, Procedure 1 does not converge, hence the system
does not admit a structural polyhedral Lyapunov function. However, if we pin
any of the nodes $X_{2}$, $X_{4}$ or $X_{5}$, the procedure converges, hence
the system is structurally stabilised by enforcing a sufficiently strong local
feedback on any of these nodes. Conversely, the procedure does not converge
even if we pin node $X_{1}$ or node $X_{3}$. This fact is explained by
noticing that pinning $X_{2}$, $X_{4}$ or $X_{5}$ cuts the loop
$X_{2}\rightarrow X_{4}\rightarrow X_{5}\rightarrow X_{2}$, and tokens
repeatedly circulating in this loop continue depleting $X_{1}$ or filling up
$X_{3}$.

Figure 5: Graph of the chemical reaction network corresponding to the system
in Example 6.

###### Example 6

The network in Fig. 5 does not admit a structural PLF: Procedure 1 does not
converge. If we pin node $X_{3}$, then the procedure converges. This node is
then the most crucial one in the network: pinning any one of the other nodes
does not yield convergence.

###### Example 7

(A translation model.) The complete translation model proposed in [20]
includes the chemical reactions:

$\begin{array}[]{c}X_{1}+X_{2}\xrightleftharpoons[\hbox{\scriptsize$g_{3}$}]{\hbox{\scriptsize$g_{12}$}}X_{3}\xrightharpoonup{\hbox{\scriptsize$g^{*}_{3}$}}X_{4}+X_{2},\\\
X_{4}\xrightharpoonup{\hbox{\scriptsize$g_{4}$}}X_{5}\xrightharpoonup{\hbox{\scriptsize$g_{5}$}}X_{6}\xrightharpoonup{\hbox{\scriptsize$g_{6}$}}X_{7}\xrightharpoonup{\hbox{\scriptsize$g_{7}$}}X_{8}+X_{9},\\\
X_{8}\xrightharpoonup{\hbox{\scriptsize$g_{8}$}}X_{1},\qquad
X_{9}\xrightharpoonup{\hbox{\scriptsize$g_{9}$}}X_{10}\xrightharpoonup{\hbox{\scriptsize$g_{10}$}}\emptyset.\end{array}$
(24)

The stability of the reduced order model [20]

$X_{1}+X_{2}\xrightleftharpoons[\hbox{\scriptsize$g_{3}$}]{\hbox{\scriptsize$g_{12}$}}X_{3}\xrightharpoonup{\hbox{\scriptsize$g^{*}_{3}$}}X_{1}+X_{2}+X_{4},~{}~{}~{}X_{4}\xrightharpoonup{\hbox{\scriptsize$g_{4}$}}\emptyset$

is considered in [3], where the system is shown to admit a piecewise-linear
Lyapunov function in rates.

We consider here the complete model associated with the reaction network (24),
corresponding to the system of equations

$\begin{array}[]{rl}&\left\\{\begin{array}[]{rcl}\dot{x}_{1}&=&-g_{12}(x_{1},x_{2})+g_{3}(x_{3})+g_{8}(x_{8})\\\
\dot{x}_{2}&=&-g_{12}(x_{1},x_{2})+g_{3}(x_{3})+g_{3}^{*}(x_{3})\\\
\dot{x}_{3}&=&g_{12}(x_{1},x_{2})-g_{3}(x_{3})-g_{3}^{*}(x_{3})\\\
\dot{x}_{4}&=&g_{3}^{*}(x_{3})-g_{4}(x_{4})\\\
\dot{x}_{5}&=&g_{4}(x_{4})-g_{5}(x_{5})\\\
\dot{x}_{6}&=&g_{5}(x_{5})-g_{6}(x_{6})\\\
\dot{x}_{7}&=&g_{6}(x_{6})-g_{7}(x_{7})\\\
\dot{x}_{8}&=&g_{7}(x_{7})-g_{8}(x_{8})\end{array}\right.\\\ &\hskip
5.69054pt\begin{array}[]{rcl}\dot{x}_{9}&=&g_{7}(x_{7})-g_{9}(x_{9})\\\
\dot{x}_{10}&=&g_{9}(x_{9})-g_{10}(x_{10})\end{array}\end{array}$

We do not consider the last two equations: if the subsystem associated with
the variables $x_{1}$-$x_{8}$ converges to an equilibrium, and in particular
$x_{7}\rightarrow\bar{x}_{7}$, then also $x_{9}$ and, in turn, $x_{10}$
converge to an equilibrium.

It is apparent that
$\dot{x}_{1}+\dot{x}_{3}+\dot{x}_{4}+\dot{x}_{5}+\dot{x}_{6}+\dot{x}_{7}+\dot{x}_{8}=0$
and $\dot{x}_{2}+\dot{x}_{3}=0,$ hence the sums of concentrations
$x_{1}+x_{3}+x_{4}+x_{5}+x_{6}+x_{7}+x_{8}\doteq w_{1}$ and $x_{2}+x_{3}\doteq
w_{2}$ remain constant, thus forming a stoichiometric compatibility class,
which is bounded because $x\geq 0$. Hence, an equilibrium indeed exists [12].
We can then replace variables $x_{1}$ and $x_{2}$ by $w_{1}$ and $w_{2}$, and
keep all the others variables. This is equivalent to considering the state
transformation $T^{-1}x=w$ and the transformed system

$\dot{w}=T^{-1}BDCTw,$

with

$T^{-1}=\left[\begin{smallmatrix}1&0&1&1&1&1&1&1\\\ 0&1&1&0&0&0&0&0\\\
0&0&1&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\ 0&0&0&0&1&0&0&0\\\ 0&0&0&0&0&1&0&0\\\
0&0&0&0&0&0&1&0\\\ 0&0&0&0&0&0&0&1\end{smallmatrix}\right]$

and then neglecting the first two rows of $B$, which are of course zero
because $w_{1}$ and $w_{2}$ are both constant, as well as the first two
columns of $C$. The resulting reduced $6$-dimensional system admits the
$BDC$-decomposition with

$B=\left[\begin{smallmatrix}1&1&-1&0&0&0&0&0\\\ 0&0&1&-1&0&0&0&0\\\
0&0&0&1&-1&0&0&0\\\ 0&0&0&0&1&-1&0&0\\\ 0&0&0&0&0&1&-1&0\\\
0&0&0&0&0&0&1&-1\end{smallmatrix}\right]$

and

$C=\left[\begin{smallmatrix}-1&-1&-1&-1&-1&-1\\\ -1&0&0&0&0&0\\\
1&0&0&0&0&0\\\ 0&1&0&0&0&0\\\ 0&0&1&0&0&0\\\ 0&0&0&1&0&0\\\ 0&0&0&0&1&0\\\
0&0&0&0&0&1\end{smallmatrix}\right],$

where $C_{i}^{\top}B_{i}=-1$, for all $i$, as expected (cf. Remark 5).
Procedure 1 converges, yielding a structural PLF with 42 vertices, and matrix
$BDC$ passes the structural non-singularity test: this proves asymptotic
stability of the complete system.

###### Example 8

(A transcription model.) The complete transcription model proposed in [20] is

$\begin{array}[]{c}X_{1}\xrightleftharpoons[\hbox{\scriptsize$g_{2}$}]{\hbox{\scriptsize$g_{1}$}}X_{2}\xrightleftharpoons[\hbox{\scriptsize$g_{3}$}]{\hbox{\scriptsize$g^{*}_{2}$}}X_{3},\\\
X_{3}+X_{4}\xrightleftharpoons[\hbox{\scriptsize$g_{5}$}]{\hbox{\scriptsize$g_{34}$}}X_{5}\xrightharpoonup{\hbox{\scriptsize$g_{5}^{*}$}}X_{6}\xrightharpoonup{\hbox{\scriptsize$g_{6}$}}X_{4}+X_{7},\\\
X_{7}\xrightharpoonup{\hbox{\scriptsize$g_{7}$}}X_{8}\xrightharpoonup{\hbox{\scriptsize$g_{8}$}}X_{9}\xrightharpoonup{\hbox{\scriptsize$g_{9}$}}X_{10}\xrightharpoonup{\hbox{\scriptsize$g_{10}$}}X_{11}+X_{12},\\\
X_{11}\xrightharpoonup{\hbox{\scriptsize$g_{11}$}}X_{1},\quad
X_{12}\xrightharpoonup{\hbox{\scriptsize$g_{12}$}}\emptyset.\end{array}$ (25)

Again, [20] also proposed a reduced-order model

$X_{1}+X_{2}\xrightleftharpoons[\hbox{\scriptsize$g_{3}$}]{\hbox{\scriptsize$g_{12}$}}x_{3}\xrightharpoonup{\hbox{\scriptsize$g^{*}_{3}$}}X_{2}+X_{4}+X_{5},~{}~{}~{}x_{3}\xrightharpoonup{\hbox{\scriptsize$g_{3}$}}X_{1},~{}~{}~{}X_{4}\xrightharpoonup{\hbox{\scriptsize$g_{4}$}}\emptyset,$

for which a piecewise-linear Lyapunov function in rates is known to exists
[3].

We consider the system of differential equations associated with the complete
reaction network (25), which is

$\left\\{\begin{array}[]{rcl}\dot{x}_{1}&=&-g_{1}(x_{1})+g_{2}(x_{2})+g_{11}(x_{11})\\\
\dot{x}_{2}&=&g_{1}(x_{1})-g_{2}(x_{2})-g^{*}_{2}(x_{2})+g_{3}(x_{3})\\\
\dot{x}_{3}&=&g_{2}^{*}(x_{2})-g_{3}(x_{3})-g_{34}(x_{3},x_{4})+g(x_{5})\\\
\dot{x}_{4}&=&-g_{34}(x_{3},x_{4})+g(x_{5})+g_{6}(x_{6})\\\
\dot{x}_{5}&=&g_{34}(x_{3},x_{4})-g(x_{5})-g_{5}^{*}(x_{5})\\\
\dot{x}_{6}&=&g_{5}^{*}(x_{5})-g_{6}(x_{6})\\\
\dot{x}_{7}&=&g_{6}(x_{6})-g_{7}(x_{7})\\\
\dot{x}_{8}&=&g_{7}(x_{7})-g_{8}(x_{8})\\\
\dot{x}_{9}&=&g_{8}(x_{8})-g_{9}(x_{9})\\\
\dot{x}_{10}&=&g_{9}(x_{9})-g_{10}(x_{10})\\\
\dot{x}_{11}&=&g_{10}(x_{10})-g_{11}(x_{11})\end{array}\right.$

We can neglect the additional equation
$\dot{x}_{12}=g_{10}(x_{10})-g_{12}(x_{12})$, because, once we prove that the
$11$-order system is stable and converges to an equilibrium, convergence of
$x_{12}$ to an equilibrium immediately follows. For this system, matrices $B$
and $C$ are

$B=\left[\setcounter{MaxMatrixCols}{14}\begin{smallmatrix}-1&1&0&0&0&0&0&0&0&0&0&0&0&1\\\
1&-1&-1&1&0&0&0&0&0&0&0&0&0&0\\\ 0&0&1&-1&-1&-1&1&0&0&0&0&0&0&0\\\
0&0&0&0&-1&-1&1&0&1&0&0&0&0&0\\\ 0&0&0&0&1&1&-1&-1&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&1&-1&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&1&-1&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&1&-1&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&1&-1&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&1&-1&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&1&-1\end{smallmatrix}\right]$

and

$C=\left[\setcounter{MaxMatrixCols}{11}\begin{smallmatrix}1&0&0&0&0&0&0&0&0&0&0\\\
0&1&0&0&0&0&0&0&0&0&0\\\ 0&1&0&0&0&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0&0&0&0\\\
0&0&1&0&0&0&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0&0&0&0\\\ 0&0&0&0&1&0&0&0&0&0&0\\\
0&0&0&0&1&0&0&0&0&0&0\\\ 0&0&0&0&0&1&0&0&0&0&0\\\ 0&0&0&0&0&0&1&0&0&0&0\\\
0&0&0&0&0&0&0&1&0&0&0\\\ 0&0&0&0&0&0&0&0&1&0&0\\\ 0&0&0&0&0&0&0&0&0&1&0\\\
0&0&0&0&0&0&0&0&0&0&1\end{smallmatrix}\right].$

Note that, given that the corresponding derivatives add up to zero, the
following sums of concentrations are constant:
$x_{1}+x_{2}+x_{3}+x_{5}+x_{6}+x_{7}+x_{8}+x_{9}+x_{10}+x_{11}\doteq w_{1}$
and $x_{4}+x_{5}+x_{6}\doteq w_{2}$ for all $t$. Along with $x\geq 0$, this
ensures boundedness of the stoichiometric compatibility class, hence the
existence of an equilibrium [12]. Again, we replace $x_{1}$ and $x_{2}$ by
$w_{1}$ and $w_{2}$, for which $\dot{w}_{1}=0$ and $\dot{w}_{2}=0$, and the
corresponding equations are removed.

For this system, Procedure 1 does not converge, hence we cannot prove
structural stability. However, according to [20], a negative regulatory action
can be present, due to a repressor signal acting on variable $x_{4}$ (the DNA
promoter). To investigate the case in which $x_{4}$ is under a feedback
action, we can pin node $x_{4}$: then, Procedure 1 converges, providing a PLF
whose unit ball has $54$ vertices (the function in the reduced space has $52$
vertices).

##  VI Concluding Discussion

The main contribution of this paper is twofold. First, we have shown that, for
dynamical systems associated with unitary chemical reaction networks, the
existence of a polyhedral Lyapunov function is equivalent to the finiteness of
the reachable set of an associated generalised Petri net, with possibly
negative token numbers. Second, we have shown that applying a pinning control
to some nodes is structurally equivalent to converting the corresponding nodes
of the generalised Petri net into _black holes_ that swallow any incoming
token.

Pinning a node means applying a strong local feedback that keeps the node
variable constant. For a biochemical reaction network, enforcing actions that
keep the concentration of a species constant seems indeed a viable control
approach, which we conjecture is actually used in natural systems to stabilise
important cellular processes. In many cases, when a species is far more
abundant than all the other chemical species involved in the reaction network,
its concentration can be regarded as constant, because it is essentially
unchanged by the process, while other concentrations are subject to ample
fluctuations: this can already be seen as an “embedded” pinning control
action.

Future research directions along these lines include considering more general
types of structural feedback laws. Another aspect we leave for the future is
how to fit this framework in a stochastic setting in which the transitions are
probabilistic, going beyond the worst case structural analysis provided here.

## References

  * [1] M. A. Al-Radhawi and D. Angeli,  “Piecewise Linear in Rates Lyapunov Functions for Complex Reaction Networks”,  _Proc. IEEE Conf. on Dec. and Control_ , 2013, pp. 4595-4560.
  * [2] M. A. Al-Radhawi and D. Angeli,  “New approach to the stability of chemical reaction networks: Piecewise linear in rates Lyapunov functions”,  _IEEE Trans. Autom. Control_ , vol. 61, no. 1, pp. 76–89, 2016.
  * [3] M. A. Al-Radhawi, D. Angeli, and E. D. Sontag,  “A computational framework for a Lyapunov-enabled analysis of biochemical reaction networks”,  _PLOS Comp. Biology_ , vol. 16, no. 2, e1007681, 2020
  * [4] M. A. Al-Radhawi, D. Angeli, and E. D. Sontag,  “A robust Lyapunov criterion for non-oscillatory behaviors in biological interaction networks”,  arXiv:2009.10702.
  * [5] D. Angeli, P. De Leenheer and E. D. Sontag,  “A Petri net approach to the study of persistence in chemical reaction networks”,  _Mathematical Biosciences_ , vol. 210, no. 2, pp. 598-618, 2007.
  * [6] M. Bartoletti, T. Cimoli and G.M. Pinna,  “Lending Petri nets”,  Science of Computer Programming, vol. 112, no. 1, pp. 75-101, 2015.
  * [7] M. Bartoletti, P. Degano, P. Di Giamberardino and R. Zunino,  “Debits and credits in Petri nets and linear logic”,  In: Logic, Rewriting, and Concurrency. Lecture Notes in Computer Science, vol. 9200, Springer, Cham, 2015.
  * [8] F. Blanchini and G. Giordano,  “Piecewise-linear Lyapunov functions for structural stability of biochemical networks”,  _Automatica_ , vol. 50, no. 10, pp. 2482–2493, 2014.
  * [9] F. Blanchini and G. Giordano,  “Polyhedral Lyapunov functions for structural stability of biochemical systems in concentration and reaction coordinates”,  _Proc. IEEE Conf. on Dec. and Control_ , pp. 3110-3115, 2015.
  * [10] F. Blanchini and G. Giordano,  “Polyhedral Lyapunov functions structurally ensure global asymptotic stability of dynamical networks iff the Jacobian is non-singular”,  _Automatica_ , 86(12):183–191, 2017.
  * [11] F. Blanchini and G. Giordano,  “Dual chemical reaction networks and implications for Lyapunov-based structural stability”,  _IEEE Control Systems Letters_ , 2021.
  * [12] F. Blanchini and G. Giordano,  “Structural analysis in biology: a control-theoretic approach”,  _Automatica_ , 126(4):109376, 2021.
  * [13] F. Blanchini and S. Miani,  “On the transient estimate for linear systems with time-varying uncertain parameters”,  _IEEE Trans. Circuits and Systems I: Fundamental Theory and Applications_ , vol. 43, no. 7, pp. 592-596, 1996.
  * [14] F. Blanchini and S. Miani,  _Set-theoretic methods in control_.  Systems & Control: Foundations & Applications (2nd ed.), Birkhäuser, Basel, 2015.
  * [15] R.K. Brayton and C.H. Tong,  “Constructive stability and asymptotic stability of dynamical systems”,  IEEE Trans. Circuits and Systems, vol. 27, no. 11, pp. 1121–1130, 1980.
  * [16] D.A. Burbano-L., G. Russo and M. di Bernardo,  “Pinning controllability of complex network systems with noise”,  _IEEE Trans. Control of Network Systems_ , vol. 6, no. 2, pp. 874-883, 2019.
  * [17] F. Chen, Z. Chen, L. Xiang, Z. Liu and Z. Yuan,  “Reaching a consensus via pinning control”,  Automatica, vol. 45, no. 5, pp. 1215-1220, 2009.
  * [18] B.L. Clarke,  _Stability of Complex Reaction Networks_.  In _Advances in Chemical Physics_ (eds I. Prigogine and S. A. Rice), 1980.
  * [19] A. Cournac and J.A. Sepulchre,  “Simple molecular networks that respond optimally to time-periodic stimulation”,  _BMC Systems Biology_ , vol. 3, pp. 29, 2009.
  * [20] D. Del Vecchio and R.M. Murray,  _Biomolecular Feedback Systems_.  Princeton University Press, 2014.
  * [21] P. DeLellis, M. di Bernardo and F. Garofalo,  “Adaptive pinning control of networks of circuits and systems in Lur’e form”,  IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 60, no. 11, pp. 3033-3042, 2013.
  * [22] F. Fabiani, G. Belgioioso, F. Blanchini, P. Colaneri and S. Grammatico,  “Convergence in uncertain linear systems”,  _Automatica_ , vol. 119, 2020.
  * [23] M. Feinberg,  _Foundations of Chemical Reaction Network Theory_ , Springer, 2019.
  * [24] D. Fiore, A. Guarino and M. di Bernardo,  “Analysis and Control of Genetic Toggle Switches Subject to Periodic Multi-Input Stimulation”,  _IEEE Control Systems Letters_ , vol. 4, pp. 278-283, 2019.
  * [25] V. C. Gerogiannis, A. D. Kameas and P. E. Pintelas,  “Comparative study and categorization of high-level Petri nets”,  Journal of Systems and Software, 1998
  * [26] G. Giordano, C. Cuba Samaniego, E. Franco and F. Blanchini,  “Computing the structural influence matrix for biological systems”,  _J. Math. Biol._ , vol. 72, no. 7, pp. 1927–1958, 2016.
  * [27] R. Jungers,  _The Joint Spectral Radius: Theory and Applications_ ,  Lecture Notes in Control and Information Sciences, vol. 385, Springer, Berlin Heidelberg, 2009.
  * [28] M. Köhler-Bußmeier and M. Kudlek,  “Linear properties of zero-safe nets with debit tokens”,  Fundamenta Informaticae, vol. 84, no. 1-4, pp. 329-342, 2008.
  * [29] X. Li, X. Wang and G. Chen,  “Pinning a complex dynamical network to its equilibrium”,  _IEEE Trans. Circuits Syst. I Reg. Papers_ , vol. 51, no. 10, pp. 2074-2087, 2004.
  * [30] H. Maeda, S. Kodama, and Y. Ohta,  “Asymptotic behavior of nonlinear compartmental systems: Nonoscillation and stability”,  _IEEE Transactions on Circuits and Systems_ , vol. 25, no. 6, pp. 372–378, 1978.
  * [31] A. P. Molchanov and E. S. Pyatnitskii,  “Lyapunov functions that define necessary and sufficient conditions for absolute stability of nonlinear nonstationary control systems”, I:  Autom. Remote Control, 47(3):344–354, 1986; II: 47(4):443–451, 1986: III: 47(5):620–630, 1986.
  * [32] T. Murata and H. Yamaguchi,  “A Petri net with negative tokens and its application automated reasoning”,  In: _Proc. 33rd Midwest Symposium on Circuits and Systems_ , 1990, pp. 762-765.
  * [33] Y. Orouskhani, M. Jalili, and X. Yu,  “Optimizing dynamical network structure for pinning control”,  Scientific Reports, vol. 6, 24252, 2016.
  * [34] M. Porfiri and M. di Bernardo,  “Criteria for global pinning-controllability of complex networks”,  Automatica, 2008.
  * [35] G. Russo, M. di Bernardo, and E. Sontag,  “Global entrainment of transcriptional systems to periodic inputs”,  _PLoS Computational Biology_ , 6:e1000739, 2010.
  * [36] D. Sokolov, I. Poliakov, and A. Yakovlev,  “ Asynchronous Data Path Models ”,  Proceedings - 7th International Conference on Application of Concurrency to System Design, 2007.
  * [37] S. Soliman,  “Invariants and other structural properties of biochemical models as a constraint satisfaction problem”,  _Algorithms for Molecular Biology_ , vol. 7, no. 1, p. 15, 2012.
  * [38] X. Wang and H. Su,  “Pinning control of complex networked systems: A decade after and beyond”,  Annual Reviews in Control, vol. 38, no. 1, pp. 103–111, 2014.

## Proof of Lemma 1

The set ${\cal S}$ is positively invariant for (17) for any large $\mu>0$.
Consider the modified set

${\cal S}_{\nu}=\\{z\in{\cal S}\colon\|z_{1}\|\leq\nu\\},$

which corresponds to the portion between the cyan planes in Fig. 6. For any
$\nu>0$ (no matter how small) there exist $\hat{\mu}$ such that ${\cal
S}_{\nu}$ becomes positively invariant for $\mu\geq\hat{\mu}$.

Figure 6: Graphical representation of Lemma 1.

Indeed,

$\dot{z}_{1}=F_{11}z_{1}-\mu z_{1}+F_{12}z_{2},$

where both $F_{12}z_{2}$ and $F_{11}z_{1}$ are bounded terms, because ${\cal
S}_{\nu}$ is a compact set. Therefore, we can write
$|z_{1}^{\top}F_{12}z_{2}|\leq\alpha$ and $|z_{1}^{\top}F_{11}z_{1}|\leq\beta$
for suitable positive coefficients $\alpha$ and $\beta$. Consider the
candidate Lyapunov function $W(z_{1})=\frac{1}{2}\|z_{1}\|^{2}$ and the ball
$\|z_{1}\|^{2}\leq\nu^{2}$. Then

$\dot{W}(z_{1})=z_{1}^{\top}F_{11}z_{1}-\mu
z_{1}^{\top}z_{1}+z_{1}^{\top}F_{12}z_{1}\leq-\mu
z_{1}^{\top}z_{1}+\alpha+\beta,$

hence $\dot{W}(z_{1})<0$ for $\|z_{1}\|^{2}>\nu^{2}$ provided that
$\mu\geq\hat{\mu}=(\alpha+\beta)/\nu^{2}$.

Therefore, assume that $\|z_{1}\|\leq\nu$, $\mu\geq\hat{\mu}$. Take the
initial condition $z_{1}(0)=0$ and $\bar{z}_{2}=z_{2}(0)$ on the boundary of
${\cal S}_{2}$.

By contradiction, assume that the solution $\bar{z}_{2}(t)$ of
$\dot{z}_{2}=F_{22}z_{2}$ leaves the set ${\cal S}_{2}$. Consider this
solution in the extended $z$ space,
$\hat{z}(t)=[0~{}\bar{z}_{2}(t)^{\top}]^{\top}$. There exists a time instant
$T$ such that $\hat{z}(T)=[0~{}\bar{z}_{2}(T)^{\top}]^{\top}$ is outside the
compact ${\cal S}$ and there exists a neighbourhood ${\cal U}$ (the ball in
Fig. 6) centred at $\bar{z}(T)$, which has no intersection with ${\cal S}$.
Note that $\hat{z}(T)$ does not depend on $\mu$.

Now, consider the solution $\tilde{z}(t)$ of the full system (which depends on
$\mu$) with the same initial condition
$\bar{z}=[0~{}\bar{z}_{2}(T)^{\top}]^{\top}$.

We complete the proof by showing that $\tilde{z}(T)$ gets arbitrarily close to
$\hat{z}(T)$ if $\mu$ is large enough. The first component satisfies
$\|\tilde{z}_{1}\|\leq\nu$. The second component satisfies
$\dot{\tilde{z}}_{2}=F_{22}\tilde{z}_{2}+F_{21}\tilde{z}_{1}.$ Then, the
difference $\tilde{z}_{2}(t)-\bar{z}_{2}(t)$ satisfies

$\frac{d}{dt}[\tilde{z}_{2}-\bar{z}_{2}]=F_{22}[\tilde{z}_{2}-\bar{z}_{2}]+F_{21}\tilde{z}_{1},~{}~{}\|\tilde{z}_{1}(t)\|\leq\nu,~{}~{}\forall
t$

with $\tilde{z}_{2}(0)-\bar{z}_{2}(0)=0$. Hence

$\displaystyle\|\tilde{z}_{2}(T)-\bar{z}_{2}(T)\|$ $\displaystyle=$
$\displaystyle\left\|\int_{0}^{T}~{}e^{F_{22}(T-t)}F_{21}\tilde{z}_{1}(t)dt\right\|$
$\displaystyle\leq$
$\displaystyle\|\nu\|\left\|\int_{0}^{T}~{}e^{F_{22}(T-t)}F_{21}dt\right\|.$

Given any small $\rho>0$, we may ensure
$\|\tilde{z}_{2}(T)-\bar{z}_{2}(T)\|\leq\rho$ by forcing a small enough $\nu$.
Since both $\nu$ and $\rho$ can be arbitrarily small, we get that
$\hat{z}(T)\in{\cal U}$, hence it is outside ${\cal S}$, against the
invariance assumption. We have reached a contradiction, which completes the
proof.

|  Franco Blanchini (Senior Member, IEEE) was born on 29 December 1959, in
Legnano (Italy). He is the Director of the Laboratory of System Dynamics at
the University of Udine. He has been involved in the organization of several
international events: in particular, he was Program Vice-Chairman of the
conference Joint CDC-ECC 2005, Seville, Spain; Program Vice-Chairman of the
Conference CDC 2008, Cancun, Mexico; Program Chairman of the Conference
ROCOND, Aalborg, Denmark, June 2012 and Program Vice-Chairman of the
Conference CDC 2013, Florence, Italy. He is co-author of the book “Set
theoretic methods in control”, Birkhäuser. He received the 2001 ASME Oil & Gas
Application Committee Best Paper Award as a co-author of the article
“Experimental evaluation of a High-Gain Control for Compressor Surge
Instability”, the 2002 IFAC prize Survey Paper Award as the author of the
article “Set Invariance in Control - a survey”, Automatica, November 1999, for
which he also received the High Impact Paper Award in 2017, and the 2017 NAHS
Best Paper Award as a co-author of the article “A switched system approach to
dynamic race modelling”, Nonlinear Analysis: Hybrid Systems, 2016. He was
nominated Senior Member of the IEEE in 2003. He has been an Associate Editor
for Automatica, from 1996 to 2006, and for IEEE Transactions on Automatic
Control, from 2012 to 2016. From 2017 to 2019 he has been an Associate Editor
for Automatica. He has been a Senior Editor for IEEE Control Systems Letters.
---|---
|  Carlos Andrés Devia (Student Member, IEEE) received his B.Sc. and M.Sc. in
electronic engineering in 2015 and 2018 respectively, from Pontificia
Universidad Javeriana, Bogotá, Colombia. From 2017 to 2019 he was an
instructor at the same institution. Since 2019 he is a PhD candidate at the
Delft Center for Systems and Control, Delft University of Technology, The
Netherlands. His research interests include networked systems with emphasis on
biological systems and opinion formation models.
---|---
|  Giulia Giordano (Member, IEEE) received the B.Sc. and M.Sc. degrees summa
cum laude in electrical engineering and the Ph.D. degree (Hons.) in systems
and control theory from the University of Udine, Italy, in 2010, 2012, and
2016, respectively. She visited the California Institute of Technology,
Pasadena (CA), USA, in 2012, and the University of Stuttgart, Germany, in
2015. She was a Research Fellow at Lund University, Sweden, from 2016 to 2017,
and an Assistant Professor at the Delft University of Technology, The
Netherlands, from 2017 to 2019. She is currently an Assistant Professor at the
University of Trento, Italy. She was recognised with the Outstanding Reviewer
Letter from the IEEE Transactions on Automatic Control in 2016 and from the
Annals of Internal Medicine in 2020. She received the EECI Ph.D. Award 2016
for her thesis “Structural Analysis and Control of Dynamical Networks”, the
NAHS Best Paper Prize 2017, as a coauthor of the article “A Switched System
Approach to Dynamic Race Modelling”, Nonlinear Analysis: Hybrid Systems, 2016,
and the SIAM Activity Group on Control and Systems Theory Prize 2021. Her main
research interests include the study of dynamical networks, the analysis of
biological systems, and the control of networked systems.
---|---


#  Is Wireless Bad for Consensus in Blockchain?

Seungmo Kim, Senior Member, IEEE S. Kim is with the Department of Electrical
and Computer Engineering, Georgia Southern University in Statesboro, GA, USA
and can be reached at<EMAIL_ADDRESS>
###### Abstract

This paper examines how wireless communication affects the performance of
various blockchain consensus mechanisms, focusing on their scalability and
decentralization. It introduces an analytical framework for quantifying these
effects, backed by extensive simulations, underscoring its broad applicability
to various consensus mechanisms despite wireless communication’s
unreliability.

###### Index Terms:

Consensus, Scalability, Decentralization, Wireless communications

##  I Introduction

###  I-1 Motivation

In the dynamic landscape of modern connectivity, the proliferation of wireless
technologies has become nothing short of permeating every facet of our daily
lives and redefining the way we connect, communicate, and consume information.
In fact, over 55% of website traffic comes from mobile devices, and 92.3% of
internet users access the internet using a mobile phone [1]. In fact, a wide
variety of wireless technologies immerse our daily lives, including Wi-Fi [2],
cellular technologies such as the fifth-generation (5G) cellular [3] and the
Long-Term Evolution (LTE) [4], to vehicle-to-everything (V2X) communications
[5] and even wearable technologies [6, 7].

This makes a strong case where increasingly more blockchains will be
established on wirelessly connected nodes. One should notice of the key
challenge here: the reliability of wireless communications is generally lower
in comparison to wired communications (e.g., ethernet), attributed to various
random factors such as noise, interference, and limited bandwidth [8]. One
should also note that this lower reliability affects the performance of
consensus in blockchain [9, 10].

All open public blockchains are based on the idea that they should be able to
reach consensus across a distributed network, even when there are conflicts,
without putting control in one place [11]. Nonetheless, the additional
dynamicity that is brought by the wireless networking has not been thoroughly
discussed in the literature as of yet. Therefore, it will be a valuable
academic attempt to formally analyze the performance of a blockchain system
established on a wireless network.

###  I-2 State of the Art

PoW [12] and PoS [13] remain the two most common mechanisms, especially in the
cryptocurrencies’ context [14]. PoW is a form of cryptographic proof in which
one party (the prover) proves to others (the verifiers) that a certain amount
of a specific computational effort has been expended [15]. The main purpose of
PoW is to deter manipulation of data by establishing large energy and
hardware-control requirements to be able to do so, which inevitably leaves the
it criticized by environmentalists for their energy consumption. Meanwhile, as
an effort to avoid such high computational cost that PoW causes, PoS is a type
of consensus mechanism for blockchains that are designed to elect validators
in proportion to their quantity of holdings in the associated cryptocurrency.

Moreover, there have been proposed a variety of techniques as an effort to
improve the performance of blockchains: namely, parallel structure [16], off-
chain [17], reinforcement learning [18], and combination of PoW and PoS [19].
Nonetheless, the blockchain trilemma sets a limit on this desire: no
blockchain can achieve improvements on all three fronts of scalability,
security, and decentralization at once. In fact, despite the large research
and experimental effort, all known approaches turn out to leave tradeoffs
[16].

This makes a compelling case where we urgently need a comprehensive framework
evaluating the performance of blockchains with wirelessly connected nodes with
respect to the blockchain trilemma. Many of the current consensus algorithms
already consider the possibilities that nodes leave, and new nodes join.
However, we claim that more is needed: due to the uncertainty in wireless
connections, a blockchain system established on a wireless network will draw a
completely different environment than a one with ethernet-connected nodes
[20].

Before our work, diverse avenues of recent research was sparked, exploring
topics such as blockchain performance in the context of IoT [21],
communication resource consumption [22], efficient consensus mechanisms among
wireless nodes [23], connectivity-adaptive contract mechanisms [24], and
Byzantine fault-tolerant (BFT) consensus among wireless nodes [25].
Specifically, the existing analytical framework (e.g., universal scalability
law [26]) overlooked the impacts of wireless connectivity, which may cause
serious imprecision these days with such a high proportion of wireless
technologies in any given network.

Furthermore, most of the previous work in the literature of blockchain and
distributed systems present experimental results, which can only be applied to
certain practical scenarios with a set of specific parameter settings.
Therefore, this paper attempts to provide a more general framework that can be
applied to any other blockchains for measurement of their consensus
performance.

###  I-3 Contributions

Addressing the aforementioned limitations of the existing literature, this
papers sheds light on precisely evaluating the performance of consensus
mechanisms when the nodes are wirelessly connected. The specific technical
contributions are summarized as follows:

  * •

It provides an analytical framework for calculation of scalability and
decentralization of a blockchain consensus process;

  * •

It draws probabilitistic analysis for formulating the impacts of wireless
connections on the scalability and decentralization.

  * •

It presents a comparitive study among PoW, PoS, vs PoC with respect to
scalability and decentralzation.

Figure 1: Example scenario of gossip protocol among 400 nodes with
$\mathsf{p}_{\text{fail}}=0.5$

##  II System Model

The key focus of this paper is the nodes being wirelessly connected. As such,
the system model is supposed in such a way that the impacts of wireless
connection on the blockchain performance can be clearly characterized.

The nodes are connected to each other via wireless communications. Each node
holds a 100 m of communication range. This paper assumes a general ad-hoc
network in which no central mediator exists; instead, the nodes are connected
to their own neighbors, e.g., device-to-device communications [27]. The
rationale is that a centralized network (e.g., cellular, Wi-Fi, etc.) will
likely employ a permissioned blockchain (or a private blockchain), whose
performance is not exactly what we are interested in measuring.

As a means to accomplish finality of data among nodes, distributed systems
often adopt a gossip protocol for propagation of information [28]. Fig. 1
shows an example of the gossip protocol operated on a wireless ad-hoc network.
Each node propagates the block to all the neighboring nodes within its
communication range. Each connection is supposed to cause a Byzantine fault at
the rate of $\mathsf{p}_{\text{fail}}$, which is shown as a red line. Such a
failed connection yields a failed reception at the receiver node, which is
drawn as a red dot.

A two-dimensional space $\mathbb{R}^{2}$ is defined as a 1 km-by-1 km square,
as illustrated in Fig. 1. We model spatial distribution of the stationary
nodes as a Poisson point process (PPP) in $\mathbb{R}^{2}$ with density
$\lambda$. Fig. 1 visualizes an example distribution of nodes with
$\lambda=400$. Notice that the $\lambda$ is the parameter that will be varied
according to the type of consensus mechanism. For instance, PoS has a smaller
$\lambda$ as compared to PoW since its consensus is operated only among the
validators who are with more than a required amount of stake.

Figure 2: Scalability R versus the number of nodes
($\mathsf{p}_{\text{fail}}=5\%$, 10 transactions per block)

Figure 3: Decentralization level G vs. node failure rate
(PoW–$\mathsf{r}_{\text{cls}}=\\{0.1,0.3,0.5\\}$;
PoS–$\mathsf{r}_{\text{v}}=0.2$; PoC–$\mathsf{r}_{\text{sfl}}=0.9$,
$\Delta_{\text{sfl}}=0$)

##  III Consensus among Wireless Nodes

Node failures occur in a “clustered” manner; the magenta-colored boxes
indicate the areas where the failures transpire. It is assumed that the node
failures take place within the clusters only; the example presented in Fig. 1
assumes the probability of each node failure to be
$\mathsf{p}_{\text{fail}}=0.5$. We will vary this $\mathsf{p}_{\text{fail}}$
as a method to evaluate the performance of the consensus mechanisms–viz., PoW,
PoS, and PoC.

The key discussion of this paper is that the wireless connection can increase
the odds of Byzantine faults. Specifically, we identify “clustered occurrence”
of faults attributed to interference among each other [29] as the key
characteristic that the wireless connection adds to Byzantine faults, which
yields that the nodes located in certain areas face a higher probability of
disconnection. We attribute such clustered disconnections to interference
among nodes. In fact, an ad-hoc wireless network is more susceptible to
interference [30], due to its decentralized network structure where no central
mediator exists to coordinate scheduling among nodes.

Specifically, we define the consensus mechanisms that we analyze in this
paper–i.e., PoW, PoS, and PoC–and identify their key characteristics that may
be affected by the wireless connections. Moreover, it is noteworthy that we
pay attention to PoC and compare it to PoW and PoS, the two most predominating
consensus mechanisms. The rationale is that PoC adopts a significantly
different principle from that of PoW and PoS, which thus distinguishes it from
the vast majority of recent consensus mechanisms built on PoW and PoS.

##  IV Analysis on Consensus Performance

Now, we present the theoretical analysis framework that this paper proposes,
which features two metrics measuring the performance of a consensus mechanism.

###  IV-1 Scalability

The literature of distributed systems often rely on the universal scalability
law [31, 32] as the key scalability measure. Nonetheless, this quantity is
usually obtained empirically, which makes it almost impossible to analytically
quantify $\alpha$ and $\beta$. Moreover, we could find an extensive empirical
analysis of the parameters for Zookeeper [26], one of the popular distributed
systems [33]. In fact, Zookeeper has appeared as on many of the blockchain
systems [34, 35]. This makes a case that we adopt the result of S found for
the Zookeeper as a scalability benchmark for this paper’s result.

Motivated from the challenge, we define a generalized metric for the
scalability in this paper, which is given by

$\displaystyle\mathsf{R}=\displaystyle\frac{\text{\\# transactions}}{\text{\\#
seconds}}=\frac{\mathsf{n}_{\text{tx}}}{\mathsf{T}_{c}}$ (1)

where $\mathsf{T}_{c}$ denotes the length of time for a single consensus.

###  IV-2 Decentralization

We propose to adopt Gini coefficient (denoted by G) to represent the level of
decentralization. To elaborate, we modify G in such a way that the
“inequality” indicates the level of participations being concentrated to a
fewer number of nodes, which is understood as “centralized.”

Building on this understanding, one can take a more formal way to write this
decentralization analysis. We define the Gini coefficient as half of the
relative mean absolute difference, which is equivalent to the definition based
on the Lorenz curve [36]. The mean absolute difference is the average absolute
difference of all pairs of items of the population, and the relative mean
absolute difference is the mean absolute difference divided by the average,
$\bar{x}$, to normalize for scale. If $x_{i}$ is the number of node $i$’s
parcipations in consensus, and there are $n$ nodes, then G is given by

$\displaystyle\textsf{G}=\frac{\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n}\left|x_{i}-x_{j}\right|}{\displaystyle
2\sum_{i=1}^{n}\sum_{j=1}^{n}x_{j}}=\frac{\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n}\left|x_{i}-x_{j}\right|}{2n^{2}\bar{x}}.$
(2)

##  V Numerical Results

Fig. 3 compares the scalability, $\mathsf{R}$, among the three consensus
mechanisms–i.e., PoW, PoS, and PoC. One can easily find that PoC is the most
scalable among the three since it is designed to select a fixed number of
witnesses, which also explains the “flatness” of its scalability versus the
number of nodes. Due to such selectivity, PoC should pose a lower level of
decentralization, as shall be presented in Fig. 3.

The scalability of PoS depends on $\mathsf{r}_{\text{v}}$, indicating the rate
of validators. We set the parameter as a Poisson random variable with a
justification that there must be a “certain amount” of tokens that the largest
number of participants hold.

Fig. 3 indicates that a higher wireless disconnection rate increases the
decentralization level. One should also notice that the three consensus
mechanisms pose the order of PoW $>$ PoC $>$ PoS in terms of the
decentralization level. For PoS, the rate of validators among all the nodes is
set to $\mathsf{r}_{\text{v}}=0.2$. For PoC, the rate of node replacement in
each round of consensus is set to $\mathsf{r}_{\text{sfl}}=0.9$, and the
number of rounds between a replacement is $\Delta_{\text{sfl}}=0$. To wit,
$\mathsf{r}_{\text{sfl}}$ gives how many of the current witnesses will be
replaced with other ones; and $\Delta_{\text{sfl}}$ means how frequently the
replacement occurs.

The most striking result is that wireless connection helps improve the
decentralization in blockchain! The reason for this result is straightforward.
The intermittent wireless connection hinders such participation. Therefore,
the level of “inequality” is decreased among nodes participating in
consensuses.

The rate of cluster (denoted by $\mathsf{r}_{\text{cls}}$) was varied with PoW
only. The rationale is that PoW is particularly susceptible to such
clusteredness of node failure transpiration. In contrast, PoS and PoC will be
less impacted since they already adopt “selective” participation of nodes in a
consensus; to wit, if failures occur at non-participating nodes, they have
almost no impact on consensus.

##  VI Future Work

We will extend the findings of this paper to analyzing how mobility [37]
affects the consensus in blockchain. Moreover, we will investigate impacts of
off-chain PoC migration (as part of the recent migration to Solana [38]) on
this paper’s result.

## References

  * [1] J. Howarth, “Internet traffic from mobile devices (Nov. 2023),” [Online]. Available: https://explodingtopics.com/blog/mobile-internet-traffic
  * [2] S. Kim, J. Choi, and C. Dietrich, “PSUN: an OFDM - pulsed radar coexistence technique with application to 3.5 GHz LTE,” Hindawi Mobile Inf. Syst. J., Jul. 2016.
  * [3] S. Kim, E. Visotsky, P. Moorut, K. Bechta, A. Ghosh, and C. Dietrich, “Coexistence of 5G with the incumbents in the 28 and 70 GHz bands,” IEEE J. Sel. Areas Commun., vol. 35, iss. 8, Aug. 2017.
  * [4] S. Kim and C. Dietrich, “Coexistence of outdoor Wi-Fi and radar at 3.5 GHz,” IEEE Wireless Commun. Lett., vol. 6, iss. 4, Aug. 2017.
  * [5] S. Kim and B. J. Kim, “Novel backoff mechanism for mitigation of congestion in DSRC broadcast,” arXiv:2005.08921, May 2020.
  * [6] S. Kim, Y. Sharif, and I. Nasim, “Human electromagnetic field exposure in wearable communications: a review,” arXiv:1912.05282, Dec. 2019.
  * [7] I. Nasim and S. Kim, “Human EMF exposure in wearable networks for internet of battlefield things,” arXiv:2003.01552, Mar. 2020.
  * [8] R. Jordan and C. Abdallah, “Wireless communications and networking: An overview,” IEEE Antennas Propag. Mag., vol. 44, no. 1, 2002.
  * [9] X. Zhang, W. Xia, Q. Cui, X. Tao, and R. Liu, “Efficient and trusted data sharing in a sharding-enabled vehicular blockchain,” IEEE Netw., Aug. 2022.
  * [10] S. Okegbile, J. Cai, and A. Alfa, “Performance analysis of blockchain-enabled data-sharing scheme in cloud-edge computing-based IoT networks,” IEEE Internet Things J., vol. 9, no. 21. Jun. 2022.
  * [11] M. Bedawala, “A deeper dive into consensus mechanisms.”
  * [12] S.Nakamoto, “Bitcoin: A peer-to-peer electronic cash system,” Oct. 2008. [Online]. Available: https://bitcoin.org/bitcoin.pdf
  * [13] Ethereum Foundation, “Proof-of-stake (PoS),” Sep. 2023. [Online]. Available: https://ethereum.org/en/developers/docs/consensus-mechanisms/pos/
  * [14] European Parliament, “Cryptocurrencies and blockchain,” PE 619.024, Jul. 2018.
  * [15] C. Dwork and M. Naor, “Pricing via processing or combatting junk mail,” LNCS Advances Cryptology, vol. 740, 1992.
  * [16] G. Monte, D. Pennino, and M. Pizzonia, “Scaling blockchains without giving up decentralization and security: A solution to the blockchain scalability trilemma,” in Proc. Wksp Cryptocurrencies Blockchains Distributed Syst. 2020.
  * [17] S. Reno and M. Haque, “Solving blockchain trilemma using off-chain storage protocol,” IET Inf. Secur., vol. 4, Jul. 2023.
  * [18] S. Kim and A. S. Ibrahim, “Byzantine-fault-tolerant consensus via reinforcement learning for permissioned blockchain implemented in a V2X network,” IEEE Trans. Intell. Veh., vol. 8, iss. 1, Jan. 2023.
  * [19] Trifecta Blockchain Team, “Trifecta: The blockchain trilemma solved,” White Paper, Oct. 2019. [Online]. Available: https://pramodv.ece.illinois.edu/pubs/Whitepaper2019-9.pdf
  * [20] T. Rathod, N. Jadav, M. Alshehri, S. Tanwar, R. Sharma, R. Felseghi, and M. Raboaca, “Blockchain for future wireless networks: A decade survey,” MDPI Sensors, vol. 22, no. 11, May 2022.
  * [21] S. Alrubei, E. Ball, J. Rigelsford, and C. Willis, “Latency and performance analyses of real-world wireless IoT-blockchain application,” IEEE Sensors J., vol. 20, no. 13, Mar. 2020.
  * [22] L. Zhang, H. Xu, O. Onireti, M. Imran, and B. Cao, “How much communication resource is needed to run a wireless blockchain network?,” IEEE Netw., vol. 36, no. 1, Sep. 2021.
  * [23] X. Xie, C. Hua, J. Hong, P. Gu, and W. Xu, “AirCon: Over-the-air consensus for wireless blockchain networks,” IEEE Trans. Mobile Comput., Jul. 2023.
  * [24] W. Wang, J. Chen, Y. Jiao, J. Kang, W. Dai, and Y. Xu, “Connectivity-aware contract for incentivizing IoT devices in complex wireless blockchain,” IEEE Internet Things J., vol. 10, no. 12, Jun. 2023.
  * [25] Y, Zou, L. Yang, G. Jing, R. Zhang, Z. Xie, H. Li, and D. Yu, “A survey of fault tolerant consensus in wireless networks,” Elsevier High-Confidence Comput., Jan. 2024.
  * [26] M. Bevilacqua-Linn, M. Byron, P. Cline, J. Moore, and S. Muir, “Sirius: Distributing and coordinating application,” in Proc. USENIX Ann. Techn. Conf. 2014.
  * [27] 3GPP, “Universal Mobile Telecommunications System (UMTS); LTE; Proximity-based services (ProSe); Stage 2,” ETSI TS 123 303, V17.1.0, Jul. 2023.
  * [28] M. Jelasity, Gossip, Springer Berlin Heidelberg, pp. 139–162, [Online]. Available: https://link.springer.com/chapter/10.1007/978-3-642-17348-6_7
  * [29] H. ElSawy, E. Hossain, and M. Haenggi, “Stochastic geometry for modeling, analysis, and design of multi-tier and cognitive cellular wireless networks: A survey,” IEEE Commun. Surveys Tut., vol. 15, no. 3, Third Quarter 2013.
  * [30] P. Cardieri, “Modeling interference in wireless ad hoc networks,” IEEE Commun. Surveys Tut., vol. 12, no. 4, Fourth Quarter 2010.
  * [31] P. Jogalekar and M. Woodside, “Evaluating the scalability of distributed systems,” IEEE Trans. Parallel Dist. Syst., vol. 11, no. 6, Jun. 2000.
  * [32] N. Lu and X. Shen, “Scaling laws for throughput capacity and delay in wireless networks–A survey,” IEEE Commun. Surveys Tut., vol. 16, no. 2, 2014.
  * [33] P. Hunt, M. Konar, F. Junqueira, and B. Reed, “ZooKeeper: Wait-free coordination for internet-scale systems,” in Proc. USENIX Ann. Technol. Conf. 2010.
  * [34] E. Bandara, D. Tosh, P. Foytik, S. Shetty, N. Ranasinghe, and K. Zoysa, “Tikiri–Towards a lightweight blockchain for IoT,” Elsevier Future Generation Comput. Syst., vol. 119, 2021.
  * [35] J. Duan, A. Karve, V. Sreedhar, and S. Zeng, “Service management of blockchain networks,” in Proc. IEEE Int. Conf. Cloud Comput. 2018.
  * [36] A. Sen, On Economic Inequality, Oxford University Press, 1977.
  * [37] M. Grossglauser and D. Tse, “Mobility increases the capacity of ad hoc wireless networks,” IEEE/ACM Trans. Netw., vol. 10, no. 4, Aug. 2002.
  * [38] Helium Foundation, “HIP-70: Scaling the Helium network,” May 2023.