Parallel One-Step Control of Parametrised Boolean Networks

Boolean network (BN) is a simple model widely used to study complex dynamic behaviour of biological systems. Nonetheless, it might be difficult to gather enough data to precisely capture the behavior of a biological system into a set of Boolean functions. These issues can be dealt with to some extent using parametrised Boolean networks (ParBNs), as it allows to leave some update functions unspecified. In this paper, we attack the control problem for ParBNs with asynchronous semantics. While there is an extensive work on controlling BNs without parameters, the problem of control for ParBNs has not been in fact addressed yet. The goal of control is to ensure the stabilisation of a system in a given state using as few interventions as possible. There are many ways to control BN dynamics. Here, we consider the one-step approach in which the system is instantaneously perturbed out of its actual state. A naive approach to handle control of ParBNs is using parameter scan and solve the control problem for each parameter valuation separately using known techniques for non-parametrised BNs. This approach is however highly inefficient as the parameter space of ParBNs grows doubly-exponentially in the worst case. In this paper, we propose a novel semi-symbolic algorithm for the one-step control problem of ParBNs, that builds on a symbolic data structures to avoid scanning individual parameters. We evaluate the performance of our approach on real biological models.


Introduction
Cell reprogramming is currently one of the most critical challenges in computational biology. The goal of cell reprogramming is to control a cell's phenotype. This ability opens many opportunities, mainly in regenerative medicine. In order to reach the desired phenotype, the correct transcription factors must be identified. That is close to impossible to be done only using in vitro biological experiments due to the very high number of possibilities how the cell might be interfered with. This is where in silico analysis and computational models of cell dynamics come into play. Formal methods and their integration provide a promising technology that allows fully automatic identification of control strategies by using computational models.
A cell can be viewed as a set of genes and their mutual regulators. The compact abstraction of these relationships can be modelled using Boolean networks (BNs). BNs are becoming very popular means for in silico experiments, as they are both simple and expressive [1]. Moreover, BNs have applications not only in molecular biology, but also in many other areas including circuit theory and computer science. BNs are composed of two essential parts. The first part is a finite set of Boolean variables representing genes or other biochemical substances. The second part is a set of The goal of the Boolean network control is to influence the behaviour of the network so that it stabilises in a particular attractor. A typical way to change the behaviour is to perturb the values of some variables. In this paper, we consider one-step state perturbations in which we apply all the perturbations of variables simultaneously and only once. After the perturbation, the system is left to behave normally. Solutions of the Boolean network control problem for a model of a cell in silico provide the basis for experimental designs allowing to reprogram the cell in vitro. Since practical realisation of particular perturbations requires non-trivial effort, the number of perturbations typically needs to be minimised. To that end, the control problem is usually enriched with some optimisation criteria.
In this paper, we focus on the source-target variant of the ParBN control problem where both the source and the target states are given in advance. To the best of our knowledge, we provide the first efficient solution to the one-step target control of ParBNs. It is worth noting that the parameters bring many new challenges to the control problem. First, the attractors might significantly change with the change in the parametrisation of the network. Second, the minimal state perturbation does not necessarily work for all parametrisations, and therefore, the notion of the optimal control strategy needs to be adapted to such a situation. Third, the parameter-space explosion (in addition to the state-space explosion) makes the problem computationally demanding. The parameter space is in the worst case, doubly-exponential [5]. Each ParBN parametrisation generates a unique BN model with a unique STG. That is why using algorithms developed for asynchronous non-parametrised BNs and computing the control set for each parametrisation distinctly (parameter scan) is not feasible, and a new approach needs to be developed.

Contribution
We propose an efficient alternative to the naïve parameter scan approach to compute the one-step target control of ParBNs. Technically, our approach relies on the integration of several formal methods. In particular, we employ a representation of ParBNs based on a symbolic edge-coloured graph using BDDs [4,6]. For this representation, we extend a well-established algorithm for asynchronous BNs control employing one-step perturbations [7]. The technique is based on the identification of attractor's strong basin. Our novel algorithm is able to compute the strong basins of all ParBN parametrisations simultaneously instead of computing them one by one (individually for every parametrisation). We show that for highly parametrised models, our approach is significantly more efficient than the naïve approach.

Related Work
Typically, two elementary types of control problems are distinguished: (i) achieving a single desired target attractor irrespective of the current state (target control) [8], (ii) achieving control between every pair of attractors (full control) [9]. Here, we focus on the target control.
The target control problem has been studied in non-parametrised BNs with both the synchronous and asynchronous update. In synchronous case, the method [8] identifies the control kernel, a minimal set of nodes, inferred from the explicit STG of the BN. In [10] a network graph aggregation approach is employed at the level of the regulatory network avoiding construction of the full STG. In [11], the regulatory network is divided to partially-dependent parts identified as strongly connected components. In the case of asynchronous BNs, the approach in [12] computes a set of relevant BN variables based on the identification of particular motifs in the regulatory network.
Traditional approaches mentioned above assume the control to be implemented by involving a permanent perturbation applied continuously for an extended period of time. However, this is not possible to be efficiently achieved in practise [13]. To that end, the concept of one-step perturbation has been introduced in [7] in the context of asynchronous BNs (the perturbation is applied in the initial state, and after that the system evolves according to its original dynamics). The algorithms in [7] address both target and full control, and they solve the existential problem as well as the optimal problem (identifying the minimal set of perturbations). The method is based on an efficient identification of a strong basin. In this paper, we employ a similar idea for the target control in the novel context of ParBNs.
It is also worth noting that in [14,15] the authors work with the concept of temporary perturbations filling the gap between one-step and permanent perturbations. Moreover, complex control strategies considering temporal sequences of perturbations are studied in [16][17][18][19]. All those results are developed for non-parametrised BNs only.
The approaches mentioned in this subsection cannot be directly lifted to work with ParBNs (the naïve solution considering iteration of existing methods is infeasible due to the combinatorial explosion of the parameter space). As already stated above, we propose the first efficient solution to the one-step target control of ParBNs.

Preliminaries
In this section, we introduce Boolean networks (BNs) and define related terms regarding long-term behaviour of BNs. Then we expand the notion of BNs by adding parameters, allowing for unspecified or unknown behaviour in the network.

Boolean Networks
Boolean networks are a simple model widely used to study complex dynamic behaviour of biological systems. That is why we define Boolean networks in a way that closely relates to regulatory networks, which represent biological processes using directed dependency graphs of biochemical entities: Definition 2.1 (Boolean network). A Boolean network is a tuple B = (V, R, F) such that: . .} is a finite set of Boolean state variables.
i.e. the subset of V regulating A.

Regulatory network
With every BN B it is possible to associate a directed graph (V, R) called a regulatory network or a dependency graph of B. This graph captures influences among the variables of B. When visualising a BN, its regulatory network is usually displayed as a directed graph (with update functions specified separately). For a BN B, one also often considers certain general properties of its regulations, which can then be depicted in the regulatory network.

Regulation types
We say that a regulation (A, B) ∈ R is observable if there exists a state such that changing the value of A also changes the value of F B , formally: Intuitively, this means that the presence of one biochemical entity has an observable influence on another entity. When a regulation is not marked observable, it can have an influence on the regulated entity, but we do not enforce it. Such regulations are drawn with a question mark.
In addition to observability, we also consider two possible monotonicity properties of a regulation: activation and inhibition. Regulation is activating if by increasing A it is not possible to decrease F B . For example, if A and C activate B, the possible functions F B are A ∨ C or A ∧ C. On the contrary, regulation is inhibiting, if by increasing A one can not increase the value of F B . For example, if A and C inhibit B, the possible functions F B are ¬A ∨ ¬C or ¬A ∧ ¬C. There can be regulations which are neither activating nor inhibiting. However, most regulations in this paper are either inhibiting or activating as this is typical for biological models. Graphically, activating regulation is depicted as a regular green arrow while inhibiting regulation is drawn as a flat red arrow (see Fig. 2a).

Dynamics
The complete dynamical behaviour of a boolean network B is captured by the directed state-transition graph (STG) G = (S, T ), where S = Π(B) and T ⊆ S × S. The definition of the transition relation T depends on the updating scheme that defines the way variables update their states along time. In this paper, we consider asynchronous updating. At a discrete time step, the system non-deterministically applies some F A ∈ F to a state s. We then obtain a transition relation → which is defined as follows: Note that for every (s, t) ∈→ the Hamming distance hdis(s, t) = 1 as during one step only one variable can change its value. For (s, t) ∈→, we simply write s → t. We use → * to denote a reflexive and transitive closure of →, also writing s → * t when (s, t) ∈→ * . Also note that the transition relation is non-deterministic. We denote Async(B) the state-transition graph of B under asynchronous updating scheme.
When studying the long-term behaviour of a BN, we typically only consider fair infinite paths in the state-transition graph. In a fair path, if a transition is enabled infinitely often, it has to be taken infinitely often. Therefore, the system cannot infinitely delay the available transitions, and it is not possible to cycle forever in a non-terminal strongly connected component.

Attractors
The long-term behaviour of BNs is captured by the notion of attractors. In biological models, we observe a phenotype in which the system eventually stabilises, whereas, in BN computational model, we observe attractors which are understood as terminal strongly connected components of the STG. In the following, we use these two terms interchangeably. such that for all s, t ∈ A, s → * t, and for all s ∈ A and t ∈ Π(B), s → t implies t ∈ A.
We denote a set of all attractors of B as A(B) or simply as A if B is clear from the context. For a fixed state s, we define Att(s) to be an attractor A ∈ A such that s ∈ A, or an empty set when s does not belong to any attractor. Furthermore, for an attractor A ∈ A(B) we also define notion of its weak and strong basins.
A strong basin is a set of states from which it is not possible to reach any other attractor than attractor A:

WB(B, A )
Notice that due to the fairness property, once a strong basin of an attractor A is reached, the system eventually stabilises in the given attractor. For better illustration, Figure 1 depicts weak and strong basins of some BN.

Parametrised Boolean Networks
Given a complex real-life system it might be very challenging to precisely determine all the update functions F of a Boolean network. Parametrised Boolean networks [3,4] provide a framework to deal with the lack of precise knowledge about the updating mechanism in a system. This extension assumes a set of logical parameters which determine the behaviour of update functions. Therefore, parametrised logical update functions either return a Boolean value (they behave normally) or a logical parameter representing the uncertainty of the consequent behaviour: • P ⊆ {0, 1} P is a subset of valid parametrisations; V} is a family of parametrised logical update functions. The signature of each F A is given as For p ∈ P , we write p(P) to denote the value of P in p and we also use the same notation p[P → k] for substitution as we used for states. The notion of the state space of a ParBN is identical to that of a BN. By fixing p ∈ P , we obtain B p = (V, R, F(p)) (a standard BN), A p (the set of attractors of B p ), and Att p (s) (the attractor of state s in the parametrisation p).  Figure 2: (a) A regulatory network of a simplified ParBN describing the DNA damage mechanism adapted from [20]. Every regulation is either activating (green) or inhibiting (red) and observable, except for (DNA, DNA), which is not necessarily observable. (b) All possible valid update functions F M2N satisfying the static constraints (monotonicity, observability). (c) Valid update functions F DNA , F M2C and F P53 satisfying the static constraints. Here, P i denote the parameters (P) of the ParBN.

Dynamics
The asynchronous semantics of a ParBN B is represented using an edge-labelled state-transition graph Async( B), where each transition s → t is labelled with a subset of parametrisations P (s, t) ⊆ P for which it is enabled. That is, p ∈ P (s, t) if and only if s → t in Async( B p ). For a fixed s, we denote successors(s) the set of all successors of s (states with Hamming distance one). In Fig. 2, we show a small example of a ParBN. Fig. 3 then presents its asynchronous semantics for selected subset of parametrisations.
In general, the size of the set of all possible parametrisations (parameter space) can be even doubly-exponential in the number of Boolean variables. In particular, the number of Boolean functions in a model with n variables is 2 2 n . It is thus critical to restrict the parameter space as much as possible. In many biological models, regulations are usually supplemented with static constraints limiting their outcomes [21,22].
We already presented observability, activation and inhibition as specific properties of regulations. In a parametrised setting, these properties can be used as constraints to restrict the parametrisation space. We assume that every regulation in a ParBN can be marked with a subset of these three constraints. Then for all p ∈ P of B, B p must adhere to these constraints, e.g. a regulation marked observable in B must be observable in B p and the same for activation and inhibition.
In Fig. 2a, a ParBN is displayed where all regulations are marked as either activating or inhibiting. Figures 2b and 2c then show the possible update functions satisfying these static constraints together with the corresponding logical parameters. Note that the fully parametrised model would have 16 parameters and 65536 parametrisations, but by applying the static constraints, only 27 parametrisations remain valid, significantly reducing the size of the associated edge-coloured graph.

Attractors
Notice that the standard notion of an attractor cannot be directly transferred to ParBNs, because a state can belong to an attractor only in certain parametrisations (for different parametrisations, the attractors do not have to overlap). We say that a subset A ⊂ Π( B) is an attractor in a parametrisation p ∈ P if A is an attractor of B p . Furthermore, given a state s, we can define Ap(s) as the subset of P , such that for each p ∈ Ap(s), Att p (s) = ∅.

Control Problem for Parametrised Boolean Networks
A computational model is controllable if we can assure that from some initial state, it reaches a desired final state in a finite amount of steps. This property was well-studied in the context of synchronous BNs and is also being pioneered for asynchronous BNs. However, the control problem for ParBNs has not been explored yet. ParBNs offer a more flexible representation of biological systems by allowing some uncertainty in the specification of the logic behind regulations. Since in reality information on regulatory mechanisms is typically ambiguous or unknown, studying the control problem for ParBNs enables new attractive, real-life applications, such as the discovery of candidate transcription factors for cell reprogramming (i.e., to change a cell's phenotype) when the model is not fully known.

Control Problem for Boolean Networks
There are multiple ways to control the behaviour of a BN, mainly differentiated between state and function perturbations. State perturbations force the system to change its current state. Alternatively, function perturbations adjust the update function(s) and therefore modify the edges of the system's state-transition graph. Typically, changing the processes in a cell (function perturbation) is more difficult then artificially adding or extracting a biochemical substance (state perturbation). For this reason, this paper focuses on state perturbations.
State perturbations can be further differentiated based on their temporal characteristics, mainly one-step, sequential, temporary, and permanent control. In one-step control, we change the values of the controlled variables once, and then the network evolves as originally defined. Sequential control identifies a sequence of perturbations that are applied at different time steps. When applying the control temporarily, there exists a finite amount of computational steps after which the control is released. Finally, the most intrusive control is permanent control of variables. However, this scenario is rather unrealistic in practice, as the given substance would need to be added or extracted forever. Moreover, the permanent perturbation might disrupt the original long-term behaviour of the network or introduce completely new behaviours.
Finally, we consider different control objectives, defined in [7] as follows: 1. Source-target control: Given a source state s ∈ Π(B), and a target attractor T ∈ A, find such control that when applied, the BN always converges from the state s to the attractor T . 2. Target control: Given a target attractor T ∈ A, find a control for every source attractor S ∈ A (such that S = T ) that when applied, the BN converges from S to T .
3. Full control: For all pairs of distinct attractors S, T ∈ A, find a control which guarantees that the BN converges from S to T . 4. All-pairs control: Given a subset of source attractors S ⊆ A and target attractors T ⊆ A, for every pair S ∈ S and T ∈ T , find a control which guarantees that the BN converges from S to T .
In some cases, one may also choose when the control is applied. Here, we consider immediate control, i.e. the control applied to the given state. Alternatively, during sequential control, the perturbations can be applied multiple times at different time points. This way we can sometimes control the system using fewer perturbations. However, this approach brings new issues, such as dealing with the non-determinism of the BN (different perturbations may be necessary in different branches of the non-deterministic network's behaviour). Moreover, we would need to be able to precisely observe current state of the BN currently. These issues can be addressed to some extent using attractor-based sequential approach [17].
In summary, the control problems for BNs differ in the following aspects: • What do we want to control; goal: What is the initial state of the BN? Where we want to end? Do we want to control only one scenario or multiple scenarios?
• What control we apply: We can either perturb states (once, temporarily, forever) or functions of variables; • When we apply control: Only once from an initial state in contrast with applying control to an arbitrary state and any amount of times.

One-step Control Set of Parametrised Boolean Networks
In this paper, we focus on state perturbations applied immediately in the initial state. This is the elementary way to perturb the system when solving the source-target control problem. The respective notion of one-step control is stated formally in the following definition.
Next we focus on establishing the control problem (when computing the one-step perturbation control) in ParBNs. Conceptually, the parametrised control problem is to some extent similar to the parameter synthesis problem. The goal of parameter synthesis is to find parametrisations that ensure some desired behaviour in the network. Such behaviour can also involve reachability of specific attractors. However, the key distinction between parameter synthesis problem and control problem is that in control, one determines perturbations that lead to a particular objective (with respect to given parametrisations), possibly resulting in behaviour that is not achievable only by tuning the network parameters. In this work, we thus treat parameters as unknown properties of the system rather than components that can be influenced.
One of the greatest challenges of ParBN control is that the attractors change based on the parametrisation. To that end, it makes sense to solve the control problem only for parametrisations which admit the given attractor. The parametrised control problem then computes a mapping that associates a potential control with the maximal set of parametrisations for which the control is applicable (the so-called controlled parametrisations). Formally, the considered parametrised control problem is stated in the following definition. Definition 3.2 (Source-target control in ParBN). Given a ParBN B, a source state s and a target state t, find a mapping Cp : C → 2 Ap(t) which assigns each control C ∈ C a maximal (possibly empty) set of parametrisations p for which when C is applied, B p converges from s to Att p (t). We refer to Cp(C) as the control enabling parametrisations for C.
Intuitively, control enabling parametrisations Cp(C) = p are the parametrisations, for which the target state is achieved by applying C in the non-parametrised case. Notice, that source-target problem in context of ParBN is aiming to drive the network into a target state of some attractor instead of an attractor itself. This is because parameterisations might contain attractors which are considered similar as all of them contain the given state. Therefore, in all parameterisations Apt the given state t is entered infinitely often by the controlled BN. If this is not the case of some ParBN and some parametersiations contain attractors which are out of the interest, the set of parameterisations AP (t) might be replaced with an arbitrary one.
It is worth noting that it is always possible to bring the system into the given target state by setting the ParBN's variables to the values adequately (with the values of the given target state). We call this control trivial. However, when controlling a system, we typically look for a control which requires the fewest interventions as possible. Therefore, we want to 8 minimise the number of the controlled variables and the trivial solution might not be optimal. In a non-parametrised setting, this is typically the only considered optimisation criterion.
In a ParBN, the situation is further complicated due to the dependence on the parameters. To reach some attractor, it is sufficient to reach its strong basin after the application of a control. Nonetheless, the strong basin of an attractor can vary according to the parametrisation, and the control thus "works" only for its control enabling parametrisations. We are interested in maximal sets of control enabling parametrisations. To that end, we consider the notion of robustness that normalizes the number of control enabling parametrisations.
Unfortunately, it is not always possible to ensure that some control is both minimal and the most robust. For example, there may be a control set C small in size, which only works for a small fraction of the parameter space. Consequently, while C is easily applicable, it may be unlikely to work in reality, as the real behaviour of the system can also follow one of the parametrisations which are not controlled by C. We discuss this issue in more detail in Section 5.

The Algorithmics
We are now ready to describe our computational framework for solving the one-step state perturbation control of ParBN. We start by introducing our approach for finding strong basins in a ParBN. Then we explain the framework for exploring ParBN's STG and for manipulation with ParBN's parametrisations. Finally, we build a concise workflow for ParBN control employing the proposed algorithms.

Semi-symbolic Parametrised Strong Basin Search Algorithm
We assume the set of parametrisations of a ParBN is represented as a reduced ordered binary decision diagram (BDD) [6]. The decision variables of the BDD are the parameters of the network (P), meaning that every path from the root to a leaf in such a BDD represents a parametrisation of the ParBN. Common logical operations on such BDDs (and, or, negation, ...) then correspond to set operations (intersection, union, complement, ...). Furthermore, static constraints (activation, inhibition, observability) can be formalised using Boolean formulae over P, and we can, therefore, create a BDD which enforces all imposed constraints and represents the set of all valid parametrisations.
Recall that we represent Async( B) as an edge-labelled state-transition graph, where each transition s → t has an associated set of parametrisations P (s, t) ⊆ P represented as a BDD. A parametrised state set is a mapping Π( B) → 2 P assigning to each state a set of parametrisations. Furthermore, we suppose that the state space is represented explicitly, meaning that all operations on states are performed element-wise (typically in parallel).
We consider parametrised reachability procedures -given a source state s and a parameter set P , these procedures compute a maximal parametrised state set of all forward/backward reachable states from the source state s (containing all reachable states t where each t is associated with a maximal set of parametrisations P t for which t is reachable from s): The chosen representation (explicit state space and symbolic parametrisations) allows to compute the reachability procedures in parallel. For the underlying implementation, we are working with internal libraries of the tool Aeon [23] which provides the most of the necessary functionality, including a convenient format for specifying ParBNs and parallel reachability procedures.
The key observation for controlling ParBNs is that an attractor (by definition) is always reached from its strong basin. From all other states, it is possible to reach also some other attractor(s); therefore, reaching the target attractor is not guaranteed. When the attractor's strong basin for all parametrisations is known, we can compute the control mapping Cp from a source state to the target attractor by considering Hamming differences between the source and the states of the strong basin.
Our algorithm is based on the fixed-point approach for strong basin computation in non-parametrised BNs [24]. The premise is that only states from which it is possible to reach the attractor (weak basin) can be part of the strong basin. The weak basin of the attractor is computed using some standard reachability algorithms, for example, BFS. Then states are iteratively removed if it is possible to leave the basin from them (and therefore not reach the given attractor). Finally, if there are no states left to remove, the fixed-point is achieved, and the strong basin is found.
For parametrised control, we first need to determine Ap(t) for the target state t, since in all other parametrisations, t is not a part of an attractor. This process is described in Algorithm 1. In Algorithm 2, we then extend the approach from [24] onto ParBN. For each state we remember under which parametrisations it is in the basin, starting with the parametrised weak basin. Then we iterate over the basin's states. For each state, we consider its successors, and we compute the parametrisations for which some successor is not in the strong basin. In these parametrisations, it is possible to leave the basin, and therefore these parametrisations are not part of the strong basin, and so we remove them for this state. We can process the states in parallel and compute for them the parameterisations which are not part of the strong basin. To avoid writing of the strong basin structure while it is intensively read from, we store the parameterisation which should be removed separately and after it is computed, we update the strong basin structure sequentially. The procedure is iterated until nothing more can be removed from the strong basin.

Control Computation Workflow
Now we can describe a complete workflow for computing the source-target control in ParBNs, depicted in Fig. 4. The workflow consists of three inputs and three computation steps, resulting in the control mapping Cp. Given an input ParBN and a target attractor state target, we start by computing valid parametrisations set Ap(target) using Algorithm 1. Then, the parametrised strong basin is computed from the ParBN with target attractor state target and its valid parametrisations Ap(target) using Algorithm 2. After that, from the strong basin and the source state, we obtain the complete control mapping. Notice that we do not need to know the source state for computing the strong basin and that we can re-use the target's strong basin for obtaining control for different sources.
To compute the control mapping, observe that viable controls correspond exactly to the Hamming differences (the variables with opposite values) between the source state and the states of the strong basin; yielding one viable control C for every state s in the strong basin SB. Any other control C does not reach the strong basin and therefore does not guarantee to reach the target attractor (for these, Cp(C ) = ∅). A control C is then viable only for parametrisations for which s appears in the strong basin, we thus set Cp(C) = SB(s).
Finally, we can compute the size and robustness of each control as we have all knowledge regarding the variables which need to be controlled and for which parametrisations the controls work. If it is desired, we can construct a witness BN for a control C (a non-parametrised BN where the given control works) by fixing some parameter from Cp(C). The prototype implementation containing all parts of the workflow is available 1 .
In highly parametrised models, it is not unusual to obtain a control with size 0, where in some parametrisations no action would be needed to control the model. This is because, in some parametrisations, the source might already be a part of the target's strong basin. If this case is considered unrealistic (since there is probably a need for control, thus these parametrisations appear to be invalid), the parametrisations Ap(target) can be replaced with a custom set of parametrisations.
The resulting set of all available controls can be then used for e.g. cell reprogramming. However, we might obtain many possible controls, and we need to decide which one should be applied. To do that, we can decide based on the size of control or its robustness. The control set can be arbitrarily filtered discarding controls with size bigger than the trivial control (which always has 100% robustness) or bigger than some set size. Similarly, the control set can be pruned based on too low robustness. The interplay of these two factors is non-trivial, and it is left to the actual application to decide which control would best suit its needs.

Evaluation
We evaluate our approach on two real-life BN models. We compare the performance of our approach using a different number of parameters implanted into the models, resulting in different size of relevant parameter space. We conducted all measurements using a machine equipped with AMD Ryzen Threadripper 2990WX 32-Core Processor and 64GB of memory.
The first model is a cell-fate decision model [25]. The model provides a high-level view of possible different cell fates such as pro-survival, necrosis or apoptosis. We used a fully-parametrised version (all update functions are completely parametrised) of this model and we selected seven biologically relevant attractors. We computed strong basins only for these attractors in several parametrised versions of the model differing in the number of unknown parameters (see Table 1).
The second model, a myeloid differentiation network [26], was designed to model a muscle tissue cell differentiation from common myeloid cell to specialised muscle cells (megakaryocytes, erythrocytes, granulocytes and monocytes). The original non-parametrised network has eleven nodes and six attractors. We derived several parametrised versions of the model by arbitrarily parametrising update functions of the model. Similarly to the previous case, we used only attractors of the original network when computing strong basins for parametrised versions of the model. The results are again shown in Table 1. Table 1: Results of strong basins computation. The values are stated as ranges because we computed strong basins of all attractors considered in given models. The second column shows the number of model's parameters. The third column shows count ranges of parametrisations which contain the given attractor. The fourth (resp. fifth) column displays ranges of the number of states in the weak (resp. strong) basins. The last column contains ranges of times needed to compute strong basins. Next, we "virtually" compare our approach to the naïve parameter-scan approach. In [7], a strong basin of (nonparametrised) asynchronous BNs is computed using a block decomposition method with 4 ms needed to finish the computation for the myeloid model. Even if the reported HW was slower than in our case and we assume that the strong basin computation for one parametrisation would last only 1 ms, the fully parametrised myeloid model contains an attractor which is present in 3.7 × 10 9 parametrisations. Therefore, the expected time for computing a strong basin for all parametrisations with 32-fold parameterisation would last more than a day compared to less than 27 seconds achieved using our parameter-based semi-symbolic approach.

Model
We evaluate the scalability of our approach on a fully parametrised myeloid model. The results are shown in Table 2.
The computation was restricted to the specified amount of CPUs. The final speed-up achieved on our machine, when using 32 CPUs compared to a non-parallel CPU usage was about 10-fold. Now let us have a look at an example of how the results might be interpreted and an suitable control selected. Suppose that we want to reprogram an erythrocyte cell (phenotype having factors EKLF=1 and GATA-2=0) into a monocyte cell (phenotype having factors cJun=1 and EgrNab=1) of the myeloid model. First, we obtain a strong basin of the monocyte attractor having 1472 states yielding us 1472 possible control sets. We can observe, that trivial; i.e., the most robust control strategy (setting variables so that we reach the attractor right after applying the control) has size 8. We can discard all controls with bigger size 8 and bigger as they are less optimal in all aspects than the trivial control. In our case, there are 1259 controls with smaller size than the trivial one.
Resulting smallest control sets have the size 1. However, the best robustness among these controls is 46%. Therefore, it is quite likely that it will not work in the practice; for a real cell. If we allow the control to have the size 2, we can 12 use the control with 76.8% robustness. Further increasing the size provides control with the size 3 and the robustness 92% and so is highly likely to reprogramm the cell's phenotype. To achieve only slightly more robustness (93.8%) we need to use a control of the size 5. In this highly parametrised model, there is no control with 100% robustness being smaller than the trivial one. It is not a rule of the thumb that with bigger size of control we obtain better robustness, for example, the control with the size 11 (the only one, with all variables, reversed compared to the original one) has the robustness of only 50%.
It can be seen that the unknown properties of ParBN make selection of one particular "best" control complicated. That is why a careful in vitro experimentation is needed to verify the correctness of the control in the reality. Nonetheless, the obtained control set still can help and highly reduce the exponential number of potential transcription factors combinations.

Conclusion
We have introduced the control problem for parametrised Boolean networks and we proposed a algorithm for solving the source-target variant of this problem using one-step perturbations. The core procedure of the algorithm is a fixed-point computation of the parametrised strong basin of the given target. The method we proposed is semi-symbolic -it relies on a unique integration of symbolic (based on BDD representation) and explicit (based on enumerative model checking) formal methods. We demonstrated that our approach is capable to control highly parametrised models in seconds.
Owing to the doubly-exponential explosion of the number of possible parametrisations, such a result cannot be achieved with the naïve parameter scan approach.
In the future work, we would like to scale up the algorithm for searching the strong basin, i.e., by using some symbolic approach for exploring ParBN's state transition graph. Moreover, we would like to consider other variants of ParBN control problems based on variants of non-parametrised control problems (see Subsection 3.1). Last but not least, we would like to incorporate the developed methods into an existing ParBN toolkit -the AEON tool [23].