# Refining Boolean models with the partial most permissive scheme

Nadine Ben Boina, Brigitte Mossé, Anaïs Baudot, Élisabeth Remy

<figure><img src="image_1.png" /></figure>

## To cite this version:

Nadine Ben Boina, Brigitte Mossé, Anaïs Baudot, Élisabeth Remy. Refining Boolean models with the partial most permissive scheme. 2025. hal-04799493

HAL Id: hal-04799493

https://hal.science/hal-04799493v1

Preprint submitted on 26 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Refining Boolean models with the partial most
permissive scheme

Nadine Ben Boina¹,², Brigitte Mossé¹, Anaïs Baudot²,³,⁴, and
Elisabeth Remy¹

¹Aix Marseille Univ, CNRS, I2M (UMR 7373), Turing Center for Living systems,
Marseille, France

²Aix Marseille Univ, INSERM, MMG, Marseille, France

³CNRS, Marseille, France

⁴Barcelona Supercomputing Centre, Barcelona, Spain

**Abstract**

**Motivation** In systems biology, modelling strategies aim to decode how molecular components interact to generate dynamical behaviour. Boolean modelling is more and more used, but the description of the dynamics generated by discrete variables with only two values may be too limited to capture certain dynamical properties. Multivalued logical models can overcome this limitation by allowing more than two levels for each component. However, multivaluing a Boolean model is challenging.

**Results** We present MRBM, a method for efficiently identifying the components of a Boolean model to be multivalued in order to capture specific fixed-point reachabilities in the asynchronous dynamics. To this goal, we defined a new updating scheme locating reachability properties in the most permissive dynamics. MRBM is supported by mathematical demonstrations and illustrated on a toy model and on two models of stem cell differentiation.

**Availability and Implementation** Data and code available at:
https://github.com/NdnBnBn/MRBM.git.

**Contact** elisabeth.remy@univ-amu.fr
# 1 Introduction

Systems biology aims to understand how the interactions between biological components such as molecules, cells, tissues, or organs, give rise to specific phenotypes and behaviours Bruggeman and Westerhoff (2007). Mathematical models, integrating data from experiments and literature to represent biological systems and facilitate their study, play a pivotal role in systems biology. They have been instrumental in studying a broad spectrum of biological systems, ranging from signal transduction Heinrich et al. (2002) to complex developmental processes Tomlin and Axelrod (2007). Importantly, mathematical models reduce the cost and time associated with *in vivo* or *in vitro* experiments.

Various mathematical frameworks are used to model biological interaction networks Chen et al. (2010); Ay and Arnosti (2011); Torres and Santos (2015). We focus here on logical modeling, a formalism that qualitatively captures the main dynamical properties that explain the overall behaviour of the system, without requiring precise parameters. In this formalism, the activity levels of the components of a biological system are represented as discrete variables. Logical models encompass both Boolean models (BMs) and multivalued models. In BMs, the activity level of the components are represented as binary variables (inactive/active) Thomas (1973). However, BMs may not be sufficient to capture the complex dynamics of biological systems. Multivalued models can address this limitation by allowing some components of the model to have more than two levels of activity, while retaining the conceptual simplicity of logical modeling Thomas (1978). When a multivalued model is derived from an existing BM, without adding new components nor modifying the effects of the regulations (activation or inhibition), we refer to this multivalued model as a multivalued refinement of the BM Paulevé et al. (2020).

The dynamics of logical models is driven by logical functions and an updating scheme. The logical functions specify the conditions of activation of each component according to the current activity levels of the components of the system, while the updating scheme dictates how the system is updated. The choice of updating scheme is crucial as it impacts the dynamics of the system. Two updating schemes are commonly used, the synchronous and the asynchronous updating schemes<fnref n="1" /> Garg et al. (2008). Under the synchronous scheme, the activity levels of all the components are updated simultaneously, resulting in a deterministic dynamics. Under the asynchronous scheme, a single component's activity level is updated at a time. The resulting dynamics is not deterministic, as some states may have several successors. The asynchronous scheme is often preferred in modeling biological systems. Indeed, the asynchronous scheme enables the emergence of complex behaviours that are not revealed with the synchronous scheme Paulevé and Sené (2022). The most permissible updating scheme (m.p. scheme), recently introduced in Paulevé et al. (2020), is a non-deterministic scheme that only applies to BMs. The m.p. scheme adds to

<footnote n="1">¹sometimes denoted fully asynchronous</footnote>
each component two intermediate activity levels which mark the transition between
the two Boolean levels. These intermediate levels represent the fact that the com-
ponent is in the process of increasing its level from 0 to 1 (level 'i') or decreasing
it from 1 to 0 (level 'd'). When a regulator of a component is at an intermediate
level, its logical function takes into account both the situation where the regulator
is present and the situation where it is absent.

Attractors, which represent the long-term behaviours of a model, are key features
of the dynamics. Once identified, they can be associated with distinct phenotypic
readouts Schwab et al. (2020). The existence of a trajectory from a given state to
a specific attractor proves the reachability of the attractor from this state. The
number of attractors present in the dynamics, and their reachability, are impacted
by the choice of the updating scheme. In particular, the dynamics of a BM under
the m.p. scheme captures all the trajectories obtained under the synchronous and
asynchronous schemes, as well as all the trajectories obtained from any multivalued
refinement of the BM. As a result, more reachability properties are satisfied in the
m.p. dynamics than in the other dynamics (see Paulevé et al. (2020)).

We introduce here MRBM (Multivalued Refinement of Boolean Model), a method
which aims at identifying components to be multivalued in a refinement of a BM
in order to provide the desired reachabilities within the asynchronous dynamics.
MRBM utilizes the partial m.p. scheme, a new updating scheme that we are propos-
ing here (Section 2.4) and that we adapted from the m.p. scheme. In the partial
m.p. , only a subset of the model’s components is updated using the m.p. scheme, the
remaining components being updated with the asynchronous scheme. The resulting
dynamics helps to pinpoint the components of the BM that need to be multival-
ued. In this manuscript, we give a detailed description of the method MRBM. We
illustrated its potential with a toy model (Section 3.1). We then use the MRBM
method to propose multivalued refinements of two BMs of stem cell differentiation.
The first model (Section 3.2) describes the aging effects on hematopoietic stem
cell differentiation. The second model (Section 3.3) explores the differentiation of
stem cells in the root of *Arabidopsis thaliana*. In both cases, we find a multivalued
refinements which display reachabilities that were not present in the BMs.

2 Method

2.1 Boolean Modeling

Set $\mathbb{B} = \{0, 1\}$ and $n$ an integer $> 0$. A Boolean Model (BM) of dimension $n$ is a map
$f = (f_1, \dots, f_n) : \mathbb{B}^n \to \mathbb{B}^n$, where each Boolean state $x = (x_1, \dots, x_n) \in \mathbb{B}^n$
specifies the activity level of $n$ components of the model $g_1, \dots, g_n$. The values of the
logical regulatory function $f_i(x) \in \mathbb{B}$ provide the target levels of the $i$-th component.
This model is associated with a *regulatory graph*, i.e. a directed signed graph where
the $n$ nodes $g_1, \dots, g_n$ of the graph are the components of the model, and the edges
between nodes represent the regulations between components. Each edge has a sign,
positive for activation, negative for inhibition (Figure 1A). Note that dual edges are allowed. Dual edges represent interactions that can be activations or inhibitions depending on the context.

<figure><img src="image_3.png" /><figcaption>Figure 1: **Toy Boolean Model (BM).** (A) Regulatory graph associated with a BM f composed of three components $g_1, g_2, g_3$, five interactions (red edges for inhibition, green edges for activation), and logical regulatory functions (operators "&", "|", "!" stand for AND, OR, NOT respectively). (B) The full State Transition Graph (STG) of the toy model in the asynchronous scheme. Each node represents a state $x = (x_1, x_2, x_3)$ of the model and each edge represents an asynchronous transition between two consecutive states. Fixed points are colored in blue. (C) STG of the toy model in the m.p. dynamics, representing all the m.p. trajectories from the initial state 111. (D) STG in the m.p. dynamics of the toy model, representing all the m.p. trajectories from the initial state 111, restricted to the Boolean states. (E) The full STG of the toy model in the m.p. dynamics, displaying only the Boolean states. (F) The $g^1$ m.p. STG (full STG in the partial m.p. dynamics applied on gene $g_1$) of the toy model, displaying only the Boolean states. The trajectories from the state 010 to the state 011 are highlighted in red (see Section 3.1.2).</figcaption></figure>
The dynamics of the model is represented by the State Transition Graph (STG). The STG is a directed graph ($\mathbb{B}^n, T$) in which $\mathbb{B}^n$ is the set of nodes that represent the states of the model, and $T \subset \mathbb{B}^n \times \mathbb{B}^n$ is the set of edges that represent the transitions between states. The full STG contains all states of the system. We can also represent a sub-graph of the STG, which contains all the states that can be reached from a given initial state. Given a BM, an updating scheme determines the stepwise evolution of the system. Within the asynchronous scheme, the activity level of only one component can be updated at a time according to the function $f$ (Figure 1B). More precisely, there is a transition from the state $x = (x_1, \dots, x_n)$ to the state $y = (y_1, \dots, y_n)$ if there exists $i_0 \in \{1, \dots, n\}$ such that $x_{i_0} \neq f_{i_0}(x)$ and

$$y_{i_0} = x_{i_0} + \text{sign}(f_{i_0}(x) - x_{i_0}), \\ y_j = x_j, \text{ for } j \neq i_0.$$

Note that the asynchronous scheme generates a non-deterministic dynamics, as each state can have several successors.

A transition from state $x$ to state $y$ using the asynchronous scheme is denoted by $x \xrightarrow{asyn} y$. A trajectory from state $x$ to state $y$ is a sequence of transitions $x \xrightarrow{asyn} x' \xrightarrow{asyn} \dots \xrightarrow{asyn} y$. A state $y$ is reachable from $x$ if there exists a trajectory from $x$ to $y$. We denote this reachability by $x \stackrel{f}{\rightsquigarrow}_{asyn} y$.

The attractors of the STG are the terminal strongly connected components (SCCs) of the graph, gathering the asymptotic states of the model. We distinguish two types of attractors, fixed points and cyclic attractors. Fixed points are terminal SCCs of size 1, whereas cyclic attractors are terminal SCCs of size $\ge 2$. Given an attractor $\omega$, the basin of attraction $\mathcal{B}_{asyn}^f(\omega)$ is the set of states from which a state of $\omega$ is reachable.

The above concepts on trajectories and attractors are directly adaptable to multivalued models and other updating schemes.

## 2.2 Multivalued Refinements of Boolean Models

A natural extension of BMs is to have multivalued variables attached to some components, enabling these components to have more than two levels of activity.

Given $m_j \in \mathbb{N}^*$ for $j \in \{1, \dots, n\}$, with at least one $m_j > 1$, let $X = \prod_{j=1}^n \{0, 1, \dots, m_j\}$. A multivalued model on $X$ is a map $h = (h_1, \dots, h_n) : X \to X$ satisfying the condition $h_j(x) - x_j \in \{-1, 0, +1\}$, for each $x = (x_1, \dots, x_n) \in X$ and $j \in \{1, \dots, n\}$.

Let us consider $f = (f_1, \dots, f_n)$ a BM of dimension $n$. Following Paulevé et al. (2020), and given $X = \prod_{j=1}^n \{0, 1, \dots, m_j\}$ as above, we introduce for $x \in X$ the set

$$\alpha(x) = \{x' \in \mathbb{B}^n ; \forall j \in \{1, \dots, n\}, (x_j = 0 \Rightarrow x'_j = 0) \text{ and } (x_j = m_j \Rightarrow x'_j = 1)\}. \\ \text{For instance, if } X = \{0, 1, 2\}^3, \alpha(021) = \{010; 011\}. \text{ A multivalued model } h \text{ on } X$$
is said to be a *refinement of* f if, for any $x = (x_1, \dots, x_n) \in X$ and $j \in \{1, \dots, n\}$,

$$h_j(x) < x_j \implies \exists x' \in \alpha(x), f_j(x') < x'_j,$$

$$h_j(x) > x_j \implies \exists x' \in \alpha(x), f_j(x') > x'_j.$$

This means that, in a multivalued refinement, the level of activity of a component may decrease (resp. increase) at a given state only if its value goes to 0 (resp. 1) under f at some Boolean state of $\alpha(x)$.

Starting from f, we use the following approach to parameterize a multivalued refinement, i.e. to define the function h. For each multivalued component $g_j$, we consider the logical formulas of its targets $g_k$ (this means that $g_k$ is regulated by $g_j$). For each conjunctive clause (i.e. conjunctive logical assertion) of these formulas, we choose a threshold s in $\{0, \dots, m_j - 1\}$. This threshold specifies that the regulation of $g_k$ occurs when $x_j \in \{s+1, \dots, m_j\}$. The threshold s is chosen according to the target and to the clause (see details in Supplementary Data Section A).

In the remaining sections of this manuscript, all the refinements we consider are built as described above. In the Results section, we describe these refinements h through the functions $\mathcal{H}_j$ of $x \in X$ with values in $\{-1, +1\}$. These functions $\mathcal{H}_j$ indicate whether gene j is called to decrease or to increase its level $x_j$ under h (when possible, i.e. remaining between 0 and $m_j$).

Note that in most cases, given a BM f, if there exists an asynchronous path from $x \in \mathbb{B}^n$ to $y \in \mathbb{B}^n$, a corresponding path also exists in the asynchronous dynamics of any refinement of f. Moreover, for any fixed point in the asynchronous dynamics of f, a corresponding fixed point exists in the asynchronous dynamics of any refinement of f belonging to $\prod_{j=1}^n \{0, m_j\}$, and conversely. There is one exception. When any multivalued component is self-inhibited, these results are wrong: reachability properties are modified, and even the number and nature of attractors may change (Supplementary Data Section B).

Finally, it should be noted that there is no need to multivalue an output (component without target) in a refinement. Indeed, the multivaluation of an output would have no consequences in terms of reachability properties.

## 2.3 Most Permissive Scheme

The *most permissive scheme* (m.p. scheme) is a non-deterministic updating scheme for BM proposed by Paulevé et al. (2020). With this updating scheme, each component can take four activity levels: 0, 1, i, and d. We set $X_{m.p.} = \{0, 1, i, d\}^n$. The increasing and decreasing levels i and d refer to intermediate levels. A component must go through these intermediate levels to reach its target level, either 1 or 0. A component at an intermediate level can be considered by its targets as being at either level 0 or level 1. We associate a set of Boolean states to each state x of $X_{m.p.}$:

$$\gamma(x) = \{x' \in \mathbb{B}^n; \forall j \in \{1, \dots, n\}, (x_j = 0 \Rightarrow x'_j = 0) \text{ and } (x_j = 1 \Rightarrow x'_j = 1)\}.$$
For instance, $\gamma(di1) = \{001; 011; 101; 111\}$. Then, given a BM $f$, there is a transition from a state $x \in X_{m.p.}$ to a state $y \in X_{m.p.}$ if there exists $j_0 \in \{1, ..., n\}$ such that one of the following occurs

$$ 
\begin{aligned}
& -x_{j_0} \in \{0, d\} \text{ and } \exists x' \in \gamma(x) \text{ such that } f_{j_0}(x') = 1 \text{ and } y_{j_0} = i, \\
& -x_{j_0} \in \{1, i\} \text{ and } \exists x' \in \gamma(x) \text{ such that } f_{j_0}(x') = 0 \text{ and } y_{j_0} = d, \\
& -x_{j_0} = i \text{ and } y_{j_0} = 1, \\
& -x_{j_0} = d \text{ and } y_{j_0} = 0,
\end{aligned}
$$

and $y_j = x_j$, for $j \neq j_0$.

All states with at least one component activity level equal to $i$ or $d$ will be referred to as m.p. states, and other states as Boolean states. For the sake of clarity, when representing the m.p. STG we do not show the m.p. states. For instance, the trajectory: $111 \to 1d1 \to dd1 \to 0d1 \to 001$ (Figure 1C) is represented by the transition $111 \to 001$ (Figure 1D). The full m.p. STG of the toy model is represented (Figure 1E).

An interesting property emerges from the m.p. scheme: given a BM, the m.p. scheme captures all the trajectories present in any of its multivalued refinements. This completeness property of the m.p. is described in Paulevé et al. (2020). Consequently, if a reachability is not achievable using the m.p. scheme, then no multivalued refinement will present this reachability.

## 2.4 Partial Most Permissive Scheme

We introduce here the partial m.p. scheme, derived from the m.p. scheme. Each of these new updating schemes is associated with the choice of a subset $J$ of $\{1, ..., n\}$. We then assign the four possible levels $0, 1, i$, and $d$ to components $g_j$ such that $j \in J$. We refer to these chosen components as "m.p. components", and we keep the Boolean levels for all other components. So we consider the set of states $X_{J_{m.p.}} = \prod_{j=1}^{n} A_j$, where $A_j = \{0, 1, i, d\}$ if $j \in J$, and $A_j = \{0, 1\}$ otherwise.

Given a BM $f$ and $J$ as above, we then define the partial $J$-m.p. scheme on $X_{J_{m.p.}}$ in the following way: there is a transition from a state $x \in X_{J_{m.p.}}$ to a state $y \in X_{J_{m.p.}}$ if there exists $j_0 \in \{1, ..., n\}$ such that one of the following occurs

$$ 
\begin{aligned}
& -j_0 \in J \text{ and } x_{j_0} \in \{0, d\} \text{ and } \exists x' \in \gamma(x) \text{ such that } f_{j_0}(x') = 1 \text{ and } y_{j_0} = i, \\
& -j_0 \in J \text{ and } x_{j_0} \in \{1, i\} \text{ and } \exists x' \in \gamma(x) \text{ such that } f_{j_0}(x') = 0 \text{ and } y_{j_0} = d, \\
& -j_0 \in J \text{ and } x_{j_0} = i \text{ and } y_{j_0} = 1, \\
& -j_0 \in J \text{ and } x_{j_0} = d \text{ and } y_{j_0} = 0, \\
& -j_0 \notin J \text{ and } x_{j_0} = 0 \text{ and } \exists x' \in \gamma(x) \text{ such that } f_{j_0}(x') = 1 \text{ and } y_{j_0} = 1, \\
& -j_0 \notin J \text{ and } x_{j_0} = 1 \text{ and } \exists x' \in \gamma(x) \text{ such that } f_{j_0}(x') = 0 \text{ and } y_{j_0} = 0,
\end{aligned}
$$

and $y_j = x_j$, for $j \neq j_0$.

Note that if the size of the set $J$ is equal to $n$ (the dimension of the model), then the partial $J$-m.p. scheme is the m.p. scheme. If $J$ is empty, the partial $J$-m.p. scheme
is the asynchronous scheme. The $^g1$ m.p. STG of the toy model is represented (Figure 1F).

The completeness property of the m.p. scheme can be generalized to the partial m.p. scheme: given a BM $f$ and a set $J$, the $^J$m.p. dynamics encompasses the dynamics of any multivalued refinement of $f$ with $\{g_j; j \in J\}$ as the set of multivalued components (see proof in Supplementary Data Section C).

We can emphasize the following properties:

* If $J_1$ and $J_2$ are subsets of $\{1, \dots, n\}$ such that $J_1$ is contained in $J_2$, then $X_{J_1 m.p.}$ is contained in $X_{J_2 m.p.}$; moreover, if there exists a trajectory between two states in the $^J1$ m.p. dynamics of $f$, the same holds in the $^J2$ m.p. dynamics. In particular, if there exists a trajectory between two Boolean states in the asynchronous dynamics of $f$, the same property holds in any partial m.p. dynamics.

* The asynchronous dynamics and the partial m.p. dynamics of $f$ share the same fixed points.

For the sake of simplicity, in the Results section, we will identify the set $J$ as the set of m.p. components.

## 2.5 Assessment of fixed points reachability property

Assessing the reachability of attractors involves determining whether specific states can reach an attractor, or assessing the entirety of states that can reach an attractor (i.e. characterizing the basin of attraction of an attractor). In this section, we will describe the methods we used to identify the existence of reachabilities in the dynamics of a logical model and to compute the basins of attraction.

The MRBM approach is based on comparing the reachability properties of fixed points across the dynamics obtained with different updating schemes for the same BM. The fixed points of a BM do not change with the updating scheme. However, the basins of attraction of these fixed points differ according to the updating scheme, reflecting changes in reachability. So, to identify these changes in reachability, we compare the basins of attraction. The state spaces are not the same depending on the choice of update: the m.p. states are not present in the Boolean asynchronous dynamics. In order to make a relevant comparison between asynchronous, partial m.p. and m.p. dynamics, we have only considered Boolean states.

In the same way, there is a one-to-one correspondence between the fixed points of a BM and those of any of its multivalued refinements, provided there is no self-inhibited multivalued component (see Section 2.2). To compare the basins of attractions in a BM with the basins in a multivalued refinement, we only took into account the states of the subspace $\prod_{j=1}^n \{0, m_j\}$. This subspace is in correspondence with the set $\mathbb{B}^n = \prod_{j=1}^n \{0, 1\}$.
### 2.5.1 Reachability Identification

The tools available to identify reachabilities do not support all the possible updating schemes. In this context, for BMs updated with the asynchronous or partial m.p. schemes, we verified reachability with model checking Baier and Katoen (2008). To do so, we used the NuSMV model checker software extension Cimatti et al. (2002) provided in the PyBoolNet Python package Klarner et al. (2020) and we expressed the attractor reachability property to be verified using computational tree logic (CTL) formulas, as described in Baier and Katoen (2008); Klarner et al. (2020). For BMs updated with the m.p. scheme, we used the Python package mpbn Paulevé et al. (2020). This package is based on Answer-Set Reprogramming Gebser et al. (2013) and on the solver Clingo Gebser et al. (2014).

### 2.5.2 Basins of Attraction Identification

Given a logical model (Boolean or any multivalued refinement) and its asynchronous dynamics, we use model checking to compute the basin of attraction of a fixed point $\omega$. This is done by identifying, for any state $x$, whether there exists a path from $x$ to $\omega$.

In the partial m.p. scheme, basins of attraction are computed in the same way. There is however a preliminary step to generate the partial m.p. dynamics of a BM using the bioLQM toolkit Naldi (2018). In the m.p. scheme, the basin of attraction are computed using the Python package mpbn Paulevé et al. (2020).

We use the size of the basins of attraction of fixed points as a measure of attractor reachability. For a logical model $f$, we define the size of the basin of attraction of $\omega$, denoted by $\#B^f(\omega)$, as the number of Boolean states (or equivalent states for a multivalued refinement) present in the basin. We report the sizes of the basins as percentages relative to the sizes of the Boolean states space ($\frac{\#B^f(\omega)}{2^n} \times 100\%$).

## 2.6 Availability and Implementation

The MRBM method and the BMs used in this manuscript are available on GitHub at: https://github.com/NdnBnBn/MRBM.git. This repository contains the scripts necessary to run the MRBM method within the MRBM directory. Three supplementary directories contain the toy model and the two stem cell differentiation models presented in the Results section. Each of these directories contains the original BM files, the output generated after running MRBM, and the multivalued refinements described in the results section. The method relies on the Colomoto Docker environment Naldi et al. (2018).
# 3 Results

## 3.1 The MRBM Method: defining Multivalued Refinements of Boolean Models

### 3.1.1 Description of the MRBM method

Let $f$ be a BM such that its asynchronous dynamics displays several fixed points and no cyclic attractor.

Suppose that a reachability property $\mathcal{P}$ is satisfied in the m.p. dynamics of the BM $f$, but not in the asynchronous dynamics. The MRBM method aims to identify multivalued refinements of $f$ such that the asynchronous dynamics exhibits the reachability property $\mathcal{P}$.

We know that $\mathcal{P}$ is satisfied when the most permissive update is applied to the whole system. The question behind the MRBM method is: can $\mathcal{P}$ be satisfied by applying the most permissive update to only a subset of the components of the system (the other components being asynchronously updated)? If yes, this subset will contain the components that need to be multivalued for the asynchronous dynamics to satisfy $\mathcal{P}$. Thus, the MRBM method is based on the use of the partial $^J$-m.p. dynamics of the BM $f$ (defined in section 2.4), where $J$ is the set of components to be identified.

Since self-inhibited and output components are excluded from the sets $J$ (see Section 2.2), we will only consider the remaining components, called the admissible components. We initiate the method by examining whether there are sets $J$ of admissible components of size 1 that satisfy the reachability property $\mathcal{P}$ in the $^J$-m.p. dynamics of $f$. If there are no such set $J$ of size 1, we increment the set size by one and reexamine. This iterative process continues until we identify set(s) of size $l \le n$ that fulfill the reachability property $\mathcal{P}$.

The next step of the MRBM method is to build a multivalued refinement with the components $g_j \in J$ multivalued and the other components remaining Boolean. Towards this goal, the logical functions of each multivalued node and its targets need to be parameterised. We favor refinements where the maximum activity of the multivalued components is as low as possible (principle of parsimony). Hence, we decided to start this parameterisation with 3 levels of activity (i.e. $m_j = 2$ if $g_j \in J$) and to increase the maximal level of activity only if necessary. Finding the multivalued functions of the refinement consists in assigning a threshold (between 1 and $m_j$) to each of the outgoing regulations of the multivalued components. Defining multivalued logical functions such that the asynchronous dynamics of the model satisfies the reachability property $\mathcal{P}$ is not an easy task and several strategies exist. One strategy, that we call exhaustive strategy, assesses all possible logical functions, i.e. assesses all possible threshold assignations to the clauses of these logical functions (see Supplementary Data Section B). The main issue with this strategy is the combinatorial explosion of possibilities. Another strategy, that we call the *ad hoc* strategy, is to identify the sequence of regulations occurring in a path of interest
in the $^J$m.p. dynamics, and deduce the thresholds of regulation from it. Importantly, in both strategies, the proposed multivalued refinements can be supported or validated by the use of available biological data. This is especially relevant when using exhaustive strategy, as this process can actually generate multiple solutions. These solutions can be filtered using biological data to retain the most biologically relevant refinements. Once the logical functions have been defined, model checking approaches are used to verify the reachability property $\mathcal{P}$ (see Section 2.5).

### 3.1.2 Illustration of the MRBM method on a toy Boolean Model

To illustrate the MRBM method, we considered the toy BM $f$ of Figure 1A, which describes regulations between three components $g_1, g_2, g_3$. Both the asynchronous and the m.p. dynamics of the BM contain two fixed points and no cyclic attractor.

First, we considered a specific reachability property $\mathcal{P}_1$, the existence of a trajectory from the state 010 to the fixed point 011. We can observe that there are trajectories from 010 to 011 in the m.p. STG (in red in Figure 1F), but not in the asynchronous STG (Figure 1B). Our aim was to identify a multivalued refinement of the toy BM for which the reachability property $\mathcal{P}_1$ is satisfied in the asynchronous dynamics. To do this, we used the MRBM method. Following the method described in Section 3.1.1, we first searched for sets of admissible components $J$ of size one for which the reachability $010 \stackrel{f}{\rightsquigarrow}_{J_{m.p.}} 011$ is satisfied in the partial $^J$m.p. dynamics of the BM. Among all the sets of size 1 $\{g_1\}, \{g_2\}$ and $\{g_3\}$, only the set $J = \{g_1\}$ fulfils the required property (Figure 1G).

<figure><img src="image_4.png" /><figcaption>Figure 2: Multivalued Refinement of the toy model. (A) Regulatory graph of the multivalued refinement of the BM of Figure 1A, and associated logical functions. The rectangle represents the multivalued component $g_1$ ($m_1 = 2$); outgoing edges from $g_1$ have labels that specify the threshold levels from which the regulations occur. (B) Asynchronous STG of the multivalued refinement $h$. Each node represents a state $(x_1, x_2, x_3)$. Trajectories representing the reachability $010 \stackrel{h}{\rightsquigarrow}_{asyn} 011$ are highlighted in red.</figcaption></figure>

Therefore, component $g_1$ has to be multivalued. To define the multivalued functions of $g_1$ and its targets, we applied the ad-hoc strategy. We investigated the shortest trajectory from state 010 to the fixed point 011 in the $g_1$ m.p. dynamics of the toy model: $010 \to 000 \to i00 \to i01 \to i11 \to d11 \to 011$ (in each state, the
number in red indicates the component that changed during the state transition). In this trajectory, we observed an increase in the activity level of $g_3$ when $g_1$ is at an intermediate level ($i00 \to i01$). Hence, in a multivalued refinement with $m_1 = 2$, we proposed a threshold 1 for the activation from $g_1$ to $g_3$. Additionally, the level of activity of $g_2$ increases when its activator $g_3$ is present and its inhibitor $g_1$ is at an intermediate level ($i01 \to i11$). We deduced that $g_1$ needs to be fully active in order to inhibit $g_2$; so we set a threshold at level 2 for the inhibition from $g_1$ to $g_2$. The resulting proposed multivalued refinement is described in Figure 2A. We verified that its asynchronous dynamics satisfies the reachability property $\mathcal{P}_1$, that is the reachability from the state 010 to the fixed point 011 (see Figure 2B).

<table><caption>Table 1: Sizes of the Basins of Attraction of the toy Model Under Different Updating Schemes. Sizes of the basins of attraction of the fixed points 101 and 011 (as a percentage of the size of the Boolean state space) obtained when using either the asynchronous, m.p., or partial $J$-m.p. schemes ($J = \{g_1\}, \{g_2\}, \text{ or } \{g_3\}$). Red values signal any increase in the size of a basin of attraction when compared to the size of the basin of attraction of the asynchronous dynamics.</caption><thead><tr><th></th><th>Asyn</th><th>m.p</th><th><sup>g</sup><sup>1</sup>m.p.</th><th><sup>g</sup><sup>2</sup>m.p.</th><th><sup>g</sup><sup>3</sup>m.p.</th></tr></thead><tbody><tr><td><span style="color:red;">#B<sup>f</sup>(101)</span></td><td>87.5</td><td>87.5</td><td>87.5</td><td>87.5</td><td>87.5</td></tr><tr><td><span style="color:red;">#B<sup>f</sup>(011)</span></td><td>50</td><td><span style="color:red;">75</span></td><td><span style="color:red;">75</span></td><td>50</td><td>50</td></tr></tbody></table>

Secondly, we considered a set of reachability properties $\mathcal{P}_2$. We observed that the size of the basin of attraction of the fixed point 101 does not change with the updating scheme, whereas the size of the basin of attraction of the fixed point 011 changes ($\#B_{m.p.}^f(011) > \#B_{asyn}^f(011)$; see Table 1). As all the asynchronous trajectories are captured in the m.p. STG (Section 2.3), this change in size can be interpreted as a loss of reachability in the asynchronous dynamics. Our aim was to build a refinement restoring all these lost reachabilities (property $\mathcal{P}_2$).

Using MRBM, we found that we get $\#B_{m.p.}^f(011) = \#B_{m.p.}^f(011)$ with the partial m.p. scheme with $J = \{g_1\}$. Hence, the asynchronous dynamics of the multivalued refinement defined above allows us to recover a basin of attraction of the same size as for the BM with the m.p. dynamics.

## 3.2 Application of MRBM to refine a BM of Early Hematopoietic Stem Cell differentiation

The BM depicted in Figure 3A is a model of Hematopoietic Stem Cell (HSC) differentiation detailed in Hérault et al. (2023). The model contains 15 components. Its dynamics displays five attractors, all of which are fixed points, each representing a distinct HSC fate (i.e. a distinct final stage of differentiation): lymphoid (pLymph), neutrophils and mastocytes (pNeuMast), erythrocytes (pEr), megakaryocytes (pMk), and an inactive state (named zeros). We excluded the fixed point
A

<figure><img src="image_5.png" /></figure>

B

<figure><img src="image_6.png" /></figure>

C

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

D

<figure><img src="image_9.png" /><figcaption>Figure 3: Boolean Model of Early Hematopoietic Stem Cell Differentiation. (A) Regulatory graph of the BM (edges in red represent inhibitions and edges in green activations). The logical rules are detailed in Hérault et al. (2023). (B) Configuration (i.e. activity level of each component) of the 5 fixed points and the 5 transient key states (in column). Each cell of the table represents the activity level of a component (0=inactive, in white, 1 = active in dark green, * = 0 or 1, in light green). (C) Reachabilities of the fixed points (in column) from transient states (in row) in m.p. and asynchronous ("asyn) dynamics. ✓ stands for "reachable", and X for "not reachable". The cells are divided if there is a difference between the m.p. and the asynchronous dynamics. (D) Reachabilities of the fixed points from transient states (in row) in the partial J m.p. dynamics. The accessibility of transient states at fixed points in the partial dynamics of J m.p. is satisfied for the J sets specified in the table cell.</figcaption></figure>

zeros from our study as it was not discussed in Hérault et al. (2023). The authors of Hérault et al. (2023) have defined (from data analysis) five transient key stages in HSC differentiation: initial HSCs (iHSC), self-renewal HSCs (srHSC), quiescent HSCs (qHSC), pre-differentiating HSCs (preDiff), and pME (denoted "Transient states" in Figure 3B). We tested the fixed points reachabilities from these transient states in the m.p. and asynchronous dynamics (Figure 3C). The two dynamics display seven differences in terms of reachability (Figure 3C). With the asynchronous scheme, pLymph is not reachable from any of the transient stages of HSC differentiation, pNeuMast is not reachable from qHSC, and pEr is not reachable from qHSCs, pre-diff, and pME. We aimed to identify a multivalued refinement that
would provide all those seven reachability properties ($\mathcal{P}$) under the asynchronous scheme.

Using the MRBM method, we identified three admissible sets $J$ of size one for which the partial m.p. dynamics captures part of the desired reachabilities (Figure 3D):

- With $J = \{Fli1\}$, pEr is reachable from pME and preDiff.

- With $J = \{Gata2\}$ or $J = \{Spi1\}$, pLymph is reachable from srHSC.

Increasing the size of the sets $J$ by 1, we identified additional sets capturing additional reachabilities of interest:

- With $J = \{Fli1, Gata1\}$, pNeuMast and pEr are reachable from qHSC.

- With $J = \{Egr1, Gata2\}$, $J = \{Egr1, Spi1\}$, $J = \{Gata1, Gata2\}$, $J = \{Gata2, Junb\}$ or $J = \{Junb, Spi1\}$, pLymph is reachable from the state iHSC.

However, pLymph is still not reachable from qHSC. We found that the set $J = \{Egr1, Fli1, Gata2, Gata1\}$ is required for the property $\mathcal{P}$ to be satisfied in the partial $J$-m.p. dynamics.

<figure><img src="image_10.png" /><figcaption>Figure 4: Subgraph of a Multivalued Refinement of the BM of Early HSC Aging. Multivalued components are represented as rectangles, while other components are depicted as ovals. The label on arrows indicates the threshold level from which the regulation by a multivalued component occurs.</figcaption></figure>

We then built a multivalued refinement of the BM that satisfies properties ($\mathcal{P}$) under the asynchronous scheme. Finding the multivalued functions of the refinement consists in assigning a threshold 1 or 2 to each of the outgoing regulations of the 4 multivalued components. First we set to 2 the maximum activity level of the four components of $J$ ($Egr1$, $Fli1$, $Gata2$, and $Gata1$). By analyzing in the partial $J$-m.p. dynamics the paths corresponding to the seven reachabilities we want to recover (ad-hoc strategy, similarly to the path analysis in the case of the toy model), we decided to set the threshold of 4 regulations at level 2: the inhibition of $Klf1$ by $Fli1$, the inhibition of $Spi1$ by $Gata2$, the activation of $Fli1$ by $Gata1$, and the activation of $Junb$ by $Egr1$ (Figure 4). The logical functions of the multivalued refinement are available in Supplementary Data Section D.

As described in (Hérault et al. (2023)), available biological data support part of this parameterisation. Indeed, starting from a pre-differentiated stage (preDiff) cells can differentiate into either neutrophils/mastocytes (pNeuMast), erythrocytes
(pEr) or megakaryocytes (pMk). The circuit composed of *Gata1*, *Fli1*, and *Klf1* controls the cell fate decision between pEr and pMk. Moreover, an analysis of the dynamics of the BM shows that the switch between the pEr fixed point and pMk fixed point from the preDiff state depends on the existence of two distinct thresholds for *Fli1*'s influence on its targets *Klf1* and *Gata1* (see Hérault et al. (2023) and references therein).

## 3.3 Application of MRBM to refine a Boolean Model of Asymmetric Stem Cell Division in *Arabidopsis Thaliana* Root

As a second case study, we considered the BM of asymmetric stem cell division in the root of *Arabidopsis thaliana* presented in García-Gómez et al. (2020) (Figure 5A). The model has six attractors, all of which are fixed points, each of them corresponding to different cell types observed within the root of the plant, including the quiescent center (QC), peripheral and central pro-vascular initials (PPI and CPI), cortex/endodermis initials (CEI), columella initials (CoI), and transition domain (TD) (Figure 5B).

By comparing the asynchronous dynamics to the m.p. dynamics of the model, we noted that the sizes of the basins of attraction of three fixed points (QC, CoI, PPI) are larger within the m.p. dynamics than within the asynchronous dynamics (Figure 5C). As already mentioned in the toy model description (Section 3.1.2), one can interpret this size difference as a loss of reachability in the asynchronous dynamics. To identify a multivalued refinement of the BM whose asynchronous dynamics provides basins of attraction as large as in the m.p. dynamics, we applied the MRBM method. Considering sets *J* of admissible m.p. components of size one, we obtained:

$$ - \text{if } J = \{JKD\} \text{ or } J = \{SCR\}, \text{ then } \#\mathcal{B}_{J_{m.p.}}^{f}(PPI) = \#\mathcal{B}_{m.p.}(PPI), $$

$$ - \text{if } J = \{SHR\}, \text{ then } \#\mathcal{B}_{J_{m.p.}}^{f}(CoI) = \#\mathcal{B}_{m.p.}(CoI). $$

A set *J* of size two is required for the basin of attraction of the fixed point QC: if *J* = {MGP, SCR}, then \#B<sub>J<sub>m.p.</sub></sub><sup>f</sup>(QC) = \#B<sub>m.p.</sub><sup>f</sup>(QC). Thus, *J* = {SHR, MGP, SCR} is sufficient to ensure identical basin of attraction sizes between the partial m.p. and the m.p. dynamics (Figure 5C).

Thus, the use of MRBM efficiently restricts the search space for determining a multivalued model. In this example, the use of MRBM enabled us to identify precisely the three components on which it was relevant to search for multi-level parameterisation, without having to test all 816 possibilities (number of 3-element parts of an 18-element set).

The next step is the specification of multivalued refinements. We started looking for multivalued refinements with the maximum value of each multivalued component (that are the components in the previously identified sets *J*) fixed at 2. We used the exhaustive strategy, testing systematically all possible multivalued functions, i.e. all possible assignments of the thresholds (1 and 2) to each of the outgoing edges of
<figure><img src="image_11.png" /><figcaption>Figure 5: Boolean model of the Asymmetric Stem Cell Division in *Arabidopsis Thaliana* Root. (A) Regulatory graph of the BM; the logical functions are available in García-Gómez et al. (2020). (B) Summary of the 6 fixed points. Each row represents a component, each column corresponds to a fixed point, and each cells to the activity level (1 for active, in dark green, 0 for inactive, in white). (C) Sizes (in percentage) of the basins of attraction of the BM under different updating schemes. The rows represent basins of attraction and the columns correspond to the updating schemes: asynchronous ("Asyn"), m.p. and partial m.p. . For each partial $^J$m.p. dynamics, the contributing set $J$ is indicated in the column headers. We grouped in the same column the partial $^J$m.p. dynamics leading to same sizes of basin of attraction. We only report the smallest set(s) of m.p. components resulting in a specific size of basin of attraction. Values in red correspond to increases in the sizes of the basin of attraction as compared to the sizes obtained in the asynchronous dynamics.</figcaption></figure>

the multivalued components. We thus provided multivalued refinements which led to expected sizes of asynchronous basin of attraction for all cases except two ($J = \{SCR, MGP\}$ and $J = \{SCR, SHR, MGP\}$). Indeed, for $J = \{SCR, MGP\}$, when the maximum value of the multivalued components is restricted to 2, only one fixed point among QC and PPI has the same size of basins of attraction in both m.p. and asynchronous dynamics (top panel of Figure 6 represents all the possible situations). We thus increased the maximum value of activity of *SCR* to 3. Doing so, we were able to find a refinement which recovers the same sizes for the two basins of attraction (PPI and QC) in both dynamics (Figure 6, bottom panels). This 3-levels parameterisation of *SCR* also makes it possible to find a refinement with $J = SCR, SHR, MGP$ with which we found the same size for the three basins of attraction (PPI, QC and CoI) in both dynamics.

The resulting multivalued refinements are available in the Supplementary Data, Section D.
<figure><img src="image_12.png" /><figcaption>Figure 6: Subgraph of Multivalued Refinements of the BM of Asymmetric Stem Cell Division in *Arabidopsis thaliana* Root. Multivalued components are depicted as rectangles, while Boolean components are shown as ovals. The label on arrows indicates the threshold level from which the regulation by a multivalued component occurs. On the top panels $m_{SCR} = 2$, and on the bottom panels $m_{SCR} = 3$. We noted below each panel if we managed to recover the expected size(s) of basin of attraction depending on the set of multivalued components (see Figure 5C). The name of the basin of attraction is followed by a checkmark for a success, or a cross for a failure.</figcaption></figure>

# 4 Discussion

Although Boolean models have proved to be a valuable and effective tool for modeling biological interaction networks, they have certain limitations due to their high level of abstraction Wang et al. (2012). In particular, they capture the salient properties of the dynamics but might overlook other finer-grained properties Chamsedine and Rejniak (2020). Multivalued logical models can be useful to overcome such limitation, as they enable the consideration of more situations (Flobak et al. (2015), Selvaggio et al. (2021)). However, the refinement of a Boolean model into a multivalued model is not straightforward. This involves finding the components that need to be multivalued, as well as parameterising these components and their targets. The possibilities are very numerous and we face a combinatorial explosion. When multivaluing a Boolean model, the aim is to capture the expected properties sparingly, by minimising the number of multivalued components and their maximum level.

In this paper, we present MRBM, an original method to refine Boolean mod-
els. The aim of these refinements of Boolean models is to encompass specific dynamical properties $\mathcal{P}$ that appear in the most permissive dynamics but not in the asynchronous Boolean dynamics. We have defined a new updating scheme, the partial most permissive dynamics. This updating scheme applies the most permissive update only to a subset of network components, while asynchronous updating is applied to the other components. This allows to identify the components of the network that control the dynamical properties $\mathcal{P}$, and that would need to be multivalued. Hence, we created a multivalued refinement by multivaluing these components in the Boolean model. Doing so, we induce the dynamical properties $\mathcal{P}$ in the asynchronous dynamics of a multivalued refinement model. Importantly, the method is supported by mathematical arguments and proofs (see Supplementary Data). The effectiveness of MRBM is demonstrated on a toy model and two published Boolean models of cell differentiation, highlighting its ability to handle complex biological systems.

The MRBM method systematically identifies the components to be multivalued but the parameterisation itself remains challenging. We propose two main strategies: 1) An exhaustive strategy, in which all possible parameterisations are tested and then checked to verify if they satisfy the desired property. This strategy can be time consuming, especially as the maximum activity levels allowed for each component increase. 2) An *ad hoc* strategy, which analyses precisely relevant paths within the selected partial dynamics and thereby provides clues for deciphering the regulation thresholds. In both strategies, available a priori biological knowledges are used to guide the definition of the logical functions.

We applied the MRBM method to two published Boolean models of differentiation: a Boolean model of mammalian early hematopoietic stem cell differentiation and a Boolean model of asymmetric stem cell division in *Arabidopsis thaliana*. In both cases, we have identified reachability properties of biological interest that appear in the most permissive dynamics but not in the asynchronous dynamics. The MRBM method allowed us to successfully identify multivalued refinements of the Boolean models exhibiting these reachability properties under the asynchronous updating scheme. In the case of the early hematopoietic stem cell differentiation Boolean model, we proposed a multivalued refinement with four multivalued components (with maximal levels of 2). Existing biological data provide support for part of the parameterisation. For the Boolean model representing asymmetric stem cell division in *Arabidopsis thaliana*, we had to multivalue three components, one of them with a maximum activity level set at 3, while the others were set at 2. The maximal level of components is a feature that can be used to identify multivalued nodes requiring more detailed regulatory changes to recover the expected sizes of basins of attraction.

There are some limitations for the use of MRBM. The method is restricted to dynamics with only stable states (and no cyclic attractors) and focuses on the reachability properties of these stable states. A generalisation to the reachability properties of cyclic attractors is not straightforward. Indeed, cyclic attractors
depend on the updating scheme. Therefore there may not be a one-to-one corre-
spondence between the cyclic attractors of the most permissive dynamics and those
of the asynchronous dynamics. This makes it difficult to carry out the global anal-
ysis of the size of the basins of attraction, but a study of specific reachabilities
between an initial state and a state belonging to a cyclic attractor remains possible.
In this work, we have forbidden the multivaluation of a self inhibited component.
Indeed, the multivaluation of an auto-inhibited component leads to a change in
dynamics. It would be necessary, when studying a particular case, to observe this
new dynamics in detail before opting for the multivaluation of such a component.
The method faces computational issues for large models. We used the python
package mpbn to check the dynamical properties. The calculation of the size of the
basins of attractions was possible for relatively small models, but when studying
larger models (of more than 15 nodes), MRBM method was limited to the study of
specific reachability properties.

Finally, we can point out an intrinsic combinatorial limitation of the method.
The number of components to be multivalued (i.e. the size of the set J) can be as
large as n, the number of components of the model. Of course, the larger the set J,
the more difficult it is to find a parameterisation. It should be noted that difficulties
in finding an appropriate refinement may indicate that the initial Boolean model is
incomplete, reflecting gaps in our understanding of the system.

MRBM overall provides a systematic and innovative approach to refine Boolean
models with multivaluation. It exploits the richness of most permissive dynamics,
extracts its significant properties, and integrates them in the asynchronous Boolean
model to refine it. Importantly, the MRBM method does not rely on external data
sources, making it useful when experimental data are limited.

Acknowledgement

The project leading to this publication has received funding from France 2030, the French Government program managed by the French National Research Agency (ANR-16-CONV-0001), from Excellence Initiative of Aix-Marseille University - A*MIDEX (Turing Centre for Living Systems). This project was supported also supported by the ANR through the CAMUDI project (ANR-21-CE45-0001).

References

Ay, A. and Arnosti, D. N. (2011). Mathematical modeling of gene expression: a guide for the perplexed biologist. *Critical Reviews in Biochemistry and Molecular Biology*, 46(2):137-151.

Baier, C. and Katoen, J.-P. (2008). *Principles of model checking*. The MIT Press, Cambridge, Mass. OCLC: ocn171152628.
Bruggeman, F. J. and Westerhoff, H. V. (2007). The nature of systems biology. *Trends in Microbiology*, 15(1):45-50.

Chamseddine, I. M. and Rejniak, K. A. (2020). Hybrid modeling frameworks of tumor development and treatment. *Wiley Interdisciplinary Reviews. Systems Biology and Medicine*, 12(1):e1461.

Chen, W. W., Niepel, M., and Sorger, P. K. (2010). Classic and contemporary approaches to modeling biochemical reactions. *Genes & Development*, 24(17):1861-1875.

Cimatti, A., Clarke, E., Giunchiglia, E., Giunchiglia, F., Pistore, M., Roveri, M., Sebastiani, R., and Tacchella, A. (2002). NuSMV 2: An OpenSource Tool for Symbolic Model Checking. In Brinksma, E. and Larsen, K. G., editors, *Computer Aided Verification*, Lecture Notes in Computer Science, pages 359-364, Berlin, Heidelberg. Springer.

Flobak, Baudot, A., Remy, E., Thommesen, L., Thieffry, D., Kuiper, M., and Lægreid, A. (2015). Discovery of Drug Synergies in Gastric Cancer Cells Predicted by Logical Modeling. *PLoS computational biology*, 11(8):e1004426.

García-Gómez, M. L., Ornelas-Ayala, D., Garay-Arroyo, A., García-Ponce, B., Sánchez, M. d. l. P., and Álvarez Buylla, E. R. (2020). A system-level mechanistic explanation for asymmetric stem cell fates: Arabidopsis thaliana root niche as a study system. *Scientific Reports*, 10(1):3525.

Garg, A., Di Cara, A., Xenarios, I., Mendoza, L., and De Micheli, G. (2008). Synchronous versus asynchronous modeling of gene regulatory networks. *Bioinformatics* (*Oxford, England*), 24(17):1917-1925.

Gebser, M., Kaminski, R., Kaufmann, B., and Schaub, T. (2013). *Answer Set Solving in Practice*. Synthesis Lectures on Artificial Intelligence and Machine Learning. Springer International Publishing, Cham.

Gebser, M., Kaminski, R., Kaufmann, B., and Schaub, T. (2014). Clingo = ASP + Control: Preliminary Report.

Heinrich, R., Neel, B. G., and Rapoport, T. A. (2002). Mathematical models of protein kinase signal transduction. *Molecular Cell*, 9(5):957-970.

Hérault, L., Poplineau, M., Duprez, E., and Remy, (2023). A novel Boolean network inference strategy to model early hematopoiesis aging. *Computational and Structural Biotechnology Journal*, 21:21-33.

Klarner, H., Heinitz, F., Nee, S., and Siebert, H. (2020). Basins of Attraction, Commitment Sets, and Phenotypes of Boolean Networks. *IEEE/ACM transactions on computational biology and bioinformatics*, 17(4):1115-1124.
Naldi, A. (2018). *BioLQM: A Java Toolkit for the Manipulation and Conversion of Logical Qualitative Models of Biological Networks*. *Frontiers in Physiology*, 9:1605.

Naldi, A., Hernandez, C., Levy, N., Stoll, G., Monteiro, P. T., Chaouiya, C., Helikar, T., Zinovyev, A., Calzone, L., Cohen-Boulakia, S., Thieffry, D., and Paulevé, L. (2018). The CoLoMoTo Interactive Notebook: Accessible and Reproducible Computational Analyses for Qualitative Biological Networks. *Frontiers in Physiology*, 9:680.

Paulevé, L., Kolčák, J., Chatain, T., and Haar, S. (2020). Reconciling qualitative, abstract, and scalable modeling of biological networks. *Nature Communications*, 11(1):4256.

Paulevé, L. and Sené, S. (2022). Boolean Networks and Their Dynamics: The Impact of Updates. In De Maria, E., editor, *Systems Biology Modelling and Analysis*, pages 173–250. Wiley, 1 edition.

Schwab, J. D., Kühlwein, S. D., Ikonomi, N., Kühl, M., and Kestler, H. A. (2020). Concepts in Boolean network modeling: What do they all mean? *Computational and Structural Biotechnology Journal*, 18:571–582.

Selvaggio, G., Chaouiya, C., and Janody, F. (2021). In Silico Logical Modelling to Uncover Cooperative Interactions in Cancer. *International Journal of Molecular Sciences*, 22(9):4897.

Thomas, R. (1973). Boolean formalization of genetic control circuits. *Journal of Theoretical Biology*, 42(3):563–585.

Thomas, R. (1978). Logical analysis of systems comprising feedback loops. *Journal of Theoretical Biology*, 73(4):631–656.

Tomlin, C. J. and Axelrod, J. D. (2007). Biology by numbers: mathematical modelling in developmental biology. *Nature Reviews. Genetics*, 8(5):331–340.

Torres, N. V. and Santos, G. (2015). The (Mathematical) Modeling Process in Biosciences. *Frontiers in Genetics*, 6:354.

Wang, R.-S., Saadatpour, A., and Albert, R. (2012). Boolean modeling in systems biology: an overview of methodology and applications. *Physical Biology*, 9(5):055001.