Marker and source-marker reprogramming of Most Permissive Boolean networks and ensembles with BoNesis

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 ﬁxed points and attractors. We consider marker properties, which specify that some components are ﬁxed to a speciﬁc value. We study 4 variants of the marker reprogramming problem: the reprogramming of ﬁxed points, of minimal trap spaces, and of ﬁxed points and minimal trap spaces reachable from a given initial conﬁguration with the most permissive update mode. The perturbations consist of ﬁxing a set of components to a ﬁxed 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 modiﬁed interactively.

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é, Kolčák, Chatain, and Haar, 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., A = 1, C = 0). After reprogramming, all the configurations in all (reachable) attractors must be compatible with these markers. Importantly, the target attractors are not necessarily attractors of the original (wild-type) BN: the reprogramming can destroy and create new attractors. In particular, if there is no attractor in the original model matching with the marker, the reprogramming will identify perturbations that will create such an attractor and ensure its reachability. This is a substantial difference with many of the methods in the literature. Moreover, the approach we present here can return all the possible solutions, possibly up to a given maximum number of perturbations to apply, and possibly avoiding uncontrollable nodes. 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 and Delaplace (2019) with the ActoNet method and by Moon, Lee, Chopra, and Kwon (2022), based on bilevel integer programming. To our knowledge, none of the other problems have been addressed completely in the literature.
The software BoNesis (https://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, Froidevaux, Paulevé, and Zinovyev, 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 (https://en.wikipedia.org/ wiki/Unate_function), 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 (https://nbviewer.org/ Loïc Paulevé 3 github/bnediction/reprogramming-with-bonesis/blob/release/paper.ipynb) and interactively executed either online (https://mybinder.org/v2/gh/bnediction/reprogramming-with-bonesis/release?urlpath=tree/paper.ipynb) using mybinder online free service, or locally, following instructions given later in this paper. A BN f is locally monotone whenever every of its local functions are unate: for each i ∈ {1, · · · , n}, there exists an ordering of components i ∈ {≤, ≥} n such that Intuitively, a BN is locally monotone whenever each of its local function can be expressed in propositional logic such that each variable appears either never or always with the same sign.
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 locallymonotone BNs.

Update modes.
Given a BN f and a configuration x, the update mode specifies how to compute the next configuration. There is a vast zoo of update modes (Paulevé and Sené, 2022), but traditionally, two modes are usually considered in biological modeling: the synchronous (or parallel) deterministic mode, where the next configuration is given by its application to f (x is succeeded by f (x)), and the fully asynchronous (often denoted only asynchronous) where the next configuration results from the application of only one local function, chosen non-deterministically.
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é, Kolčák, Chatain, and Haar, 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" ( ) and "decrease" ( ): when a component can be activated, it will first go through the "increase" state where it can be interpreted as either 0 or 1 by the other components, until eventually reaching the Boolean 1 state; and symmetrically for deactivation. A formal definition of MP dynamics is given later in this section.

Perturbations.
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, P = {2 → 0, 4 → 1}. Given a perturbation P, the perturbed BN f /P is given by, for each component i ∈ {1, ... , n}:
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 Q 1 · · · Q m but ignoring repetitions: an ∃∃∀∀∀∃∃∀-expression has the same decision complexity as an ∃∀∃∀-expression.
Computational complexity (Papadimitriou, 1995) 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 ∃-expressions is the infamous (Boolean) SAT(isfiability) problem, which is NP-complete: it can be solved by a non-deterministic polynomial time algorithm, and it is among the hardest problems in this class: any problem in NP can be (efficiently) transformed into a SAT problem. In practice, our computers being deterministic, the resolution of the SAT problems employs algorithms running in worst-case time and space exponential with the number of variables in the Boolean expression. However, modern SAT solvers can approach expressions with thousands to millions of variables.
The decision of satisfiability of ∀-expressions can be seen as a complementary problem to ∃expression: ∀X φ(X ) is satisfiable if and only if ∃X ¬φ(X ) is not satisfiable: it is a coNP-complete problem. It is not known whether coNP = NP.
Then, the alternation of quantifiers makes the problem climbing into the so-called polynomial hierarchy (see https://en.wikipedia.org/wiki/Polynomial_hierarchy). The ∃ ...-expressions are Σ P kcomplete problems, where k is the number of alternating quantifiers (starting with ∃), while ∀ ...expressions are Π P k -complete (Σ P 1 =NP and Π P 1 =coNP). It is not known yet whether all these complexity classes are equal, but in practice, algorithms of resolution scale rapidly down with their height in the polynomial hierarchy. Each of these complexity classes are included in PSPACE, the class of problems solvable in polynomial space. PSPACE-complete problems, such as the verification of properties of asynchronous BNs, are known to be difficult to tackle in practice (currently limited to a couple of hundreds of variables in the case of BNs).
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.

Elementary dynamical properties and their complexity
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é, Kolčák, Chatain, and Haar, 2020 is also one of its vertices (h is closed by f ). In particular, the (sub-)hypercube * n is always a trap space.
A sub-hypercube h is smaller than a sub-hypercube h , denoted by h h whenever c(h) ⊆ c(h ). Equivalently, this means that each non-free component of h is fixed to the same value in 1.3.2. MP attractors are minimal trap spaces. The attractors of MP dynamics are the minimal trap spaces of the Boolean function f (Paulevé, Kolčák, Chatain, and Haar, 2020), i.e., the trap spaces which do not include strictly smaller trap spaces. Thus, we denote MP attractors by subhypercubes, i.e., an MP attractor A is a vector in {0, 1, * } n . Therefore, a component with a * value in an MP attractor A indicates that the component that can always oscillate between 0 and 1 in the (cyclic) attractor.
The computational complexity of decision problems related to minimal trap spaces has been extensively addressed in (Moon, Lee, and Paulevé, 2022) with different representations of BNs 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 coNP coNP -complete problem, i.e., equivalent to the decision of satisfiability of ∀∃ expressions. In the case of locally-monotone BNs, deciding whether a subhypercube is a trap space is in PTIME, whereas deciding whether it is a minimal trap spaces is a coNP-complete problem, i.e., equivalent to the decision of satisfiability of ∀-expressions.

MP reachability of attractors. Given a configuration x and an MP attractor
there is an MP trajectory from x to any configuration y ∈ A if and only if A is smaller than the smallest trap space containing x. We write reach(x, y ) the existence of such a trajectory.
Let us denote by TS ( 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 PTIME NP in general (this problem is actually in NP when allowing a number of variables quadratic with n), and PTIME in the locallymonotone 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.

Belonging to an MP attractor.
In the following, we will consider the problem of deciding whether a given configuration x belongs to an MP attractor of f . We write IN-ATTRACTOR(x) such a property. This can be verified in two steps: (1) compute the smallest trap spaces containing x, noted TS(x), and (2) verify whether TS(x) is a minimal trap space. This later property is true if and only if for any vertex y of TS(x), the minimal trap space containing y is equal to TS(x):

IN-ATTRACTOR(x) ≡ ∀y ∈ c(TS(x)), TS(y ) = TS(x) .
Finally, given a set of perturbations P, we write TS P (x) for the small trap space of perturbed BN (f /P) containing x, and IN-ATTRACTOR P (x) the property of x belonging to an attractor of the perturbed BN (f /P).

BoNesis
BoNesis (Paulevé, 2023a) 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 F, and (2) a set of Boolean dynamical properties φ, and can enumerate the BNs f ∈ F such that f |= φ, i.e., f verifies the properties φ.
Currently, the domain of BNs F can be one of the following: 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 F = {f 1 , · · · , f m }.
• Any locally-monotone BN matching with a given influence graph G: An influence graph is a signed digraph between components of the form • Any locally-monotone BN matching with a partially-defined BN following the AEON framework (Beneš, Brim, Pastva, and Šafránek, 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 (https://potassco.org/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 ∃and ∃∀-expressions, thus having an expressiveness higher than classical SAT.
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.

Results
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 https://doi.org/10.5281/ zenodo.7733095. Its execution requires the Jupyter notebook system, Python, and the Python package bonesis to be installed. Alternatively, the notebook can be executed within the CoLo-MoTo Docker distribution (Naldi et al., 2018)  Alternatively, the computation of reprogramming perturbations from single BNs can be performed using the command line program bonesis-reprogramming, provided alongside the bonesis Python package. We detail its usage in each case.

Marker reprogramming of Boolean networks
We first consider the reprogramming of a single BN f of dimension n. In the framework of BoNesis, this means the domain of BNs is the singleton F = {f }.
A marker M is a map associating a subset of components of f to a Boolean value. For instance, Given k ∈ N, we denote by M ≤k the sets of maps associating at most k components among {1, · · · , n} to a Boolean value.
The objective of marker-reprogramming is to identify perturbations so that all the attractors of the perturbed f match with the marker M. The source-marker reprogramming then focuses on the attractors reachable from a given initial configuration only, thus potentially requiring fewer perturbations.
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: 1. Marker reprogramming of fixed points (P1); 2. Source-marker reprogramming of fixed points (P2); 3. Marker reprogramming of attractors (P3); 4. 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: is the command line bonesis-reprogramming model.bnet M k, with model.bnet the path to a file in BooleanNet format (http://colomoto.org/biolqm/doc/format-bnet.html), 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 P of at most k components so that all the fixed points of f /P match with the given marker M. The associated decision problem can be expressed as the following ∃∀-expression, hence being at most in Σ P 2 : (1) There exists a perturbation being a map of at most k components to a Boolean value, such that for all configurations x ∈ B n , if x is a fixed point of the perturbed BN (f /P), then x matches with the marker M".
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: 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 n): 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 Boolean-Net 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.  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: The resulting graphics is reproduced in Figure 1.
With the (fully) asynchronous update mode, the system has a single attractor, consisting of all the configurations of the network.
Recall that the fixed points are identical in asynchronous and MP. We use mpbn (https://github .com/bnediction/mpbn) to analyze the dynamical properties with the MP update mode: 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: 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. 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.

Source-marker reprogramming of fixed points (P2).
Given an initial configuration z, we identify the perturbations P of at most k components so that all the fixed points of f /P that are reachable from z in f /P match with the given marker M. The associated decision problem can be expressed as the following ∃∀-expression, hence being at most in Σ P 2 : (3) There exists a perturbation P such that for any configuration x ∈ B n , if x is a fixed point of the perturbed BN (f /P), and x is reachable from z in (f /P), then x must match with M".
As explained in the Method section, the reachability property boils down to computing the smallest trap space containing z: if it contains the fixed point x, then x is reachable from z with the MP update mode.
As with the previous case, in practice we may also want that there exists at least one fixed point reachable from z.
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 n): 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 z which fulfil the specified dynamical properties.
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:  The resulting graphics is reproduced in Figure 3. 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: In [15]: z = {"A": 1, "B": 1, "C": 0, "D": 0} f.dynamics("asynchronous", init=z) The resulting graphics is reproduced in Figure 4. 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: The program bonesis-reprogramming can perform P2 by specifying the --rechable-from option giving the initial configuration in JSON format: In

Marker reprogramming of attractors (P3).
We identify the perturbations P of at most k components so that the configurations of the all the attractors of f /P match with the given marker M (i.e., in each attractor, the specified markers cannot oscillate). The associated decision problem can be expressed as follows: ("There exists a perturbation P of at most k components, such that for all configurations x, if x belongs to an attractor of the perturbed BN f /P, then x matches with the specified markers M") By restricting the range of the universal part of the equation to the configurations which do not match with the marker M, we obtain: The IN-ATTRACTOR property being itself a quantified Boolean expression, we obtain the following ∃∀∃-expression: ∃P ∈ M ≤k , ∀x ∈ B n : x |= M, ∃y ∈ B n , y ∈ TS P (x), TS P (y ) = TS P (x) The problem of satisfiability of this quantified Boolean expression is beyond the expressiveness power of ASP which is limited to ∃∀-expressions. Nevertheless, we can approach this problem by its complementary: the existence of perturbations of size k such that at least one configuration belonging to an attractor does not match with the marker M. This complementary problem can be expressed with this following expression which is an ∃∀-expression: (9) ∃P ∈ M ≤k , ∃x ∈ B n , x |= M ∧ ∀y ∈ B n , y ∈ TS P (x) =⇒ TS P (y ) = TS P (x) .
Because the domain of candidate perturbations M ≤k is finite, one can first resolve the complementary problem, giving all bad perturbations, and returns its complement.
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 n): 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: The resulting graphics is reproduced in Figure 5.
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: In [24]: tabulate(list(f.attractors())) Out [24]: A B C D E 0 0 1 0 * * 1 1 0 * * * 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: 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. The same results can be obtained using the command line as follows.
In [28]: with open("example3.bnet", "w") as fp: fp.write(f.source()) In [ 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: The resulting graphics is reproduced in Figure 6. Unperturbed, this network has a single cyclic attractor, as A and B are in a sustained oscillation.
In [31]: tabulate(g.attractors()) Out[31]: A B C D 0 * * * * 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.

Soure-marker reprograming of attractors (P4).
Given an initial configuration z, we identify the perturbations P of at most k components so that the configurations of the all the attractors of f /P that are reachable from z match with the given marker M (i.e., in each reachable attractor, the specified markers cannot oscillate). Thus, P4 is the same problem as P3, except that we focus only on attractors reachable from z, therefore potentially requiring fewer perturbations. The associated decision problem can be expressed as follows: (10) ∃P ∈ M ≤k , ∀x ∈ B n , x |= M ∨ x / ∈ρ f /P mp (z) ∨ ¬ IN-ATTRACTOR P (x) ("There exists a perturbation P of at most k components, such that for all configurations x, either x matches with the marker M, or x does not belong to an attractor, or x is not reachable from z").

Loïc Paulevé 17
Peer Community Journal, Vol. 3 (2023), article e30 https://doi.org/10.24072/pcjournal.255 By integrating the definition of the IN-ATTRACTOR property, we obtain the following ∃∀∃expression: 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 x must be reachable from z, leading to the ∃∀-expression: (12) ∃P ∈ M ≤k , ∃x ∈ B n , x ∈ TS P (z) ∧ x |= M ∧ ∀y ∈ B n , y ∈ TS P (x) =⇒ TS P (y ) = TS P (x) .
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 n): 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. The same results can be obtained using the command line as follows. In

Reprogramming of ensembles of Boolean networks
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 F, BoNesis returns a solution whenever at least one BN of this domain satisfies the given properties. Intuitively, this means that the logic satisfiability problem is of the form ∃f ∈ F, Φ(f ). As detailed in (Chevalier, Froidevaux, Paulevé, and Zinovyev, 2019) in the scope of locally-monotone BNs, the size of the "f ∈ F" formula is, in general, exponential (binomial coefficient) with the in-degree of nodes in the influence graph. This complexity is due to the maximum number of clauses a Boolean function can have. Our encoding in BoNesis allows specifying an upper bound to this number, which enables tackling very large scale instances although giving access only to a subset of F. The encoding of F in BoNesis also supports enforcing a canonic representation of BNs in order to offer a non-redundant enumeration of the BNs, at the price of a quadratic size of the formula. However, in our case, as we are only interested in enumerating the perturbations, the canonic form of BNs is not needed.
In the case of our implementation of marker reprogramming of fixed points (P1 and P2), the expression becomes of the form: ∃f ∈ F, ∃P ∈ M ≤k , · · · Therefore, a perturbation P is returned as soon as it is a reprogramming solution for at least one BN of the input domain: P may not work on every BN in F, but at least one.
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: ∃f ∈ F, ∃coP ∈ M ≤k , · · · Therefore, a bad perturbation coP is returned as soon as it is a bad perturbation for at least one BN of the input domain. By taking the complement of these perturbations (in M ≤k ), we obtain that the returned perturbations are reprogramming solutions for all the BNs in F. Let us illustrate the ensemble reprogramming with the following example. First, let us define an influence graph to delimit the domain of admissible BNs: 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: 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): 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:  As BNs 2 and 3 have no fixed point, they fulfill the criteria "all the fixed points match with marker M".

Scalability
In order to evaluate the scalability on realistic BNs, we use the benchmark constituted by Moon, Lee, Chopra, and Kwon (2022) 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 (Paulevé, 2023b).
For these 10 models, we applied the P1 and P3 reprogramming for different maximum number of simultaneous perturbations (denoted k in the previous sections), up to k = 6. In each case, we measured the time for the first solution, for listing up to 100 solutions, and for listing all the solutions, with a timeout of 10 minutes. The experiments have been performed on an Intel(R) Xeon(R) processor at 3.3Ghz with 16BG of RAM. In the case of P1, with k = 6, it took around 1s to get at least one solution for each of the 10 models; up to 100 solutions have been listed in the same timing, except for one model which took 8s. The full listing of solutions of 3 of the larger models have timed out, the rest necessitating between 1 and 18s. With k = 4, BoNesis was able to list all the solutions for all models (up to 5min for one of the larger model).
In the case of P3, with k ≥ 4, 3 of the 10 models could not find a single solution in the given time limit; for most of the other models, a first solution was found in around 1s, a couple of models took around 1-2min. The enumeration of the first 100 solutions took a similar time with k = 4, but timed out with k = 6 for all but the 4 smallest models. With k ≤ 2, BoNesis has found all the solutions to P3 for all the 10 models in a few seconds maximum.
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, Froidevaux, Paulevé, and Zinovyev, 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, Lee, Chopra, and Kwon (2022) show that their bilevel integer programmingbased method systematically outperforms the ASP implementation of pyActoNet (Biane and Delaplace, 2019). On the same benchmark, BoNesis performed either similarly or in shorter time, albeit limited to locally-monotone BNs only.

Discussion
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 ∃∀∃expression. This problem may be relevant for modeling cases where cyclic attractors do not make sense. The programming interface of BoNesis do not permit an efficient encoding of this problem at the moment. 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, Haar, and Paulevé, 2017;Pardo, Ivanov, and Delaplace, 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 m, the reprogramming problems consists in identifying m sets of perturbations which will be applied in sequence, and their application may be restricted to attractors only (Mandon, Su, Haar, Pang, and Paulevé, 2019). Interestingly in that case, having fixed the number of reprogramming steps, the computational complexity remains identical to the one-step reprogramming with locally-monotone BNs, due to the PTIME complexity of the reachability. For instance, the 2steps reprogramming of BNs along fixed points only with the MP update mode can be expressed as the following ∃∀-expression, as P1: ∃P, Q ∈ M ≤k , ∀x, y ∈ B n , f P (x) = x ⇒ reach Q (x, y ) ⇒ f Q (y ) = y ⇒ y |= M "There exist two sets of perturbations P and Q, such that for any configuration x and any configuration y , if x is a fixed point of under the perturbation P, then if y is reachable from x under the perturbation Q, then if y is a fixed point under the perturbation Q, then it must match with the marker M". The more general 2-steps reprogramming along attractors can be expressed as follows: 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.