jupytext | kernelspec | ||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
+++ {"latex_maketitle": {"corrAuthor": "Lo\u00efc Paulev", "corrEmail": "[email protected]", "firstAuthorLast": "Paulev\u00e9", "runningTitle": "Reprogramming of Most Permissive Boolean networks with BoNesis"}}
Marker and source-marker reprogramming of Most Permissive Boolean networks and ensembles with BoNesis
Loïc Paulevé
Univ. Bordeaux, CNRS, Bordeaux INP, LaBRI, UMR 5800, F-33400 Talence, France
Boolean networks (BNs) are discrete dynamical systems with applications to the modeling of cellular behaviors. In this paper, we demonstrate how the software BoNesis can be employed to exhaustively identify combinations of perturbations which enforce properties on their fixed points and attractors. We consider marker properties, which specify that some components are fixed to a specific value. We study 4 variants of the marker reprogramming problem: the reprogramming of fixed points, of minimal trap spaces, and of fixed points and minimal trap spaces reachable from a given initial configuration with the most permissive update mode. The perturbations consist of fixing a set of components to a fixed value. They can destroy and create new attractors. In each case, we give an upper bound on their theoretical computational complexity, and give an implementation of the resolution using the BoNesis Python framework. Finally, we lift the reprogramming problems to ensembles of BNs, as supported by BoNesis, bringing insight on possible and universal reprogramming strategies. This paper can be executed and modified interactively.
Keywords: Discrete dynamical systems, Control, Minimal trap spaces, Attractors, Reachability
+++
Boolean networks (BNs) are formal discrete dynamical systems with pertinent applications for modeling cellular differentiation and fate decision processes (Saez-Rodriguez et al., 2009; Cohen et al., 2015; Schwab et al., 2021; Zañudo et al., 2021; Montagud et al., 2022). In these applications, BNs aim at capturing the stable behaviors (attractors) and the transient dynamics (trajectories) of the cell. From this perspective, BNs offer a formal framework for predicting perturbations that destabilize the system and drive it towards a desired new stable behavior. BN control or BN reprogramming, in reference to cellular reprogramming which aims at converting cell types, is thus receiving a lot of interest from the computational systems biology community (Zañudo et al., 2015; Yang et al., 2018; Biane et al., 2018; Mandon et al., 2019; Cifuentes Fontanals et al., 2020; Su and Pang, 2020; Rozum et al., 2021).
The reprogramming of BNs led to a range of methods and tools addressing different instantiations of this problem: with different type of perturbations (instantaneous, temporary, permanent), different temporal modalities (one-step, sequential), different scopes (global reprogramming or from a given initial condition), different restrictions on the target attractor (fixed points only, attractors of the original "wild-type" BN). On top of that, the update mode of the BN, which determines how the trajectories are computed, can play an important role on the predictions.
+++
In this paper, we address the BN reprogramming with the Most Permissive (MP) update mode, where attractors are the minimal trap spaces of the BN (Paulevé et al., 2020). The problems we tackle are related to marker reprogramming: the desired target attractors are specified by a set of markers, associating a subset of nodes of the network to fixed values (e.g.,
We address the following BN reprogramming problems in the scope of the MP update mode:
- P1: Marker reprogramming of fixed points: after reprogramming, all the fixed points of the BN match with the given markers; optionally, we can also ensure that at least one fixed point exists.
- P2: Source-marker reprogramming of fixed points: after reprogramming, all the fixed points that are reachable from the given initial configuration match with the given markers.
- P3: Marker reprogramming of attractors: after reprogramming, all the configurations of all the MP attractors (the minimal trap spaces) of the BN match with the given markers.
- P4: Source-marker reprogramming of attractors: after reprogramming, all the configurations of all the attractors that are reachable from the given initial configuration match with the given markers.
MP fixed points match with the fixed points of the global Boolean map of the BN and are thus identical to the fixed points of the (a)synchronous update modes. MP attractors match with so-called minimal trap spaces of the BN, which are the smallest sub-hypercubes closed by the Boolean map. Problem P1 has been already addressed in the literature, notably by Biane et al. (2018) with the ActoNet method and by Moon et al. (2022a), based on bilevel integer programming. To our knowledge, none of the other problems have been addressed completely in the literature.
The software BoNesis (github.com/bnediction/bonesis) provides a generic environment for the automated construction of BNs with MP update mode from specified structural and dynamical properties. The properties are translated into a logic satisfiability problem, expressed in Answer-Set Programming (ASP). Initially, BoNesis has been designed for performing BN synthesis (Chevalier et al., 2019): solutions of the logic model correspond to BNs that possess the specified structural and dynamical properties. Leveraging this generic declarative specification of properties, BoNesis is a versatile tool for reasoning on BNs in general, with the MP update mode: besides synthesis, it can be used to do model checking, identify fixed points and attractors in ensemble of BNs, and identifying reprogramming strategies.
In this paper, we show how the software BoNesis can be employed to solve P1, P2, P3, and P4 in the scope of locally-monotone BNs, where each local function is unate, i.e., where each local function never depends on both positively and negatively from a same component. Locally-monotone BNs cover all the models where it assumed that a node cannot be both an activator and inhibitor of a same other node, which is a common assumption when modeling biological system.
BoNesis enables reasoning on ensembles of BNs: one of its basic input is the domain of BNs to consider. This domain could be reduced to a singleton BN: in that case, the reasoning is similar to standard model checking and reprogramming. In general, the domain is specified from an influence graph, possibly with additional constraints on the underlying logical functions. For BN synthesis, this domain is used to delimit symbolically the set of candidate models: BoNesis will output the subset of them that verify the desired dynamical properties. We show how problems P1 to P4 can be partly lifted to ensembles of BNs using this approach.
+++
The paper is structured as follows. The Methods section gives the necessary background on BNs and MP update mode and formulation of elementary dynamical properties as satisfaction problems, as well as main principles of the BoNesis environment. The Results section details how the different reprogramming problems P1-P4 can be addressed using BoNesis and shows some experiments to assess their scalability. Finally, the Discussion section sketches possible extensions of the addressed problems and underlines current challenges for their resolution.
This paper is executable: it contains snippets of Python code employing BoNesis to demonstrate its usage on small examples, including command line usage. Instructions for its execution are given at the beginning of the Results section. It can be visualized online at nbviewer.org/github/bnediction/reprogramming-with-bonesis/blob/release/paper.ipynb and interactively executed either online at mybinder.org/v2/gh/bnediction/reprogramming-with-bonesis/release using mybinder online free service, or locally, following instructions given later in this paper.
+++
A Boolean network (BN) of dimension
A BN
Example. The BN
Locally monotone BNs should not be confused with monotone BNs where a component appears in all local functions with the same sign. Monotone BNs are a particular case of locally-monotone BNs.
Given a BN
However, (a)synchronous update modes do not lead to a complete qualitative abstraction of quantitative systems and preclude the prediction of trajectories that are actually feasible when considering time scales or concentration scales. The Most Permissive (MP) (Paulevé et al., 2020; Paulevé and Sené, 2021) is a recently-introduced update mode which brings the formal guarantee to capture any trajectory that is feasible by any quantitative system compatible with the Boolean network (see (Paulevé et al., 2020) for details).
The main idea behind the MP update mode is to systematically consider a potential delay when a component changes state, and consider any additional transitions that could occur if the changing component is in an intermediate state. It can be modeled as additional dynamic states "increase" (
In this paper we will consider BN perturbations that modify the local functions of some components so they become a constant function. Perturbations model mutations, where a gene is silenced or constitutively activated. Mathematically, a perturbation is a map associating a set of components to a Boolean value, for instance,
+++
A Boolean expression is a logic formula composed of Boolean variables and propositional logic operators (negation
Deciding whether such an expression is true (satisfiable) is a fundamental problem in computer science. The complexity of this problem actually depends on the alternation of quantifiers. Thus, in the following we will classify the quantified Boolean expressions by their sequence of quantifiers
Computational complexity (Papadimitriou, 1993) is a fundamental theory of computer science to classify decision problems: a (decision) problem is in class C whenever there exists an algorithm of worst-case complexity C, C referring to either a time or space complexity. For instance, the class PTIME gathers all the problems that can be decided in time polynomial with the size of the input (e.g., the length of the Boolean expression).
The decision of satisfiability of
The decision of satisfiability of
Then, the alternation of quantifiers makes the problem climbing into the so-called polynomial hierarchy1.
The reader should keep in mind that the length of the expression is a crucial parameter for the decision complexity. When variables have a finite domain, one can rewrite quantified Boolean expression in a universal-free one. However, the length of the obtained expression will be exponentially larger.
In the rest of the paper, for the sake of simplicity, we will not fully detail the size of the quantified Boolean expression we derive, and are expected to be of length linear or polynomial with the size of the BN.
+++
We present the formal aspects of the MP dynamics that are employed in the rest of the paper, i.e., related to attractors and the reachability of attractors. The proofs and full MP definition and properties can be found in (Paulevé et al., 2020).
A sub-hypercube specifies for each dimension of the BN if it is either fixed to a Boolean value, or free: it can be characterized by a vector
A sub-hypercube
A sub-hypercube
The attractors of MP dynamics are the minimal trap spaces of the Boolean function
The computational complexity of decision problems related to minimal trap spaces has been extensively addressed in (Moon et al., 2022b) with different representations of Boolean networks.
For the case of local functions represented with propositional logic, as we consider here, deciding whether a sub-hypercube is a trap space is coNP-complete problem, whereas deciding whether it is a minimal trap space is a coNPcoNP-complete problem, i.e., equivalent to the decision of satisfiability of
Given a configuration
Let us denote by SAT(h, f[i] = -x[i])
is true if and only if there exists a configuration
Algorithm TS(x: configuration)
Returns sub-hypercube h
--
h := x
repeat
changed := false
for i in 1..n:
if h[i] != * and SAT(h, f[i] = -x[i]):
h[i] := *
changed := true
while changed
In the worst case, this algorithm makes a quadratic number of calls to the SAT
problem.
Therefore, the decision of MP reachability of attractors is in PTIMENP in general2, and PTIME in the locally-monotone case.
Note that the general MP reachability property is not addressed here, but its overall complexity is identical. With (a)synchronous update modes, it is a PSPACE-complete problem.
In the following, we will consider the problem of deciding whether a given configuration
Finally, given a set of perturbations
+++
BoNesis (github.com/bnediction/bonesis) is a Python library which has been primarily designed for identifying BNs satisfying user-given dynamical properties among a given domain of BNs and with the MP update mode.
It takes as input (1) a domain of BNs
Currently, the domain of BNs
- A singleton locally-monotone BN
$\mathbb F={f}$ . In that case, BoNesis can be employed as a model checker to verify that$f$ has the specified dynamical properties. In this paper, this is the main setting we will consider, in order to predict perturbations to reprogram the attractors of$f$ . - An explicit ensemble of locally-monotone BNs
$\mathbb F={ f^1,\cdots, f^m }$ . - Any locally-monotone BN matching with a given influence graph
$\mathcal G$ :$\mathbb F = { f\mid G(f)\subseteq \mathcal G}$ . An influence graph is a signed digraph between components, i.e., of the form$({1,\cdots,n},V)$ with$V\subseteq {1,\cdots,n}\times {+1,-1}\times {1,\cdots n}$ . The influence graph of a BN$f$ , denoted by$G(f)$ has an edge$i\xrightarrow{s} j$ if and only there exists a configuration$x\in\mathbb B^n$ such that$f_j(x_1, \ldots, x_{i-1}, 1, x_{i+1},\ldots, x_n) - f_j(x_1, \ldots, x_{i-1}, 0, x_{i+1},\ldots, x_n) = s$ . - Any locally-monotone BN matching with a partially-defined BN following the AEON framework (Beneš et al., 2021).
BoNesis offers a Python programming interface to declare the dynamical properties over BNs, including reachability, fixed points and trap spaces. BoNesis relies on Answer-Set Programming (ASP) and the ASP solver clingo for the enumeration of solutions. ASP is a declarative logic programming framework for expressing combinatorial decision problems and enumerate their solutions, possibly with optimizations. ASP can be employed for efficiently solving
We emphasize that BoNesis is currently restricted to locally-monotone BNs only for which efficient logical encoding of domains of models are possible. Whereas it is a common assumption when modeling of biological systems (a node cannot be both an activator and inhibitor of a same other node), non-monotone BNs are also employed, and cannot be addressed with the current implementation.
The usage of BoNesis Python programming interface and command line will be explained along with the code snippets provided in the next sections.
+++
We show how the general declarative approach of BoNesis can be instantiated to compute the complete solutions to the P1, P2, P3, and P4 reprogramming problems on BNs, and also extend the reasoning to ensembles of BNs. Importantly, note that BoNesis currently supports only locally-monotone BNs.
This is an executable paper which demonstrates the use of BoNesis for the reprogramming of BNs. The corresponding notebook can be downloaded from nbviewer.org/github/bnediction/reprogramming-with-bonesis/blob/release/paper.ipynb. Its execution requires the Jupyter notebook system, Python, and the Python package bonesis
to be installed, see github.com/bnediction/bonesis for instructions.
Alternatively, the notebook can be executed within the CoLoMoTo Docker distribution (Naldi et al., 2018), using the Docker image colomoto/colomoto-docker:2023-03-01
, which can be launched as follows:
pip install colomoto-docker
colomoto-docker -V 2023-03-01
then open http://127.0.0.1:8888 and upload the notebook from the Jupyter interface. See github.com/bnediction/reprogramming-with-bonesis for further help.
#!pip install --user bonesis # uncomment to install bonesis
import bonesis
from colomoto_jupyter import tabulate # for display
import pandas as pd # for display
import mpbn # for analyzing individual Boolean networks with MP update mode
from colomoto.minibn import BooleanNetwork
Alternatively, the computation of reprogramming perturbations from single Boolean networks can be performed using the command line programbonesis-reprogramming
, provided alongside the bonesis
Python package. We detail its usage in each case.
+++
We first consider the reprogramming of a single BN
A marker
The objective of marker-reprogramming is to identify perturbations so that all the fixed points/attractors of the perturbed
A very important aspect of marker reprogramming is that it accounts for the creation and deletion of attractors due to the perturbation. Thus, in general, the attractors of the reprogrammed BN are different from the attractors of the input (wild-type) BN.
In this section, we tackle the following instantiations of the reprogramming problem:
- Marker reprogramming of fixed points (P1);
- Source-marker reprogramming of fixed points (P2);
- Marker reprogramming of attractors (P3);
- Source-marker reprogramming of attractors (P4).
In each case, we briefly study the complexity of the associated decision problem (existence of a perturbation given the desired reprogramming property), and give the Python and command line recipe to identify the perturbations with BoNesis. The following table summarizes the results, with the complexity in the locally-monotone case and command line usage:
Problem | Complexity | Command line |
---|---|---|
P1 | [base] --fixpoints |
|
P2 | [base] --fixpoints --reachable-from z |
|
P3 | [base] |
|
P4 | [base] --reachable-from z |
where [base]
is the command line bonesis-reprogramming model.bnet M k
, with model.bnet
the path to a file in BooleanNet format, M
specifies the marker as a JSON map, k
is the maximum number of simultaneous perturbations, and z
is the initial configuration as a JSON map.
For instance,
bonesis-reprogramming model.bnet '{"A": 0, "C": 1}' 3 \
--reachable-from '{"A":1, "B":0, "C":0,"D":0}'
+++
We identify the perturbations
\begin{equation} \exists P\in \mathbb M^{\leq k}, \forall x\in\mathbb B^n, (f/P)(x) = x \Rightarrow x\models M \end{equation}
"There exists a perturbation being a map of at most
Remark that any BN having no fixed point verify the above equation with an empty perturbation. Thus, in practice, one may also expect that the perturbed BN possesses at least one fixed point:
\begin{equation} \exists P\in \mathbb M^{\leq k}, \exists y\in\mathbb B^n, (f/P)(y) = y, \forall x\in\mathbb B^n, (f/P)(x) = x \Rightarrow x\models M\enspace. \end{equation}
With the BoNesis Python interface, this reprogramming property can be declared as follows, where f
is a BN, M
the marker (specified as Python dictionary associating a subset of components to a Boolean value), and k
the maximum number of components that can be perturbed (at most
def marker_reprogramming_fixpoints(f: BooleanNetwork,
M: dict[str,bool],
k: int, ensure_exists:bool=True):
# f: Boolean network; M: marker; k: maximum number of components to perturb
bo = bonesis.BoNesis(f)
P = bo.Some(max_size=k) # perturbations to identify
with bo.mutant(P):
# impose that all the fixed points of the perturbed BN match with M
bo.all_fixpoints(bo.obs(M))
if ensure_exists:
# impose the existence of at least one fixed point matching with M
bo.fixed(~bo.obs(M))
return P.assignments()
The line bo = bonesis.BoNesis(f)
instantiates a BoNesis object bo
with the domain f
. In this section, we assume that f
is a single Boolean network. It can be either the path to a BooleanNet file, a minibn
object, for instance as returned by the biolqm.to_minibn
function for importing models from various formats, or a Python dictionary, associating components to a Boolean function. Examples are given below.
Then, the line P = bo.Some(max_size=k)
declares a variable that will consist of a map associating at most k
components to a Boolean value (the perturbation to be identified).
The instruction with bo.mutant(P)
opens a context where dynamical properties will be verified against the BN f
with the perturbation P
applied.
Within this mutant context, we declare with bo.all_fixpoints(bo.obs(M))
that each fixed point of the perturbed model matches with M
. Moreover, whenever ensure_exists
is true, the constraint bo.fixed(~bo.obs(M))
imposes that the configuration associated with M
(~bo.obs(M)
) is a fixed point.
Finally, P.assignments()
returns an iterator over all the possible assignments of P
that fulfill the above dynamical properties.
The corresponding command line is of the form
bonesis-reprogramming model.bnet M k --fixpoints
where model.bnet
is a BN encoded in BooleanNet format, M
specifies the marker as a JSON map, and k
is the maximum number of perturbations.
The existence of at least one fixed point can be lifted with the option --allow-no-fixpoint
.
+++
Example. We illustrate the marker-reprogramming of fixed points on a small toy BN example, which has no fixed point. In the following, we use colomoto.minibn.BooleanNetwork
to define this BN.
f = BooleanNetwork({
"A": "B",
"B": "!A",
"C": "!A & B"
})
f
This example BN has two components in negative feedback: they will oscillate forever. The state of the third component C
is then determined by the state of the oscillating components. The following command returns its influence graph:
---
nbconvert_latex:
nbconvert.figure.caption: Influence graph of the example BN for P1
---
f.influence_graph()
With the (fully) asynchronous update mode, the system has a single attractor, consisting of all the configurations of the network.
---
nbconvert_latex:
nbconvert.figure.caption: Fully-asynchronous dynamics of the example BN for P1
---
f.dynamics("asynchronous")
Recall that the fixed points are identical in asynchronous and MP. We use mpbn
to analyze the dynamical properties with the MP update mode:
mf = mpbn.MPBooleanNetwork(f)
list(mf.fixedpoints())
list(mf.attractors())
Indeed, the network has no fixed points, and its attractor is the full hypercube of dimension 3.
Using the marker_reprogramming_fixpoints
snippet defined above, we identify all perturbations of at most 2 components which ensure that (1) all the fixed points have C
active, and (2) at least one fixed point exists:
list(marker_reprogramming_fixpoints(f, {"C": 1}, 2))
Indeed, fixing A
to 0 breaks the negative feedback, and make B
converge to 1. There, C
converges to state 1.
Then, remark that fixing C
to 1 is not enough to fulfill the property, as A
and B
still oscillate. Thus, one of these 2 must be fixed as well, to any value. The solution {'A': 0, 'C': 1}
is not returned as {'A': 0}
is sufficient to acquire the desired dynamical property.
In our BoNesis code snippet defined above, by default we ensure that the perturbed BN possesses at least one fixed point. When relaxing this constraint, we obtain that the empty perturbation is the (unique) minimal solution, as f
has no fixed point.
list(marker_reprogramming_fixpoints(f, {"C": 1}, 2, ensure_exists=False))
In the following, we demonstrate how to perform the same computation with the command line. By default, the reprogramming of fixed points adds the constraint that at least one fixed point must exist.
:nbconvert.display: true
with open("example1.bnet", "w") as fp:
fp.write(f.source())
%cat example1.bnet
:nbconvert.display: true
!bonesis-reprogramming example1.bnet '{"C": 1}' 2 --fixpoints
Adding the option --allow-no-fixpoint
would return an empty perturbation as unique minimal solution.
+++
Given an initial configuration
\begin{equation} \exists P\in\mathbb M^{\leq k}, \forall x\in\mathbb B^n, ((f/P)(x)=x \wedge \operatorname{reach}_P(z,x))\implies x\models M \end{equation}
"There exists a perturbation
As explained in the Method section, the reachability property boils down to computing the smallest trap space containing
\begin{equation} \exists P\in\mathbb M^{\leq k}, \forall x\in\mathbb B^n, ((f/P)(x)=x \wedge x\in\operatorname{TS}_P(z))\implies x\models M\enspace. \end{equation}
As with the previous case, in practice we may also want that there exists at least one fixed point reachable from
With the BoNesis Python interface, this reprogramming property is declared as follows, where f
is a BN, z
the initial configuration (Python dictionary), M
the marker, and k
the maximum number of components that can be perturbed (at most
def source_marker_reprogramming_fixpoints(f: BooleanNetwork,
z: dict[str,bool],
M: dict[str,bool],
k: int):
# f: Boolean network; z: initial configuration;
# M: marker; k: maximum number of components to perturb
bo = bonesis.BoNesis(f)
P = bo.Some(max_size=k) # perturbation to identify
with bo.mutant(P):
# all the fixed points reachable from z match with M
~bo.obs(z) >> "fixpoints" ^ {bo.obs(M)}
# at least one fixed point matching with M is reachable from z
~bo.obs(z) >= bo.fixed(~bo.obs(M))
return P.assignments()
Compared to the previous code snippet, this function relies on specific operators to restrict the properties to the fixed point reachable from z
.
The instruction ~bo.obs(z)
refers to a specific configuration matching with z
;
~bo.obs(z) >> "fixpoints" ^ {bo.obs(M)}
specifies that all the fixed points reachable from such configuration have to match with at least one configuration given in the set {bo.obs(M)}
, i.e., M
in this case.
This property is satisfied whenever no fixed point are reachable. Thus, the next line ensures that at least one fixed point is reachable from the configuration associated with z
: bo.fixed(~bo.obs(M))
refers to one configuration which is a fixed point (in the perturbed BN), and which matches with M
. Then, the binary operator >=
declares the existence of a trajectory from its left to its right configuration.
Notice that with this formulation, in the case whenever z
is only partially defined (some components are not associated to a Boolean value), a perturbation is returned as long as there exists at least one fully-defined configuration matching with
The corresponding command line is of the form
bonesis-reprogramming model.bnet M k --fixpoints --reachable-from z
where model.bnet
is a BN encoded in BooleanNet format, M
specifies the marker as a JSON map, k
is the maximum number of perturbations, and z
is the initial configuration as a JSON map.
The existence of at least one fixed point can be lifted with the option --allow-no-fixpoint
.
+++
Example. Let us consider the following toy BN with two positive feedback cycles:
---
nbconvert_latex:
nbconvert.figure.caption: Influence graph of the example BN for P2
---
f = BooleanNetwork({
"A": "B",
"B": "A",
"C": "!D & (A|B)",
"D": "!C"
})
f.influence_graph()
This BN has 3 fixed points, 2 of which are reachable from the configuration where A
and B
are active, and C
and D
inactive:
---
nbconvert_latex:
nbconvert.figure.caption: Fully-asynchronous dynamics of the example BN for P2 from
the configuration $1100$
---
z = {"A": 1, "B": 1, "C": 0, "D": 0}
f.dynamics("asynchronous", init=z)
list(mpbn.MPBooleanNetwork(f).fixedpoints())
list(mpbn.MPBooleanNetwork(f).fixedpoints(reachable_from=z))
Let us compare the results of the global marker-reprogramming of fixed points (P1) with the source-marker reprogramming of fixed points (P2), the objective being to have fixed points having C
active.
In the first case, putting aside the perturbation of C
, this necessitates to act on either A
or B
to prevent the existence of the fixed points where A
, B
and C
are inactive:
list(marker_reprogramming_fixpoints(f, {"C": 1}, 2))
Considering only the fixed points reachable from the configuration z
, there is no need to act on A
or B
:
list(source_marker_reprogramming_fixpoints(f, z, {"C": 1}, 2))
The command line program bonesis-reprogramming
can perform P2 by specifying the --rechable-from
option giving the initial configuration in JSON format:
with open("example2.bnet", "w") as fp:
fp.write(f.source())
:nbconvert.display: true
!bonesis-reprogramming example2.bnet '{"C": 1}' 2 --fixpoints \
--reachable-from '{"A": 1, "B": 1, "C": 0, "D": 0}'
We identify the perturbations
\begin{equation} \exists P\in\mathbb M^{\leq k}, \forall x\in\mathbb B^n, \operatorname{IN-ATTRACTOR}_P(x) \implies x\models M \end{equation}
("There exists a perturbation
By restricting the range of the universal part of the equation to the configurations which do not match with the marker
The
\begin{equation} \exists P\in\mathbb M^{\leq k}, \forall x\in\mathbb B^n: x\not\models M, \exists y\in\mathbb B^n, y\in \operatorname{TS}_P(x), \operatorname{TS}_P(y) \neq \operatorname{TS}_P(x) \end{equation}
The problem of satisfiability of this quantified Boolean expression is beyond the expressiveness power of ASP which is limited to
Because the domain of candidate perturbations
Notice that this approach is highly combinatorics, and is likely limited to identifying perturbations of small size. However, to our knowledge, this is the first implemented method which addresses the complement reprogramming of MP attractors, i.e., of the minimal trap spaces of the BN.
With the BoNesis Python interface, this reprogramming property is declared as follows, where f
is a BN, M
the marker, and k
the maximum number of components that can be perturbed (at most
def marker_reprogramming(f: BooleanNetwork,
M: dict[str,bool],
k: int):
bo = bonesis.BoNesis(f)
coP = bo.Some(max_size=k)
with bo.mutant(coP):
x = bo.cfg()
bo.in_attractor(x)
x != bo.obs(M)
return coP.complementary_assignments()
The idea of the above code snippet is to declare the property of being a bad perturbation (coP
). In the context of this perturbation, we declare with x=bo.cfg()
a configuration. Then, bo.in_attractor(x)
imposes that x
belongs to an attractor (minimal trap space) of the perturbed BN; and the x != bo.obs(M)
instruction adds the constraint that x
must not match with M
.
The .complementary_assignments()
method performs the computation of the complement of solutions. It solves the above problem with perturbation sizes ranging from 0 to k
, and returns complement minimal perturbations only.
The corresponding command line is of the form
bonesis-reprogramming model.bnet M k
where model.bnet
is a BN encoded in BooleanNet format, M
specifies the marker as a JSON map, k
is the maximum number of perturbations.
+++
Example. Let us consider the following BN:
---
nbconvert_latex:
nbconvert.figure.caption: Influence graph of the example BN for P3
---
f = mpbn.MPBooleanNetwork({
"A": "!B",
"B": "!A",
"C": "A & !B & !D",
"D": "C | E",
"E": "!C & !E",
})
f.influence_graph()
Essentially, A
and B
always stabilize to opposite states. Whenever A
is active (and B
inactive) then C
will oscillate, otherwise it stabilizes to 0. In each case D
and E
oscillate.
This lead to the following MP attractors:
tabulate(list(f.attractors()))
Let us say that our objective is to reprogram the BN such that all the attractors of the component C
fixed to 1.
The reprogramming of fixed points (P1) gives the following solutions:
list(marker_reprogramming_fixpoints(f, {"C": 1}, 3))
Putting aside the trivial solution of perturbing C
, let us analyze the BN perturbed with the D
forced to 0:
pf = f.copy()
pf["D"] = 0
tabulate(pf.attractors())
The (only) fixed point of the network indeed has C
active. However, it possesses another (cyclic) attractor, where C
is inactive.
This example points out that focusing on fixed point reprogramming may lead to predicting perturbations which are not sufficient to ensure that all the attractors show the desired marker.
The complete attractor reprogramming returns that the perturbation of D
must be coupled with a perturbation of A
or B
, in this case to destroy the cyclic attractor.
list(marker_reprogramming(f, {"C": 1}, 3))
The same results can be obtained using the command line as follows.
with open("example3.bnet", "w") as fp:
fp.write(f.source())
:nbconvert.display: true
!bonesis-reprogramming example3.bnet '{"C": 1}' 3
In other cases, the reprogramming of attractors may also lead to fewer required perturbations than focusing on fixed points, provided we enforce the existence of at least one fixed point. This can be illustrated with the following example:
---
nbconvert_latex:
nbconvert.figure.caption: Influence graph of the example BN $g$ for P3
---
g = mpbn.MPBooleanNetwork({
"A": "!B",
"B": "A",
"C": "A & B",
"D": "C"
})
g.influence_graph()
Unperturbed, this network has a single cyclic attractor, as A
and B
are in a sustained oscillation.
tabulate(g.attractors())
Enforcing that all attractors have D
fixed to 1 can be achieved either by perturbing A
, or by perturbing only C
, letting A
and B
oscillate.
list(marker_reprogramming(g, {"D": 1}, 2))
However, besides the forced activation of A
, ensuring that all and at least one fixed point have D
active always requires perturbing at least two components.
list(marker_reprogramming_fixpoints(g, {"D": 1}, 2))
Given an initial configuration
The associated decision problem can be expressed as follows:
\begin{equation} \exists P\in\mathbb M^{\leq k}, \forall x\in\mathbb B^n, x\models M \vee x\notin\bar\rho^{f/P}_{\mathrm{mp}}(z) \vee \neg\operatorname{IN-ATTRACTOR}_P(x) \end{equation}
("There exists a perturbation
By integrating the definition of the
The resolution of the problem can be approached in a very similar way to P3, i.e., by solving its complement. The equation is almost the same, with the addition that
With the BoNesis Python interface, this reprogramming property is declared as follows, where f
is a BN, z
the initial configuration, M
the marker, and k
the maximum number of components that can be perturbed (at most
def source_marker_reprogramming(f: BooleanNetwork,
z: dict[str,bool],
M: dict[str,bool],
k: int):
bo = bonesis.BoNesis(f)
coP = bo.Some(max_size=k)
with bo.mutant(coP):
x = bo.cfg()
bo.in_attractor(x)
x != bo.obs(M)
~bo.obs(z) >= x
return coP.complementary_assignments()
The above code snippet is very similar to the previous marker_reprogramming
, with the addition of the ~bo.obs(z) >= x
instruction which declares that x
, a configuration which belongs to an attractor of the perturbed BN and which does not match with M
, is reachable from z
.
+++
The corresponding command line is of the form
bonesis-reprogramming model.bnet M k --reachable-from z
where model.bnet
is a BN encoded in BooleanNet format, M
specifies the marker as a JSON map, k
is the maximum number of perturbations, and z
is the initial configuration as a JSON map.
+++
Example. Let us consider again the BN f
analyzed in the previous section. By focusing only on attractors reachable from the configuration where A
is fixed to 1 and other nodes to 0, the reprogramming required to make all attractors have C
fixed to 1 consists only of fixing D
to 0. Note that in the specific example, the reprogramming of reachable fixed point would give an equivalent result.
z = f.zero()
z["A"] = 1
list(source_marker_reprogramming(f, z, {"C": 1}, 3))
The same results can be obtained using the command line as follows.
:nbconvert.display: true
!bonesis-reprogramming example3.bnet '{"C": 1}' 3 \
--reachable-from '{"A": 1, "B": 0, "C": 0, "D": 0}'
In the previous section, the reprogramming was performed on a single BN, by giving to BoNesis the singleton domain of BNs to consider. As described in the Method section, BoNesis can reason on ensembles of BNs, either specified explicitly, or implicitly by an influence graph. The functions defined above can then be directly applied to such ensembles of BNs. In this section, we briefly discuss how the resulting reprogramming solutions should then be interpreted with respect to these ensembles.
Given a domain of BNs
In the case of our implementation of marker reprogramming of fixed points (P1 and P2), the expression becomes of the form:
In the case of our implementation of marker reprogramming of attractors (P3 and P4), because we tackle the complementary problem, the expression becomes of the form:
+++
Let us illustrate the ensemble reprogramming with the following example. First, let us define an influence graph to delimit the domain of admissible BNs:
---
nbconvert_latex:
nbconvert.figure.caption: Influence graph delimiting the ensembles of BNs in the
example
---
dom = bonesis.InfluenceGraph([
("C", "B", {"sign": 1}),
("A", "C", {"sign": 1}),
("B", "C", {"sign": -1}),
("C", "D", {"sign": 1}),
], exact=True, canonic=False) # we disable canonic encoding
dom
This domain encloses all the BNs having exactly (exact=True
) the specified influence graph, 4 distinct BNs in this case:
dom.canonic = True # we set canonic encoding for enumerating BNs
F = list(bonesis.BoNesis(dom).boolean_networks())
dom.canonic = False
pd.DataFrame(F)
Let us explore the attractors of each individual BNs:
:nbconvert.display: true
for i, f in enumerate(F):
print(f"Attractors of BN {i}:", list(f.attractors()))
In this example, we focus on reprogramming the attractors so that the component D
is fixed to 1.
On the one hand, when reprogramming fixed points only, because one BN already verifies this property, the empty perturbation is a solution:
list(marker_reprogramming_fixpoints(dom, {"D": 1}, 2))
On the other hand, the reprogramming of attractors returns solution that work on every BN:
list(marker_reprogramming(dom, {"D": 1}, 2))
Indeed, fixed C
to 1, ensures in each case that D
is fixed to 1.
+++
The computation of universal solutions for the reprogramming of fixed points can be tackled by following a similar encoding than the reprogramming of attractors, i.e., by identifying perturbations which do not fulfill the property for at least one BN in the domain (the complement results in perturbations working for all the BNs):
def universal_marker_reprogramming_fixpoints(f: BooleanNetwork,
M: dict[str,bool],
k: int):
bo = bonesis.BoNesis(f)
coP = bo.Some(max_size=k)
with bo.mutant(coP):
x = bo.cfg()
bo.fixed(x) # x is a fixed point
x != bo.obs(M) # x does not match with M
return coP.complementary_assignments()
list(universal_marker_reprogramming_fixpoints(dom, {"D": 1}, 2))
Note that in this implementation, we do not ensure the existence of a fixed point after reprogramming. This is why the perturbation fixing only A
to 1 is considered as a solution in our example:
:nbconvert.display: true
for i, f in enumerate(F):
f["A"] = 1
print(f"Attractors of BN {i} after fixing A to 1:", list(f.attractors()))
As BNs 2 and 3 have no fixed point, they fulfill the criteria "all the fixed points match with marker M
".
+++
In order to evaluate the scalability on realistic BNs, we use the benchmark constituted by Moon et al. (2022a) to evaluate the reprogramming of fixed points (P1). Their benchmark gathers 10 locally-monotone BNs and 1 non-monotone one, that BoNesis cannot address. The dimension of the 10 BNs are respectively 14, 17, 18, 20, 28, 32, 53, 59, 66, and 75. For each of these models, a target marker for reprogramming has been defined from the corresponding published studies. Moreover, a subset of nodes has been declared as uncontrollable to avoid trivial solutions. We used this benchmark to evaluate the scalability of the P1 and P3 implementation we propose in this paper.
For these 10 models, we applied the P1 and P3 reprogramming for different maximum number of simultaneous perturbations (denoted
In the case of P1, with
In the case of P3, with
These experiments testify of the difference of complexity between P1 and P3, and more precisely on the resolution approach taken for P3: the computation of the complementary sets of solutions becomes rapidly intractable for large combinations of perturbation. Indeed, in practice, there are many bad perturbations, thus their enumeration, which is necessary to compute their complement, is an important bottleneck.
Evaluating the scalability of P2 and P4 would require defining initial configurations which are meaningful for the different models, which are not available in the selected benchmark. Nevertheless, because the source constraint does not change the complexity classes, we can conjecture that their scalability should be comparable to P1 and P3 respectively. Moreover, having benchmarks at larger scale would be insightful, but none of them are available to the best of our knowledge. It should be noted BoNesis has been applied to do BN synthesis for models with 1,000 nodes (Chevalier et al., 2019), suggesting a potential applicability of BoNesis for the reprogramming of large BNs.
As stressed in the introduction, there exists only tools addressing P1 to compare with. The experiments of Moon et al. (2022a) show that their bilevel integer programming-based method systematically outperforms the ASP implementation of pyActoNet (Biane et al., 2018). On the same benchmark, BoNesis performed either similarly or in shorter time, albeit limited to locally-monotone BNs only.
+++
In this paper, we demonstrated how the BoNesis Python library can be employed to fully characterize permanent perturbations which guarantee that all the fixed points or all minimal trap spaces of the perturbed BN have a subset of components fixed to desired values (marker). We focused on reprogramming for achieving elementary dynamical properties, that are the fixed points or attractors, optionally reachable from a given configuration. Nevertheless, the snippets shown in this paper can be extended to account for more complex or specific dynamical properties after mutation, e.g., existence of additional trajectories, considering multiple initial configurations.
It should be noted that the candidate combinations of perturbations are computed solely based on the Boolean dynamics, and do not account for experimental feasibility, e.g., in the scope of models of biological systems.
Future work may consider optimization or prioritization of perturbations based on such extra information.
Currently, BoNesis enables specifying uncontrollable components which must not be perturbed (exclude
option for the Some
object, or --exclude
for the command line, taking a list of components which should be excluded from the candidate perturbations).
We considered for problems, referred to as P1, P2, P3, P4, where P1-P2 relate to the reprogramming of fixed points, and P3-P4 to the reprogramming of MP attractors (i.e., minimal trap spaces). The computational complexity of P1-P2 allows an efficient implementation using Answer-Set Programming (ASP), whereas the one of P3-P4 necessitate working around complementary solution to fit into the expressiveness of ASP, limiting their scalability. Future work may explore alternative implementations using different logic frameworks.
The identified perturbations may destroy and create new fixed points and attractors for the BN. This is a significant difference with most of the methods developed in the literature where many focus on the control towards attractors of the unperturbed BNs only. Whereas the problem P1 which has been already addressed with different methods, we are not aware of any other approach tackling P2, P3, and P4.
Besides the four reprogramming problems tackled in this paper, an additional variant would be the marker-reprogramming of fixed points which also ensures the absence of cyclic attractors. Note that its complexity is equivalent to the one of P3/P4, i.e., it can be expressed as a
This paper focused on permanent perturbations, i.e., enforcing the value of one or several components constantly over time, independently of the state of the system.
Sequential reprogramming (Mandon et al., 2017,Pardo et al., 2021) consists in applying sets of perturbations at different time. This can lead to reducing the overall number of component to perturb, by taking advantage of the natural transient dynamics of the system.
Sequential reprogramming brings the BN reprogramming settings closer to classical control theory, as the
control can depend both on time and state of the system.
Having fixed a number of steps, say
"There exist two sets of perturbations
Finally, we demonstrated how BoNesis can be employed to reason on the reprogramming of BNs, leading to either solutions that work for at least one BN of the ensemble, or working on each of them (universal reprogramming). We believe that reasoning on ensemble of models is a promising direction to address the robustness of predictions in the scope of applications in systems biology.
+++
The author would like to thank Kyungduk Moon and Kangbok Lee for stimulating discussions and for providing the model and configuration files for the benchmarks of Moon et al. (2022a).
The software BoNesis is available at github.com/bnediction/bonesis. The code of this paper uses version 0.4.93 (archived at doi:10.5281/zenodo.7657487). The notebook, models, and instructions for reproducing the results are available at github.com/bnediction/reprogramming-with-bonesis (archived at doi:10.5281/zenodo.7733095).
This work has been supported by the French Agence Nationale pour la Recherche (ANR) in the scope of the project BNeDiction (grant number ANR-20-CE45-0001).
Footnotes
-
this problem is actually in NP when allowing a number of variables quadratic with $n$ ↩