# Data-driven inference of Boolean networks from transcriptomes to predict cellular differentiation and reprogramming

Stéphanie Chevalier, Julia Becker, Yujuan Gui, Vincent Noël, Cui Su, Sascha Jung, Laurence Calzone, Andrei Zinovyev, Antonio del Sol, Jun Pang, et al.

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

► To cite this version:

Stéphanie Chevalier, Julia Becker, Yujuan Gui, Vincent Noël, Cui Su, et al. Data-driven inference of Boolean networks from transcriptomes to predict cellular differentiation and reprogramming. npj Systems Biology and Applications, 2025, 11, pp.105. 10.1038/s41540-025-00569-z. hal-04753079v2

HAL Id: hal-04753079

https://hal.science/hal-04753079v2

Submitted on 26 Sep 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>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Data-driven inference of Boolean networks from transcriptomes to predict cellular differentiation and reprogramming

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

Check for updates

Stéphanie Chevalier Ⓔ<fnref n="1" />,², Julia Becker³, Yujuan Gui Ⓔ³, Vincent Noël Ⓔ⁴,⁵,⁶, Cui Su⁷, Sascha Jung Ⓔ⁸, Laurence Calzone Ⓔ⁴,⁵,⁶, Andrei Zinovyev Ⓔ⁹, Antonio del Sol⁸,¹⁰,¹¹, Jun Pang Ⓔ⁷, Lasse Sinkkonen Ⓔ³, Thomas Sauter Ⓔ³ & Loïc Paulevé Ⓔ¹²✉

Boolean networks provide robust, explainable, and predictive models of cellular dynamics, especially for cellular differentiation and fate decision processes. Yet, the construction of such models is extremely challenging, as it requires integrating prior knowledge with experimental observation of the transcriptome, potentially relating thousands of genes. We present a general methodology for integrating transcriptome data and prior knowledge on the underlying gene regulatory network in order to generate automatically ensembles of Boolean networks able to reproduce the modeled qualitative behavior. Our methodology builds on the software BoNesis, which implements the automatic construction of Boolean networks from a specification of their expected structural and dynamical properties. We show how to transform transcriptome data into such a qualitative specification, and then how to exploit the generated ensembles of Boolean networks for identifying families of candidate models, and for predicting robust cellular reprogramming targets. We illustrate the scalability and versatility of our overall approach with two applications: the modeling of hematopoiesis from single-cell RNA-Seq data, and modeling the differentiation of bone marrow stromal cells into adipocytes and osteoblasts from bulk RNA-seq time series data. For this latter case, we took advantage of ensemble modeling to predict combinations of reprogramming factors for trans-differentiation that are robust to model uncertainties due to variations in experimental replicates and choice of binarization method. Moreover, we performed an in silico assessment of the fidelity and efficiency of the reprogramming and conducted preliminary experimental validation.

Mathematical models have demonstrated their utility in elucidating experimental findings that might challenge intuitive comprehension. Through meticulous depiction of interactions within intricate signaling pathways and by contextualizing the dynamics of gene expression, these models provide a systematic approach to unveil the regulatory mechanisms that control cellular processes and their dysregulation in diseases. Among the various mathematical formalisms, the Boolean network (Boolean network) is a simple but expressive formalism that relies on pragmatic rules to

qualitatively simulate essential systems' features. It is notably valuable in poorly understood large-scale systems, as it can be employed for systems with hundreds of components and as the inference of Boolean network models, contrary to quantitative models (typically ordinary differential equation (ODE)-based models), does not require kinetic parameters derived from in-depth and often unavailable knowledge.

Consequently, Boolean networks are increasingly used to capture the interaction dynamics within complex signaling pathways and regulatory

<footnote n="1">¹Translational Medicine, Servier, Suresnes, France. ²LISN, Univ. Paris-Saclay, CNRS, Paris, France. ³Univ. Luxembourg, Department of Life Sciences and Medicine, 6 avenue du Swing, L-4367 Belvaux, Luxembourg, France. ⁴Institut Curie, Université PSL, F-75005 Paris, France. ⁵INSERM, U900, F-75005 Paris, France. ⁶Mines ParisTech, Université PSL, F-75005 Paris, France. ⁷Univ. Luxembourg, Department of Computer Science, 6 avenue de la Fonte, L-4364 Esch-sur-Alzette, Luxembourg. ⁸Computational Biology Group, Luxembourg Centre for Systems Biomedicine (LCSB), University10 of Luxembourg, L-4362 Esch-sur-Alzette, Luxembourg. ⁹In silico R&D, Evotec, Toulouse, France. ¹⁰Computational Biology Group, CIC bioGUNE-BRTA (Basque Research and Technology7 Alliance), Bizkaia Technology Park, Derio, Spain. ¹¹Ikerbasque, Basque Foundation for Science, Bilbao, Bizkaia, 48012, Spain. ¹²Univ. Bordeaux, CNRS, Bordeaux INP, LaBRI, UMR 5800, F-33400 Talence, France. ✉e-mail: loic.pauleve@labri.fr</footnote>
mechanisms governing cellular behaviors. They have been inferred from high-throughput data for modeling a range of biologically meaningful phenomena such as the mammalian cell cycle¹, cell differentiation and specifications²⁻⁶, stress/aging-related cell behaviors⁷,⁸, cell apoptosis⁹, and cancer cell functions¹⁰⁻¹³. Recent endeavors have focused on enhancing the quantitative interpretation of the resulting models, incorporating probabilistic approaches to effectively simulate heterogeneous cell populations and dynamical interacting populations¹⁴,¹⁵, and introducing a semantics that offers the formal guarantee of completely capturing any behavior achievable by any quantitative model (multilevel or ODE) following the same logic¹⁶.

A Boolean network consists of logical rules that are associated with each variable and that predict how its state evolves with time. The regulators of a variable are combined with logical connectors or, and, and not and define conditions for a variable to be at 0 (false/inactive/absent) or (1/true/active/ present). Through structural and dynamic analyses along with simulations under perturbations, Boolean networks provide versatile opportunities for exploring mechanisms underlying biological phenomena. This versatility allows them to serve as standalone informative tools but also, for specific modeling contexts and needs, such as in systems pharmacology, to lay the groundwork for more detailed pharmacokinetic/pharmacodynamic and quantitative systems pharmacology (QSP) models using ODEs¹⁷⁻²⁰.

One the prominent challenges for applications of Boolean networks in biology is the design of their logical rules. Indeed, besides the inference of the underlying gene regulatory network, the search for logical rules faces a double combinatorial explosion: the number of possible logical rules for a single variable is exponential with the number of its regulators, and checking whether a candidate Boolean network possesses the desired dynamical features (steady states, trajectories, ...) can involve analysis that take time and memory exponential with the number of variables. Thus, in most applications of the literature, the Boolean networks have been manually designed from expert knowledge and by re-utilizing previously published models. Nevertheless, there has been recent progress on the inference of Boolean models from data²¹⁻²⁸. These methods address the inference problem with different restrictions, either on the type of data and their interpretation in order to obtain simpler dynamical properties, or on restricting the set of logical rules to those having a particular structure, in order to reduce the combinatorics of candidate models. However, these methods remain difficult to scale above hundreds of variables, and typically enforce a rather specific way of interpreting the experimental data in terms of Boolean properties.

In this paper, we present a general methodology for the inference of Boolean networks from knowledge, data, and expert interpretation of data, and the computation of predictions from the resulting models. The methodology, summarized in Fig. 1, builds on the following steps:

1. The modeling of the knowledge, essentially in terms of an admissible structure for the models.

2. The qualitative modeling of the data in terms of expected dynamical properties of the model. This step depends on the biological expertise of the system, and relies on data analysis, notably to classify gene expression into binary values.

3. The tool BoNesis²⁹, which integrates (1) and (2) and, using logic programming and combinatorial optimization algorithms, infers ensembles of Boolean networks that are compatible with modeled static and dynamical properties.

4. The analysis of sampled ensembles of models to perform predictions, including key genes and reprogramming mutations.

The modeling steps 1 and 2, that define the inference problem, allow a versatile pipeline, that is not tied to specific type of data, or a specific interpretation of them. Indeed, our approach aims at moving the modeling effort from the design of Boolean rules to the specification of the expected features of the model. Then, BoNesis relies on state-of-the-art symbolic artificial intelligence technologies to automatically construct models that satisfy the desired properties.

We showcase the implementation of this methodology on two case studies: the inference of Boolean networks from scRNA-seq data of hematopoiesis, with the identification of key genes and the analysis of families of candidate models; and the prediction of reprogramming targets for adipocyte to osteoblast conversion from bulk RNA-seq time series data.

Compared to the current state of the art (see “Methods”), our methodology enables a scalable data-driven approach from different types of experimental datasets, including single-cell or bulk RNA sequencing. We build on existing software bricks for data analysis, trajectory reconstruction, gene activity classification and generic Boolean network inference from qualitative specification. For the first case study, we leverage trajectory reconstruction and scRNA-seq binarization methods to translate scRNA-seq data into Boolean trajectory specification. The resulting ensemble of models shows a substantial intersection with manually designed models. Moreover, we show how Boolean networks can be clustered in order to identify subfamilies of compatible models that could be distinguished by the differences in the Boolean function complexity of a few genes. This demonstrated that our methodology enables the data-driven automatic identification of key genes in the hematopoietic process, as well as the ability to access the diversity and subfamilies of compatible Boolean networks. For the second case study, we show how we can take advantage of several bulk RNA-seq binarization methods to generate diverse ensembles of Boolean networks that reproduce differentiation processes observed in time series. Additionally, we propose new methods and metrics for computing and

<figure><img src="image_5.png" /><figcaption>**Fig. 1 | General principle of the conducted inference of ensembles of Boolean networks.** By offering a generic modeling language, BoNesis enables integrating prior knowledge on regulation mechanisms with different types of experimental data, after qualitative interpretation, which may depend on biological hypotheses and experimental systems. These inputs specify what an admissible model is. Then, employing logic programming, BoNesis can generate ensembles of Boolean networks that fulfill the structural and dynamical properties and, by combining with combinatorial optimization technologies, enables the prediction of important genes and derive predictions for attractor reprogramming.</figcaption></figure>
ranking combinations of mutations for triggering attractor change from such heterogeneous ensembles, aiming at predicting robust cellular reprogramming targets. Finally, we performed new scRNA-seq experiments to validate these predictions.

Importantly, these case studies demonstrate that our approach is scalable to TF-scale networks, by starting from complete TF–TF regulatory networks and, through the inference of Boolean networks comprising thousands of nodes, we are able to automatically identify subnetworks that are sufficient to explain the observed dynamics. This scalability to TF-scale network while accounting for complex dynamical properties is another advantage of our approach compared to other state-of-the-art methods, often limited to hundreds of genes.

## Results

### Case study 1: ensemble modeling of hematopoiesis from scRNA-seq data

We applied our inference pipeline (Fig. 2A) to the identification of key genes and Boolean rules that can explain the hematopoiesis observed in a mouse sample using scRNA-seq data from Nestorowa et al.³⁰. Employing trajectory reconstruction and binarization methods, we derived a logical specification of the differentiation dynamics. Then, using BoNesis, we considered any Boolean network employing TF regulations referenced in the DoRothEA database, and automatically identified the sparsest among them that are able to reproduce the differentiation dynamics. We compared the selected genes of importance with an export model of the literature, showing a substantial overlap. Finally, we highlighted the advantage of the ensemble modeling by analyzing the variability of Boolean models compatible with the input data. We notably performed clustering of sampled models, resulting in clear 3 subfamilies of models that can be distinguished on specific features of Boolean rules.

Hematopoiesis is a crucial differentiation process of blood cells for immune system regeneration. It has been extensively studied, including with mathematical and logic dynamical models³⁻⁵,³¹.

We focused on data-driven modeling of the early differentiation of mouse hematopoietic stem cells (HSCs) from scRNA-seq data of Nestorowa et al.³⁰. The data show heterogeneity of cells during HSCs differentiation, including lympho-myeloid primed progenitors (LMPPs) and common myeloid progenitors (CMPs), further differentiated into granulocyte-monocyte progenitors (GMPs) and megakaryocyte-erythrocyte progenitors (MEPs). We performed hyper-variable gene selection and trajectory reconstruction using STREAM³². The resulting trajectory has the shape of a tree with two bifurcations visible in Fig. 2B, having as root the endpoint that concentrates the hematopoietic stem cells.

To transform obtained trajectories into properties over Boolean states, we considered *six states* that must correspond to the start and end of branches (named S0–S5 in Fig. 2B). Moreover, in order to reduce the sensitivity bias of single-cell observations, we chose to consider observations formed by the union of several cells. This resulted in six clusters of a few tens to hundreds of cells, corresponding to initiation (root), two bifurcation points, and three leaves, which we considered to be steady states of the Boolean model. We classified the activity of each gene of each cluster using PROFILE³³ on individual cells and an aggregation by majority of value among 0, 1, and ND (not determined).

Then, we specified the expected dynamical properties of a Boolean network corresponding to the data as follows. There must exist trajectories linking with states following the STREAM trajectories: e.g., there must exist a trajectory from a Boolean state corresponding to the root S1 to the Boolean state matching with S0, then from that later Boolean state to a Boolean state matching with S2, and so forth. Moreover, we requested that the Boolean states corresponding to the leaves (S2, S4, S6) must be a steady state of the Boolean model, and that any steady state reachable from S3 must match with S5 and S4, and any steady state reachable from the root state must match with S2, S4, and S5. Fig. 2C summarizes the declared properties.

Besides the expected dynamical properties of the model, BoNesis requires an a priori GRN from which it will reconstruct Boolean rules that

(1) employ only genes and regulations referenced in the input GRN, and (2) form a Boolean network that possesses the expressed dynamical properties. Moreover, BoNesis is able to identify the genes that can have a constant binary expression towards the whole dynamics, and can thus be ignored in the final model. In the scope of this case study, we took as a priori GRN the full DOROTHEA TF regulation database³⁴ with regulation p to confidence level C, comprising 2777 regulations among 1001 genes, 849 of which have an expression measurement in our dataset. We performed a multi-stage combinatorial optimization procedure to identify the largest number of TF genes that cannot be considered as constant, and for each, there exists a Boolean network having the desired dynamical properties. It resulted in the selection of 39 TF genes and 137 TF–TF interactions shown in Supplementary Fig. 1. A gene set enrichment analysis is performed with METASCAPE³⁵, showing a clear enrichment of terms related to hematopoiesis (Fig. 2D), with the top term being hematopoiesis, followed by others representing more specific biological processes included in hematopoiesis.

Moreover, we compared the data-driven selected genes with three expert art logical models of hematopoiesis by Hamey et al.⁵, Collombet et al.⁴ and Moignard et al.³. These three models, composed of 20–31 genes, include a total of 53 genes, and it is worth noticing that there is no consensus on the interactions implied in the regulation of this process since only 2 genes are common to these three models. The Venn diagram of Fig. 2E shows the intersection between the components we have automatically selected thanks to BoNesis and the three state-of-the-art hematopoiesis models. Each of the state-of-the-art models shares 6 genes with our selection (for a total of 10 distinct genes in common):

* With Hamey et al.⁵: FLI1, GATA1, GFI1B, IKZF1, MYB, RUNX1;

* with Moignard et al.³: FLI1, GATA1, GFI1B, IKZF1, MYB, SPI1;

* with Collombet et al.⁴: CEBPA, EFB1, IKZF1, MEF2C, RUNX1, SPI1.

### Case study 1: model variability analysis

From the TF–TF subnetwork extracted in the previous step (Supplementary Fig. 1), we employed BoNesis to sample 1000 distinct Boolean networks that all respect the qualitative dynamical properties described previously. The sampling has been performed using heuristics to range over models with diverse logical rules. Each of the sampled Boolean networks is able to reproduce the qualitative differentiation dynamics. Moreover, the universal constraints on the reachable steady states, ensured by design, ensure that all trajectories from the root state end in one of the three observed differentiated types. Furthermore, we verified a posteriori that none of the sampled networks possess cyclic attractors.

For each selected gene, we analyzed how many Boolean functions have been assigned to it in the sampled ensemble of 1000 Boolean networks (Supplementary Fig. 2). It resulted that 12 genes always received the same Boolean function, 12 other genes received only 2–3 distinct Boolean functions. This suggests that a large part of the logic rules are highly preserved in all compatible models. Most of the diversity in the sampled ensembles is essentially focused on gene FOS (682 different Boolean functions) and TRP53 (44 Boolean functions).

Clustering the models according to the similarity of their functions can highlight different possible pathways for process regulation and also point out some biologically irrelevant groups. We performed a multidimensional scaling clustering (MDS) of sampled Boolean networks, using distance based on the inequality of Boolean functions: given two Boolean networks *f* and *g* of size *n*, $d(f,g) = \sum_{i=1}^{n} |f_i \neq g_i|$. MDS results highlight three groups of models (Supplementary Fig. 3). It should be kept in mind that the number of models in each group does not reflect their biological relevance, as the group size may result from underlying combinatorial aspects of compatible models, and may also relate to which some close solutions may be easier to find. The variability within each cluster in terms of the number of Boolean functions and influence graph is summarized in Fig. 3. We remark that cluster A (the one in red in Supplementary Fig. 3C) shows a much sparser influence graph, with few regulators per node, and thus has simple Boolean functions with little variability. On the other hand,
<figure><img src="image_6.png" /></figure>

**D** Histogram of the 20 most enriched similar terms sets from the list of the 39 genes automatically selected by BoNesis, colored according to 3 p-value thresholds and with, per set, the representative term of the set of terms.

<figure><img src="image_7.png" /><figcaption>**Fig. 2 |** Overview of the case study on scRNA-seq data-driven modeling of hematopoiesis. A Summary of the pipeline going from scRNA-seq data and DoRothEA interaction database to node selection and Boolean network variability analysis and clustering. B STREAM plot of scRNA-data after pseudo-trajectory reconstruction. The circle illustrate the selected start and end of branches from</figcaption></figure>

**E** Venn diagram showing the number of genes in common between our automatic selection with BoNesis and the state-of-the-art models (16,31,60)

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

which cells where extracted and binarized. C Excerpt of the logical specification of Boolean network dynamical properties, expressed in the BoNesis language. D Result of gene ontology enrichment analysis on nodes selected by BoNesis. E Venn diagram showing the overlap of nodes of 3 Boolean networks of the literature with the nodes selected by BoNesis.

clusters B and C show a rather dense influence graph, witnessing much more complex Boolean functions. Cluster B intra-variability resides mostly in the diversity of functions of the FOS and TRP53 genes, while cluster C also has variability of the Boolean functions of the MYC gene. Detailed variability comparison across the different clusters are provided in Supplementary Figs. 4–6.

## Case study 2: ensemble modeling of adipocyte to osteoblast conversion from bulk RNA-seq time series data

This case study demonstrates a full pipeline going from experimental bulk RNA-seq time series data and background knowledge on TF-TF networks to the prediction of genetic mutations for trans-differentiation and preliminary experimental validation. The predictions have been obtained by
<figure><img src="image_9.png" /><figcaption>**Fig. 3** | Variability of the influence graph with each of the Boolean network clusters. The label of each node includes the number of different local Boolean functions in the sampled Boolean network set. Orange indicates a unique shared function, yellow only 2 different functions. Edge is labeled with the number of</figcaption></figure>

Boolean networks that utilize the influence (over the 1000 sampled), with its thickness scaled accordingly. Graphics A refers to the red cluster, B to the yellow cluster, and C to the green cluster depicted in Supplementary Fig. 3.

<figure><img src="image_10.png" /><figcaption>**Fig. 4** | Illustration of the pipeline employed for case study 2. We combined prior knowledge TF-TF interactions extracted from METACORE(r) database with variants of Boolean dynamical properties extracted from time series bulk RNA-seq data. We employed BoNesis to sample Boolean networks fulfilling these properties. Then, we performed prediction of reprogramming determinants using CABEAN on a</figcaption></figure>

subset of these sampled models. The identified reprogramming determinants are combinations of gene knock-outs and constitute activations. The predictions are then assessed on the other sampled Boolean models using EnsembleMaBoSS simulation, leading to a scoring in terms of fidelity and efficiency, and thus a ranking of the most robust predictions.

combining inference of Boolean network ensembles with formal methods for control and ensemble simulations for scoring (Fig. 4).

Bone marrow stromal progenitor cells (MSCs) are multipotent cells capable of differentiating either into osteoblasts to form bone tissue or into bone marrow adipocytes, which play an important role in the hormonal homeostasis of the bone marrow³⁶. These cells are reciprocal in their differentiation, and a correct balance is important for bone health, with increased adipogenesis observed in obesity and during aging. Previous

studies have sought to understand the gene regulatory networks underlying MSC differentiation, and an improved understanding of these networks could allow for the identification of efficient reprogramming targets and better control of bone marrow cell composition³⁷. We have previously performed transcriptomic profiling of ST2 cells, a mouse MSC cell line, across multiple time points of parallel differentiations into adipocytes and osteoblasts using RNA-seq³⁸. Both differentiations were performed for 15 days, with RNA samples collected from undifferentiated cells (ST2D0) and
at 5 different time points (D1, D3, D5, D9, and D15) of each differentiation. Osteoblastogenesis (OD1-OD15) was performed with one composition of differentiation medium for the entire duration of 15 days, while for adipogenesis (AD1-AD15), two different media were used, with media composition changed on the third day of differentiation (AD3). In vitro differentiations are known to vary between experiments³⁹ and therefore three independent differentiation experiments were performed to capture the robustly reproducible gene expression changes. Please see “Methods” for further details. The obtained dataset formed a rich resource of dynamic gene expression profiles towards two different trajectories from the same starting point.

We employed two different binarization methods by classifying gene activity with respect to background RNA-seq data on a range of tissues, either by bootstrapping parametric distributions (RefBool), or by a simpler statistical procedure applying gene-specific cutoffs, referred to as MUQ in the following. Each gene of each time point is thus assigned to a Boolean or undefined value for each binarization method and each replicate. In total, 1560 genes received a binary value in at least one binarization method and one replicate. For a fixed binarization method, we observed opposite classifications between replicates for up to 29 genes, while up to only 2 genes when comparing across binarization methods. This can be explained by the fact that, in general, the RefBool classifies fewer genes than MUQ, and genes classified by RefBool are classified by MUQ similarly. In order to account for this variance in binarizations, we considered each of the 6 profiles (2 binarization methods times 3 replicates) as alternate model specifications. Our rationale was to constitute ensembles of models for each of these profiles and study the robustness of reprogramming prediction across them.

Our main modeling hypothesis was that the Boolean model must be able to reproduce the observed maturation trajectories from fixed cellular environments. For instance, due to the change of treatment between days 1 and 3 of adipocyte culture, we did not impose the existence of a trajectory from AD1 to AD3. From the experimental protocol, this resulted in the specification of two trajectories: the maturation of adipocytes, modeled as the existence of a Boolean trajectory going through AD3 → AD5 → AD9 → AD15; and the maturation of osteoblasts, modeled as the existence of a Boolean trajectory going through OD1 → OD3 → OD5 → OD9 → OD15.

Moreover, we assumed that ST2D0, AD15, and OD15 are observations of cells in steady states. At last, we modeled the observed cellular differentiation process by denying the existence of trajectories across the two branches, nor reverting to the precursor state: there must not exist a trajectory from AD3 to OD15, from AD3 to ST2D0, from OD1 to AD15, nor from OD1 to ST2D0.

We performed selection of genes and sampling of diverse ensembles of Boolean networks from the qualitative modeling of RNA-seq data. We followed a workflow similar to use case 1, that we repeated for each replicate and for each binarization method. Then, as the gene selection resulted in several optimal solutions, the sampling of Boolean networks has been repeated on each of them.

The prior GRN consisting of TF-TF interactions extracted from the METACORE™ database from all the TFs of the RNA-seq dataset. It resulted in a signed digraph comprising 1027 genes with 11,159 regulations, 169 of which were with an undetermined sign. This prior GRN served to define the set of candidate Boolean networks. Because several genes have more than one hundred referenced regulators, we restricted ourselves to the Boolean networks whose activation functions can be expressed with at most 32 disjunctive clauses, without any limit on the size of the clauses. Thus, while we permit a gene to depend on all of its potential regulators, we limited the number of activation contexts.

On the 1027 genes of the prior GRN, one can expect that only a fraction of them are involved in the observed differentiation process. As in case study 1, we employed BONESIS to perform gene selection by identifying Boolean networks reproducing the dynamics of as many genes as possible, while assigning as many activation functions as possible to a constant value. Depending on the replicate and binarization method, the optimization results in different optimal sets of genes to preserve, ranging from 49 to 79

genes. Finally, we performed the diverse sampling of 264 distinct Boolean networks for each of the six profiles, resulting in 1584 Boolean networks verifying the qualitative dynamical properties corresponding to at least one replicate and one binarization method, and using one of the optimal sets of genes.

### Case study 2: prediction of reprogramming targets with high fidelity and efficiency

Our objective was to predict combinations of perturbations of gene expression to trigger a transdifferentiation of adipocytes into osteoblasts. In terms of the Boolean network, this corresponded to identifying control strategies to enforce the reachability of the OD15 state from the AD15 state. In order to account for candidate model heterogeneity, our approach was to compute reprogramming targets on individual Boolean networks from a subset of the sampled ensemble and evaluate them on the full ensemble. In the end, we aimed at selecting the perturbations predicted to be most effective on a range of models reconstructed from different binarization and replicates.

We selected the CABEAN tool⁴⁰ to compute combinations of temporary gene knock-out and constitutive activations enforcing the convergence to OD15 state from AD15 state for a given individual Boolean network. For each of the 6 qualitative profiles, we applied CABEAN on 24 of the 264 Boolean networks sampled in the inference part. CABEAN failed on 15 of these 144 Boolean networks due to a memory issue. On the remaining 129 models, CABEAN identified combinations of up to 5 simultaneous perturbations leading to a reprogramming from AD15 to OD15. Because we are interested in perturbations that can be effective on as many models as possible, we kept combinations of perturbations that have been identified in at least 10% of the individual models given to CABEAN. This short list of candidates contained 34 different combinations of 2-4 simultaneous perturbations.

The reprogramming perturbations computed by CABEAN are guaranteed to be effective on their input individual Boolean network. Our objective was to assess the robustness of these (combination of) perturbations, the Boolean network ensembles inferred from different qualitative interpretations of the data. To do so, we extended the MABOSS stochastic Boolean network simulator to sample trajectories from the Boolean network ensembles: each Boolean network of the ensemble is simulated k times from the corresponding state corresponding to AD15 while enforcing the given reprogramming perturbation. For each candidate combination of perturbations, we obtain an estimation of the distribution of the steady states of Boolean networks of the ensembles after enforcing the perturbation from the AD15 state.

We defined scores to evaluate reprogramming candidates from their simulation on Boolean network ensembles for their ability to reprogram to the osteoblast phenotype, as observed at OD15. Inspired by usual cellular reprogramming assessments, we considered two measures: the efficiency relates to the proportion of cells (models) that show all the prior knowledge adipocyte gene markers (ADIPOQ, FABP4, CEBPA, LPL) and none of the osteoblast prior knowledge marker gene (ALPL, HEY1, SP7). Then, the fidelity relates to the similarity of those cells to the full OD15 state.

For a given reprogramming perturbation, we write S the set of states resulting from the ENSEMBLEMABOSS simulation, and for each state s ∈ S, we denote by p_s its estimated steady state probability. Moreover, we define 1_{Ost}(s) as being equal to 1 whenever each osteoblast marker gene is active in s and each adipocyte marker gene is inactive in s, otherwise, 1_{Ost}(s) is equal to 0.

The reprogramming efficiency is computed from the estimated steady state distribution after perturbation of Boolean network ensembles as the fraction of simulations ending in a state having all the osteoblast marker genes active and all the adipocyte marker genes inactive:

$$ efficiency(S) = \sum_{s \in S} p_s \cdot 1_{Ost}(s) \quad (1) $$
<figure><img src="image_11.png" /><figcaption>Fig. 5 | Specific TFs occur more often and with higher efficiency in the predicted reprogramming determinants. Frequency (y-axis) and maximal efficiency (x-axis) of the TFs included in the predicted reprogramming determinants applying permanent stimulation are reported. Names of the most promising TFs are given.</figcaption></figure>

The reprogramming fidelity employs a similarity measure between a full Boolean state $s \in S$ and the binarized state $\beta_{OD15}$ of the corresponding qualitative interpretation. Because the binary state of some genes of $\beta_{OD15}$ can be non-determined, the similarity is defined as the proportion of binarized genes in $\beta_{OD15}$ that have the same value in $s$. In the end, the fidelity is the expected similarity of the Osteoblast steady states of Boolean networks in the ensemble after reprogramming from AD15:

$$fidelity(S) = \left( \sum_{s \in S} p_s \cdot 1_{Ost}(s) \cdot similarity(s, \beta_{OD15}) \right) / efficiency(S) \quad (2)$$

### Case study 2: validation of node selection and reprogramming targets with literature

The TFs included in the predicted reprogramming determinants were ranked by frequency and maximal predicted efficiency (Fig. 5). Analysis of the literature related to the predicted top TFs revealed their existing association with the regulation of osteoblastogenesis, thereby providing indirect validation for our approach. For example, HEY1, a well-known target of Notch signaling⁴¹ and a regulator of osteogenesis⁴² was among the factors with the highest predicted efficiency score. Similarly, CEBPA, SP1, and TRP63 have been associated with adipogenesis, osteogenesis, and bone formation⁴³⁻⁴⁵, respectively. The repression of CEBPA is likely to reverse the adipogenic gene expression program. Likewise, the repression of NR2F2 or GATA3 was predicted to lead to high efficiency of conversion, something that is also supported by literature⁴⁶⁻⁴⁹. One of the top TFs included TCF7L2, a mediator of WNT signaling (Fig. 5). Consistently, WNT signaling is known to positively regulate osteoblastogenesis³⁰,⁵¹, with the effect mediated by TCF7L2⁵². Finally, MYC genes have been previously shown to contribute to the reprogramming of fibroblasts towards osteoblasts⁵³,⁵⁴. However, some of the well-established regulators of osteoblastogenesis, such as SP7 (also known as Osterix) or RUNX2⁵³,⁵⁴, were not included among the predicted reprogramming determinants. This could be related to, at least in the case of RUNX2, its already abundant expression in undifferentiated ST2 cells³⁸.

Based on the predicted combinatorial efficiencies, and the existing literature evidence, three combinations of three determinants each were selected for experimental validation. Reprogramming experiment 1 (RE1) involved upregulation of TRP63 and downregulation of CEBPA and SP1. RE2 involved upregulation of TCF7L2 and TRP63 and downregulation of

SP1. And RE3 involved upregulation of TCF7L2 and downregulation of CEBPA and SP1.

### Case study 2: preliminary experimental validation

In order to experimentally test the different combinations of reprogramming determinants selected for RE1, RE2 and RE3, we took advantage of lentiviral overexpression in combination with siRNA-mediated gene repression. Overexpression constructs were synthesized at Sirion Biotech and were designed to coexpress either green fluorescent protein (GFP) and red fluorescent protein (RFP) to confirm a successful lentiviral transduction and to allow separation between different constructs (specifically, co-transduction of TCF7L2 and TRP63). Day 9 adipocytes were transduced with the lentivirus of interest, and 24 h later, were transfected with the relevant siRNAs to achieve the respective combinations for each RE. In parallel, control cells were transduced with a negative control virus over-expressing only GFP and subsequently transfected with a scrambled control siRNA (please see Methods for additional details). Using reverse transcription real-time quantitative PCR (RT-qPCR) TCF7L2 and TRP63 were confirmed to be 8- and 15-fold overexpressed in undifferentiated cells, respectively, while siRNAs against SP1 and CEBPA led to 35–68% decrease in their mRNA expression (Supplementary Fig. 7). Following transfections, the adipocytes were cultured in osteoblast medium for 9 days before collection of the cells for single cell RNA-sequencing (scRNA-seq). After the sequencing and data processing, transduced cells were identified based on the presence of the mRNA sequence encoding for GFP and/or RFP and the relative changes in the transcriptomes were analyzed. Differentially expressed genes (DEGs) were obtained in the single-cell expression data based on the Wilcoxon rank-sum test with FDR < 0.05. Directionality of the observed change was determined with the common language effect size, with a >0.5 indicating up-regulation (see Supplementary Data 1). The precision for the upregulated single-cell DEGs compared against the bulk DEGs (Day 15: Adipocyte vs Osteoblast, adjusted p-value < 0.05) is 0.57, 0.51, and 0.33 for the three different reprogramming experiments RE1–RE3. This means that for RE1 57% of the genes which were called upregulated in the single-cell experiment were also upregulated in osteoblasts compared to adipocytes in the bulk data³⁸, indicating an initial consistent effect of the applied reprogramming determinant. The respective hypergeometric p-values for enrichment of the true positives in the single-cell DEGs are 4.4e-14, 1.1e-12, and 0.057, respectively. This promising tendency is
also confirmed by visual inspection of the bulk RNA-seq data of DEGs obtained in the single-cell experiment (for RE1 see Supplementary Figs. 8 and 9). The very low recall (<0.02) in all experiments, however, indicates that only a few of the relevant genes responded to this perturbation in the relatively short time frame of the experiment.

# Discussion

## Overview of the methodology

We demonstrated our methodology on two complementary case studies for a combination of data-driven and expert-driven inference of Boolean models from scRNA-seq and bulk RNA-seq data. The essence of our approach is the explicit modeling of the inference problem, with the specification of the prior knowledge, typically extracted from databases of TF-TF interactions, and the specification of the expected properties of the model. These latter must reflect the expert and qualitative interpretation of the experimental data. We relied on the BONESIS engine, which enables specifying and combining a broad range of dynamical properties, notably related to trajectories and steady states. Then, ensembles of compatible models can be sampled and further analyzed. We took advantage of ensemble modeling to analyze candidate model variability with model clustering (case study 1), as well as to provide robust reprogramming predictions (case study 2).

The flexibility of the workflow facilitates the comparison of different modeling choices and hypotheses. For instance, in case study 2, we considered different binarization methods and different replicates. We took advantage of the variability to generate an ensemble of models spanning the different hypotheses and identify consensus predictions. The choice of the input GRN can also impact the results, and it could be assessed similarly. Such an analysis is out of the scope of this paper, as our objective was not in the benchmark of prior knowledge data nor binarization methods.

Remarkably, the inference pipeline was tractable on TF-scale networks. By employing logic programming and relying on recent complexity breakthroughs in Boolean network analysis¹⁶, we have been able to fully account for thousands of transcription factors. Then, we took advantage of combinatorial optimization technologies to automatically prune non-necessary variables and identify subnetworks that drive the observed dynamics.

## Case study 1

With this case study, we leveraged ensemble modeling to analyze the diversity of models that can explain the observed differentiation process. The experimental scRNA-seq data has been processed with trajectory inference methods from which we extracted both clusters of cells corresponding to initiation, bifurcation, and differentiated states, and dynamical properties related to existence of trajectories linking those states.

From a full TF-TF interaction database, we have been able to identify a subset of core TFs from which Boolean networks can be drawn reproducing the qualitative differentiation, which present some intersection with expert models from the literature. By analyzing the ensemble of sampled models, we identified three families of models that diverge by the complexity of their logical rules. The variability analysis also emphasized model patterns that are preserved across the ensemble. Such pinpointed model features could be challenged with expert knowledge, and further experimental studies are needed in order to discriminate among subfamilies of models.

As we performed a partial enumeration, we are not assured of not missing models with characteristics different from those of the three highlighted groups. However, partial enumerations of 250 and 500 models already highlight the three groups and suggest that increasing the number of models only leads to an increase of the size of the 3 groups, without highlighting any new type of models. This motivates the chosen number of 1000 models. MDS done with 250 and 500 models are presented in Supplementary Fig. 3A and Supplementary Fig. 3B.

This modeling of hematopoiesis from a causal network automatically inferred from the dynamics of the data already provides avenues for further exploration of the mechanisms of its regulation, with the possibility of ensemble simulation of models with MaBoSS¹⁵.

## Case study 2

The applied pipeline yielded promising predictions for the reprogramming of adipocytes to osteoblasts, as confirmed by our literature validation, highlighting its utility in identifying novel targets from time-series RNA-seq data. Notably, transcription factors such as SP1, HEY1, and TRP63 frequently appeared in both temporary and permanent perturbation lists, suggesting central roles for them in adipocyte-to-osteoblast conversion. The successful prediction included factors like CEBPA and SP1, known to be crucial for adipocyte and osteoblast functions as modulators of the activity of the respective master regulators of these lineages, PPARG and SP7 (also known as Osterix)⁴,⁵. These TFs' selection as reprogramming targets is consistent with their involvement in these regulatory networks, supported by their high occurrence in prediction lists and established biological roles in these processes⁴,⁵.

The reliability of the predictive model is evidenced by its consistency with existing literature on the differentiation of mesenchymal cells. Factors such as ATF4 and NR2F2, previously implicated in positive and negative control of osteoblast differentiation, respectively, are found among our predictions⁴,⁷,⁵. Similarly, the positive effect of HEY1 and TRP63 on bone formation in the literature⁵ highlights the utility of our approach. The predictions also drew attention to the role of less commonly studied factors like AHR and RARA, which are implicated in sensing the cellular microenvironment crucial for differentiation. AHR, predicted for upregulation in two combinations, is known to shift the gene regulatory network toward a less specialized cell state, potentially facilitating the reprogramming process⁵. Similarly, RARA, noted for its negative regulatory role in adipogenesis, is part of the reprogramming predictions⁶, suggesting that its modulation could inhibit adipocyte traits while promoting osteoblast characteristics.

Follow-up studies of our findings could benefit the understanding of bone-related diseases like osteoporosis. By utilizing knowledge of specific gene targets that promote osteoblast differentiation, new therapeutic avenues may be explored that enhance bone regeneration and repair⁶,¹². Additionally, manipulating adipocyte-to-osteoblast conversion could benefit obesity-related diseases where abnormal adipogenesis is prevalent⁶.

Nevertheless, for now, our validation experiments face limitations. While gene expression changes upon reprogramming perturbations were enriched for genes involved in osteoblastogenesis, as observed by scRNA-seq, the overall efficiency of gene delivery and knockdown would need to be improved to be able to better address the quality of our predictions. Usage of stable and inducible systems that would be independent of chemicals such as doxycycline, which might perturb cellular systems beyond intended effects, should optimally be used. Additionally, the complexity of gene and protein interactions, including the necessity of considering interaction partners like beta-catenin with TCF7L2 or post-translational modifications like phosphorylation of factors like JUN and FOS, suggests a greater number of (also non-genetic) perturbations may be required for effective reprogramming than initially predicted. Integrating multi-omics data with transcriptomic profiles could provide a more comprehensive understanding of the reprogramming processes, allowing for the examination of post-translational modifications and protein activity states crucial in the transcriptional state changes.

## General recommendations for addressing a transcriptome dataset

The presented methodology converts quantitative transcriptomic data into Boolean states and trajectories through statistical analysis of the data. Thus, the chosen binarization and trajectory inference method are sensitive parameters, as a slight change in the qualitative specification may impact the sets of compatible models. The case studies employed different methods, which can be guided by the type of data. As shown in the second case study, one may even consider modeling generated from different binarization methods. In general, the binarization of a gene is performed with respect to a reference dataset. This reference dataset can be limited to the studied dataset of the experiment, in which case the binarization reflects qualitative changes
within the studied dataset. Otherwise, the reference dataset can be con-
stituted of data from multiple tissues, in which case, the binarization reflects
a more absolute qualitative state of the gene with regard to other cellular
types. Another important parameter of binarization methods is their
binarization thresholds. Usually, the methods will classify a gene as either
active, inactive, or undetermined. Having too large a proportion of unde-
termined genes may leave too few constraints for building a relevant model.
On the other hand, too large a proportion of determined binary values may
be too sensitive to the employed method and make the inference problem
non-satisfiable⁶⁴. Currently, BoNesis only identifies models that fully satisfy
the given binary observations. Future work may consider relaxing the
compatibility criteria, allowing for and minimizing mismatches between
inferred Boolean networks and input qualitative data.

The Boolean networks generated by BoNesis draw their regulation from the given prior knowledge GRN, typically extracted from a database. Thus, the quality of the results may be impacted when important regulations are missing: either there can be no Boolean network reproducing the dynamics of certain genes, or there can be more complex solutions that bypass the missing regulation through other indirect regulations. Currently, it is impossible to infer missing regulations from BoNesis. BoNesis offers the possibility to identify genes whose behavior is impossible to model with the given PKN, and ignore these genes in order to access the solutions for the other part of the network. Another direction could be to complete the prior knowledge GRN with the GRN inference method applied to the input dataset.

As underlined in the two case studies, and because the experimental data typically covers a tiny fraction of the cellular dynamics, the solution space of Boolean networks reproducing the observations is typically intractable. To address this issue, we leverage the sampling capability of BoNesis, which is based on dynamically tuning solving heuristics of the underlying logical solver, to constitute ensembles of diverse Boolean networks. Nevertheless, we emphasize that there is no theoretical guarantee that the sampled ensemble recapitulates the whole diversity of compatible models. To partly circumvent this issue, future work may consider extending BoNesis with the possibility of prioritizing certain families of Boolean functions, such as those being canalizing or nested canalizing⁶⁵.

Accounting for the previously described limitations, the application of
the presented methodology to another transcriptomic dataset can follow a
similar workflow as case study 1 for the case of single-cell data, or case study
for the case of bulk data. Having chosen a qualitative interpretation of the
data and a source for the prior knowledge GRN, a first control point is the
gene enrichment and functional analysis of selected components (Auto-
mated component selection with BoNesis section), from which the
ensembles of Boolean networks will be generated. They constitute
the minimal and sufficient genes identified by BoNesis for reproducing the
specified cellular behavior. Another important control point is the com-
pletion of gene activities that were left undetermined in the binarization
process of selected cellular states. Indeed, whenever a gene is not categorized
as either active or inactive in a Boolean state, we consider that it can be
assigned to any value. Such an analysis can result in further constraints for
BoNesis, such as imposing value for marker genes missed by the binariza-
tion method, or imposing that certain genes should be at an identical
(but undetermined) value in several states, as used in case study 2.

Finally, regarding computational complexity and cost, the heaviest step
is the generation of ensembles using BoNesis. This is due to the complexity
of the tackled Boolean synthesis problem, which can be NP-complete or
even Σ₂⁻¹-complete decision problems⁶⁶ depending on the requested dyna-
mical properties⁶⁷,⁶⁸. Moreover, the employed logical solver CLINGO⁶⁹ has
few capabilities to exploit parallel computing, and thus has no particular
advantage in being run on a computing cluster. The runtime is also very
dependent on the specification, and can range from a few seconds to several
days. The component selection part of the inference process for case studies
1 and 2 was the most computationally extensive, requiring 3 days for case
study 1, and around 10 h for case study 2 on a 2 GHz CPU with
256 Gb of RAM.

Methods

Knowledge-based (bottom-up) and data-based (top-down)
Boolean network inference methods

The manual approach is a bottom-up modeling (also called forward modeling) that designs a model through expert and literature-based knowledge to determine the Boolean functions from known and suspected interactions<sup>1,2,4,6-12,70</sup>. The resulting model is validated or refined in an iterative process according to the fitting between its dynamic features and observations of the biological phenomenon (for example, its attractors correspond to known phenotypes). This time-consuming approach requires a deep understanding of the biological system and does not ensure that all possible regulations leading to the observed behavior are explored. In contrast, methods have been developed to propose a top-down (reverse) modeling approach that derives from experimental data both the topology of the network, namely the causalities between the biological components, and the logical rules of node activation, constituting the Boolean network<sup>3,5,71-77</sup>. This purely data-based inference approach can suffer from overfitting to a dataset and circular reasoning: a dataset generates a network, which then predicts the dataset. It also confronts the vastness of the solution space without considering prior knowledge about the phenomena. Consequently, other methods have been developed midway between these previously mentioned top-down and bottom-up modeling approaches<sup>21-28,67</sup>, which encompasses the method presented here. The purpose is to leverage both experimental observations and prior knowledge on interactions related to the biological phenomenon to model. Each tool combining data- and knowledge-based modeling considers as input both experimental data, typically expression profiles, and a gene regulatory network (GRN) potentially extended with some features. The GRN, also named influence graph or regulation graph, structures the prior knowledge about interacting components according to the following definition: a directed graph that represents regulatory interactions between biological components, interactions categorized in the simplest case as activations and inhibitions<sup>78</sup>. A high number of tools have been developed to address these issues because they perform the Boolean network inference through different algorithms and are designed for different types of experimental data and different biological interpretations of it. Exploring the space of possible Boolean networks and their dynamics is a challenge because there can be millions of combinations of logical rules that can be formulated. This is why, depending on the algorithm implemented by the tool, the output models can be either Boolean networks that exactly reproduce the desired properties or Boolean networks optimized to best match these properties with no guarantee of optimality. In the first case, the exploration of solutions can be exhaustive. It should also be taken into account that this search for models (exact or by optimization) can be limited to a subspace of solutions if Boolean function patterns are imposed. If we focus on the way solutions are explored, we can distinguish between optimization-based methods<sup>21,25,28,73,74,76,79</sup> (heuristic, identify a Boolean network that best aligns with the observations and prior knowledge, based on a defined set of criteria or objective functions) and decision-based methods on satisfiability problems<sup>3,22,27,28,67</sup> (exact, determine whether there exists a Boolean network consistent with the observations and prior knowledge—based on constraints on its topology and dynamics—and find such network if it exists). The methodology we present in this paper, based on previously published works<sup>67,68</sup>, belongs to this latter category. According to the chosen approach, some methods output (near-)optimal solutions with no guarantee of reaching a global optimum, while others can theoretically explore the solution space exhaustively to output the whole set of Boolean networks that comply exactly with the biological constraints (in practice often limited to a subset of these Boolean networks as the network size increases and the number of solutions becomes tremendous). This latter ability of inferring an ensemble of models is highly informative. It enables variability analysis across the models to identify common patterns of interactions, supporting the formulation of hypotheses that benefit both the modeling process and the biological exploration of key components.

Beyond this distinction between categories of algorithms, limitations
on the formulation of the Boolean rules that are considered by the methods
also impact which solutions are explored. It can be a fixed limit on the number of regulators (e.g, a maximum of six regulators including four activators and two inhibitors⁵,⁷⁷) or a certain pattern of rules (e.g, a component is expressed if all its activators are expressed and none of its inhibitors, i.e., and logic gate between activators and or between inhibitors⁷¹,⁷⁷). Because of the combinatorial explosion, methods designed to consider all possible Boolean functions can also offer a limit on the size of the Boolean function to enable dealing with dozens of components and complex regulations. This limit can be on the number of regulators or, like the choice we made, on the number of different combinations of regulators (i.e., a limit on the number of clauses that can compose the Boolean functions). This latter choice seems to be biologically consistent: the number of direct regulators for a gene or protein can vary significantly, often being numerous. Whereas the number of functional states, reflected by specific combinations of these regulators, is typically more constrained, ensuring system resilience and controlled responses to a manageable range of conditions.

## Qualitative interpretations of experimental data

Regardless the implemented algorithm, the Boolean network inference tools can also be judiciously compared by focusing on the types of data that they can consider for the modeling: bulk or single-cell expression profiles, at cellular steady states or covering cellular evolutions, linear evolutions, or with bifurcations (e.g., lineage differentiation), observed with time series limited to two time points or longer, with/out perturbations on components, etc. This is directly linked to the dynamical features they are able to model. Measurements acquired at assumed stable cellular states are commonly interpreted as attractors in the Boolean network dynamics. On this basis, some inference tools such as CellNOptR²¹, Griffin²³, and BONITA-RD²⁵ are specifically designed for inferring Boolean networks that respect constraints on their attractors (focusing mostly on fixed points), constraints that are derived from a set of steady-state data.

A range of tools extends the modeling scope by addressing time series of expression measurements, representing a succession of cellular states ensuing one after the other. Among them is CellNOptR-dt²¹, which can fit time course data using a synchronous updating scheme. The modeling of this data is widened to asynchronous semantics with ASKEed²⁷ and IQCELL²⁷ that model time series as trajectories in which the succession of measurements constitutes a succession of transitions (reachability in a single dynamic step). Caspo-TS²², also asynchronous, relaxes the constraint of observing each transition as it considers the reachability property between the observed cellular states. It is worth noting that none of these methods takes into account the modeling of bifurcations/branching in the trajectories. Among the methods specially designed for inferring Boolean networks from single-cell expression data, we can distinguish different interpretations of their dynamical features. SCNS³, BTR²⁹ and SgpNet⁷⁶ interpret the set of single-cell measurements as a state space reached from an initial state in an asynchronous Boolean network. These tools seek, through optimization, a Boolean network for which the set of states contained in its transition graph (the “model state space”) is close to the set of binarized expression data (the “data state space”). However, the strategy presented in ref. 75 is based on another interpretation: the state of each single-cell measurement is assumed to be a potential predecessor or successor of the state of any other single-cell measurement coming from the same sample. Consequently, a solution is here a synchronous Boolean network whose dynamics matches a random set of two-time-point time series, namely, couples of single-cell measurements. The previously mentioned IQCELL²⁷ proposes a third interpretation of this data, considering a pseudo-time ordering of the cells seen as a time series and searching for an asynchronous Boolean network whose dynamics includes a sequence of states compatible with the cell ordering.

A tool that can consider a diversity of dynamical features can exploit different types of data and address the modeling of varied and complex cellular behaviors. Three methods stand out for the complexity of the behaviors that can be modeled. The tool BRE:IN²⁴ allows for more nuanced and complex dynamic behaviors than the previously mentioned tools, as it

supports complex temporal logic specifications (CTL and LTL), both in synchronous and asynchronous semantics. However, the high computational cost of checking these traces prevents scaling to networks of biological interactions of dozens of components: the larger networks for benchmarks²⁴ had 16 nodes with synchronous semantics and 11 with asynchronous semantics. Similarly, Boolean network sketches⁸⁰ also performs the inference with the help of model-checking methods, in asynchronous semantics, but with the richer logic HCTL (hybrid extension of the branching-time temporal logic CTL), allowing more expressive specifications. Unlike BRE:IN, the exploration of the Boolean functions is not limited by patterns, and this method was tested in a network with more than 300 nodes.

## Binarization of transcriptome data in case study 1

From the full scRNA-seq dataset, we employed PROFILE³³ to classify in each cell each gene to either 0, 1, or undetermined, by comparing its expression to the distribution in all cells. Then, each cluster identified using STREAM (see main text) has been binarized using a majority rule: a gene is classified following the majority of values the gene has been classified to in the individual cells of the cluster. Among the 4768 hyper-variable genes selected by STREAM, 1519 have been classified with different binary values in at least two different clusters, and 1369 are classified as binary values in all clusters.

## Binarization of transcriptome data in case study 2

We previously performed experiments for acquiring RNA sequencing data in populations of ST2 cells at different stages of differentiation towards adipocytes and osteoblasts¹⁸. Adipocyte differentiation was induced using a medium of isobutylmethylxanthine, dexamethasone, and insulin for 2 days, followed by a change of medium with rosiglitazone insulin (Sigma-Aldrich, I9278) until 9 days. Differentiation towards osteoblasts was induced using bone morphogenetic protein-4 for 15 days. Each experiment has been replicated 3 times. At days 0, 1, 3, 5, 9, 15, a subpopulation of the cells has been sequenced genome-wide. In the scope of this project, we focused on the activity of transcription factors (TFs) in the different stages of ST2 differentiation. To enable building of qualitative models for the differentiation process, an automated method was applied for transforming the quantitative RNA-seq measurements of TF activities into a qualitative assessment of their activity: active (1), inactive (0), or undetermined (intermediate). Two different methods of binarization were employed and compared: RefBool⁸¹ for binarization with respect to background RNA-seq data collected in a range of different cell types; and a statistical analysis developed in the scope of the project, for binarization with respect to the collected RNA-seq data. For this statistical analysis, representative background distributions of gene expressions in 63 mouse tissues were generated from ArchS⁴⁵² data [https://maayanlab.cloud/archs4/, Kallisto raw read counts, retrieved 2019]. Independent vertex sampling was performed per tissue to remove correlated samples. Samples were further filtered for overall read counts between 10 and 100 Mio and a median ≥ 1 to remove unusual distributions and outliers. This background data was merged with our own data (after Kallisto alignment) and then quantile normalized and converted to TPM (gene length normalization and TPM scaling). The gene-specific background distributions were then applied for the discretization of the own data as follows: Gene expressed in own data below the median of the background was discretized to 0, and above the upper quartile to 1. In addition, genes with TPM < 1 are discretized to 0. Furthermore, genes with large expression differences over all samples and time-points were identified via k-means clustering (2 clusters), with a minimum threefold-change of centroid locations, at least three data points per cluster, and fulfilling a *t*-test²² (MATLAB®) between the two clusters. Thirty-six genes were found accordingly. The upper cluster was discretized to 1 and the lower cluster to 0.

## Prior knowledge of gene regulation

For case study 1, we extracted TF–TF and TF–gene interactions referenced in the DoRothEA database⁸³ with confidence levels A–C. It resulted in a signed directed graph with 5186 components and 12,895 regulations. We
filtered the components to keep only genes with classified binary values, reducing to 1001 components and 2777 regulations.

For case study 2, we extracted from the METACORE© database (in 2019) all known interactions (Transcription regulation, Influence on expression, and Binding) between TFs⁸⁴ and seven known marker genes which show clear expression changes from adipocytes to osteoblasts in our own expression data: Four high in adipocytes (ADIPOQ, FABP4, CEBPA, LPL) and three high in osteoblasts (ALPL, HEY1, SP7). Only measured nodes and their interactions were kept. This resulted in a network of 1027 nodes (almost only TFs) with 11,000 pairwise interactions.

## Inference of Boolean networks from structural and dynamical properties

The strategy behind BoNesis is to describe, in the form of a logical problem to be solved, the search for a Boolean network compatible with a network architecture and with given dynamic properties. To do so, BoNesis integrates prior knowledge influence graph and observations within the same logic program, so that any solution of this program is a Boolean network made up of possible interactions given the PKN whose dynamic properties are compatible with the behavior of the observations.

The logic program is written in Answer-Set Programming. It describes the biological data (the prior influence graph as well as the observations and the dynamical properties linking them) and the modeling formalism (the Boolean network and the computation of its dynamics in Most Permissive Semantics) via predicates and constraints. The constraints are necessary and sufficient conditions that guarantee that any solution of the problem is a Boolean network compatible with the biological data, i.e., a Boolean network included in the architecture of the prior influence graph and whose dynamic properties are compatible with the behavior of the observations. Solutions are obtained using an answer-set solver, *clasp*⁸⁵: the models are the answer-sets satisfying the logic program.

## Automated component selection with BoNesis

Building a prior knowledge influence graph specific to a biological process is a complex and delicate expert task. Yet, determining which interactions are to be considered when building a model is an essential step in the preparation of data before synthesizing models. It is important not to miss components that are essential to the regulation mechanism in order to be able to explain the observations, but also not to consider components that do not play any role and hence penalize the construction and understanding of the system. If the interaction graph contains components that are not involved in the regulation of the observed biological behavior, many different functions can be attributed to these components without impacting the compatibility of the Boolean network with the observations. These components, without importance for the behavior, then strongly increase the number of solutions without bringing any information. Conversely, if no combination of functions can reproduce the data because key components and interactions of the process are missing, no model can be found. The complexity of this task currently limits the use of modeling.

With BoNesis, we propose a way to assist in the design of a relevant interaction graph with regard to observations. It offers to select, among a large interaction graph (as it can be extracted from a public interaction database such as DoRothEA³⁴ or Signor⁸⁶), the components that can be included in a model to explain the observations.

To this end, we set two optimization criteria. Firstly, we want the solver to search for a compatible Boolean network that maximizes the number of components in the models. An optimal solution is then the largest Boolean network, composed of components coming from the large original interaction graph, that can reproduce the observations. Yet, we also want the solution to maximize the number of components of a particular type, called *strong constants*. A strong constant is a component to which a constant function is assigned, and whose value remains constant to reproduce the observations within the dynamics of the Boolean network. Thus, within a Boolean network compatible with biological data, a node *A* is a strong constant if and only if $f(A) = v$ with $v \in \{-1, 1\}$ and that within the dynamics

of the Boolean network it is possible to reproduce the data with the node *A* always equal to *v* (all the configurations *x* associated with the observations have $x_A = v$). In other words, it is a component having neither activator nor inhibitor and whose value can remain unchanged without preventing the Boolean network from being compatible with the observations. The particularity of these strong constants is that they can be removed from the domain without impacting compatibility with data behavior. Once the domain of interactions has been reduced to components that can explain the dynamics of the observations and that are not strong constants, we can limit ourselves to the maximum strongly connected component of this graph, which is particularly interesting to focus on the interactions that regulate the observed process.

## Automated model synthesis with BoNesis

BoNesis searches for all the Boolean networks of the same size as the input interaction network and whose dynamics exactly reproduce those of the data (via the defined dynamical constraints). Hence, when BoNesis is used for model synthesis without optimization criteria, all the Boolean networks output from BoNesis are models of the same relevance with respect to the dynamics to be modeled. For biological applications, it is frequent that an exhaustive enumeration of models is not judicious, given that biological observations are rarely constraining enough considering a large-sized interaction network. Indeed, it is sufficient that a few components have little dynamic information to have an explosion in the number of models.

For case study 1, we first performed component selection (see previous section) from the 1001 genes prior knowledge influence graph and the dynamical constraints of Fig. 2C. We extracted the largest strongly connected component of the resulting graph, leading to selecting 39 genes and 137 regulations. Then, from this 39-gene influence graph, we performed a diverse sampling of 1000 Boolean networks fulfilling the dynamical constraints. For case study 2, as described in the main text, the dynamical properties consisted of (1) the existence of trajectories, (2) the stability properties, and (3) the absence of trajectories. Moreover, we add prior knowledge markers of the two differentiated cellular types: ADIPOQ, FABP4, CEBPA, LPL for adipocytes, and ALPL, HEY1, SP7 for osteoblasts. In a first stage, we performed gene selection on the METACORE(r) prior knowledge influence graph by identifying Boolean networks that maximize, by decreasing priority, the inclusion of a prior knowledge markers, the number of genes whose dynamics can be explained, and the number of strong constants, that will be removed (see previous section). For this first stage, we ignored constraints on the absence of a trajectory for complexity reasons. We performed this gene selection for each replicate and each binarization method. We then analyzed the binary values inferred by BoNesis for the selected genes that have not been classified by binarization methods in the stable phenotypes. We identified 19 genes (ATF2, CLOCK, CTCF, CUX1, GATA6, NFYB, REL, RELA, SMAD1, SMAD2, SMAD4, STAT3, TGF1, TRP53, TTF1, USF1, USF2, YY1, ZFP148) that have been inferred to have different binary states, while the data show no clear statistical ground for supporting these different qualitative states. We again performed the gene selection stage with the additional constraint that the binary value of these 19 genes must be equal in all phenotypes. Depending on the binarization method and replicates, different sets of genes have been selected, with, in some cases, multiple optimal solutions. Then, for each optimal solution, we extracted the sub-GRN consisting of the identified non-constant genes, and performed diverse sampling of Boolean networks taking into account the whole dynamical constraints, including the absence of trajectories.

## Boolean network clustering and complexity analysis (Case study 1)

To explore the three highlighted groups, we study whether the complexity of the functions of a model is characteristic of the group to which the model belongs (Supplementary Figs. 4 and 5). Specifically, Supplementary Fig. 4 shows, for each group, the percentage of the functions of its models according to the number of clauses constituting the functions.
Supplementary Fig. 5 shows, for each group, the percentage of clauses constituting its model functions according to the number of components in the clauses. We observe that red models are significantly less complex than those of the other two groups, with functions composed of a smaller number of clauses, themselves composed of a smaller number of components, compared to the orange and green group models.

Among the functions constituting the models, some are invariant between models of the same group, beyond the 12 invariable functions which are common to all groups. We therefore have three sets of invariable functions that seem to indicate three key patterns of interactions to reproduce the data dynamics. Supplementary Fig. 6 shows, for each group, the percentage of its invariable functions according to their number of clauses. Again, a clear difference in complexity between the red group and the others is highlighted. All the functions shared by the red models have a single clause, whereas they represent for green and orange models, respectively, less than half and one third of the functions (some of their invariable functions having even, respectively, up to 7 and 8 clauses). Hence, the red group models appear to be more parsimonious explanations of the regulation mechanism of hematopoiesis than the other two groups.

## Computational prediction of reprogramming determinants (case study 2)

The prediction of the combinations of perturbations for case study 2 has been performed via the software CABEAN⁴⁰. CABEAN implements several methods for the source-target control of asynchronous Boolean networks, and these methods can be used to identify the minimal and exact control sets of perturbations that ensure the inevitable reachability of the target attractor from a given source attractor. Based on the application time of the perturbations, CABEAN supports several types of controls: instantaneous control applies perturbations instantaneously; temporary control applies perturbations for a sufficient time and then releases them to retrieve the original dynamics; and permanent control applies the control for all the following time steps. All these control methods⁸⁷ are based on the computation of strong and weak basins of attractors, which explore both the structure and dynamical properties of asynchronous Boolean networks. More specifically, an instantaneous control drives the network dynamics from the source attractor to a state in the strong basin of the target attractor, from which there only exist paths to the target attractor. On the other hand, temporary control and permanent control can make use of the spontaneous evolutions of the network dynamics by moving into the weak basin of the target attractor, from which there exist paths to the target attractor. To ensure the inevitable reachability of the target attractor, a temporary control should drive the network dynamics to a state in the strong basin of the target attractor at the end of control, while a permanent control stirs the network from the source state to a state in the strong basin of the target attractor in the resulting transition system under control. More recently, CABEAN⁸⁸ has been extended with target control methods for asynchronous Boolean networks⁸⁹, which can be used to identify perturbations that can drive the dynamics of a Boolean network from any initial state to the desired target attractor.

In this way, CABEAN fits well with the objective for case study 2, i.e., enforcing the convergence to OD15 state from AD15 state. CABEAN was applied to each inferred Boolean network from BoNesis to compute combinations of temporary perturbations (gene knock-out), and it identified combinations of up to 5 simultaneous perturbations leading to a reprogramming from AD15 to OD15. We kept combinations of perturbations that have been identified in at least 10% of the individual models given to CABEAN. In the end, we have obtained 34 different combinations of 2 to 4 simultaneous perturbations, which are given as input for the next analysis steps.

## Stochastic simulations of ensembles of Boolean networks (case study 2)

To perform simulations on our ensembles of models, we relied on the software Ensemble-MaBoSS⁶⁸, which extends the stochastic Boolean

simulator MaBoSS¹⁵ to ensembles. MaBoSS usually takes a single Boolean model and performs numerous simulations to obtain many stochastic trajectories. It later uses this set of trajectories to compute time-dependent probabilities for every visited Boolean state. During these simulations, it also stores every fixed point observed and produces the set of fixed points observed in all these simulations. To perform the simulation of an ensemble of models, we need to compute an equal set of trajectories for every model of the ensemble. We can then use the consolidated set of trajectories from all the models and compute the time-dependent probabilities specific to the whole ensemble of models. We can also produce the same list. To analyze the composition of the ensemble, we also compute the time-dependent probabilities for each model. This way, we get individual results that can then be used for ensemble analysis, such as model clustering.

## Case study 2: cell culture

The mouse bone marrow-derived stroma cell line ST2 was established from Whitlock-Witte type long-term bone marrow culture of BC8 mice³⁰. The ST2 cells were cultured in growth medium: Roswell Park Memorial Institute (RPMI) 1640 medium (Gibco, 32404-014) supplemented with 10% fetal bovine serum (FBS) (Gibco, 10270-106) and 1% L-Glutamine (Lonza, BE17-605E) at 37 °C, 5% CO₂. All experiments were carried out with cells less than 10 passages. For differentiation of adipocytes, ST2 cells were seeded 2 days before differentiation. The cells reached 100% confluency after 24 h of culture, and were further maintained in growth medium for 24 h (Day 0). Adipogenic differentiation was initiated on Day 0 by culturing in adipogenic differentiation medium I consisting of growth medium, 0.5 mM isobutylmethylxanthine (IBMX) (Sigma-Aldrich, I5879), 0.25 µM dexamethasone (DEXA) (Sigma-Aldrich, D4902) and 5 µg/mL insulin (Sigma-Aldrich, I9278). From Day 2, the cells were cultured in adipogenic differentiation medium II consisting of growth medium, 500 nM rosiglitazone (RGZ) (Sigma-Aldrich, R2408) and 5 µg/mL insulin (Sigma-Aldrich, I9278). The adipogenic differentiation medium II was replaced every second day.

For differentiation of osteoblasts, ST2 cells were seeded 2 days before differentiation, reached 100% confluency after 24 h of culture, and were further maintained in growth medium for 24 h (Day 0). Osteogenic differentiation was initiated on Day 0 by culturing in osteogenic differentiation medium consisting of growth medium and 100 ng/mL bone morphogenetic protein-4 (BMP4) (PeproTech, 315-27). The osteogenic differentiation medium was replaced every second day.

## Case study 2: RNA extraction and cDNA synthesis

Total RNA was extracted from cells using Quick-RNATM MiniPrep (Zymo Research, R1055). RNA concentration was measured by Nanodrop 2000c (Thermo Fisher Scientific, E597). cDNA was synthesized with the following cocktail: 1 µg total RNA, 0.5 mM dNTPs (Thermo Fisher Scientific, R0181), 2.5 µM oligo dT-primer (Eurofins GmbH), 1 U/µL Ribolock RNase inhibitor (Thermo Fisher Scientific, E00381), and 1 U/µL RevertAid Reverse transcriptase (Thermo Fisher Scientific, EP0352). The cocktail was maintained at 42 °C for 1 h, followed by 70 °C for 10 min to stop the reaction.

## Case study 2: reverse transcription real-time quantitative PCR (RT-qPCR)

RT-qPCR was performed to measure the RNA expression using the Applied Biosystems 7500 Fast Real-Time PCR System. Each reaction contained the following cocktail: 5 µL of cDNA, 5 µL of primer mix (forward and reverse primers, both in 2 µM concentration), and 10 µL of Absolute Blue qPCR SYBR Green Low ROX Mix (Thermo Fisher Scientific, AB4322B). The PCR reaction was the following: 95 °C for 15 min and repeating 40 cycles of 95 °C for 15 s, 55 °C for 15 s, followed by 72 °C for 30 s. The gene expression level was calculated using the 2^(-ΔΔCt) method⁹¹. The ΔΔCt refers to (ΔCt(target gene) - ΔCt(housekeeping gene)) from the treatment - (ΔCt(target gene) - ΔCt(housekeeping gene)) from the control. Rpl13a was used as the housekeeping gene, and the primer sequences can be found in the Supplementary Data 2.
## Case study 2: viral transduction

Day 9 adipocytes were transduced with lentivirus (Sirion Biotech GmbH; details can be found in the Supplementary Data 2). For a 48-well plate, the estimated cell number was 200,000. The amount of lentivirus was calculated to achieve the expected multiplicity of infection (MOI). The transduction cocktail was as follows: lentivirus stock diluted with RPKM and 8 µg/mL polybrene transfection reagent (EMD Millipore, TR-1003-G) to achieve a final volume of 100 µL. The culture medium was removed, and the cells were washed once with 1xDPBS (Gibco, 14190-094), then supplemented with the transduction cocktail. The cells were kept in the transduction cocktail for 6 h. After the 6-h transduction, the transduction cocktail was removed, and the cells were supplemented with the growth medium. The transduction efficiency was controlled by observing GFP and RFP levels by microscopy.

## Case study 2: RNA interference

On Day 1 post-transduction (Day 1 PT), the cells were transfected with siRNA according to the manufacturer's recommendation (Horizon Discovery; details can be found in the Supplementary Data 2). In brief, siRNA was diluted to 5 µM solution in DNase/RNase free water (Invitrogen, 10977-035). In separate tubes, siRNA and DharmaFECT 1 transfection reagent (Horizon Discovery, T-2001-03) were diluted with RPMI. To prepare the transfection cocktail for 1 well of a 48-well plate, in Tube 1, 25 µL of the siRNA in serum-free medium was prepared by adding 1.25 µL of 5 µM siRNA to 23.75 µL of RPMI. In Tube 2, 1 µL of the DharmaFECT 1 in serum-free medium was prepared by adding 1 µL of DharmaFECT 1-24 µL of RPMI. The mixture was incubated for 5 min, then the two tubes were mixed following incubation for 20 min. After incubation, 200 µL of growth medium was added to achieve a final volume of 250 µL. The culture medium was removed, and the cells were supplemented with the transfection medium for 24 h. For single and double siRNA transfection, the final concentration was maintained at 100 nM.

## Case study 2: single-cell RNA-seq

The single-cell RNA-seq (scRNA-seq) was performed according to Chrominum Next GEM Single Cell 3' Reagent Kits v3.1 Review D (CG000204). The ChromiumTM Next GEM Single Cell 3' GEM, Library Gel Bead Kit 3.1 (1000128) consisted the following: Chrominum Next GEM Single Cell 3' Gel Bead Kit v3.1, 4 (1000129), Chromium Next GEM Single Cell 3' GEM Kit v3.1 (1000130), Chromium Next GEM Single Cell 3' Library Kit v3.1 (1000158), Single Index Kit T SetA (1000213), DynabeadsTM MyOneTM SILANE (2000048). To achieve a single cell suspension, the cells were first treated with 1 mg/mL Collagenase A (Roche, 10103586001) for 15 min. After Collagenase A treatment, the cells were trypsinized with 150 µL Trypsin (Lonza, BE17-161E) for 5 min. The trypsin was quenched with 500 µL growth medium. Wide-bore tips (Thermo Scientific ART 1000G Self-Sealing Barrier Pipet Tips, 2079G) were used to pipette up and down. The suspension was centrifuged at 1000 rpm for 5 min, and the supernatant was removed. Single cell suspension was achieved by adding 400 µL growth medium and pipetting up and down with wide-bore tips. The cell number was determined with C-CHIP (NanoEntek, DHC-N01). The suspension was centrifuged again at 1000 rpm for 5 min, and the supernatant was removed and replaced with 1xDPBS + 0.04% BSA to achieve the target cell number. The suspension was filtered through a 40 µM Flowmi Cell strainer (Merck, BAH136800040-50EA). The cell number was determined again with C-CHIP, which would be used to determine the input for scRNA-seq. In brief, cells were loaded with a targeted recovery rate of 10,000 cells per sample. scRNA-seq library quality was assessed by Agilent DNA High sensitivity Bioanalyzer chip (Agilent, 5067-4626) and further sequenced on a 150 cycles High Output Kit using Illumina NextSeqTM 500 with targeted sequencing depth of 20,000 read pairs per cell.

## Data availability

The software BoNesis, employed for Boolean network inference, is available at https://bnediction.github.io/bonesis under the GPLv3-compatible free software license CeCiLL. The software Cabean, employed for reprogramming

prediction, is available at https://satoss.uni.lu/software/CABEAN. Both tools are available through the CoLoMoTo Docker distribution ⁹² at https://colomoto.github.io/colomoto-docker/. Case study 1: data and scripts are available at https://doi.org/10.5281/zenodo.16990754. Case study 2: The scRNA-seq data have been deposited in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB82479 (https://www.ebi.ac.uk/ena/browser/view/PRJEB82479). Matlab code for statistical analysis and binarization of bulk RNA-seq data, as well as data and scripts for Boolean network inference and control predictions are available at https://doi.org/10.5281/zenodo.16990770.

Received: 21 October 2024; Accepted: 17 July 2025;

Published online: 26 September 2025

## References

1. Fauré, A., Naldi, A., Chaouiya, C. & Thieffry, D. Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle. *Bioinformatics* **22**, e124-e131 (2006).

2. Martínez-Sosa, P. & Mendoza, L. The regulatory network that controls the differentiation of T lymphocytes. *Biosystems* **113**, 96-103 (2013).

3. Moignard, V. et al. Decoding the regulatory network of early blood development from single-cell gene expression measurements. *Nat. Biotechnol.* **33**, 269-276 (2015).

4. Collombet, S. et al. Logical modeling of lymphoid and myeloid cell specification and transdifferentiation. *Proc. Natl Acad. Sci. USA* **114**, 5792-5799 (2017).

5. Hamey, F. K. et al. Reconstructing blood stem cell regulatory network models from single-cell molecular profiles. *Proc. Natl Acad. Sci. USA* **114**, 5822-5829 (2017).

6. Ikonomi, N., Kühlwein, S. D., Schwab, J. D. & Kestler, H. A. Awakening the HSC: dynamic modeling of HSC maintenance unravels regulation of the TP53 pathway and quiescence. *Front. Physiol.* **11**, 848 (2020).

7. Meyer, P. et al. A model of the onset of the senescence associated secretory phenotype after DNA damage induced senescence. *PLOS Comput. Biol.* **13**, e1005741 (2017).

8. Siegle, L. et al. A Boolean network of the crosstalk between IGF and Wnt signaling in aging satellite cells. *PLoS ONE* **13**, e0195126 (2018).

9. Mai, Z. & Liu, H. Boolean network-based analysis of the apoptosis network: Irreversible apoptosis and stable surviving. *J. Theor. Biol.* **259**, 760-769 (2009).

10. Zhang, R. et al. Network model of survival signaling in large granular lymphocyte leukemia. *Proc. Natl Acad. Sci. USA* **105**, 16308-16313 (2008).

11. Cohen, D. P. A. et al. Mathematical modelling of molecular pathways enabling tumour cell invasion and migration. *PLOS Comput. Biol.* **11**, e1004571 (2015).

12. Werle, S. D. et al. Unraveling the molecular tumor-promoting regulation of cofilin-1 in pancreatic cancer. *Cancers* **13**, 725 (2021).

13. Montagud, A. et al. Patient-specific boolean models of signalling networks guide personalised treatments. *eLife* **11**, e72626 (2022).

14. Shmulevich, I., Dougherty, E. R., Kim, S. & Zhang, W. Probabilistic Boolean networks: a rule-based uncertainty model for gene regulatory networks. *Bioinformatics* **18**, 261-274 (2002).

15. Stoll, G. et al. MaBoSS 2.0: an environment for stochastic Boolean modeling. *Bioinformatics* **33**, 2226-2228 (2017).

16. Paulevé, L., Kolčák, J., Chatain, T. & Haar, S. Reconciling qualitative, abstract, and scalable modeling of biological networks. *Nat. Commun.* **11**, 4256 (2020).

17. Krumsiek, J., Pölsterl, S., Wittmann, D. M. & Theis, F. J. Odefy—from discrete to continuous models. *BMC Bioinforma.* **11**, 233 (2010).

18. Chudasama, V. L., Ovacik, M. A., Abernethy, D. R. & Mager, D. E. Logic-based and cellular pharmacodynamic modeling of bortezomib responses in U266 human myeloma cells. *J. Pharmacol. Exp. Therap.* **354**, 448-458 (2015).
19. Bloomingdale, P., Nguyen, V. A., Niu, J. & Mager, D. E. Boolean network modeling in systems pharmacology. *J. Pharmacokinet.* **45**, 159–180 (2018).

20. Putnins, M. & Androulakis, I. P. Boolean modeling in quantitative systems pharmacology: challenges and opportunities. *Crit. Rev. Biomed. Eng.* **47**, 473–488 (2019).

21. Terfve, C. et al. CellNOptR: a flexible toolkit to train protein signaling networks to data using multiple logic formalisms. *BMC Syst. Biol.* **6**, 133 (2012).

22. Ostrowski, M., Paulevé, L., Schaub, T., Siegel, A. & Guziolowski, C. Boolean network identification from perturbation time series data combining dynamics abstraction and logic programming. *Biosystems* **149**, 139–153 (2016).

23. Muñoz, S., Carrillo, M., Azpeitia, E. & Rosenblueth, D. A. Griffin: a tool for symbolic inference of synchronous boolean molecular networks. *Front. Genet.* **9**, 39 (2018).

24. Goldfeder, J. & Kugler, H. Bre:in—a backend for reasoning about interaction networks with temporal logic. In (eds L. Bortolussi and G. Sanguinetti) *Computational Methods in Systems Biology*, pp 289–295 (Springer International Publishing, Cham, 2019).

25. Palli, R., Palshikar, M. G. & Thakar, J. Executable pathway analysis using ensemble discrete-state modeling for large-scale data. *PLoS Comput. Biol.* **15**, e1007317 (2019).

26. Gjerga, E. et al. Converting networks to predictive logic models from perturbation signalling data with CellNOpt. *Bioinformatics* **36**, 4523–4524 (2020).

27. Vaginay, A., Boukhobza, T. & Smail-Tabbone, M. Automatic synthesis of Boolean networks from biological knowledge and data. In (eds B. Dorronsoro, L. Amodeo, M. Pavone & P. Ruiz) *Optimization and Learning*, pp 156–170 (Springer International Publishing, Cham, 2021).

28. Aghamiri, S. S. & Delaplace, F. Taboo Boolean network synthesis based on tabu search. *IEEE/ACM Trans. Comput. Biol. Bioinf.* https://doi.org/10.1109/TCBB.2021.3063817 (2021).

29. Chevalier, S., Boyenval, D., Magaña-López, G., Roncalli, T., Vaginay, A. & Paulevé, L. *BoNesis: a Python-Based Declarative Environment for the Verification, Reprogramming, and Synthesis of Most Permissive Boolean Networks*, pp 71–79. (Springer Nature Switzerland, 2024) https://doi.org/10.1007/978-3-031-71671-3_6.

30. Nestorowa, S. et al. A single-cell resolution map of mouse hematopoietic stem and progenitor cell differentiation. *Blood* **128**, e20–e31 (2016).

31. Hérault, L., Poplineau, M., Duprez, E. & Remy, É. A novel Boolean network inference strategy to model early hematopoiesis aging. *Comput. Struct. Biotechnol. J.* **21**, 21–33 (2023).

32. Chen, H. et al. Single-cell trajectories reconstruction, exploration and mapping of omics data with stream. *Nat. Commun.* **10**, 1903 (2019).

33. Béal, J., Montagud, A., Traynard, P., Barillot, E. & Calzone, L. Personalization of logical models with multi-omics data allows clinical stratification of patients. *Front. Physiol.* **9**, 1965 (2019).

34. Holland, C. H., Szalai, B. & Saez-Rodriguez, J. Transfer of regulatory knowledge from human to mouse for functional genomics analysis. *Biochim. Biophys. Acta* **1863**, 194431 (2020). Transcriptional Profiles and Regulatory Gene Networks.

35. Zhou, Y. et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. *Nat. Commun.* **10**, 1523 (2019).

36. Ambrosi, T. H. et al. Adipocyte accumulation in the bone marrow during obesity and aging impairs stem cell-based hematopoietic and bone regeneration. *Cell Stem Cell* **20**, 771–784.e6 (2017).

37. Meyer, M. B., Benkusky, N. A., Sen, B., Rubin, J. & Pike, J. W. Epigenetic plasticity drives adipogenic and osteogenic differentiation of marrow-derived mesenchymal stem cells. *J. Biol. Chem.* **291**, 17829–17847 (2016).

38. Gérard, D. et al. Temporal enhancer profiling of parallel lineages identifies AHR and GLIS1 as regulators of mesenchymal multipotency. *Nucleic Acids Res.* **47**, 1141–1163 (2019).

39. Grancharova, T. A comprehensive analysis of gene expression changes in a high replicate and open-source dataset of differentiating hiPSC-derived cardiomyocytes. *Sci. Rep.* **11**, 15845 (2021).

40. Su, C. & Pang, J. CABEAN: a software for the control of asynchronous Boolean networks. *Bioinformatics* **37**, 879–881 (2020).

41. Tu, X. et al. Physiological notch signaling maintains bone homeostasis via RBPjk and hey upstream of NFATc1. *PLoS Genet.* **8**, e1002577 (2012).

42. Sharff, K. A. et al. Hey1 basic helix-loop-helix protein plays an important role in mediating BMP9-induced osteogenic differentiation of mesenchymal progenitor cells. *J. Biol. Chem.* **284**, 649–659 (2009).

43. Lin, F. T. & Lane, M. D. CCAAT/enhancer binding protein alpha is sufficient to initiate the 3T3-L1 adipocyte differentiation program. *Proc. Natl Acad. Sci. USA* **91**, 8757–8761 (1994).

44. Niger, C. et al. The transcriptional activity of osterix requires the recruitment of Sp1 to the osteocalcin proximal promoter. *Bone* **49**, 683–692 (2011).

45. Gu, J. et al. Mouse p63 variants and chondrogenesis. *Int. J. Clin. Exp. Pathol.* **6**, 2872–2879 (2013).

46. Xie, X., Qin, J., Lin, S.-H., Tsai, S. Y. & Tsai, M.-J. Nuclear receptor chicken ovalbumin upstream promoter-transcription factor II (COUP-TFII) modulates mesenchymal cell commitment and differentiation. *Proc. Natl Acad. Sci. USA* **108**, 14843–14848 (2011).

47. Lee, K.-N. et al. Orphan nuclear receptor chicken ovalbumin upstream promoter-transcription factor II (COUP-TFII) protein negatively regulates bone morphogenetic protein 2-induced osteoblast differentiation through suppressing runt-related gene 2 (runx2) activity. *J. Biol. Chem.* **287**, 18888–18899 (2012).

48. Lee, M. N. et al. FGF2 stimulates COUP-TFII expression via the MEK1/2 pathway to inhibit osteoblast differentiation in C3H10T1/2 cells. *PLoS ONE* **11**, e0159234 (2016).

49. Wang, L. & Di, L.-J. Wnt/β-catenin mediates AICAR effect to increase GATA3 expression and inhibit adipogenesis. *J. Biol. Chem.* **290**, 19458–19468 (2015).

50. Day, T. F., Guo, X., Garrett-Beal, L. & Yang, Y. Wnt/beta-catenin signaling in mesenchymal progenitors controls osteoblast and chondrocyte differentiation during vertebrate skeletogenesis. *Dev. Cell* **8**, 739–750 (2005).

51. Baron, R. & Kneissel, M. WNT signaling in bone homeostasis and disease: from human mutations to treatments. *Nat. Med.* **19**, 179–192 (2013).

52. Yin, C. et al. Transcription factor 7-like 2 promotes osteogenic differentiation and boron-induced bone repair via lipocalin 2. *Mater. Sci. Eng. C Mater. Biol. Appl.* **110**, 110671 (2020).

53. Yamamoto, K. et al. Direct conversion of human fibroblasts into functional osteoblasts by defined factors. *Proc. Natl Acad. Sci. USA* **112**, 6152–6157 (2015).

54. Mizoshiri, N. et al. Transduction of Oct6 or Oct9 gene concomitant with myc family gene induced osteoblast-like phenotypic conversion in normal human fibroblasts. *Biochem. Biophys. Res. Commun.* **467**, 1110–1116 (2015).

55. Madsen, M. S., Siersbæk, R., Boergesen, M., Nielsen, R. & Mandrup, S. Peroxisome proliferator-activated receptor γ and c/ebpa synergistically activate key metabolic adipocyte genes by assisted loading. *Mol. Cell. Biol.* **34**, 939–954 (2014).

56. Bernot, D. et al. Down-regulation of tissue inhibitor of metalloproteinase-3 (TIMP-3) expression is necessary for adipocyte differentiation. *J. Biol. Chem.* **285**, 6508–6514 (2010).

57. Makowski, A. J. et al. The loss of activating transcription factor 4 (atf4) reduces bone toughness and fracture toughness. *Bone* **62**, 1–9 (2014).
58. Gu, J. et al. Mouse p63 variants and chondrogenesis. *Int. J. Clin. Exp. Pathol.* **6**, 2872 (2013).

59. Gérard, D. et al. Temporal enhancer profiling of parallel lineages identifies AHR and GLIS1 as regulators of mesenchymal multipotency. *Nucleic Acids Res.* **47**, 1141–1163 (2018).

60. Muenzner, M. et al. Retinol-binding protein 4 and its membrane receptor Stra6 control adipogenesis by regulating cellular retinoid homeostasis and retinoic acid receptor α activity. *Mol. Cell. Biol.* **33**, 4068–4082 (2013).

61. Oh, W.-T. et al. Wnt-modulating gene silencers as a gene therapy for osteoporosis, bone fracture, and critical-sized bone defects. *Mol. Ther.* **31**, 435–453 (2023).

62. Clabaut, A. et al. Adipocyte-induced transdifferentiation of osteoblasts and its potential role in age-related bone loss. *PLoS ONE* **16**, e0245014 (2021).

63. Chen, Q. Fate decision of mesenchymal stem cells: adipocytes or osteoblasts? *Cell Death Differ.* **23**, 1128–1139 (2016).

64. Magaña López, G., Calzone, L., Zinovyev, A. & Paulevé, L. scBoolSeq: linking scRNA-Seq Statistics and Boolean Dynamics. *PLOS Comput. Biol.* **20**, 1–25 (2024).

65. Kadelka, C. et al. A meta-analysis of Boolean network models reveals design principles of gene regulatory networks. *Sci. Adv.* https://doi.org/10.1126/sciadv.adj0822 (2024).

66. Papadimitriou, C. H. *Computational Complexity*. (Addison-Wesley, 1995).

67. Chevalier, S., Froidevaux, C., Paulevé, L. & Zinovyev, A. Synthesis of Boolean Networks From Biological Dynamical Constraints Using Answer-set Programming. In *2019 IEEE 31st International Conference on Tools with Artificial Intelligence (ICTAI)*, pp 34–41 https://doi.org/10.1109/ICTAI.2019.00014 (2019).

68. Chevalier, S., Noël, V., Calzone, L., Zinovyev, A. & Paulevé, L. Synthesis and simulation of ensembles of boolean networks for cell fate decision. In (eds A. Abate, T. Petrov & V. Wolf) *Computational Methods in Systems Biology*, pp 193–209 (Springer International Publishing, Cham, 2020).

69. Gebser, M., Kaminski, R., Kaufmann, B. & Schaub, T. Clingo = ASP + control: preliminary report. *CoRR*, abs/1405.3694 (2014).

70. Thakar, J., Pilione, M., Kirimanjeswara, G., Harvill, E. T. & Albert, R. Modeling systems-level regulation of host immune responses. *PLOS Comput. Biol.* **3**, e109 (2007).

71. Dunn, S.-J., Martello, G., Yordanov, B., Emmott, S. & Smith, A. G. Defining an essential transcription factor program for naïve pluripotency. *Science* **344**, 1156–1160 (2014).

72. Udyavar, A. R. Novel hybrid phenotype revealed in small cell lung cancer by a transcription factor network model that can explain tumor heterogeneity. *Cancer Res.* **77**, 1063–1074 (2017).

73. Liu, X., Wang, Y., Shi, N., Ji, Z. & He, S. Gapeore: Boolean network inference using a genetic algorithm with novel polynomial representation and encoding scheme. *Knowl. Based Syst.* **228**, 107277 (2021).

74. Hung-Cuong, T. & Yung-Keun, K. Cga-bni: a novel constrained genetic algorithm-based boolean network inference method from steady-state gene expression data. *Bioinformatics* **37**, i383–i391 (2021).

75. Schwab, J. D. et al. Reconstructing Boolean network ensembles from single-cell data for unraveling dynamics in the aging of human hematopoietic stem cells. *Comput. Struct. Biotechnol. J.* **19**, 5321–5332 (2021).

76. Gao, S., Sun, C., Xiang, C., Qin, K. & Lee, T. H. Learning asynchronous boolean networks from single-cell data using multiobjective cooperative genetic programming. *IEEE Trans. Cybern.* **52**, 2916–2930 (2022).

77. Heydari, T. et al. Iqcell: A platform for predicting the effect of gene perturbations on developmental trajectories using single-cell RNA-seq data. *PLOS Comput. Biol.* **18**, e1009907 (2022).

78. Marku, M. & Pancaldi, V. From time-series transcriptomics to gene regulatory networks: a review on inference methods. *PLOS Comput. Biol.* **19**, e1011254 (2023).

79. Lim, C. Y. et al. BTR: training asynchronous Boolean models using single-cell expression data. *BMC Bioinforma.* **17**, 355 (2016).

80. Beneš, N., Brim, L., Huvar, O., Pastva, S. & Šafránek, D. Boolean network sketches: a unifying framework for logical model inference. *Bioinformatics* **39**, btad158 (2023).

81. Jung, S., Hartmann, A. & Del Sol, A. RefBool: a reference-based algorithm for discretizing gene expression data. *Bioinformatics* **33**, 1953–1962 (2017).

82. Lachmann, A. et al. Massive mining of publicly available RNA-seq data from human and mouse. *Nat. Commun.* **9**, 1366 (2018).

83. Garcia-Alonso, I., Holland, C. H., Ibrahim, M. M., Turei, D. & Saez-Rodriguez, J. Benchmark and integration of resources for the estimation of human transcription factor activities. *Genome Res.* **29**, 1363–1375 (2019).

84. Heinäniemi, M. et al. Gene-pair expression signatures reveal lineage control. *Nat. Methods* **10**, 577–583 (2013).

85. Gebser, M., Kaufmann, B. & Schaub, T. Conflict-driven answer set solving: from theory to practice. *Artif. Intell.* **187**, 52–89 (2012).

86. Licata, L. et al. SIGNOR 2.0, the SIGNaling Network Open Resource 2.0: 2019 update. *Nucleic Acids Res.* **48**, D504–D510 (2019).

87. Su, C., Paul, S. & Pang, J. Controlling large Boolean networks with temporary and permanent perturbations. In *Proceedings of the 23rd International Symposium on Formal Methods*, volume 11800 of *Lecture Notes in Computer Science*, pp 707–724. (Springer-Verlag, 2019).

88. Su, C. & Pang, J. CABEAN 2.0: Efficient and efficacious control of asynchronous Boolean networks. In *Proceedings of the 24th International Symposium on Formal Methods*, volume 13047 of *Lecture Notes in Computer Science*, pp 581–598. (Springer-Verlag, 2021).

89. Su, C. & Pang, J. Target control of asynchronous Boolean networks. *IEEE/ACM Trans. Comput. Biol. Bioinforma.* **20**, 707–719 (2023).

90. Ogawa, M. et al. B cell ontogeny in murine embryo studied by a culture system with the monolayer of a stromal cell clone, st2: B cell progenitor develops first in the embryonal body rather than in the yolk sac. *EMBO J.* **7**, 1337–1343 (1988).

91. Livak, K. J. & Schmittgen, T. D. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-delta delta C(T)) method. *Methods* **25**, 402–408 (2001).

92. Naldi, A. et al. The CoLoMoTo interactive notebook: accessible and reproducible computational analyses for qualitative biological networks. *Front. Physiol.* https://doi.org/10.3389/fphys.2018.00680 (2018).

## Acknowledgements

This research was primarily supported by the bilateral French Agence Nationale pour la Recherche (ANR) and Luxembourg Fonds pour la Recherche Nationale (FNR) project AlgoReCell (ANR-16-CE12-0034; FNR INTER/ANR/15/11191283). S.C. also acknowledges support from ITMO Cancer. A.Z. and L.C. also acknowledge support from ANR as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). L.P. acknowledges support from ANR in the scope of the project BNeDiction (grant number: ANR-20-CE45-0001) and RD2BOOL (grant number: ANR-23-CE45-0014), and in the scope of the France 2030 project AI4scMED operated by ANR (grant number ANR-22-PESN-0002).

## Author contributions

L.C., A.Z., J.P., L.S., T.S. and L.P. designed the methodological principles; S.C., Y.G., V.N., C.S., J.P., T.S. and L.P. implemented and performed computational analyses; J.B., Y.B., S.J., A.S., L.S. and T.S. provided biological expertise and assessment of the predictions; J.B., Y.G. and L.S. performed the experimental work. S.C., J.P., L.S., T.S., and L.P. wrote the initial version of the paper. All authors reviewed the paper.
**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41540-025-00569-z.

**Correspondence** and requests for materials should be addressed to Loïc Paulevé.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2025