# Deep Reinforcement Learning for Controlled Traversing of the Attractor
Landscape of Boolean Models in the Context of Cellular Reprogramming

Andrzej Mizera  University of Warsaw  IDEAS NCBR  Jakub Zarzycki  University
of Warsaw  IDEAS NCBR

###### Abstract

Cellular reprogramming can be used for both the prevention and cure of
different diseases. However, the efficiency of discovering reprogramming
strategies with classical wet-lab experiments is hindered by lengthy time
commitments and high costs. In this study, we develop a novel computational
framework based on deep reinforcement learning that facilitates the
identification of reprogramming strategies. For this aim, we formulate a
control problem in the context of cellular reprogramming for the frameworks of
BNs and PBNs under the asynchronous update mode. Furthermore, we introduce the
notion of a pseudo-attractor and a procedure for identification of pseudo-
attractor state during training. Finally, we devise a computational framework
for solving the control problem, which we test on a number of different
models.

##  1 Introduction

Complex diseases pose a great challenge largely because genes and gene
products operate within a complex system – the gene regulatory network (GRN).
There is an inherent dynamic behaviour emerging from the structural wiring of
a GRN: gene expression profiles, i.e., states of a GRN, evolve in time to
finally reach stable states referred to as attractors. Attractors correspond
to cell types or cell fates [13]. During normal development of a multi-
cellular organism, not all attractors are manifested. Some of the ‘abnormal
attractors’, associated with diseases, become accessible by disturbance of the
GRN’s dynamics. This is seldom a consequence of a disruption in a single gene,
but rather arises as an aftermath of GRN perturbations [3]. This could be
cured by guiding cells to desired ‘healthy’ attractors with experimental
techniques of cellular reprogramming, i.e., the artificial changing of cell
fate. Unfortunately, finding effective interventions that trigger desired
changes using solely wet-lab experiments is difficult, costly, and requires
lengthy time commitments. This motivates us to consider _in-silico_
approaches.

Although various computational frameworks are commonly used to model GRNs, the
formalism of Boolean networks (BNs) and its extension, i.e., probabilistic
Boolean networks (PBNs), have the advantage of being simple yet capable of
capturing the important dynamic properties of the system under study. As such,
they facilitate the modelling of large biological systems. This is especially
relevant in the context of _in-silico_ identification of effective cellular
reprogramming strategies, which requires large GRNs to be modelled.

Identification of cellular reprogramming strategies can be stated as a control
problem of BN and PBN models of GRNs. Although many BN/PBN control methods
exist in the literature, the existing structure- and dynamics-based state-of-
the-art computational techniques are limited to small and mid-size networks,
i.e., of up to a hundred of genes or so, usually requiring the systems to be
decomposed in some way. This is often insufficient for cellular reprogramming
considerations.

The issue of scalability can be addressed by devising new methods based on
deep reinforcement learning (DRL) techniques, which have proved very
successful in decision problems characterised by huge state-action spaces. To
contribute to the realisation of this idea, we formulate a control problem in
the context of cellular reprogramming for the frameworks of BNs and PBNs under
the asynchronous update mode. Furthermore, we introduce the notion of a
pseudo-attractor and a procedure for identifying pseudo-attractor states
during DRL agent training. Finally, these contributions allow us to devise a
DRL-based framework for solving the control problem. We consider our
contributions as a relevant step towards achieving scalable control methods
for large Boolean models of GRNs for identifying effective and efficient
cellular reprogramming strategies.

The paper is structured as follows. Related work is discussed in Sec. 2.
Preliminaries are provided in Sec. 3. We formulate our control problem in the
context of cellular reprogramming in Sec. 4 and devise our DRL-based control
framework in Sec. 5. The experiments performed to evaluate our framework and
the obtained results are presented in Sec.6 and Sec. 7, respectively. Finally,
we conclude our study in Sec. 8.

##  2 Related work

###  2.1 Dynamics-based approaches to GRN control

Identification of proper control strategies for non-linear systems requires
both network structure and dynamics [11]. Thus, we focus on dynamic-based and
DRL-based methods for BN/PBN control. An efficient method based on the ‘divide
and conquer’ strategy was proposed in [18] to solve the minimal one-step
source-target control problem by using instantaneous, temporary, and permanent
gene perturbations. The minimal sequential source-target control and the
target control problems of BNs were considered in [23] and [22], respectively.
All these methods were implemented as a software tool CABEAN [24]. Recently,
semi-symbolic algorithms were proposed in [4] to stabilise partially specified
asynchronous BNs in states exhibiting specific traits. In [19], the control
problem for the most permissive BN update mode in the context of fixed points
and minimal trap spaces is considered.

###  2.2 DRL-based approches to GRN control

The application of reinforcement learning for controlling GRNs was pioneered
by in [21] with focus on how to control GRNs by avoiding undesirable states in
terms of steady state probabilities of PBNs. The main idea was to treat the
time series gene expression samples as a sequence of experience tuples and use
a batch version of Q-Learning to produce an approximated policy over the
experience tuples. Later, the BOAFQI-Sarsa method that does not require time
series samples was devised in [16]. A batch reinforcement learning method,
mSFQI, was proposed in [15] for control based on probabilities of gene
activity profiles. Recently, the study of [1] used a Deep Q-Network with
prioritised experience replay, for control of synchronous PBNs to drive the
networks from a particular state towards a more desirable one. Finally, a DRL-
based approximate solution to the control problem in synchronous PBNs was
proposed in [14]. The proposed method finds a control strategy from any
network state to a specified target attractor using a Double Deep Q-Network
model.

##  3 Preliminaries

###  3.1 Boolean and probabilistic Boolean networks

Boolean networks is a well established framework for the modelling of GRNs. A
Boolean network consists of nodes, that can be in one of two states, and
functions describing how the individual nodes interact with each other. PBNs
are an extension of the formalism of BNs.

######  Definition 3.1.

(Boolean Network) A  _Boolean Network_ is defined as a pair $(V,F)$, where
$V=\\{x_{1},x_{2},\ldots,x_{n}\\}$ is a set of binary-valued nodes (also
referred to as genes) and $F=\\{f_{1},f_{2},\ldots,f_{n}\\}$ is a set of
Boolean predictor functions, where
$f_{i}(x_{i_{1}},x_{i_{2}},\ldots,x_{i_{k}})$ defines the value of node
$x_{i}$ depending on the values of the $k\leq n$ parent nodes
$x_{i_{1}},x_{i_{2}},\ldots,x_{i_{k}}$ with $i_{j}\in[1..n]$ for $j\leq k$.

Since interactions in biology are usually more complex we need a more general
model of a GRN. We achieve it by allowing for each node to have multiple
Boolean functions. Formally, probabilistic Boolean networks are defined as
follows:

######  Definition 3.2.

(Probabilistic Boolean Network) A  _probabilistic Boolean network_ is defined
as a pair $(V,\mathcal{F})$, where $V=\\{x_{1},x_{2},\ldots,x_{n}\\}$ is a set
of binary-valued nodes (also referred to as genes) and
$\mathcal{F}=(F_{1},F_{2},\ldots,F_{n})$ is a list of sets. Each node
$x_{i}\in V$, $i=1,2,\ldots,n$, has associated a set $F_{i}\in\mathcal{F}$ of
Boolean predictor functions:
$F_{i}=\\{f^{i}_{1},f^{i}_{2},\ldots,f^{i}_{l(i)}\\}$, where $l(i)$ is the
number of predictor functions of node $x_{i}$. Each $f^{i}_{j}\in F_{i}$ is a
Boolean function defined with respect to a subset of $V$ referred to as parent
nodes for $f^{i}_{j}$ and denoted $\textrm{Pa}(f^{i}_{j})$. For each node
$x_{i}\in V$ there is a probability distribution
$\mathbf{c}^{i}=(c^{i}_{1},c^{i}_{2},\ldots,c^{i}_{l(i)})$ on $F_{i}$, where
each predictor function $f^{i}_{j}\in F_{i}$ has an associated selection
probability denoted $c^{i}_{j}$; it holds that $\sum_{j=1}^{l(i)}c^{i}_{j}=1$.

A PBN in which each node only admits only one Boolean function is a  _Boolean
network_.

###  3.2 Network dynamics

We define a  _state_ of a BN/PBN as an $n$-dimensional vector
$\mathbf{s}\in\\{0,1\\}^{n}$, where the $i$-th element represents the state of
gene $x_{i}$ for $i\in[1..n]$. A BN/PBN evolves in discrete time steps. It
starts in an initial state $\mathbf{s}_{0}$ and its state gets updated in
every time step in accordance to the predictor functions. In this study, we
focus on the asynchronous updating, which is preferable in the context of GRN
modelling. Under the asynchronous scheme, a single gene $x_{i}$ is selected
and updated in accordance with its predictor function $f_{i}$ (BNs) or one
randomly selected from $F_{i}$ in accordance with $\mathbf{c}^{i}$. The
network dynamics can be depicted in the form of a  _state transition graph_.
Based on this concept of, we can introduce the notion of a BN/PBN attractor.

######  Definition 3.3.

(State Transition Graph (STG)) A state transition graph of a BN/PBN of $n$
genes under the asynchronous update mode is a graph $G(S,\rightarrow)$, where
$S=\\{0,1\\}^{n}$ is the set of all possible states and $\rightarrow$ is the
set of directed edges such that a directed edge from $s$ to $s^{\prime}$,
denoted $s\rightarrow s^{\prime}$, is in $\rightarrow$ if and only if
$s^{\prime}$ can be obtained from $s$ by a single asynchronous update.

######  Definition 3.4.

(Attractor) An  _attractor_ of a BN/PBN is a bottom strongly connected
component in the STG of the network. A  _fixed-point attractor_ and a  _multi-
state attractor_ are bottom strongly connected components consisting of a
single state or more than one state, respectively.

######  Example 3.5.

We consider a PBN of 4 genes $V=\\{x_{0},x_{1},x_{2},x_{3}\\}$ regulated in
accordance with the following Boolean functions:

$f_{0}^{1}(x_{0})=x_{0}$
---
$f_{0}^{2}(x_{0},x_{1},x_{2},x_{3})=x_{0}\&\neg(x_{0}\&\neg x_{1}\&\neg
x_{2}\&x_{3})$
$f_{1}^{1}(x_{0},x_{1})=\neg x_{0}\&x_{1}$
$f_{1}^{2}(x_{0},x_{1},x_{2},x_{3})=\neg x_{0}\&(x_{1}|(x_{2}\&x_{3}))$
$f_{2}^{1}(x_{0},x_{1},x_{2},x_{3})=\neg x_{0}\&(x_{1}\&x_{2}\&x_{3})$
$f_{2}^{2}(x_{0},x_{1},x_{2},x_{3})=x_{0}\&(\neg x_{1}\&\neg x_{2}\&\neg
x_{3})$
$f_{3}^{1}(x_{0},x_{1},x_{2},x_{3})=\neg x_{0}\&(x_{1}|x_{2}|x_{3})$
$f_{3}^{2}(x_{0},x_{1},x_{2},x_{3})=\neg x_{0}\&(x_{1}|x_{2}|x_{3})$

Under the asynchronous update mode, the dynamics of the PBN is governed by the
STG depicted in Fig. 1.

###  3.3 Reinforcement Learning

The main task of reinforcement learning (RL) is to solve sequential decision
problems by optimising a cumulative reward function. A  _policy_ is a strategy
that determines which action to take and an  _optimal policy_ is one
determined by selecting the actions that maximise the future cumulative
reward. It be obtained by solving the _the Bellman equation_ , which expresses
the relationship between the value of a state and the expected future rewards:

$V(s)=\max_{a}[R_{a}(s,s^{\prime})+\gamma\sum_{s^{\prime}}P(s^{\prime}\mid
s,a)V(s^{\prime})],$

where $V(s)$ is the value of state $s$, $R_{a}(s,s^{\prime})$ is the immediate
reward, $P(s^{\prime}\mid s,a)$ is the transition probability to the next
state $s^{\prime}$, and $\gamma$ is the discount factor. The equation guides
the RL agent’s decision-making by considering both immediate rewards and the
discounted value of future states, forming the basis for reinforcement
learning algorithms. To find an approximate solution to the Bellman equation,
the $Q$ function is considered which is defined as the total discounted reward
received after taking action $a$ in state $s$:

$Q(s,a)=R_{a}(s,s^{\prime})+\gamma\sum_{s^{\prime}}P(s^{\prime}\mid
s,a)\max_{a^{\prime}}Q(s^{\prime},a^{\prime}).$

###  3.4 Q function approximations

In the case of large state-action spaces, the $Q$ function values often cannot
be determined, therefore they are approximated using DRL. It was shown that as
the agent explores the environment this approximation converges to the true
values of $Q$ [27]. Under the assumption that the environment is stationary,
i.e., the reward function and transition probabilities do not change in time,
one can keep evaluating the agent on new states without affecting its ability
to train itself [17].

(0,0,0,1)(0,0,1,0)(0,0,1,1)(0,1,0,0)(0,1,0,1)(1,0,0,0)(0,1,1,1)(0,1,1,0)(1,0,0,0)(1,0,1,0)(1,0,1,1)(1,1,0,0)(1,1,0,1)(1,1,1,0)(1,1,1,1)(0,0,0,0)

Figure 1: STG of the PBN defined in Example 3.5 under the asynchronous update
mode. Shaded states are the attractor states of the three attractors, i.e.,
two fixed-point attractors $A_{1}=\\{(0,0,0,0)\\}$ and
$A_{2}=\\{(0,1,0,1)\\}$, and one multi-state attractor
$A_{3}=\\{(1,0,0,0),(1,0,1,0)\\}$.

####  3.4.1 Branching Dueling Q-Network

Different DRL-based approaches can be used for the approximation of the $Q$
function. In this study we will focus on the Branching Dueling Q-Network (BDQ)
approach introduced in [25] as an extension of another well-known approach,
i.e., the Dueling Double Deep Q-Network (DDQN) [10].

BDQ deep neural network structures are designed to address complex and high-
dimensional action spaces. Instead of using a single output layer for all
actions, BDQ has multiple branches, each responsible for a specific subset of
actions. BDQ aims to enhance the scalability and sample efficiency of
reinforcement learning algorithms in complex scenarios. It can be thought of
as an adaptation of the dueling network into the action branching
architecture. The dueling architecture uses two separate artificial neural
networks, i.e., the _target network_ for evaluation and the _controller
network_ for selection of actions. Its main benefit is that it avoids
overestimating q-values, can more rapidly identify action redundancies, and
generalises more efficiently by learning a general q-value that is shared
across many similar actions.

##  4 Formulation of the control problem

###  4.1 Pseudo-attractors

Unfortunately, obtaining the attractor landscape for a large BN/PBN network,
i.e., the family of all its attractors, is a challenging problem by itself and
one cannot expect to be in possession of this information in advance. Because
our aim is to devise a scalable computational framework for the control of
large network models based on DRL, we need to be able to identify the BN/PBN
attractors during training, i.e., exploration of the DRL environment. For this
purpose, we first introduce the notion of a  _pseudo-attractor_. Then, we
proceed to define the problem of _source-target attractor control_.

In general, identifying attractors of a large Boolean network is a
computationally demanding task. Finding an attractor with the shortest period
is an NP-hard problem [2]. Moreover, in the case of classical PBNs, the fix-
point and limit cycle attractors correspond to the irreducible sets of states
in the underlying Markov chain [20]. For large-size PBNs with different
predictors for numerous individual genes, the limit cycle attractors may be
large, i.e., they may consist of many states. Nevertheless, usually states of
an irreducible set are not revisited with the same probability. From the point
of view of the control problem in the context of cellular reprogramming, only
the frequently revisited states of an attractor are the relevant ones since
they correspond to phenotypical cellular states that are observable in the
lab. This makes these states ‘recognisable’ for the application of cellular
reprogramming interventions in practice in accordance with the control
strategy suggested by our computational framework. We refer to the subset of
frequently revisited states of an attractor as a  _pseudo-attractor_
associated with the attractor and define it formally as follows.

######  Definition 4.1 (Pseudo-attractor).

Let $A$ be an attractor of a PBN in the classical formulation, i.e., an
irreducible set of states of the Markov chain underlying the PBN. Let $n:=|A|$
be the size of the attractor $A$ and let $\mathbb{P}_{A}$ be the unique
stationary probability distribution on $A$. The _pseudo-attractor_ associated
with $A$ is the maximal subset $PA\subseteq A$ such that for all $s\in PA$ it
holds that $\mathbb{P}_{A}(s)\geq\frac{1}{n}$. The states of a pseudo
attractor are referred to as _pseudo-attractor states_.

The correctness of the definition is guaranteed by the fact that the state
space of the underlying Markov chain of a PBN is finite and that the Markov
chain restricted to the attractor states is irreducible. It is a well known
fact that all states of a finite and irreducible Markov chain are positive
recurrent. In consequence, the attractor restricted Markov chain has a unique
stationary distribution. Furthermore, for any PBN attractor there exists a
non-empty pseudo-attractor as stated by the following lemma.

######  Observation 4.2.

Let $A$ be an attractor of a PBN. Then there exists a pseudo-attractor
$PA\subseteq A$ such that $|PA|\geq 1$.

###### Proof.

Let $n$ be the size of the attractor $A$, i.e., $n:=|A|$. Since the underlying
Markov chain of the PBN restricted to $A$ is irreducible and positive
recurrent, it has a unique stationary distribution, which we denote
$\mathbb{P}_{A}$. We proceed to show that there exists at least one state
$s^{\prime}\in A$ such that $\mathbb{P}_{A}(s^{\prime})\geq\frac{1}{n}$. For
this, let us assume that no such state exists. Then, we have that

$\sum_{s\in A}\mathbb{P}_{A}(s)<\sum_{s\in A}\frac{1}{n}=n\cdot\frac{1}{n}=1.$

The left-hand side of the above inequality is strictly less than 1 and hence
$\mathbb{P}_{A}$ is not a probability distribution on $A$, which leads to a
contradiction. In consequence, $|PA|\geq 1$. ∎

######  Observation 4.3.

In the case of the uniform stationary distribution on an attractor, the
associated pseudo-attractor is equal to the attractor: Let $A$ be an attractor
of a PBN such that the unique stationary distribution of the underlying Markov
chain of the PBN restricted to $A$ is uniform. Then, for the pseudo-attractor
$PA$ associated with $A$ it holds that $PA=A$.

###### Proof.

Let $n$ be the size of the attractor $A$. By the assumption of uniformity of
the stationary distribution $\mathbb{P}_{A}$ it holds that
$\mathbb{P}_{A}(s)=\frac{1}{n}$ for each $s\in A$. Since the pseudo-attractor
$PA$ is the maximal subset of $A$ such that $\mathbb{P}_{A}(s)\geq\frac{1}{n}$
for each $s\in A$, it folows that $PA=A$. ∎

Finally, we argue that Def. 4.1 of the pseudo-attractor straightforwardly
extends to BNs under the asynchronous update mode and that Obs. 4.2 and Obs.
4.3 remain valid in this case. Indeed, the asynchronous dynamics of the BN
restricted to a multi-state attractor of the network is a finite and
irreducible Markov chain. Therefore, in the continuation, we use the notion of
the pseudo-attractor both in the context of PBNs and BNs.

###  4.2 Source-target attractor control

With the biological context of cellular reprogramming in mind, we proceed to
define our control problem for BN and PBN models of GRNs. We start with
providing the definition of an  _attractor-based control strategy_ , also
referred to as _control strategy_ for short. Then, we define _Source-Target
Attractor Control_ , and immediately follow with an example. Note that in Def.
4.4 pseudo-attractor states are considered and not pseudo-attractors. This is
due to the fact that the procedure that will be introduced in Sec. 5.1
identifies pseudo-attractor states but does not assign them to individual
pseudo-attractors. Note that our definition of the source-target attractor
control is a generalisation of the ‘attractor-based sequential instantaneous
control (ASI)’ problem for BNs defined in [23] as our formulation of the
control problem extends to the formalism of PBNs and pseudo-attractor states.
An exact ‘divide-and-conquer’-type algorithm for solving the ASI problem for
BNs was provided in [23], and implemented in the software tool CABEAN [24].

######  Definition 4.4.

(Attractor-based Control Strategy) Given a BN/ PBN and a pair of its source-
target (pseudo-)attractors, a attractor-based control strategy is a sequence
of interventions which drives the network dynamics from the source to the
target (pseudo-)attractor. Interventions are understood as simultaneous flips
(perturbations) of values for a subset of genes in a particular network state
and their application is limited to (pseudo-)attractor states. We will denote
simultanious inverventions as sets, e.g. $\\{x_{1},x_{3},x_{7}\\}$ and
strategies as lists of sets, e.g.
$[\\{x_{1}x_{7}\\},\\{x_{2}\\}\\{x_{2},x_{4}\\}]$. Furthermore, the _length of
a control strategy_ is defined as the number of interventions in the control
sequence. We refer to an attractor-based control strategy of the shortest
length as the _minimal attractor-based control strategy_.

######  Definition 4.5 (Source-Target Attractor Control).

Given a BN/ PBN and a pair of source-target attractors or pseudo-attractor
states, find a minimal attractor-based control strategy.

######  Example 4.6.

The PBN from Example 3.5 may be controlled from state $(1,0,1,0)$ to
$(0,0,0,0)$ by intervening on $x_{1}$ and allowing the PBN to evolve in
accordance with its original dynamics:

$(1,0,1,0)\xrightarrow{i=1}(0,0,1,0)\xrightarrow{\text{evolution}}(0,0,0,0).$

However, the evolution is non-deterministic and the PBN may evolve to another
attractor, see Fig. 1:

$(0,0,1,0)\xrightarrow{\text{evolution}}(0,1,0,1).$

The only way to be sure to move to $(0,0,0,0)$ is to flip genes
$\\{x_{1},x_{3}\\}$ either simultaneously, which gives a strategy of length
one, or one-by-one, which gives a strategy of length two, i.e.,
$[\\{x_{3}\\},\\{x_{1}\\}]$.

##  5 DRL-based framework for the source-target attractor control

We propose a DRL-based computational framework, i.e., pbn-STAC, for solving
the source-target attractor control problem. Since our control problem is to
some extent similar to the one considered in [14] and our implementation is
based on the implementation therein, we compare our framework assumptions and
solutions to theirs during the presentation of pbn-STAC. In contrast to the
synchronous PBN update mode in [14], we consider the asynchronous update mode,
which is commonly considered as more appropriate for the modelling of
biological systems. The approach of [14] allows DRL agents to apply control
actions in any state of the PBN environment. Since our focus is on the
modelling of cellular reprogramming, we believe that this approach may be hard
to apply in experimental practice. It would require the ability to discern
virtually all cellular states, including the transient ones, which is
impossible with currently available experimental techniques. Since attractors
correspond to cellular types or, more generally, to cellular phenotypic
functional states, which are more easily observable in experimental practice,
we allow our DRL agent to intervene only in (pseudo-) attractor states in
consistency with the control problem formulation in Sec. 4.

In the control framework of [14], an action of the DRL agent can perturb at
most one gene at a time. However, for our formulation of the control problem
this is too restrictive. We have encountered examples of source-target
attractor pair where no control strategy consisting of such actions exists.
Therefore, we need to relax this restriction. However, we do not want to
intervene on too many genes at once as it would be rather pointless – in the
extreme case of allowing all genes to be perturbed at once, one could simply
flip all of the unmatched gene values. Furthermore, such an intervention would
also be hard to realise or even be unworkable in real biological scenarios –
it is expensive and sometimes even impossible to intervene on many genes at
once in the lab. Hence, we introduce a parameter which value defines an upper
limit for the number of genes that can be simultaneously perturbed. Based on
experiments (data not shown), we set this value to three. This setting is
sufficient for obtaining successful control strategies for all of our case
studies, yet low enough not to trivialise the control problem. Of course, the
value can be tuned to meet particular needs.

The DRL agent in [14] learns how to drive the network dynamics from any state
to the specified target attractor. With the context of cellular reprogramming
in mind, we consider in our framework only attractors as control sources and
targets with both of them specified. This models the process of transforming a
cell from one type into another. To be able to solve the source-target
attractor control problem, we define the reward function $R_{a}(s,s^{\prime})$
as:

$R_{a}(s,s^{\prime})=1000*\mathbbm{1}_{TA}(s^{\prime})-|a|,$

where $\mathbbm{1}_{TA}$ is an indicator function of target attractor, and
$|a|$ is the number of genes perturbed by applying action $a$. The loss
function is defined as the Mean Squared Error (MSE) between the predicted
Q-values and the target Q-values, calculated using the Bellman equation.

To train our DRL agent in each episode, we randomly choose a source-target
attractor pair and terminate each episode after 20 unsuccessful steps. This
approach however requires all the attractors to be known prior to training.
For networks with small numbers of nodes, the attractors can be computed.
However, as already mentioned, obtaining the list of all attractors for large
networks is a challenging problem by itself and one cannot expect the list to
be available in advance. To address this issue, we have introduced the notion
of a pseudo-attractor in Def. 4.1. Now we proceed to present a procedure for
detecting pseudo-attractor states which is exploited by our framework for
solving the control problem for large networks, i.e., ones for which
information on attractors is missing.

###  5.1 Pseudo-attractor states identification procedure

Identification of pseudo-attractors is hindered in large-size PBN models.
Nevertheless, pseudo-attractor states can be identified with simulations due
their property of being frequently revisited. We propose the following Pseudo-
Attractor States Identification Procedure (PASIP) which consists of two steps
executed in two phases: Step I during PBN environment pre-processing and Step
II with two cases, referred to as Step II-1 and Step II-2, during DRL agent
training.

##### Pseudo-Attractor States Identification Procedure

  1. I

During environment pre-processing, a pool of $k$ randomly selected initial
states is considered, from which PBN simulations are started. Each PBN
simulation is run for initial $n_{0}=200$ time steps, which are discarded,
i.e., the so-called burn-in period. Then, the simulation continues for
$n_{1}=1000$ time steps during which the visits to individual states are
counted. All states in which at least $5\%$ of the simulation time $n_{1}$ is
spent are added to the list of pseudo-attractor states.

  2. II

During training, the procedure discerns two cases:

     * II-1

The simulation of the PBN environment may enter a fix-point attractor not
detected in Step I. If the simulation gets stuck in a particular state for
$n_{2}=1000$ steps, the state is added to the list of pseudo-attractor states.

     * II-2

During training, the simulation of the PBN environment may enter a multi-state
attractor that has not been detected in Step I. For this reason, a history of
the most recently visited states is kept. When the history buffer reaches the
size of $n_{3}=10000$ items, revisits for each state are counted and states
revisited more than 15% of times are added to the list of pseudo-attractor
states. If no such state exists, the history buffer is cleared and the
procedure continues for another $n_{3}$ time steps. The new pseudo-attractor
states are added provided no known pseudo-attractor state was reached.
Otherwise, the history information is discarded.

Notice that the procedure allows us to identify the pseudo-attractor states,
but does not allow us to assign them to individual pseudo-attractors.
Therefore, when training a DRL agent with the use of pseudo-attractor states,
we consider the control strategies between all source-target pairs of pseudo-
attractor states. This is why our formulation of the control problem, the DRL
agent is restricted to apply its actions only in PBN (pseudo-)attractor
states. Therefore, in the case of large networks where no information on
attractors is available, the environment pre-processing phase is important as
it provides an initial pool of pseudo-attractor states for the training.

When identifying pseudo-attractors, we do not know the size of the PBN
attractor with which the pseudo-atractor is associated. Therefore, we cannot
determine the exact probability threshold of Def. 4.1 for identifying
individual pseudo-attractors. The proposed procedure addresses this issue as
follows. In Step I, the chosen $5\%$ identification threshold enables the
identification of a pseudo-attractor being a complete attractor, i.e., all its
states, of size up to $20$ states, which follows from the following
observation.

######  Observation 5.1.

For any PBN attractor $A$, the size of the associated pseudo-attractor $PA$
found by Step I of the pseudo-attractor identification procedure with $k\%$
identification threshold is exactly upper bounded by

$|PA|\leq\begin{cases}\frac{100}{k}-1,&100\text{ mod }k=0\text{ and
}|A|>\frac{100}{k}\\\
\left\lfloor\frac{100}{k}\right\rfloor,&\text{otherwise}.\end{cases}$

###### Proof.

Let $S=\frac{100}{k}-1$ if $k\mid 100$ and
$|A|=\left\lfloor\frac{100}{k}\right\rfloor$ and assume that we have an
attractor $A$ of size strictly greater than $S$. Then by the pigeonhole
principle one of the states has to be visited less then $\frac{100}{S}\%<k\%$
of times. So it will not be fully recovered by the procedure. Hence
$|PA|<S+1$.

Contrary, if we have an attractor $A$ of size $S$, which has uniform
distribution then $PA$ will equal exactly $|A|$, so the upper bound for the
size of $|PA|$ is at least $S$. ∎

In light of Def. 4.1 and Obs. 4.3, the associated pseudo-attractor of an
attractor of size $20$ can be identified only if the stationary distribution
on the attractor is uniform. If the attractor size is less than 20, it is
still possible to include all attractor states in the pseudo-attractor even if
the distribution is non-uniform. Notice that with decreasing size of the
unknown attractor, our procedure allows more and more pronounced deviations
from the uniform distribution while preserving the complete attractor
detection capability provided the stationary probabilities of all attractor
states are above the threshold.

If an attractor is of size larger than $20$ states, Step I of our procedure
with $5\%$ identification threshold will identify the associated pseudo-
attractor only if the stationary distribution is non-uniform and the pseudo-
attractor will contain only the most frequently revisited states. The maximum
possible size of the identified pseudo-attractor in this case is $19$, which
follows from Obs. 5.1. This is a desired property of our procedure as it keeps
the number of pseudo-attractor states manageable which has significant
positive influence on stabilising the model training as will be discussed
below.

The environment pre-processing phase provides an initial set of pseudo-
attractor states. The initial set is expanded in Step II during the model
training phase. Step II-1 allows to identify plausible fix-point attractors.
Step II-2 enables the identification of plausible multi-state attractors.
However, here the focus is on smaller attractors than in the case of Step 1:
we classify states as pseudo-attractor states if they are revisited at least
$15\%$ of time, which corresponds to attractors of size 6. This is to restrict
the number of spurious pseudo-attractor states in order to stabilise model
training as explained next.

We have encoutered an issue related to late discovery of pseudo-attractor
states during training. As can be observed in Fig. 2(a), the procedure may
detect a new pseudo-attractor at any point in time which destabilises
training: a new state is detected at around 90000 steps, which causes abrupt,
significant increase of the average episode length. We propose a remedy to
this problem in Sec. 5.2. Our experiments with small networks, i.e., ones for
which exact attractors could be computed, revealed that it is beneficial to
underestimate the set of attractor states in Step I of the procedure as the
missed ones are usually discovered later during the training phase.

For big networks, e.g., with hundreds of nodes, the set of pseudo-attractors
may take a long time to stabilise. Yet this approach provides us with the
ability to process networks too big to be handled by traditional methods. The
computations of pseudo-attractors can be parallelised in a rather
straightforward way to speed up the detection. Furthermore, the notion of
pseudo-attractors can easily be generalised to other types of GRN models,
e.g., PBNs with perturbation, which is yet another well-established GRN
modelling framework.

###  5.2 Exploration probability boost

The approach of [14] implements the $\varepsilon$-greedy policy in order to
balance exploitation and exploration of the DRL agent during training. The
$\varepsilon$-greedy policy introduces the _exploration probability_
$\varepsilon$ and with probability $1-\varepsilon$ follows the greedy policy,
i.e., it selects the action $a^{*}=\arg\max_{a\in\mathcal{A}}Q(s,a)$, or with
the $\varepsilon$ exploration probability selects a random action. We set the
initial $\varepsilon$ value to 1 and linearly decrease it to 0.05 over the
initial 3000 steps of training.

Combining the original $\varepsilon$-greedy policy with online pseudo-
attractor states identification gives rise to unstable training. When trying
to train the DRL agent for our control problem while keeping identifying
pseudo-attractor states during training, stability issues discussed in Sec.
4.1 were observed. To alleviate this negative influence on training, we
introduce the _exploration probability boost_ (EPB) to the
$\varepsilon$-greedy policy. The idea of EPB is to increase the exploration
probability $\varepsilon$ after each discovery of a new pseudo-attractor to
$\max(\varepsilon,0.3)$ if the current value of $\varepsilon$ is less than
0.3. After the increase, the linear decrease to 0.05 follows with the rate of
the initial decrease. As revealed by our experiments, this simple technique
makes learning much more stable. This is illustrated in Fig. 2(b), where the
agent discovers new pseudo-attractor states at around the 150000-th training
step and the use of the improved $\varepsilon$-greedy policy allowed us to
reduce the increase of the average episode length in a significant way and
resulted in a quick return to the previously trained low value of the average
episode length.

((a)) Training without EPB.

((b)) Training with EPB.

Figure 2: Examples of average episode lengths during training run with and
without EPB. New pseudo-attractor states are being identified during training.

###  5.3 pbn-STAC implementation

We implement pbn-STAC as a fork of gym-PBN [8], an environment for modelling
of PBNs, and pbn-rl [9], a suite of DRL experiments for a different PBN
control problem formulated in [14]. In pbn-STAC, we have adapted the original
code of gym-PBN and pbn-rl to our formulation of the PBN control problem,
i.e., the source-target attractor control. First, we extend gym-PBN by adding
the asynchronous PBN environment to it. Second, to allow for simultaneous
perturbation of a combination of genes within a DRL action, we replace the
original DDQN architecture with the BDQ architecture [25], which, contrary to
DDQN, scales linearly with the dimension of the action space. The architecture
of our BDQ network is depicted in Fig. 3.

Figure 3: Schematic illustration of the BDQ network architecture

Third, we implement the pseudo-attractor states identification procedure and
the exploration probability boost technique. Finally, the framework takes as
input a source-target pair of attractors or pseudo-attractor states. In the
case of a multi-state target attractor, a training episode is regarded as
successful if any of the target attractor states is reached. For a multi-state
source attractor, we uniformly sample one of its states and set it as the
initial state. In this way, different source attractor states are considered
as initial during DRL agent training. Our DRL-based framework for the source-
target attrator control is made available via the dedicated pbn-STAC GitHub
repository [28].

##  6 Experiments

###  6.1 BN and PBN models of GRNs

Melanoma models. We infer BN and PBN environments of various sizes for the
melanoma GRN using the gene expression data provided by Bittner _et al._ in
[6]. This is a well-known dataset on melanoma, which is extensively studied in
the literature, see, e.g., [14, 5, 7]. To infer the BN/PBN structures, we
follow the approach of [14] implemented in gym-PBN [14]. It is based on the
coefficient of determination (COD), which is a measure of how well the
dependent variable can be predicted by a model, a perceptron in the case of
[14].

The original dataset of Bittner _et al._ is quantised by the method of [14].
Then, the BN and PBN models of sizes 7, 10, and 28 are obtained from these
data. The models are denoted as BN-x or PBN-x, respectively, where x is the
number of genes. To infer the predictors for the models, we set the number of
predictors for each gene to 1 for BN models and to 3 for PBN models. For each
gene, the algorithm selects the Boolean functions with the maximum COD values.
For more details on the inference method, we refer to [14].

Case study of B. bronchiseptica. We test our DRL-based control framework on an
existing model of a real biological system, i.e., the network of immune
response against infection with the respiratory bacterium _Bordetella
bronchiseptica_ , which was originally proposed and verified against empirical
findings in [26]. The computational model, denoted IRBB-33, is an asynchronous
BN consisting of 33 genes.

###  6.2 Performance evaluation methodology

We evaluate the performance of pbn-STAC in solving the control problem
formulated in Sec. 4 on BN and PBN models of melanoma of various sizes, i.e.,
incorporating 7, 10, and 28 genes. Moreover, we consider IRBB-33, the 33-genes
BN model. The dynamics under the asynchronous update mode is considered for
all models. The evaluation consists of the agent interacting with the
environment by taking actions, where an action consists of flipping the values
of a particular subset of genes in an attractor or pseudo-attractor state. We
recover a control strategy for a given source-target pair learned by a trained
DRL agent byinitialising the BN/PBN environment with the source and target and
letting it run while applying the actions suggested by the DRL agent in the
source and all the intermediate (pseudo-)attractor states encountered on the
path from the source to the target. To evaluate the performance of pbn-STAC on
a particular BN/PBN model, we recover control strategies for all possible
ordered source-target pairs of the model’s attractors or pseudo-attractor
states. For all the BN models of melanoma and IRBB-33, we are able to compute
all their attractors and optimal control strategies for all pairs of
attractors using the CABEAN software tool with the attractor-based sequential
instantaneous source-target control (ASI) method. We use the information on
the exact attractors and optimal control strategies for BN models as ground
truth for the evaluation of pbn-STAC.

For PBN-7 and PBN-10, we compute the attractors with the NetworkX package
[12], which facilitates the analysis of complex networks in Python.
Unfortunately, due to very large memory requirements, we are unable to obtain
the attractors of the 28-genes PBN model of melanoma with this approach, so we
consider pseudo-attractor states. The optimal-length control strategies are
obtained for PBN-7 and PBN-10 models by exhaustive search.

Notice that due to the nondeterministic nature of our environments, i.e., the
asynchronous update mode, the results may vary between runs. Therefore, for
each source-target pair, we repeat the run 10 times. For each recovered
control strategy, we count its length, and record the information whether the
target attractor is reached. For a given BN/PBN model, we report the
percentage of successful control strategies found and the average length of
the successful control strategies.

##  7 Results

###  7.1 Identification of pseudo-attractor states

We evaluate the performance of PASIP proposed in Sec. 5.1. For this purpose,
we run pbn-STAC with PASIP on the considered BN and PBN models. We present the
obtained results in Tab. 1. For each model, except the melanoma PBN-28 for
which the exact attractors could not be obtained, we provide the information
on the number of exact attractors, the total number of attractor states, and
the total numbers of identified pseudo-attractor states with our procedure. We
measure the precision of our approach defined as
$\textrm{TP}/(\textrm{TP}+\textrm{FP})$, where TP is the number of true
positives, i.e., the number of pseudo-attractor states that are attractor
states, and FP is the number of false positives, i.e., the number of states
identified as pseudo-attractor states which are not part of any of the
network’s attractor. We can conclude that for all cases in which the exact
attractors are known, our procedure does not introduce any FPs. Moreover, it
can identify the attractor states with $100\%$ precision in all but one case,
i.e., the PBN-28 network which has 2412 fix-point attractors and for which our
procedure correctly identifies 1053 of them. This justifies our strong belief
that running our procedure for longer time would result in definitely higher
precision also in the case of BN-28. In summary, the presented results show
that PASIP is reliable.

Model | #Attr. | #Attr. states | #PA-states | Precision
---|---|---|---|---
BN-7 | 6 | 6 | 6 | 100%
BN-10 | 26 | 26 | 26 | 100%
BN-28 | 2412 | 2412 | 1053 | 43.65%
IRBB-33 | 3 | 3 | 3 | 100%
PBN-7 | 4 | 4 | 4 | 100%
PBN-10 | 6 | 6 | 6 | 100%
PBN-28 | unknown | unknown | 14 | N/A
Table 1: Comparison of the number of exact attractor states and pseudo-attractor states identified by PASIP for various BN and PBN models. The fact that we were unable to obtain the exact attractors for the PBN-28 model is indicated with ‘unknown’. Attr. is short for Attracor and PA stands for Pseudo-attractor. Model | #Attractors | Optimal Strategy | pbn-STAC
---|---|---|---
BN-7 | 6 | 1.0 | 3.98
BN-10 | 26 | 1.1 | 2.14
BN-28 | 2412 | 1.1 | -
IRBB-33 | 3 | 1.0 | 9.2
PBN-7 | 4 | 1.1 | 5.5
PBN-10 | 6 | 1.2 | 15.2
PBN-28 | unknown | unknown | 60.7
Table 2: Average lengths of pbn-STAC control strategies and optimal control
strategies obtained with CABEAN (BNs) or exhaustive search (PBNs) for all
source-target pairs of individual models. The fact that we were unable to
obtain the optimal strategy for the PBN-28 model is indicated with ‘unknown’.

###  7.2 Control of BN models of melanoma

We evaluate the ability of pbn-STAC to solve the control problem by comparing
the obtained results to the optimal ASI control strategies computed with
CABEAN. As can be seen in Tab. 2, the strategies obtained with pbn-STAC for
larger BN models tend to be longer on average compared to the optimal ones.
However the overhead is rather stable across different models. We investigate
the issue of longer control strategies further by computing a histogram of
control strategy lengths for the BN-7 model provided in Fig. 5(a). It is
apparent that in most cases the control strategies are short and close to the
optimal ones. Nevertheless, there are a few cases of longer control strategies
that give rise to the higher average values. The longer strategies are present
due to the fact that the interventions suggested by the trained DRL-agent
often place the system in a so-called _weak basin of attraction_ of an
attractor, i.e., a set of states from which the attractor is reachable, but
not necessarily – the dynamics can still lead the system to another attractor
from these states due to non-determinism arising from the asynchronous update
mode. The strategies computed by CABEAN are optimal since they are obtained by
considering the so-called _strong basins of attraction_ , i.e., states from
which only a single attractor can be reached. Nevertheless, determining strong
basins is challenging, not to say impossible, for large networks (see [18] for
details). In light of this and the fact that pbn-STAC can handle larger
networks, the obtained results can be seen as reasonable and acceptable.

Unfortunately, due to the huge number of attractors of the BN-28 model, the
training of pbn-STAC on this model needs to be run for much longer time and we
did not manage to finish it within our time limits. Notice that our training
procedure considers all ordered pairs of the attractors. Handling such cases
requires further research.

((a)) Usual reward

((b)) Improved reward

Figure 4: Training of the DRL agent on the IRBB-33 environment with different
reward schemes.

###  7.3 Control of the IRBB-33 model

In the case of the IRBB-33 network, we have to modify the reward scheme. As
can be seen in Fig. 4(a), the reward scheme introduced in Sec. 5, referred to
as the _mixed reward_ , does not lead to any improvement of the average
episode length during training over 200 000 training steps. After trying
different reward schemes for this network (data not shown), we found that the
following scheme

$\displaystyle R_{a}(s,s^{\prime})=-|a|+100*(\mathbbm{1}_{TA}(s^{\prime})-1),$

It improves the training of the DRL agent significantly, as can be seen in
Fig. 4(b), where the convergence is achieved in tens of thousands of steps.

((a)) BN-7

((b)) BN-10

((c)) PBN-28

((d)) BN-33

Figure 5: Histogram of the control strategy lengths for the BN-7 model.

The average control strategy length obtained with pbn-STAC is $9.2$, as
presented in Tab. 2. The length is again larger than in the optimal case, but
the overhead is stable with respect to the results obtained for the BN models
of melanoma. Again, as can be observed in Fig. 5(d), in the majority of cases
the strategies are of length one, which perfectly corresponds with the optimal
strategy. Unfortunately, there are a few very long ones, which give rise to
the higher average value.

###  7.4 Control of PBN models of melanoma

We run the pbn-STAC control framework on the three PBN models of melanoma. For
PBN-28, we are not able to compute the set of exact attractors, but we
identify 14 pseudo-attractor states. Unfortunately, we can not obtain the
optimal control strategies for this network with exhaustive search.

As shown in Tab. 2, the control strategies found by pbn-STAC are on average
longer than the optimal ones. Moreover, their lengths seem to increase with
the size of the network faster than in the case of BN models. Unfortunately,
the optimal result is not available for the PBN-28 model to make a comparison.
Although the average length for the this network is high, the distribution is
heavily skewed with a long tail of longer control strategies as can be seen in
Fig. 5(c). Nevertheless, once again the majority of the source-target pairs
are controllable with very few interventions. This characteristic of the
control strategies obtained with pbn-STAC is consistent across different
models and theirs types.

##  8 Conclusions

In this study we formulated a control problem for the BN and PBN frameworks
under the asynchronous update mode that corresponds to the problem of
identifying effective cellular reprogramming strategies. We have developed and
implemented a computational framework, i.e., pbn-STAC, based on DRL that
solves the control problem. It allows to find proper control strategies that
drive a network from the source to the target attractor by intervening only in
other attractor states that correspond to phenotypical functional cellular
states that can be observed in the lab. Since identifying attractors of large
BNs/PBNs is a challenging problem by itself and we consider our framework as a
contribution towards developing scalable control methods for large networks,
we introduced the notion of a pseudo-attractor and developed a procedure that
identifies pseudo-attractor states during DRL agent training. We evaluate the
performance of pbn-SEC on a number of networks of various sizes and a
biological case study and compare its solutions with the exact, optimal ones
wherever possible.

The obtained results show the potential of the framework in terms of
effectiveness and at the same time reveal some bottlenecks that need to be
overcome to improve the performance. The major identified issue is related to
the long tails of the distributions of the lengths of strategies identified by
pbn-SEC, i.e., there are many strategies of lengths close to the optimal ones,
and a few which are very long. This negatively influences the average value.
Addressing this problem would allow us to significantly improve the
performance of pbn-STAC and make it effective on large models. We consider
these developments and the evaluations of pbn-STAC on models of large sizes as
future work.

Finally, we perceive our framework as rather straightforwardly adaptable to
other types of PBNs, such as PBNs with perturbations, or Probabilistic Boolean
Control Networks.

## References

  * [1] A. Acernese and et al.  Double Deep-Q Learning-Based Output Tracking of Probabilistic Boolean Control Networks.  IEEE Access, 2020.
  * [2] T. Akutsu, S. Kuhara, O. Maruyama, and S. Miyano.  Identification of genetic networks by strategic gene disruptions and gene overexpressions under a boolean model.  Theoretical Computer Science, 298(1), 2003.
  * [3] A.-L. Barabási, N. Gulbahce, and J. Loscalzo.  Network medicine: a network-based approach to human disease.  Nature Reviews Genetics, 12(1):56–68, 2011.
  * [4] N. Beneš and et al.  Phenotype Control of Partially Specified Boolean Networks.  In Proc. 21st International Conference on Computational Methods in Systems Biology (CMSB’23). Springer-Verlag, 2023.
  * [5] N. Beneš, L. Brim, O. Huvar, S. Pastva, and D. Šafránek.  Boolean network sketches: a unifying framework for logical model inference.  Bioinformatics, 39(4), Apr. 2023.
  * [6] M. Bittner and et al.  Molecular classification of cutaneous malignant melanoma by gene expression profiling.  Nature, 406(6795):536–540, Aug. 2000.
  * [7] Q. Du, Y. Lin, C. Ding, L. Wu, Y. Xu, and Q. Feng.  Pharmacological activity of matrine in inhibiting colon cancer cells vm formation, proliferation, and invasion by downregulating claudin-9 mediated emt process and mapk signaling pathway.  Drug Design, Development and Therapy, Volume 17:2787–2804, Sept. 2023.
  * [8] C. Evangelos.  gym-pbn.  https://github.com/UoS-PLCCN/gym-PBN.
  * [9] C. Evangelos.  pbn-rl.  https://github.com/UoS-PLCCN/pbn-rl.
  * [10] S. Fujimoto, H. van Hoof, and D. Meger.  Addressing function approximation error in actor-critic methods, 2018.
  * [11] A. J. Gates and L. M. Rocha.  Control of complex networks requires both structure and dynamics.  Scientific Reports, 6:Article 24456, 2016.
  * [12] A. A. Hagberg, D. A. Schult, and P. J. Swart.  Exploring network structure, dynamics, and function using NetworkX.  In Proceedings of the 7th Python in Science Conference (SciPy2008), 2008.
  * [13] S. Huang, G. Eichler, Y. Bar-Yam, and D. E. Ingber.  Cell fates as high-dimensional attractor states of a complex gene regulatory network.  Physical Review Letters, 2005.
  * [14] S. Moschoyiannis, E. Chatzaroulas, V. Šliogeris, and Y. Wu.  Deep Reinforcement Learning for Stabilization of Large-scale Probabilistic Boolean Networks.  IEEE Transactions on Control of Network Systems, 10(3):1412–1423, 2022.
  * [15] C. E. H. Nishida, R. A. C. Bianchi, and A. H. R. Costa.  A framework to shift basins of attraction of gene regulatory networks through batch reinforcement learnin.  Artificial Intelligence in Medicine, 107, 2020.
  * [16] C. E. H. Nishida, A. H. R. Costa, and R. A. C. Bianchi.  Control of gene regulatory networks basin of attractions with batch reinforcement learning.  In Proc. 7th Brazilian Conference on Intelligent Systems. IEEE CS, 2018.
  * [17] S. Padakandla.  A survey of reinforcement learning algorithms for dynamically varying environments.  ACM Computing Surveys, 54(6):1–25, July 2021.
  * [18] S. Paul, C. Su, J. Pang, and A. Mizera.  An efficient approach towards the source-target control of Boolean networks.  IEEE/ACM Transactions on Computational Biology and Bioinformatics, 17(6):1932–1945, 2020.
  * [19] L. Paulevé.  Marker and source-marker reprogramming of Most Permissive Boolean networks and ensembles with BoNesis.  Peer Community Journal, 3:Article e30, 2023.
  * [20] I. Shmulevich and et al.  Probabilistic Boolean networks: a rule-based uncertainty model for gene regulatory networks.  Bioinformatics, 18(2), 2002.
  * [21] U. Sirin, F. Polat, and R. Alhajj.  Employing batch reinforcement learning to control gene regulation without explicitly constructing gene regulatory networks.  In Proc. 23rd International Joint Conference on Artificial Intelligence, pages 2042–2048. AAAI Press, 2013.
  * [22] C. Su and J. Pang.  A dynamics-based approach for the target control of Boolean networks.  In Proc. 11th ACM Conference on Bioinformatics, Computational Biology, and Health Informatics, pages 50:1–50:8. ACM Press, 2020.
  * [23] C. Su and J. Pang.  Sequential temporary and permanent control of Boolean networks.  In Proc. 18th International Conference on Computational Methods in Systems Biology, volume 12314 of Lecture Notes in Computer Science, pages 234–251. Springer-Verlag, 2020.
  * [24] C. Su and J. Pang.  CABEAN: A software for the control of asynchronous Boolean networks.  Bioinformatics, 37(6):879–881, 2021.
  * [25] A. Tavakoli, F. Pardo, and P. Kormushev.  Action branching architectures for deep reinforcement learning.  In AAAI Conference on Artificial Intelligence, 2018.
  * [26] J. Thakar, A. K. Pathak, L. Murphy, R. Albert, and I. M. Cattadori.  Network model of immune responses reveals key effectors to single and co-infection dynamics by a respiratory bacterium and a gastrointestinal helminth.  PLoS Computational Biology, 8(1):e1002345, Jan. 2012.
  * [27] C. J. C. H. Watkins and P. Dayan.  Q-learning.  Machine Learning, 8(3–4):279–292, May 1992.
  * [28] J. Zarzycki.  https://github.com/jakub-zarzycki2022/gym-pbn-stac, 2023.