# Robust Control of Partially Specified Boolean Networks

Luboš Brim
Faculty of Informatics
Masaryk University
<EMAIL_ADDRESS>
&Samuel Pastva
Faculty of Informatics
Masaryk University
<EMAIL_ADDRESS>
&David Šafránek
Faculty of Informatics
Masaryk University
<EMAIL_ADDRESS>
&Eva Šmijáková
Faculty of Informatics
Masaryk University
<EMAIL_ADDRESS>

###### Abstract

Regulatory networks (RNs) are a well-accepted modelling formalism in
computational systems biology. The control of RNs is currently receiving a lot
of attention because it provides a theoretical and computational basis for
cell reprogramming – an attractive technology developed in regenerative
medicine. By solving the control problem, we learn which parts of a biological
system should be perturbed in order for it to stabilise the system in the
desired phenotype.

In this paper, we use _Boolean networks (BNs) with asynchronous update_ to
represent the dynamics of an RN as a discrete finite-state system.
Furthermore, we allow the specification of the Boolean model representing a
given RN to be incomplete (partial). This is crucial in cases where the exact
Boolean update functions of the BN are not fully known (e.g. due to a lack of
measured data). To that end, we utilise the formalism of _partially specified
Boolean networks_ which allows us to cover every possible behaviour reflecting
the unspecified parts of the system. Such an approach inevitably causes a
significant state explosion. This problem is efficiently addressed by using
_symbolic methods_ to represent both the unspecified model behaviour as well
as all possible perturbations of the system. Our framework supports several
ways to control the system by employing one-step, temporary, and permanent
perturbations.

Additionally, to make the control design efficient and practically applicable,
the optimal control should be minimal in terms of the size – the number of
perturbed system components. Moreover, in a partially specified model, a
control may achieve the desired stabilisation only for a subset of the
possible fully specified instantiations of the model. To address these aspects
in the control design, we utilise several quantitative measures. In
particular, apart from the _size_ of perturbation necessary to achieve
control, we also examine its _robustness_ – a portion of fully specified
models for which the control is applicable.

We show that the proposed symbolic methods solving the control problem for
partially specified BNs are efficient and scale well with the number of
unspecified model elements. The provided experiments demonstrate that our
algorithms can handle relatively large BNs. We also evaluate the robustness
metrics in cases of all three studied control types. The robustness metric
tells us how big a proportion of fully defined systems the given perturbation
works. Our experiments support the hypothesis that one-step perturbations may
be less robust than temporary or permanent perturbations.

This is a full version of a paper that is submitted to a journal.

_K_ eywords Boolean Network $\cdot$ Robustness $\cdot$ Perturbation $\cdot$
Temporary Control $\cdot$ Permanent Control $\cdot$ Symbolic Algorithm

##  1 Introduction

The control theory is on the rise in computational systems biology thanks to
its practical potential [1, 2, 3]. In particular, the solution to the control
problem for an _in silico_ model of a cell can provide the fundamental basis
for experimental design for reprogramming of the cell _in vitro_ [4]. The
common applications of cell reprogramming are differentiation of T-cells [5,
6, 7], reprogramming of embryonic stem cells [8, 9, 10, 11, 12], curing
obesity by converting fat cells to the brown type [13], and recovering
malfunctioning cells [14] or even organs [15].

_Regulatory network_ (RN) is a widely used model for _in silico_ experiments
because it allows capturing numerous relationships in complex biochemical
systems [16, 17]. Nonetheless, RNs represent a highly abstract level of view
on a system, and more specific (mechanistic) model adaptations are usually
employed for realising computational tasks such as a simulation of the system
dynamics. Boolean models of RNs [18, 19], or Boolean networks (BNs) for short,
are often employed to study the RN dynamics. In BNs, the dynamical behaviour
is specified by a set of Boolean update functions defining the evolution of
variables over time. In this paper, we focus on BNs with the asynchronous
update that reflects well the concurrent dynamics of individual components
(genes) in the RN [20, 21]. This updating scheme results in a non-
deterministic state-transition graph (STG), as opposed to a simpler,
deterministic graph which arises under the less realistic synchronous scheme
[22].

A critical problem of BN modelling is to fully determine the update functions
from data or literature [23]. In many cases, the knowledge of update functions
is incomplete, and the modellers must address this problem [24]. Other
possible situations when multiple model adaptations should be considered at
the same time are differences between intrinsic and extrinsic observations
[25], or structural changes in RNs occurring over time [26]. We address these
issues by targeting analysis of _partially specified_ Boolean models. We
represent the partial specification in BN by allowing the use of
_uninterpreted functions_ as part of BN’s update functions. Intuitively, such
an uninterpreted function stands for some unknown but fixed part of the
network’s dynamics. By fixing the behaviour of these uninterpreted functions,
we can then obtain a wide range of standard (i.e. fully specified) BNs, which
we refer to as _instantiations_. Technically, in our framework a partially
specified BN in represented as a BN with inputs.

The goal of _Boolean network control_ is to influence the behaviour of the
network so that it stabilises in a given _attractor_ – a minimal component of
the system, which once reached, cannot be escaped. A typical way to influence
(change) the BN’s behaviour is to perturb the values of some of its variables.
Since the practical realisation of perturbations requires non-trivial effort,
we typically require the number of perturbations to be as small as possible.
To that end, the control is usually optimised by minimising the number of
controlled variables [27].

However, the second factor that needs consideration when the network is not
fully known is that a perturbation may succeed only for some fraction of
instantiations of the partially specified network. Meanwhile, for other
instances, the target attractor may not be reached. We quantify this aspect
using a _robustness_ metric – a portion of fully specified model
instantiations for which the resulting perturbation can control the respective
network to achieve the desired attractor.

There are multiple possible ways a perturbation can be applied to a biological
system. In this paper, we focus on _temporary_ and _permanent_ perturbations.
As such, we consider that the value of perturbed variables changes
simultaneously and is held fixed for multiple time-steps. In the temporary
case, the system is eventually left to behave as unperturbed. In contrast, the
variables are fixed ad infinitum during permanent perturbations.

There is a large body of existing work on the BN control problem. Most
commonly, the _source-target_ variant of the BN control problem is considered,
where both the source and the target state (or attractor) are specified. Other
types of control include target control [28] (the controlled system must reach
the target attractor regardless of the initial state) and full-control [29] (a
control strategy between all attractor pairs is desired).

In the case of BNs utilising the synchronous update, the problem is far less
complicated, with efficient approaches already shown in [28, 30, 31]. The main
reason is that in a synchronous BN, due to the determinism of the network, its
attractors are simpler (only single-state and cyclic attractors appear), and
from each initial state, only one attractor is reachable. Meanwhile, an
asynchronous network can contain larger, more complex attractors, and multiple
attractors can be reachable from any network state. Owing to the more detailed
mechanistic view of the concurrent regulatory events [21], control of
asynchronous BNs is currently receiving a lot of attention. For example, the
approach developed in [7, 32, 33] computes a set of relevant BN variables
based on the identification of particular motifs in the RN. However,
perturbations based on such stable motifs control the BNs only when applied
permanently.

In some cases, such permanent perturbations may not be feasible [34]. To that
end, the concept of a _one-step_ perturbation has been developed in [27]. This
type of perturbation is applied in the initial state, and after that, the
system evolves according to its original dynamics. This approach has motivated
a branch of methods based on the identification of the _strong basin_ of the
target attractor (i.e. the states which _must_ eventually reach said
attractor). Subsequently, the concept of temporary perturbations was
introduced and studied in [35, 36, 37], filling the gap between one-step and
permanent perturbations. Temporary perturbations were shown to be very
efficient: While the perturbation does not need to be enforced indefinitely,
the number of perturbed variables necessary to control the network is usually
minimal compared to a one-step perturbation. Even more complex control
strategies, considering temporal perturbation sequences are studied in [38,
39, 40, 41]. The tool CABEAN comprises these methods [42].

Another line of research addressing temporary perturbations is based on the
so-called trap-spaces [43]. In this approach, the target attractor is also
only partially known. The developed method allows the specification of the
target via a phenotype – a subset of desired variable values, instead of an
attractor. Such phenotype control is explored in more detail in [44], where a
method based on model checking is developed. However, the scalability of this
method may be limited by the necessity to scan through the available
perturbations.

All results mentioned so far have been developed only for fully specified BNs.
It is worth noting that the incomplete specification is far more
computationally demanding. Each possible full model specification generates a
new BN with a unique STG. Therefore, the _state-space explosion_ (state-space
size is exponential in the number of BN variables) is multiplied by the number
of possible BN instances. In the worst case, this number is doubly-exponential
with respect to the number of inputs of the network update functions [45]. The
size of the resulting multi-graph makes control algorithms developed for fully
specified asynchronous BNs (by computing the control for each BN instance
distinctly) intractable. Consequently, addressing the control problem in
partially specified BNs requires substantially novel algorithms.

In [46], the control of partially specified BNs have been addressed for the
first time. In particular, the paper introduced a method for one-step source-
target control of BNs with logical parameters (encoding the partial
specification of BNs). On the technical side, the method has been supported by
a semi-symbolic algorithm. In particular, the unspecified parts are
represented symbolically using binary decision diagrams (BDDs) [47] while the
STG is explicit [48]. The algorithm is based on the identification of the
attractor’s strong basin for all parametrisations simultaneously.

In this paper, the methodology and the algorithmics for the control of
partially specified BNs are significantly advanced. First, we extend the semi-
symbolic approach to an entirely symbolic graph representation. This allows
processing models with a higher number of variables. Second, we develop new
algorithms for more permissive types of perturbations (in terms of minimal
size and robustness) – temporary and permanent. Last but not least, the
novelty of our approach lies in a new method we use to simulate the perturbed
BN. Instead of constructing a new graph for each perturbation, we model all
perturbations in a single transition-labelled graph. We show that our approach
scales well for highly unspecific models and outperforms the semi-symbolic
approach from [46] primarily in terms of processable model size. On
biologically relevant models, we also demonstrate that temporary and permanent
perturbations are more robust compared to available one-step perturbations.

##  2 Theory

This section presents a formal introduction to the topics discussed in this
paper. Mainly, it introduces the modelling framework of partially specified
Boolean networks and the control problems for this type of model.

At first, let us define $\mathbb{B}=\\{0,1\\}$ to be the set of Boolean
values, and $\mathbb{B}_{*}=\\{0,1,*\\}$ to be an extension of $\mathbb{B}$
which also admits a _free_ value $*$ (i.e. neither $\mathit{true}$ nor
$\mathit{false}$). Additionally, we write $\mathbb{B}^{n}$ to denote the set
of all $n$-element vectors over $\mathbb{B}$. For each vector
$x\in\mathbb{B}^{n}$, by convention, $x_{i}$ is the $i$-th element and
$x[i\mapsto b]$ is the copy of $x$ with the $i$-th element set to
$b\in\mathbb{B}$.

Asynchronous Boolean networks: As the name suggests, a Boolean network (BN)
assumes $n$ Boolean _variables_ , whose state evolves over time. The state of
such a network is therefore an $n$-element vector $s\in\mathbb{B}^{n}$ (we say
that $\mathbb{B}^{n}$ is the _state space_ of the BN). Additionally, we assume
that the network can have $m$ Boolean _inputs_ (sometimes also referred to as
_constants_), which have arbitrary but fixed values. Consequently, we say that
$\mathbb{B}^{m}$ is the input space of the BN and $c\in\mathbb{B}^{m}$ is an
input valuation. Finally, the Boolean network $B$ itself is a collection of
update functions $B=(b_{1},\ldots,b_{n})$, such that each
$b_{i}:\mathbb{B}^{n}\times\mathbb{B}^{m}\to\mathbb{B}$ takes the current
network state and input valuation, and outputs a new value for the $i$-th
variable. Typically, if function $b_{i}$ depends on the value of variable $j$,
we say that $j$ _regulates_ $i$. These relationships are prescribed by a
directed _influence graph_ (also called a _regulatory network_) which can
serve as a visual representation of the BN structure. Note that the inputs and
variables are not explicitly distinguished in some literature since an input
can be simply seen as a variable with identity as its update function.

We assume that the variables of the network are updated _asynchronously_ ,
i.e. each state-transition updates exactly one variable based on the output of
$b_{i}$. Consequently, we can describe the dynamics of a BN with inputs using
a _coloured state-transition graph_. A coloured graph
$\textsc{STG}(B)=(C,V,E)$ consists of three elements. A finite set of
_colours_ $C=\mathbb{B}^{m}$ corresponding to the input valuations of the
network. A finite set of _vertices_ $V=\mathbb{B}^{n}$ corresponding to the
state space of the network. And finally, the edges $E\subseteq V\times C\times
V$ are constructed based on the update functions using the asynchronicity
assumption explained above:

$\displaystyle(s,c,t)\in E\Leftrightarrow\exists i\in[1,n].(t=s[i\mapsto
b_{i}(s,c)]\land s\not=t)$

Additionally, to guarantee that $E$ is total for every $c$, we assume that a
self-loop $(s,c,s)\in E$ is present for a state $s$ that has no other outgoing
edges. We write $s\xrightarrow{c}t$ when $(s,c,t)\in E$, and
$s\xrightarrow{c}^{*}t$ when $(s,c,t)\in E^{*}$ (i.e. the reflexive and
transitive closure of $E$). Intuitively, the coloured graph $\textsc{STG}(B)$
can be seen as a family of standard directed graphs over a common set of
vertices, where the choice of input valuation determines the edges of the
graph. Because all variables and inputs are Boolean, this structure has a
fairly straightforward symbolic representation in terms of binary decision
diagrams, as we later demonstrate.

Additionally, for any coloured graph we can define a _run_ $\pi$ under colour
$c$ as a maximal sequence of states $s,s_{1},\ldots$ such that
$(s_{i},c,s_{i+1})\in E$ for every $i$. Note that such $\pi$ is necessarily
infinite due to $E$ being total. State $s_{0}$ is the _initial state_ of such
a run. State $s_{0}$ is the _initial state_ of such a run. We write $\pi(i)$
to denote the $i$-th state $s_{i}$ of $\pi$, and $inf(\pi)$ to denote all
states which appear in $\pi$ infinitely often.

Since Boolean networks primarily function as models of biological systems
which have some inherent (but in our case unquantified) stochasticity, we
typically assume that the system can only exhibit _fair_ runs. That is,
available transitions cannot be delayed indefinitely [49]. To capture this
property, we write $\Pi_{c}(s_{0})$ to denote the set of all fair runs under
colour $c$ which originate in the state $s_{0}$.

Partially specified Boolean networks: A Boolean network with inputs allows us
to easily encode a wide range of biochemical systems in a machine-friendly
format. However, for systems with a high degree of uncertainty, it often fails
to capture this uncertainty in a way that is understandable to a human reader.

To address this issue, the actual input of our method are _partially
specified_ Boolean networks, which explicitly allow parts of the update
functions to be designated as unknown. Assume that $f_{1}^{(a_{1})}$,
$f_{2}^{(a_{2})}$, $\ldots$ are symbols standing in for some uninterpreted
(fixed but arbitrary) Boolean functions ($a_{i}$ being the function arity). A
partially specified BN then assumes $n$ variables and $p$ uninterpreted
functions. In this kind of network, every update function $b^{\prime}_{i}$ is
specified as a Boolean expression that can contain the function symbols
$f_{1},\ldots,f_{p}$.

This formalism is often easier to comprehend, as the uncertainty in dynamics
is tied to the update functions instead of inputs. However, it is also not
immediately clear how such a network should be represented symbolically. It is
relatively easy to translate a partially specified network into a BN with
inputs, though. Any uninterpreted function $f_{i}^{(a)}$ can be encoded in
terms of $2^{a}$ fresh Boolean inputs if we consider that $j$-th input denotes
the output of $f_{i}^{(a)}$ in the $j$-th row of its function table. Formally,
this translation can be achieved using a repeated application of the following
expansion rule:

$\displaystyle f(\alpha_{1},\ldots,\alpha_{a})\equiv(\alpha_{1}\Rightarrow
f^{\prime}_{1}(\alpha_{2},\ldots,\alpha_{a}))\land(\neg\alpha_{1}\Rightarrow
f^{\prime}_{2}(\alpha_{2},\ldots,\alpha_{a}))$

Here, $f^{\prime}_{1}$ and $f^{\prime}_{2}$ are fresh uninterpreted functions
of arity $a-1$, and $\alpha_{i}$ are arbitrary Boolean expressions. Using this
rule, we can always convert a partially specified network to a Boolean network
with inputs. Nonetheless, the number of inputs will be exponential with
respect to the arity of the employed uninterpreted functions (since each
application of the rule replaces one uninterpreted function with two fresh
ones).

Network attractors and source-target control: To understand the dynamics of a
BN (or a partially specified BN translated to a BN with inputs), we have to be
familiar with the notion of _attractor_. Informally, an attractor is the
smallest set of states that cannot be escaped and hence the system eventually
converges into it. Formally, an attractor $A\subseteq V$ for colour $c\in C$
is the smallest set such that for every $s\in A$, if $s\xrightarrow{c}t$ for
some $t\in V$, then $t\in A$. Note that for every fair run $\pi$, $inf(\pi)$
is always some attractor of the network. In Boolean networks, different
attractors typically correspond to different real-world phenotypes which the
system can exhibit.

Therefore, for a real-world Boolean network, one of the fundamental questions
is whether (and if so, how, or even how efficiently) the network can be
controlled into a specific attractor. We assume that the network is controlled
by a _perturbation_ which fixes the values of certain variables to constants.
Formally, a perturbation is a vector $q\in\mathbb{B}_{*}^{n}$, where for every
$q_{i}=*$, we say that the $i$-th variable is unperturbed, whereas for
$q_{i}=0$ or $q_{i}=1$, we say that the variable is perturbed to either $0$ or
$1$. For a state $s\in\mathbb{B}^{n}$, we then write $q(s)$ to denote a copy
of $s$ where the perturbation $q$ has been applied. Formally, $q(s)_{i}=s_{i}$
when $q_{i}=*$, and $q(s)_{i}=q_{i}$ otherwise. Finally, we overload the
notation $\Pi_{c}^{q}(s)$ to denote the set of fair runs under colour $c$
starting in state $s$, where the runs are restricted to the edges of
$\textsc{STG}(B)$ that only update the _unperturbed_ variables in $q$. For the
sake of maintaining the property that every fair run is infinite, we assume
that this restriction admits a self-loop on every state that cannot update any
unperturbed variable (i.e. any state that would otherwise have no outgoing
transitions). That is, intuitively, along the run in $\Pi_{c}^{q}(s)$, the
perturbed variables are never updated, and if a state from which no outgoing
edges using unperturbed variable exist, the state is repeated forever.

While the notion of a perturbation should now be clear, the exact temporal
nature of the perturbation can be still different based on the biological
assumptions of the perturbed system and the real-world perturbation method. In
every instance, a system in state $s$ is perturbed into the state $q(s)$ by
applying the perturbation $q$. However, we recognise three types of state
perturbations based on the subsequent evolution of the system: (a) _one-step_
perturbation: after perturbation, all system variables again evolve freely;
(b) _permanent_ perturbation: once perturbed, only unperturbed variables can
evolve. Perturbed variables are fixed to their respective constant values; (c)
_temporary_ perturbation: after perturbation, the perturbed variables remain
constant for some arbitrary but finite number of time-steps but are eventually
allowed to evolve freely again.

Based on these three types of perturbations, we can then formalise three
control problems. In every case, we assume a source state
$source\in\mathbb{B}^{n}$ together with a target state
$target\in\mathbb{B}^{n}$, and a fixed input valuation $c\in\mathbb{B}^{m}$:

  * •

One-step control: Identify every perturbation $q\in\mathbb{B}_{*}^{n}$ s.t.
for every $\pi\in\Pi_{c}(q(source))$, $target\in inf(\pi)$. That is, after
performing the perturbation, $target$ is always visited infinitely often
without further restricting the network in any way.

  * •

Permanent control: Identify every perturbation $q\in\mathbb{B}_{*}^{n}$ s.t.
for every $\pi\in\Pi_{c}^{q}(q(source))$, $target\in inf(\pi)$. That is,
$target$ is always visited infinitely often assuming the perturbed variables
remain constant.

  * •

Temporary control: Identify every perturbation $q\in\mathbb{B}_{*}^{n}$ s.t.
for every $\pi\in\Pi_{c}^{q}(q(source))$, there is some $j\in\mathbb{N}$ where
for every $\pi^{\prime}\in\Pi_{c}(\pi(j))$, $target\in inf(\pi^{\prime})$.
Intuitively, the system can evolve for an arbitrary number of time steps while
assuming the perturbed variables are constant, but eventually, every such run
reaches a state from which the system can evolve into $target$ as if
unperturbed.

In the case of one-step and temporary control, such a perturbation only exists
when $target$ is a member of some attractor of $B$ (or perturbed $B$ in case
of permanent control). Otherwise, no fair run will visit it infinitely often,
regardless of applied perturbations. It is therefore typically reasonable to
only consider input valuations for which $target$ is indeed an attractor state
[46]. On the contrary, the permanent perturbation may induce the appearance of
new attractors, including attractors containing the $target$ state.
Nonetheless, this scenario is solving a much harder biological task -
introduction of a new phenotype (phenotypic variation) [50]. We leave this
line of research for future work and consider only colours containing an
attractor with the $target$ state.

Ideally, we would also like to select for perturbations that are biologically
feasible (typically, due to physical limitations, not every variable can be
perturbed) and that perturb the smallest number of variables (to minimise the
cost and complexity of a perturbation). However, the validity of a
perturbation is also tied to the associated input valuation $c$, which may
describe some unknown parts of the system that we cannot influence. As such,
most perturbations will correctly control the network only for some subset of
input valuations. Depending on how large this subset is, we can then talk
about _robustness_ $\rho$ of a one-step, permanent, or temporary perturbation
$q$ with respect to the network inputs [46]:

$\displaystyle\centering\rho(q)=\frac{|\\{c\in\mathbb{B}^{m}\mid q\text{
controls }B\text{ under }c\\}|}{|\\{c\in\mathbb{B}^{m}\mid target\text{ is an
attractor state}\\}|}\@add@centering$

Robustness is a critical factor when selecting viable perturbations in BNs
with inputs. To compute this property, our goal is to compute a _control
relation_ $\mathcal{C}_{?}\subseteq\mathbb{B}^{n}_{*}\times\mathbb{B}^{m}$
(here, $?$ is either $O$ for one-step, $P$ for permanent, or $T$ for
temporary). This relation includes all pairs
$(q,c)\in\mathbb{B}^{n}_{*}\times\mathbb{B}^{m}$ such that $q$ controls $B$
under input valuation $c$ from $source$ to $target$ assuming one-step
($\mathcal{C}_{O}$), permanent ($\mathcal{C}_{P}$), or temporary
($\mathcal{C}_{T}$) control.

Our robustness measure gives us the ability to compare the “quality” of
particular perturbations. In order to also evaluate the “quality” of different
control relations (or control approaches), we extend the notion of robustness
to be measurable on a set of perturbations. Specifically, given a set of
perturbations $Q$, we define _maximal robustness_ $\rho_{max}(Q)$ to be the
value of maximal robustness within the given set, and _union robustness_
$\rho_{union}(Q)$ to be the union of all input valuations controlled by the
given set with respect to the relevant colours:

$\rho_{max}(Q)=max\\{\rho(q)\mid q\in Q\\}$
$\rho_{union}(Q)=\frac{|\bigcup_{q\in Q}\\{c\in\mathbb{B}^{m}\mid q\text{
controls }B\text{ under }c\\}|}{|\\{c\in\mathbb{B}^{m}\mid target\text{ is an
attractor state}\\}|}$

Intuitively, the maximal robustness $\rho_{max}$ gives us the best robustness
that can be achieved within a particular set using a single perturbation.
Meanwhile, union robustness $\rho_{union}$ is the measure of how well a system
can be controlled _overall_ using perturbations from $Q$. Both measures are
worth exploring, as each represents a slightly different optimisation goal:
While a set $Q_{1}$ with a high $\rho_{max}$ is generally more likely to
provide a practically viable control strategy, a set $Q_{2}$ with a high
$\rho_{union}$ may achieve control even for cases that are not feasible using
$Q_{1}$.

##  3 Methods

Real-world Boolean networks can have a large number of inputs and variables
(and consequently, many admissible perturbations). It is not possible to
explore the STGs of such networks explicitly. We propose several symbolic
algorithms to address this issue. These are based on the existing symbolic
approaches to control of BNs _without_ inputs (i.e. networks with complete
knowledge). However, the addition of inputs introduces extra complexity, which
our methods need to address. In particular, identifying a single minimal
perturbation is not sufficient, as the robustness of such a perturbation can
be very low, making it unlikely to work in practice. To provide a reasonably
reliable method, one needs to consider a large set of perturbations.
Consequently, for performance reasons, our algorithms cannot test only a
single perturbation at a time but need to explore the whole perturbation space
symbolically.

Symbolic computation model: We employ symbolic representation using binary
decision diagrams (BDDs) [47]. In this representation, a BDD corresponding to
a Boolean formula $\varphi$ is used to represent a set (or a relation) $X$ of
Boolean vectors where the elements of $X$ are exactly the satisfying
valuations of $\varphi$. Consequently, any subset of vertices or colours of
$\textsc{STG}(B)$, i.e. $V=\mathbb{B}^{n}$ and $C=\mathbb{B}^{m}$, or any
relation $R\subseteq V\times C$ can be represented in this way. To represent
sets of perturbed variables, we will also refer to a set $P=\mathbb{B}^{n}$.
But note that the members of this set are not perturbations per se (i.e.
members of $\mathbb{B}_{*}^{n}$). For each $p\in P$, value $p_{i}=1$ indicates
only that the variable $i$ is perturbed, not the perturbation value itself.

The reasoning behind this choice of encoding is that after the network is
perturbed, the perturbation values are contained within the state of the
network variables. To get the full picture of the network’s state under
multiple different perturbations and input valuations, we can thus use a
relation $X\subseteq V\times C\times P$. For $(s,c,p)\in X$, $p_{i}=1$
indicates that, assuming input valuation $c$, $i$-th variable is perturbed,
and $s_{i}$ gives the actual value of the perturbation (for $p_{i}=0$, $s_{i}$
is the normal mutable state of the network). Finally, we use
$\mathcal{S}=V\times C\times P$ as a shorthand for this extended set of
perturbable states and input valuations. We also write that a pair $(s,p)\in
V\times P$ is _equivalent_ to a perturbation $q\in\mathbb{B}_{*}^{n}$,
$(s,p)\equiv q$, if $q_{i}=*$ when $p_{i}=0$, and $q_{i}=s_{i}$ when
$p_{i}=1$. Observe that this approach is more efficient than simply
incorporating the set of all possible perturbations $\mathbb{B}^{n}_{*}$
directly into the encoding. Each $\mathbb{B}_{*}$ would require more than one
bit per network variable to encode, wherease our approach only adds one bit on
top of the existing coloured state space.

On these symbolic sets and relations, we can perform standard set operations
($\cap,\cup,\setminus,\times$, etc.) as these correspond to logical operations
on BDDs. Additionally, we use $\exists_{P}(X\subseteq\mathcal{S})$ to denote
existential projection over the perturbed variables, i.e.
$\exists_{P}(X)=\\{(s,c)\mid\exists p\in P.(s,c,p)\in X\\}$. This can also be
implemented using a single BDD operation. Finally, we need a mechanism to
symbolically explore the edges of the STG $\textsc{STG}(B)$, but we also need
to consider variable perturbations in this process. To do so, we define the
two following symbolic operations:

$\displaystyle\textsc{Pre}(X\subseteq\mathcal{S})$
$\displaystyle=\\{(s,c,p)\in\mathcal{S}\mid\exists t\in V.\ \exists
i\in[1,n].\ t=s[i\mapsto b_{i}(s,c)]\land s\not=t\land p_{i}=0\\}$
$\displaystyle\textsc{Post}(X\subseteq\mathcal{S})$
$\displaystyle=\\{(t,c,p)\in\mathcal{S}\mid\exists s\in V.\ \exists
i\in[1,n].\ t=s[i\mapsto b_{i}(s,c)]\land s\not=t\land p_{i}=0\\}$

Recall that this is almost exactly the condition under which $(s,c,t)\in E$ of
$\textsc{STG}(B)$. As such, Pre and Post compute the sets of predecessors and
successors of the states in $X$ under their respective input valuations.
However, for each computed transition, we also require that $p_{i}=0$, i.e.
the $i$-th variable is _not_ perturbed. Hence values of the perturbed
variables remain constant, and in case $p_{i}=0$ for all $i$, we get exactly
the behaviour of $\textsc{STG}(B)$ without any perturbations. These two
operations can also be implemented using BDDs, since $b_{i}$ is a Boolean
function (which itself trivially corresponds to a BDD), and the rest are
either logical operations or quantifications that can be realised within the
domain of our BDD representation.

Utility algorithms: Our algorithms rely on the utility procedures defined in
Algorithm 1. Here, Bwd is a standard backward reachability procedure, albeit
over members of our extended state space $\mathcal{S}$. Meanwhile, the result
of Trap is a so-called _trap set_ [43]: a subset of $X$ from which one cannot
escape $X$. Note that if $X$ contains some attractors of $B$,
$\textsc{Trap}(X)$ will contain exactly those attractors and states from $X$
that can only reach these attractors. Meanwhile, if $X$ does not contain any
attractor, $\textsc{Trap}(X)$ will be empty.

Finally, note that while our Post and Pre procedures _respect_ the setting of
perturbed variables within subsets of $\mathcal{S}$, they do not actually
perform any perturbations (i.e. they do not _force_ perturbed variables to
their perturbed values). To work around this issue, Algorithm 1 also describes
a CanPerturb function. This function takes a _source_ state, and a set of
_target_ states, and computes the subset of _target_ states that can be
reached from the _source_ via a perturbation:

To compute this information, we notice that the elements of targets not
belonging to $\textsc{CanPerturb}(\texttt{source},\texttt{target})$ are based
on exactly those combinations of $(s,p)\in V\times P$ for which there is at
least one variable $i$ where $s_{i}$ does not match $\texttt{source}_{i}$, but
$p_{i}=0$ (variable is not perturbed). That is, a “perturbation step” cannot
change a value of an unperturbed variable: all values in source are either
perturbed, in which case a perturbation will update their value to match $s$,
or unperturbed, in which case they must already agree with $s$. We thus go
through all variables and remove all cases where this situation occurs using
symbolic operations.

Fn  _Bwd $(X\subseteq\mathcal{S})$_

repeat $X\leftarrow X\cup\textsc{Pre}(X)$ until  _$X\textbf{ reaches
fixpoint}$_ ;

return _$X$_ ;

Fn  _Trap $(X\subseteq\mathcal{S})$_

/* Eliminate elements that escape from $X$ in one step. */

repeat $X\leftarrow X\setminus\textsc{Pre}(\textsc{Post}(X)\setminus X)$ until
_$X\textbf{ reaches fixpoint}$_ ;

return _$X$_ ;

Fn  _CanPerturb $(\texttt{source}\in V,\texttt{targets}\subseteq\mathcal{S})$_

for  _$i\in\\{1\ldots n\\}$_ do

$\texttt{diff}_{i}\leftarrow\\{s\in V\mid s_{i}\not=\texttt{source}_{i}\\}$;

$\texttt{zero}_{i}\leftarrow\\{p\in P\mid p_{i}=0\\}$;

$\texttt{targets}\leftarrow\texttt{targets}\setminus(\texttt{diff}_{i}\times
C\times\texttt{zero}_{i}))$;

return _targets_ ;

Algorithm 1 Basic graph algorithms.

Fn  _Permanent $(\texttt{source}\in V,\texttt{target}\in V)$_

$\texttt{basin}\leftarrow\textsc{Trap}(\textsc{Bwd}(\\{\texttt{target}\\}\times
C\times P))$;

$\texttt{map}\leftarrow\textsc{CanPerturb}(\texttt{source},\texttt{basin})$;

return
_$(q,c)\in\mathcal{C}_{P}\Leftrightarrow(s,c,p)\in\texttt{map}\land(s,p)\equiv
q$_ ;

Fn  _OneStep $(\texttt{source}\in V,\texttt{target}\in V)$_

$\texttt{basin}\leftarrow\textsc{Trap}(\textsc{Bwd}(\\{\texttt{target}\\}\times
C\times\\{0\\}^{n}))$;

$\texttt{basin}\leftarrow\exists_{P}(\texttt{basin})\times P$;

$\texttt{map}\leftarrow\textsc{CanPerturb}(\texttt{source},\texttt{basin})$;

return
_$(q,c)\in\mathcal{C}_{O}\Leftrightarrow(s,c,p)\in\texttt{map}\land(s,p)\equiv
q$_ ;

Fn  _Temporary $(\texttt{source}\in V,\texttt{target}\in V)$_

$\texttt{basin}\leftarrow\textsc{Trap}(\textsc{Bwd}(\\{\texttt{target}\\}\times
C\times\\{0\\}^{n}))$;

$\texttt{basin}\leftarrow\exists_{P}(\texttt{basin})\times P$;

$\texttt{basin}\leftarrow\textsc{Trap}(\textsc{Bwd}(\texttt{basin}))$;

$\texttt{map}\leftarrow\textsc{CanPerturb}(\texttt{source},\texttt{basin})$;

return
_$(q,c)\in\mathcal{C}_{T}\Leftrightarrow(s,c,p)\in\texttt{map}\land(s,p)\equiv
q$_ ;

Algorithm 2 Control procedures.

Control procedures: Algorithm 2 then describes bottom-up procedures which
compute the relations $\mathcal{C}_{O}$, $\mathcal{C}_{P}$, and
$\mathcal{C}_{T}$. In the case of Permanent control, the variable basin
contains all combinations of states, input valuations and perturbed variables
from which target is visited infinitely often on fair paths. Specifically, Bwd
finds all states that can reach target using some fair run. Out of these, Trap
eliminates the ones that may not always reach target but have the opportunity
to also escape into other attractors. In particular, this also completely
eliminates input valuations where target is not a part of any attractor.
Finally, CanPerturb is used to only retain combinations that can be reached
from source via a perturbation.

OneStep relies on a very similar principle, however, the initial computation
of basin is happening without any perturbations enabled: $P$ is exchanged for
a singleton set $\\{0\\}^{n}$. This adheres to the intuition that in the one-
step control problem, the evolution after the perturbation is not constrained.
However, once basin is computed, we use $\exists_{P}$ to substitute
$\\{0\\}^{n}$ for $P$. This effectively “enables” all possible perturbations
and procedure CanPerturb can then assesses which of such perturbations are
feasible from the source state.

Finally, Temporary is conceptually a combination of both approaches: first, a
basin is computed with all perturbations disabled. Then, $\\{0\\}^{n}$ is
exchanged for $P$, which enables perturbations. Afterwards, the combination of
Bwd and Trap procedures is repeated, yielding a larger set of states that
eventually always reach target only if the corresponding perturbation is
applied temporarily. Finally, CanPerturb again limits this set to
perturbations that are actually feasible from source.

One last consideration is regarding the format of relations $\mathcal{C}_{?}$.
Notice that in our chosen symbolic representation (BDDs), all elements have
Boolean domains. Consequently, we cannot represent a control relation
$\mathcal{C}\subseteq\mathbb{B}_{*}^{n}\times\mathbb{B}^{m}$ directly
($\mathbb{B}_{*}^{n}$ is not a set of Boolean vectors). However, the result of
CanPerturb already contains all information contained within $\mathcal{C}_{?}$
– in the return statements, we thus simply describe a way of translating the
information between $\mathcal{C}$ and the existing coloured vertex set map. As
such, the map variable essentially represents one of the possible Boolean
encodings of a non-Boolean relation
$\mathcal{C}\subseteq\mathbb{B}_{*}^{n}\times\mathbb{B}^{m}$.

##  4 Results

To evaluate our methods, we have implemented algorithms from the previous
section by using internal libraries of the tool AEON [51]. The prototype
implementation is available as an open-source Github
repository111https://github.com/sybila/biodivine-pbn-control. Our experiments
were run using a computer with AMD Ryzen Threadripper 2990WX 32-Core Processor
and 64GB of memory. The results of experiments are discussed in Section 5.

Our first experiment measures the performance of methods on five biological
models with varying state space size (induced by the number of variables), a
number of inputs, and a possible number of perturbations. All these three
components influence the size of the coloured graph under consideration.
Whereas the number of vertices is equal to the number of states, the number of
all STG instances is equivalent to the product of colours and possible
perturbations.

The first model is a _myeloid_ differentiation network [52] which models a
bone marrow tissue cell differentiation from common myeloid cell to
specialised blood cells (megakaryocytes, erythrocytes, granulocytes and
monocytes). The second model depicts a heart development done by _cardiac_
kernel transcription factors [53]. The third model explains regulatory
interactions, which may cause _ERBB_ kinase over-expression – a marker of
breast cancer [54]. The fourth model focuses on specific conditions which lead
to a metastatic _tumour_ [55]. The last model, the Mitogen-Activated Protein
Kinase (_MAPK_) network, consists of signalling pathways involved in diverse
cellular processes, such as cell cycle, survival, apoptosis and
differentiation [24]. We derived partially specified versions of these models
by substituting a portion of the network’s update functions with uninterpreted
functions of the same arity. This emulates the model development process,
during which the general structure of the system is already determined, but
the exact network dynamics are not yet fully established. All studied
partially specified models can be found in the Github repository.

We have selected two or three attractors from every fully specified model
version obtained from the Cell Collective platform [56]. We then use partially
specified versions of these models. All algorithms are computed for all pairs
of source and target attractors. For every model, using the introduced
symbolic algorithms, we compute all types of control – one-step, permanent,
and temporary. For every computed attractor, we also state how big is the
relevant colour space, as it is expected to have an influence on the
computation time. The results are grouped by the target attractor, as the
times for these results are similar given the nature of the algorithm (the
most demanding part – computation of target trap space depends solely on the
choice of the target). The results can be found in Table 1. This experiment
demonstrates that our framework can handle even relatively large models (20
and more variables) in a convenient range of times.

Table 1: Performance of algorithms achieved on the biological models. Sub-
tables correspond to particular models. The sub-table header displays the name
of the model, reference to its source, and factors influencing the size of the
underlying graph – the number of states (induced by the number of the model
variables), the approximate number of model colours, and employed types of
control. The results are computed for some two or three arbitrary attractors.
Therefore, for every target attractor, we list the number of relevant colours
where the attractor is present and the average of times needed to compute one-
step, permanent and temporary control from the source attractors to the said
target attractor. Computation times marked as N/A did not finish in 4 hours.
Myeloid [52] – $2^{11}$ states, $\sim 2^{17}$ colours, $2^{11}$ perturbations
---
Control target $t$ |  Attract. colours |  One-step |  Permanent |  Temporary
$A_{1}$ |  $\sim 2^{10}$ |  8 ms |  15 ms |  10 ms
$A_{2}$ |  $\sim 2^{14}$ |  38 ms |  107 ms |  128 ms
$A_{3}$ |  $\sim 2^{16}$ |  62 ms |  124 ms |  607 ms
Cardiac [53] – $2^{15}$ states, $\sim 2^{28}$ colours, $2^{15}$ perturbations
Control target $t$ |  Attract. colours |  One-step |  Permanent |  Temporary
$A_{1}$ |  $\sim 2^{27}$ |  174 ms |  350 ms |  2 s
$A_{2}$ |  $\sim 2^{27}$ |  520 ms |  2058 ms |  13 s
$A_{3}$ |  $\sim 2^{27}$ |  453 ms |  1050 ms |  77 s
ERBB [54] – $2^{20}$ states, $\sim 2^{28}$ colours, $2^{20}$ perturbations
Control target $t$ |  Attract. colours |  One-step |  Permanent |  Temporary
$A_{1}$ |  $\sim 2^{27}$ |  905 ms |  35 s |  8.2 min
$A_{2}$ |  $\sim 2^{28}$ |  7141 ms |  40 s |  2 hrs
Tumour [55] – $2^{32}$ states, $\sim 2^{16}$ colours, $2^{32}$ perturbations
Control target $t$ |  Attract. colours |  One-step |  Permanent |  Temporary
$A_{1}$ |  $\sim 2^{14}$ |  1049 ms |  35 s |  8 min
$A_{2}$ |  $\sim 2^{14}$ |  887 ms |  56 s |  15 min
$A_{3}$ |  $\sim 2^{14}$ |  947 ms |  7.5 s |  20 min
MAPK [24] – $2^{53}$ states, $\sim 2^{16}$ colours, $2^{53}$ perturbations
Control target $t$ |  Attract. colours |  One-step |  Permanent |  Temporary
$A_{1}$ |  $\sim 2^{15}$ |  6.4 s |  3 hrs |  N/A
$A_{2}$ |  $\sim 2^{15}$ |  5.5 s |  1.2 hrs |  N/A

The second demonstration inquires about the scalability of our algorithms with
respect to the number of colours. The results are shown in Table 2 for the
tumour model [55]. It can be seen that the times needed to compute the results
grow much slower than the number of model colours. In particular, a jump from
$2^{14}$ to $2^{21}$ colours, which represents a $128\times$ increase in the
size of the colour space, results in an approx. $3-6\times$ increase in
runtime.

Table 2: Algorithms scalability with respect to colours demonstrated on the tumour model [55]. We selected one arbitrary target attractor in the model for all following computations. In the header, the number of colours containing an attractor with the computed target is listed. The shown time is the worst case of source-target control in all considered source attractors. Rel. colours |  $\mathbf{2^{0}}$ |  $\mathbf{\sim 2^{4}}$ |  $\mathbf{\sim 2^{10}}$ |  $\mathbf{\sim 2^{14}}$ |  $\mathbf{\sim 2^{21}}$
---|---|---|---|---|---
One-step |  120 ms |  170 ms |  546 ms |  1053 ms |  3324 ms
Permanent |  3.2 s |  13.5 s |  14.2 s |  57.2 s |  262 s
Temporary |  0.5 min |  1.3 min |  4.4 min |  20.4 min |  130 min

Next, we compare our one-step approach to the semi-symbolic one-step control
method described in [46]. In the semi-symbolic method, only the BN inputs were
specified symbolically, using BDDs. In contrast, using our fully symbolic
approach presented here, the perturbations, as well as the sets of states, are
represented symbolically, all by using BDDs. In Table 3, we show a comparison
of these two techniques.

Table 3: Results of one-step control – comparison of the fully symbolic approach to semi-symbolic approach presented in [46]. The values are stated as ranges because we computed one step-control on a set of attractors. The first column contains name of model and number of its states. The second column shows the number of model’s colours. The third column shows count ranges of colours which contain the given target attractor. The fourth and fifth columns display ranges of times to compute one-step controls (strong basins) using semi-symbolic respectively symbolic methods. The computation time marked as N/A was not feasible to compute due to insufficient memory. Model | Colours | Attract. colours | Semi-symbolic | Symbolic
---|---|---|---|---
Cell-Fate [57] $2^{19}$ states | 1 | 1 | 4.4 – 9.19 s | 5 – 21 ms
$\sim 2^{10}$ | 1 – 4 | 4.63 – 13.42 s | 103 – 317 ms
$\sim 2^{24}$ | 7 – 56 | 5.82 – 26.59 s | 3 – 17.4 s
Myeloid [52] $2^{11}$ states | 1 | 1 | 8 – 30 ms |  $\leq$ 1 ms
$\sim 2^{26}$ | 63 – 2,052 | 14 – 214 ms | 43 – 460ms
$\sim 2^{49}$ |  $\sim 2^{16}$ – $\sim 2^{24}$ | 147 – 1717 ms | 1.4 – 14 min
| MAPK [24]
---
$2^{53}$ states
$\sim 2^{16}$ | $\sim 2^{15}$ | N/A | 5.5 – 6.4 s
Figure 1: Maximal robustness $\rho_{max}$ per set of perturbations with a
given size – comparison between different types of controls. The $x$-axis
shows the size of perturbations in the set and the $y$-axis displays the
maximal robustness. Each chart is showing robustness for a different model.

$1$$2$$3$$4$$5$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$SizeMaximal
robustnessa) Cardiac (15 variables)One-stepPermanentTemporary

$1$$2$$3$$4$$5$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$SizeMaximal
robustnessb) ERBB (20 variables)One-stepPermanentTemporary

$1$$2$$3$$4$$5$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$SizeMaximal
robustnessc) Tumour (32 variables)One-stepPermanentTemporary

Figure 2: Union robustness $\rho_{max}$ per set of perturbations with a given
size – comparison between different types of controls. The $x$-axis shows the
size of perturbations in the set and the $y$-axis displays the maximal
robustness. Each chart is showing robustness for a different model.

$1$$2$$3$$4$$5$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$SizeUnion
robustnessa) Cardiac (15 variables)One-stepPermanentTemporary

$1$$2$$3$$4$$5$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$SizeUnion
robustnessb) ERBB (20 variables)One-stepPermanentTemporary

$1$$2$$3$$4$$5$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$SizeUnion
robustnessc) Tumour (32 variables)One-stepPermanentTemporary

Finally, we conduct experiments evaluating the presented robustness metrics
across different types of controls. In particular, we use maximal robustness
$\rho_{max}$ and union robustness $\rho_{union}$ metrics which we measure on
sets of fixed perturbation sizes. The results are shown in Fig. 1 and Fig. 2.
The models are the same for both experiments and the same as in our first
experiment (see Table 1). For demonstration purposes, an arbitrary but fixed
source and target attractor have been selected among the ones considered in
Table 1.

##  5 Discussion

In the presented work, we show how to employ symbolic methods to solve the
control of partially specified BNs. The symbolic representation is used to
encode three types of entities: a set of the STG vertices, a set of colours
(input valuations) which admit an edge in the STG, and a set of perturbations
that permit the given edge. When we increase the number of colours (unknown
parts of the system), only one dimension of the resulting symbolic
representation grows. In opposition, when we increase the number of system’s
variables, both the size of the STG, and the available perturbations grow.
Therefore, it is expected that increasing the model size has a bigger impact
on the performance than increasing the number of model’s parameters. This
hypothesis is to some extent supported by our two experiments captured in
Table 1 and Table 2. We can observe that when the colour space increases in
isolation, the model still remains processable and the computation time grows
much slower than the absolute number of colours. On the other hand, increasing
the size of the state space can result in analysis being infeasible (see the
MAPK case).

Nonetheless, it must be highlighted that in both of these cases, the impact of
scaling up colours versus variables is to some extent unpredictable. This is
given by the intrinsic mechanics of the introduced STGs. As the nature of our
method relies on fixed-point computation, such that each iteration eliminates
only the direct predecessors of a particular set, in cases where e.g. a long
path in the graph arises, our method will not be able to scale well: the
algorithm will need many iterations to converge, each eliminating only a small
fraction of the graph vertices. Another unforeseeable effect on the method’s
performance have the BDDs themselves. The ordering of the BDD diagram has a
significant impact on the effectiveness of all used operations and no single
ordering strategy can perform best in all cases [58].

Another insight that can be obtained from Table 1 is how one-step, temporary
and permanent controls differ in their complexity in terms of computational
time. We can notice that one-step control is significantly over-performing
other methods. This is because, in the case of one-step control, we do not
need to consider perturbation state-space when computing the trap space of the
target attractor. In contrast, the temporary control is by far the most
computationally demanding compared to the other types control. This is again
given by the fact that we need to compute two distinct trap spaces. On the one
hand, the first computed trap space is the same as in the case of one-step
control and should be thus easy to compute (the perturbation space is
irrelevant). However, on the other hand, the seed for the second computed trap
space is a much bigger vertex set, and therefore the resulting second trap
space (which includes perturbations) of the first trap space (without
perturbations) can generate a much more complex symbolic structure.

Apart from the time-wise performance, we also compared the robustness of
different types of control (see Figure 1 and Figure 2). At first, let us look
at what particular types of robustness can tell us about the inspected
systems. Under the assumption that the real biological system corresponds to
one of the instantiations of our partially specified BN model, what the
maximal robustness tells us is how high is the chance that _a particular
perturbation_ will work in the real model. Alternatively, we use union
robustness to evaluate a chance that there is _some perturbation_ of the given
size that will work. For example, in the ERBB model (Figure 1b and Figure 2b),
we can notice that there is a 60% chance that the most robust perturbation
(temporary or permanent) of size two will work for the real, fully specified
model. However, even if this particular perturbation fails, we know that 80%
of instantiations can be influenced by a perturbation of this size. This means
that out of the 40% of instantiations that cannot be influenced by the
perturbation with the maximal robustness, a half is still controllable by
perturbing only two variables. For all three models, it is thus highly likely
that a temporary perturbation of size no bigger than 3 exists which controls
the network. That is a very small portion of the total model variables which
could be possibly perturbed.

As the one-step control is a special case of the temporary control, it is
expected that temporary perturbations are achieving the same or higher
robustness than the one-step case. This is also demonstrated in the
experiments, where we find no instance in which the one-step control could
achieve a higher maximal or union robustness. We can also observe that
temporary and permanent controls mostly perform similarly well. This is in
line with [36] where the minimal temporary and permanent controls also
performed the same in the most of the cases. In our case, only in the cardiac
model, the temporary control outperforms the permanent one in the union
robustness (Figure 2a). This means there are some perturbations which work
when applied temporarily as opposed when applied permanently. There are two
admissible explanations for this. The first may be a simple fact that the trap
space computed for temporary control is expected to be bigger, since the seed
of the trap space is not only the given target attractor, but its trap space
in an unperturbed case. The second explanations for this is a scenario when a
permanent perturbation disrupts the BN dynamics in such a way that the target
attractor completely ceases to exist.

Another insight of our experiments is a comparison of the novel symbolic
approach to the semi-symbolic approach from [46]. The experiment captured in
Table 3 shows that the symbolic approach is performing better on larger models
with fewer colours. Nonetheless, if the model is small but has many colours,
the semi-symbolic approach can outperform the fully symbolic one. This is
understandable due to the already discussed heuristic nature of BDDs: in some
cases, the ordering of BDD nodes might not be optimal, resulting in large
underlying memory consumption and longer computation. Nonetheless, the
greatest benefit of the novel symbolic approach lies in the ability to process
models with a much bigger state space. This is the case of the MAPK model,
where the explicit state-space of the semi-symbolic approach completely fails
to fit into the memory.

Finally, we compare our results with the related work discussed in Section 1.
It is important to note that the previously published methods other than [46]
do not support partially specified BNs. Therefore, we deem our method
incomparable to the previous line of research in terms of performance.
Hypothetically, the other methods could be used to iteratively process all
instantiations of a partially specified BN, resulting in a naïve parameter-
scan approach. Nonetheless, as the number of instantiations grows
exponentially, this approach cannot scale in the same way and thus quickly
become infeasible. More than that, our approach allows easy evaluation of
robustness, which would otherwise also be more complicated when the results
are computed naïvely. Considering other aspects of BN control, our method
supports all types of static perturbations (one-step, temporary, permanent)
[41]. However, it does not support more complex, dynamical application of
sequential perturbations [39, 41]. With that type of perturbations it is
possible to achieve the minimal perturbation of a very small size. However, we
consider such approach too computationally demanding for the case of partially
specified BNs at the moment – the combinatorial space of available
perturbations, colours, and attractors (each present only for a subset of
colours) would be enormous. Finally, our method supports only source-target
control. We leave other types of control, e.g., the target control [37] and
the phenotype control [44] for future work.

##  6 Conclusion

In this work, we have developed methods and efficient algorithms to solve the
source-target control problem for partially specified BNs using one-step,
temporary and permanent state perturbations. The methods we proposed are
entirely symbolic and are based on BDDs. We demonstrated that our approach is
capable of controlling models with a high degree of uncertainty in a
conveniently short time. Such a result cannot be achieved with the naïve scan
over all BN inputs due to the doubly exponential explosion of the number of
the inputs representing the fully specified BN instances. We also showed that
our novel methods are superior to the existing semi-symbolic one-step control
approach [46] in terms of model size we can process as well as the robustness
of temporary and permanent perturbations. We were able to evaluate robustness
between different types of control thanks to the metrics we introduced –
maximal and union robustness.

In future work, we would like to improve the algorithms for computing the
strong basin, e.g., by employing suitable network decomposition procedures
[59, 60]. We would also like to allow target control or only a partial
specification of the state space into which the system is controlled (i.e.
trap spaces or phenotypes). Another possible extension to all algorithms would
be to implement a restriction on the control size or on variables that may be
perturbed. Next, we aim to explore practical applications of the developed
methods further. If these inquiries show that permanent perturbations are
practically applicable, we might study further this type of perturbations and
its potential application for phenotypes variance. Last but not least, we plan
to incorporate the control methods into our toolkit for partially specified
BNs – the AEON tool [51].

## References

  * [1] Hiroaki Kitano.  Computational systems biology.  Nature, 420(6912):206–210, 2002.
  * [2] Sarah S Richardson and Hallam Stevens.  Postgenomics: Perspectives on biology after the genome.  Duke University Press, 2015.
  * [3] Enrico Borriello and Bryan Daniels.  The basis of easy controllability in boolean networks.  Nature Communications, 12:5227, 09 2021.
  * [4] Anne B. C. Cherry and George Q. Daley.  Reprogramming cellular identity for regenerative medicine.  Cell, 148(6):1110–1122, 2012.
  * [5] Natasa Miskov-Zivanov, Michael S. Turner, Lawrence P. Kane, Penelope A. Morel, and James R. Faeder.  The duration of T cell stimulation is a critical determinant of cell fate and plasticity.  Science Signaling, 6(300), 2013.
  * [6] Assieh Saadatpour, Rui-Sheng Wang, Aijun Liao, Xin Liu, Thomas P. Loughran, István Albert, and Réka Albert.  Dynamical and structural analysis of a t cell survival network identifies novel candidate therapeutic targets for large granular lymphocyte leukemia.  PLOS Computational Biology, 7(11):1–15, 2011.
  * [7] Jorge G. T. Zañudo and Réka Albert.  Cell fate reprogramming by control of intracellular network dynamics.  PLOS Computational Biology, 11(4):1–24, 2015.
  * [8] Carla A. Herberts, Marcel SG Kwa, and Harm PH Hermsen.  Risk factors in the development of stem cell therapy.  Journal of Translational Medicine, 9(1), 2011.
  * [9] Amar M. Singh and Stephen Dalton.  The cell cycle and myc intersect with mechanisms that regulate pluripotency and reprogramming.  Cell Stem Cell, 5(2):141–149, 2009.
  * [10] Kazutoshi Takahashi and Shinya Yamanaka.  Induction of pluripotent stem cells from mouse embryonic and adult fibroblast cultures by defined factors.  Cell, 126(4):663–676, 2006.
  * [11] Richard A. Young.  Control of the embryonic stem cell state.  Cell, 144(6):940–954, 2011.
  * [12] Marius Wernig, Alexander Meissner, Ruth Foreman, Tobias Brambrink, Manching Ku, Konrad Hochedlinger, Bradley E. Bernstein, and Rudolf Jaenisch.  In vitro reprogramming of fibroblasts into a pluripotent ES-cell-like state.  Nature, 448(7151):318–324, 2007.
  * [13] Yanbei Zhu, Rongze Yang, John McLenithan, Daozhan Yu, Hong Wang, Yaping Wang, Devinder Singh, John Olson, Carole Sztalryd, Dalong Zhu, et al.  Direct conversion of human myoblasts into brown-like adipocytes by engineered super-active PPAR$\gamma$.  Obesity, 23(5):1014–1021, 2015.
  * [14] Adilson E Motter, Natali Gulbahce, Eivind Almaas, and Albert-László Barabási.  Predicting synthetic rescues in metabolic networks.  Molecular Systems Biology, 4(1), 2008.
  * [15] Michael S Goligorsky.  New trends in regenerative medicine: reprogramming and reconditioning.  Journal of the American Society of Nephrology, 30(11):2047–2051, 2019.
  * [16] Eric Davidson and Michael Levin.  Gene regulatory networks.  Proceedings of the National Academy of Sciences, 102(14), 2005.
  * [17] Michael Hecker, Sandro Lambeck, Susanne Toepfer, Eugene van Someren, and Reinhard Guthke.  Gene regulatory network inference: Data integration in dynamic models—a review.  Biosystems, 96(1):86–103, 2009.
  * [18] Roberto Barbuti, Roberta Gori, Paolo Milazzo, and Lucia Nasti.  A survey of gene regulatory networks modelling methods: from differential equations, to boolean and qualitative bioinspired models.  Journal of Membrane Computing, pages 1–20, 2020.
  * [19] Réka Albert.  Boolean modelingof genetic regulatory networks.  In Complex networks, pages 459–481. Springer, 2004.
  * [20] René Thomas, Denis Thieffry, and Marcelle Kaufman.  Dynamical behaviour of biological regulatory networks—i. biological role of feedback loops and practical use of the concept of the loop-characteristic state.  Bulletin of mathematical biology, 57(2):247–276, 1995.
  * [21] Julian D. Schwab, Silke D. Kühlwein, Nensi Ikonomi, Michael Kühl, and Hans A. Kestler.  Concepts in Boolean network modeling: What do they all mean?  Computational and Structural Biotechnology Journal, 18:571–582, 2020.
  * [22] Julio Aracena, Eric Goles, Andrés Moreira, and Luis Salinas.  On the robustness of update schedules in boolean networks.  Biosystems, 97(1):1–8, 2009.
  * [23] Natalie Berestovsky and Luay Nakhleh.  An evaluation of methods for inferring boolean networks from time-series data.  PloS one, 8(6):e66031, 2013.
  * [24] Luca Grieco, Laurence Calzone, Isabelle Bernard-Pierrot, François Radvanyi, Brigitte Kahn-Perles, and Denis Thieffry.  Integrative modelling of the influence of mapk network on cancer cell fate decision.  PLoS Comput Biol, 9(10):e1003286, 2013.
  * [25] Liesbet Geris and David Gomez-Cabrero.  An Introduction to Uncertainty in the Development of Computational Models of Biological Processes, pages 3–11.  Springer, 2016.  Generated from Scopus record by KAUST IRTS on 2021-02-16.
  * [26] Alberto JM Martin, Calixto Dominguez, Sebastián Contreras-Riquelme, David S Holmes, and Tomas Perez-Acle.  Graphlet based metrics for the comparison of gene regulatory networks.  PloS one, 11(10), 2016.
  * [27] Alexis Baudin, Soumya Paul, Cui Su, and Jun Pang.  Controlling large Boolean networks with single-step perturbations.  Bioinformatics, 35(14):i558–i567, 2019.
  * [28] Junil Kim, Sang-Min Park, and Kwang-Hyun Cho.  Discovery of a kernel for controlling biomolecular regulatory networks.  Nature Scientific Reports, 3(2223), 2013.
  * [29] Bernold Fiedler et al.  Dynamics and control at feedback vertex sets. i: Informative and determining nodes in regulatory networks.  J. Dyn. Differ. Equ., 25(3):563–604, 2013.
  * [30] Yin Zhao, Bijoy K Ghosh, and Daizhan Cheng.  Control of large-scale Boolean networks via network aggregation.  IEEE Trans. on Neural Networks and Learning Systems, 27(7):1527–1536, 2015.
  * [31] Mohammad Moradi, Sama Goliaei, and Mohammad-Hadi Foroughmand-Araabi.  A Boolean network control algorithm guided by forward dynamic programming.  PLOS ONE, 14:1–21, 2019.
  * [32] Jordan C Rozum, Jorge Gómez Tejeda Zañudo, Xiao Gan, Dávid Deritei, and Réka Albert.  Parity and time reversal elucidate both decision-making in empirical models and attractor scaling in critical boolean networks.  Science Advances, 7(29), 2021.
  * [33] Jordan C Rozum, Dávid Deritei, Kyu Hyong Park, Jorge Gómez Tejeda Zañudo, and Réka Albert.  pystablemotifs: Python library for attractor identification and control in Boolean networks.  Bioinformatics, 12 2021.
  * [34] Sean P Cornelius, William L Kath, and Adilson E Motter.  Realistic control of network dynamics.  Nature Communications, 4(1):1–9, 2013.
  * [35] Hugues Mandon, Stefan Haar, and Loïc Paulevé.  Temporal reprogramming of boolean networks.  In CMSB, volume 10545 of LNCS, pages 179–195. Springer, 2017\.
  * [36] Cui Su, Soumya Paul, and Jun Pang.  Controlling large Boolean networks with temporary and permanent perturbations.  In International Symposium on Formal Methods, pages 707–724. Springer, 2019.
  * [37] Cui Su and Jun Pang.  A dynamics-based approach for the target control of boolean networks.  CoRR, abs/2006.02304, 2020.
  * [38] Wassim Abou-Jaoudé, Pauline Traynard, Pedro T. Monteiro, Julio Saez-Rodriguez, Tomáš Helikar, Denis Thieffry, and Claudine Chaouiya.  Logical modeling and dynamical analysis of cellular networks.  Frontiers in Genetics, 7:94, 2016.
  * [39] Hugues Mandon, Cui Su, Stefan Haar, Jun Pang, and Loïc Paulevé.  Sequential reprogramming of Boolean networks made practical.  In Computational Methods in Systems Biology, pages 3–19. Springer, 2019.
  * [40] Jérémie Pardo, Sergiu Ivanov, and Franck Delaplace.  Sequential reprogramming of biological network fate.  In CMSB, pages 20–41. Springer, 2019.
  * [41] Cui Su and Jun Pang.  Sequential temporary and permanent control of Boolean networks.  In Computational Methods in Systems Biology, pages 234–251. Springer, 2020.
  * [42] Cui Su and Jun Pang.  CABEAN: a software for the control of asynchronous Boolean networks.  Bioinformatics, 2020.
  * [43] Laura Cifuentes Fontanals, Elisa Tonello, and Heike Siebert.  Control strategy identification via trap spaces in Boolean networks.  In Computational Methods in Systems Biology, pages 159–175. Springer, 2020.
  * [44] Laura Cifuentes-Fontanals, Elisa Tonello, and Heike Siebert.  Control in boolean networks with model checking.  arXiv preprint arXiv:2112.10477, 2021.
  * [45] Rui-Sheng Wang, Assieh Saadatpour, and Reka Albert.  Boolean modeling in systems biology: an overview of methodology and applications.  Physical biology, 9(5), 2012.
  * [46] Luboš Brim, Samuel Pastva, David Šafránek, and Eva Šmijáková.  Parallel one-step control of parametrised boolean networks.  Mathematics, 9(5), 2021.
  * [47] Randal E. Bryant.  Graph-based algorithms for boolean function manipulation.  IEEE Trans. Computers, 35(8):677–691, 1986.
  * [48] Nikola Beneš, Luboš Brim, Samuel Pastva, Jakub Poláček, and David Šafránek.  Formal analysis of qualitative long-term behaviour in parametrised Boolean networks.  In Formal Methods and Software Engineering, pages 353–369. Springer, 2019.
  * [49] Martin Abadi and Leslie Lamport.  Conjoining specifications.  ACM Transactions on Programming Languages and Systems, 17(3):507–534, 1995.
  * [50] Anders Forsman.  Effects of genotypic and phenotypic variation on establishment are important for conservation, invasion, and infection biology.  Proceedings of the National Academy of Sciences, 111(1):302–307, 2014.
  * [51] Nikola Beneš, Luboš Brim, Samuel Pastva, and David Šafránek.  AEON: attractor bifurcation analysis of parametrised boolean networks.  In CAV 2020, volume 12224 of Lecture Notes in Computer Science, Cham, 2020. Springer International Publishing.
  * [52] Jan Krumsiek, Carsten Marr, Timm Schroeder, and Fabian J. Theis.  Hierarchical differentiation of myeloid progenitors is encoded in the transcription factor network.  PLOS ONE, 6(8):1–10, 2011.
  * [53] Ashley J Waardenberg, Mirana Ramialison, Romaric Bouveret, and Richard P Harvey.  Genetic networks governing heart development.  Cold Spring Harbor perspectives in medicine, 4(11):a013839, 2014\.
  * [54] Özgür Sahin, Holger Fröhlich, Christian Löbke, Ulrike Korf, Sara Burmester, Meher Majety, Jens Mattern, Ingo Schupp, Claudine Chaouiya, Denis Thieffry, et al.  Modeling erbb receptor-regulated g1/s transition to find novel targets for de novo trastuzumab resistance.  BMC systems biology, 3(1):1–20, 2009.
  * [55] David PA Cohen, Loredana Martignetti, Sylvie Robine, Emmanuel Barillot, Andrei Zinovyev, and Laurence Calzone.  Mathematical modelling of molecular pathways enabling tumour cell invasion and migration.  PLoS Comput Biol, 11(11):e1004571, 2015.
  * [56] Tomáš Helikar et al.  The cell collective: toward an open and collaborative approach to systems biology.  BMC systems biology, 6(1):1–14, 2012.
  * [57] Laurence Calzone, Laurent Tournier, Simon Fourquet, Denis Thieffry, Boris Zhivotovsky, Emmanuel Barillot, and Andrei Zinovyev.  Mathematical modelling of cell-fate decision in response to death receptor engagement.  PLOS Computational Biology, 6(3):1–15, 2010.
  * [58] Justin E Harlow III and Franc Brglez.  Design of experiments and evaluation of bdd ordering heuristics.  International Journal on Software Tools for Technology Transfer, 3(2):193–206, 2001.
  * [59] Andrzej Mizera, Jun Pang, Hongyang Qu, and Qixia Yuan.  Taming asynchrony for attractor detection in large boolean networks.  IEEE/ACM transactions on computational biology and bioinformatics, 16(1):31–42, 2018.
  * [60] Cui Su, Jun Pang, and Soumya Paul.  Towards optimal decomposition of boolean networks.  IEEE/ACM Transactions on Computational Biology and Bioinformatics, 2019.