# Refining Boolean models with the partial most permissive scheme

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

► To cite this version:

Nadine Ben Boina, Brigitte Mossé, Anaïs Baudot, Elisabeth Remy. Refining Boolean models with the partial most permissive scheme. Bioinformatics, 2025, 41 (4), pp.btaf123. 10.1093/bioinformatics/btaf123. hal-05226720

HAL Id: hal-05226720

https://amu.hal.science/hal-05226720v1

Submitted on 30 Jan 2026

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Systems biology

## Refining Boolean models with the partial most permissive scheme

Nadine Ben Boina<fnref n="1" />,<fnref n="2" />, Brigitte Mossé<fnref n="1" />, Anaïs Baudot<fnref n="2" />,<fnref n="3" />,<fnref n="4" />✉, Elisabeth Remy<fnref n="1" />,*

### Abstract

**Motivation:** In systems biology, modeling strategies aim to decode how molecular components interact to generate dynamical behavior. Boolean modeling 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:** The MRBM method and the BMs used in this article are available on GitHub at: https://github.com/NdnBnBn/MRBM, and archived in Zenodo (doi: 10.5281/ZENODO.14979798).

### 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 behaviors (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 behavior 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 is 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—sometimes denoted fully asynchronous (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 behaviors that are not revealed with the synchronous scheme (Paulevé and Sené 2022). The most permissive updating

<footnote n="1">¹Aix Marseille Univ, CNRS, I2M (UMR 7373), Turing Center for Living systems, Marseille, France</footnote>

<footnote n="2">²Aix Marseille Univ, INSERM, MMG, Marseille, France</footnote>

<footnote n="3">³CNRS, Marseille, France</footnote>

<footnote n="4">⁴Barcelona Supercomputing Centre, Barcelona, Spain</footnote>

<footnote>*Corresponding author. Aix Marseille Univ, CNRS, I2M (UMR 7373), Turing Center for Living systems, Marseille, France. E-mail: elisabeth.remy@univ-amu.fr.</footnote>

<footnote>Associate Editor: Macha Nikolski</footnote>
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 each component two intermediate activity levels which mark the transition between the two Boolean levels. These intermediate levels represent the fact that the component 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 behaviors 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 Multivalued Refinement of Boolean Model (MRBM), 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 proposing 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 multivalued. In this article, 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 Materials and methods

## 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, ..., f_n) : \mathbb{B}^n \to \mathbb{B}^n$, where each Boolean state $x = (x_1, ..., x_n) \in \mathbb{B}^n$ specifies the activity level of $n$ components of the model $g_1, ..., 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, ..., 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 (Fig. 1A). Note that dual edges are allowed. Dual edges

<figure><img src="image_2.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$ 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>

represent interactions that can be activations or inhibitions depending on the context.

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$ (Fig. 1B). More precisely, there is a transition from the state $x = (x_1, ..., x_n)$ to the state $y = (y_1, ..., y_n)$ if there exists $i_0 \in \{1, ..., n\}$ such that $x_{i_0} \neq f_{i_0}(x)$ and
$$ 
\begin{aligned}
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.
\end{aligned}
$$

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{f_{asyn}} y$. A trajectory from state $x$ to state $y$ is a sequence of transitions $x \xrightarrow{f_{asyn}} x' \xrightarrow{f_{asyn}} \dots \xrightarrow{f_{asyn}} y$. A state $y$ is reachable from $x$ if there exists a trajectory from $x$ to $y$. We denote this reachability by $x \xrightarrow{f_{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 $B_{asyn}(\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_{i=1}^n \{0, 1, \dots, m_j\}$. A multivalued model on $X$ is a map $b = (b_1, \dots, b_n) : X \to X$ satisfying the condition $b_i(x) - x_i \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_{i=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)\}$. For instance, if $X = \{0, 1, 2\}^3$, $\alpha(021) = \{010; 011\}$. A multivalued model $b$ 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\}$,

$$ 
\begin{aligned}
h_j(x) < x_j &\Rightarrow \exists x' \in \alpha(x), f_i(x') < x'_i, \\
h_j(x) > x_j &\Rightarrow \exists x' \in \alpha(x), f_i(x') > x'_i.
\end{aligned}
$$

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 $b$. 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 Section SA).

In the remaining sections of this article, all the refinements we consider are built as described above. In the Section 3, we describe these refinements $b$ through the functions $H_j$ of $x \in X$ with values in $\{-1, +1\}$. These functions $H_j$ indicate

whether gene $j$ is called to decrease or to increase its level $x_j$ under $b$ (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_{i=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 Section SB).

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.}$:

$$ 
\begin{aligned}
\gamma(x) &= \{x' \in \mathbb{B}^n; \forall j \in \{1, \dots, n\}, \\
&\quad (x_j = 0 \Rightarrow x'_j = 0) \text{ and } (x_j = 1 \Rightarrow x'_j = 1)\}
\end{aligned}
$$

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, \dots, n\}$ such that one of the following occurs:

$$ 
\begin{aligned}
-x_{j_0} &\in \{0, d\} \text{ and } \exists x' \in \gamma(x) |f_{j_0}(x') = 1 \text{ and } y_{j_0} = i, \\
-x_{j_0} &\in \{1, i\} \text{ and } \exists x' \in \gamma(x) |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 \xrightarrow{d} 1 \xrightarrow{0} 1 \xrightarrow{0} 1 \xrightarrow{0} 01$ (Fig. 1C) is represented by the transition $111 \to 001$ (Fig. 1D). The full m.p. STG of the toy model is represented (Fig. 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, \dots, 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_{l,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_{l,m,p.}$ in the following way: there is a transition from a state $x \in X_{l,m,p.}$ to a state $y \in X_{l,m,p.}$ if there exists $j_0 \in \{1, \dots, n\}$ such that one of the following occurs:

$$ 
\begin{align*}
& -j_0 \in J \text{ and } x_{j_0} \in \{0, d\} \text{ and } \exists x' \in \gamma(x) \mid 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) \mid 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) \mid 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) \mid f_{j_0}(x') = 0 \text{ and } y_{j_0} = 0,
\end{align*}
$$

and $y_i = x_i$, 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 $g_1$-m.p. STG of the toy model is represented (Fig. 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 Section SC).

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_{l,m,p.}$ is contained in $X_{l_2,m,p.}$; moreover, if there exists a trajectory between two states in the $J_1$-m.p. dynamics of $f$, the same holds in the $J_2$-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) and 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 $f$ 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'(\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'(\omega)}{2^n} \times 100\%$).

## 2.6 Availability and implementation

The MRBM method and the BMs used in this article are available on GitHub at: https://github.com/NdnBnBn/MRBM, and archived in Zenodo (doi: 10.5281/ZENODO.14979798). 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 Section 3. 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 dynam- ics. The MRBM method aims to identify multivalued refine- ments 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 $^1$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 reachabil- ity property $\mathcal{P}$ in the $^1$m.p. dynamics of $f$. If there are no such set $J$ of size 1, we increment the set size by one and reex- amine. 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 multival- ued refinement with the components $g_i \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 parameterized. 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 parameterization with 3 levels of activity (i.e. $m_i = 2$ if $g_i \in J$) and to increase the maximal level of activity only if necessary. Finding the multivalued functions of the refine- ment consists in assigning a threshold (between 1 and $m_i$) to each of the outgoing regulations of the multivalued compo- nents. 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 Section SB). The main issue with this strategy is the combinatorial explosion of possibilities. Another strat- egy, that we call the *ad hoc* strategy, is to identify the se- quence of regulations occurring in a path of interest in the $^1$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 avail- able biological data. This is especially relevant when using ex- haustive strategy, as this process can actually generate multiple solutions. These solutions can be filtered using bio- logical data to retain the most biologically relevant refine- ments. Once the logical functions have been defined, model checking approaches are used to verify the reachability prop- erty $\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 Fig. 1A, which describes regulations between three compo- nents $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 the fixed point 011. We can observe that there are trajectories from 010 to 011 in the m.p. STG STG (in red in Fig. 1F), but not in the asynchronous STG (Fig. 1B). Our aim was to iden- tify a multivalued refinement of the toy BM for which the reachability property $\mathcal{P}_1$ is satisfied in the asynchronous dy- namics. 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 $0.10 \to_{J,m.p.}^f$ 0.11 is satisfied in the partial $^1$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\}$ fulfills the required prop- erty (Fig. 1G).

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^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 ob- served an increase in the activity level of $g_3$ when $g_1$ is at an intermediate level ($i00 \to i01$). Hence, in a multivalued refine- ment with $m_1 = 2$, we proposed a threshold 1 for the activa- tion 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 Fig. 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 Fig. 2B).

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 ($\#\mathcal{B}_{m.p.}^f(011) > \#\mathcal{B}_{asyn}^f(011)$; see Table 1). As all

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Multivalued refinement of the toy model. (A) Regulatory graph of the multivalued refinement of the BM of Fig. 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 \xrightarrow{h}_{asyn} 011$ are highlighted in red.</figcaption></figure>
<table><caption>**Table 1.** Sizes of the basins of attraction of the toy model under different updating schemes.ª</caption><thead><tr><th></th><th>Asyn</th><th>m.p</th><th>g<sup>1</sup>m.p</th><th>g<sup>2</sup>m.p</th><th>g<sup>3</sup>m.p</th></tr></thead><tbody><tr><td>#B<sup>f</sup> (101)</td><td>87.5</td><td>87.5</td><td>87.5</td><td>87.5</td><td>87.5</td></tr><tr><td>#B<sup>f</sup> (011)</td><td>50</td><td>75</td><td>75</td><td>50</td><td>50</td></tr></tbody></table>

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 $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 Fig. 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 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 Fig. 3B). We tested the fixed points reachabilities from these transient states in the m.p. and asynchronous dynamics (Fig. 3C). The two dynamics display seven differences in terms of reachability (Fig. 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 ($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 (Fig. 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 $P$ to be satisfied in the partial ¹m. p. dynamics.

We then built a multivalued refinement of the BM that satisfies properties ($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 ¹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$ (Fig. 4). The logical functions of the multivalued refinement are available in Supplementary Section SD.

As described in Hérault et al. (2023), available biological data support part of this parameterization. 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 Garcia-Gómez et al. (2020) (Fig. 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/endo dermis initials (CEI), columella initials (CoI), and transition domain (TD) (Fig. 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 (Fig. 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:

* if $J = \{JKD\}$ or $J = \{SCR\}$, then $\#B_{m.p.}^f(PPI) = \#B_{m.p.}^f(PPI)$,

* if $J = \{SHR\}$, then $\#B_{m.p.}^f(CoI) = \#B_{m.p.}^f(CoI)$.

<footnote>ª 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 ¹m.p. schemes ($J = \{g_1\}, \{g_2\}$, 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.</footnote>
<figure><img src="image_4.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 five fixed points and the five 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 $\sqrt{m.p.}$ dynamics. The accessibility of transient states at fixed points in the partial dynamics of $\sqrt{m.p.}$ is satisfied for the $J$ sets specified in the table cell.</figcaption></figure>

A set $J$ of size two is required for the basin of attraction of the fixed point QC: if $J = \{MGP, SCR\}$, then $\#B_{m.p.}^f(QC) = \#B_{m.p.}^f(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 (Fig. 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 parameterization, 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 the multivalued components. We thus provided multivalued refinements which led to expected sizes of asynchronous basins 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 asyn-
chronous dynamics (top panel of Fig. 6 represents all the pos-
sible situations). We thus increased the maximum value of
activity of SCR to 3. Doing so, we were able to find a refine-
ment which recovers the same sizes for the two basins of

attraction (PPI and QC) in both dynamics (Fig. 6, bottom panels). This 3-levels parameterization of SCR also makes it possible to find a refinement with $J = SCR, SHR, MGP$ 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 Section SD*.

<figure><img src="image_5.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>

## 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 ab-
straction (Wang *et al.* 2012). In particular, they capture the
salient properties of the dynamics but might overlook other
finer-grained properties (Chamseddine 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 re-
finement of a Boolean model into a multivalued model is not
straightforward. This involves finding the components that
need to be multivalued, as well as parameterizing these com-
ponents and their targets. The possibilities are very numerous

<figure><img src="image_6.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 six 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 m.p. dynamics, the contributing set J is indicated in the column headers. We grouped in the same column the partial 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>
<figure><img src="image_7.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 Fig. 5C). The name of the basin of attraction is followed by a checkmark for a success, or a cross for a failure.</figcaption></figure>

and we face a combinatorial explosion. When multivaluing a Boolean model, the aim is to capture the expected properties sparingly, by minimizing the number of multivalued components and their maximum level.

In this paper, we present MRBM, an original method to refine Boolean models. 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 parameterization itself

remains challenging. We propose two main strategies: (i) An exhaustive strategy, in which all possible parameterizations 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 components increase. (ii) 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 parameterization. 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 generalization 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 correspondence between the cyclic attractors of the most permissive dynamics and those of the asynchronous dynamics. This makes it difficult to carry out the global analysis 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 multi-valued (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 parameterization. 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.

## Author contributions

Nadine Ben Boina (Conceptualization [equal], Data curation [equal], Investigation [equal], Methodology [equal], Software [lead], Writing—original draft [equal], Writing—review & editing [equal]), Brigitte Mossé (Supervision [equal], Writing—review & editing [equal]), Anaïs Baudot (Conceptualization [equal], Data curation [equal], Formal analysis [equal], Funding acquisition [equal], Investigation [equal], Methodology [equal], Project administration [equal], Supervision [equal], Writing—original draft [equal], Writing—review & editing [equal]), and Elisabeth Remy (Conceptualization [equal], Data curation [equal], Formal analysis [equal], Funding acquisition [equal], Investigation [equal], Methodology [equal], Project administration [equal], Supervision [equal], Writing—original draft [equal], Writing—review & editing [equal])

## Supplementary data

Supplementary data are available at *Bioinformatics* online.

Conflict of interest: None declared.

## Funding

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 also supported by the ANR through the CAMUDI project [ANR-21-CE45-0001].

## References

Ay A, Arnosti DN. Mathematical modeling of gene expression: a guide for the perplexed biologist. *Crit Rev Biochem Mol Biol* 2011; 46:137–51.

Baier C, Katoen J-P. *Principles of Model Checking*. Cambridge, MA: The MIT Press, 2008.

Bruggeman FJ, Westerhoff HV. The nature of systems biology. *Trends Microbiol* 2007;15:45–50.

Chamseddine IM, Rejniak KA. Hybrid modeling frameworks of tumor development and treatment. *Wiley Interdiscip Rev Syst Biol Med* 2020;12:e1461.

Chen WW, Niepel M, Sorger PK. Classic and contemporary approaches to modeling biochemical reactions. *Genes Dev* 2010;24:1861–75.

Cimatti A, Clarke E, Giunchiglia E et al. NuSMV 2: an OpenSource tool for symbolic model checking. In: Brinksma E, Larsen KG (eds.), *Computer Aided Verification, Lecture Notes in Computer Science*. Berlin, Heidelberg: Springer, 2002, 359–64.

Flobak A, Baudot A, Remy E et al. Discovery of drug synergies in gastric cancer cells predicted by logical modeling. *PLoS Comput Biol* 2015; 11:e1004426.

García-Gómez ML, Ornelas-Ayala D, Garay-Arroyo A et al. A system-level mechanistic explanation for asymmetric stem cell fates: *Arabidopsis thaliana* root niche as a study system. *Sci Rep* 2020; 10:3525.

Garg A, Di Cara A, Xenarios I et al. Synchronous versus asynchronous modeling of gene regulatory networks. *Bioinformatics* 2008; 24:1917–25.

Gebser M, Kaminski R, Kaufmann B et al. Answer Set Solving in Practice. *Synthesis Lectures on Artificial Intelligence and Machine Learning*. Cham: Springer International Publishing, 2013.

Gebser M, Kaminski R, Kaufmann B et al. *Clingo = ASP + Control: Preliminary Report*. 2014. https://doi.org/10.48550/arXiv.1405.3694

Heinrich R, Neel BG, Rapoport TA. Mathematical models of protein kinase signal transduction. *Mol Cell* 2002;9:957–70.

Hérault L, Poplineau M, Duprez E et al. A novel Boolean network inference strategy to model early hematopoiesis aging. *Comput Struct Biotechnol* 2023;21:21–33.

Klarner H, Heinitz F, Nee S et al. Basins of attraction, commitment sets, and phenotypes of boolean networks. *IEEE/ACM Trans Comput Biol Bioinform* 2020;17:1115–24.

Naldi A. BioLQM: a java toolkit for the manipulation and conversion of logical qualitative models of biological networks. *Front Physiol* 2018;9:1605.

Naldi A, Hernandez C, Levy N et al. The CoLoMoTo interactive notebook: accessible and reproducible computational analyses for qualitative biological networks. *Front Physiol* 2018;9:680.

Paulevé L, Kolčák J, Chatain T et al. Reconciling qualitative, abstract, and scalable modeling of biological networks. *Nat Commun* 2020; 11:4256.
Paulevé L, Sené S. Boolean networks and their dynamics: the impact of updates. In: De Maria E (ed.), *Systems Biology Modelling and Analysis*, 1 edn. Wiley, 2022, 173–250.

Schwab JD, Kühlwein SD, Ikonomi N *et al*. Concepts in Boolean network modeling: what do they all mean? *Comput Struct Biotechnol J* 2020;18:571–82.

Selvaggio G, Chaouiya C, Janody F. In silico logical modelling to uncover cooperative interactions in cancer. *Int J Mol Sci* 2021;22:4897.

Thomas R. Boolean formalization of genetic control circuits. *J Theor Biol* 1973;42:563–85.

Thomas R. Logical analysis of systems comprising feedback loops. *J Theor Biol* 1978;73:631–56.

Tomlin CJ, Axelrod JD. Biology by numbers: mathematical modelling in developmental biology. *Nat Rev Genet* 2007; 8:331–40.

Torres NV, Santos G. The (mathematical) modeling process in biosciences. *Front Genet* 2015;6:354.

Wang R-S, Saadatpour A, Albert R. Boolean modeling in systems biology: an overview of methodology and applications. *Phys Biol* 2012;9:055001.