Boolean models are a well-established framework to model developmental gene regulatory networks (DGRN) for acquisition of cellular identity. During the reconstruction of Boolean DGRNs, even if the network structure is given, there is generally a very large number of combinations of Boolean functions (BFs) that will reproduce the different cell fates (biological attractors). Here we leverage the developmental landscape to enable model selection on such ensembles using the relative stability of the attractors. First we show that 5 previously proposed measures of relative stability are strongly correlated and we stress the usefulness of the one that captures best the cell state transitions via the mean first passage time (MFPT) as it also allows the construction of a cellular lineage tree. A property of great computational convenience is the relative insensitivity of the different measures to changes in noise intensities. That allows us to use stochastic approaches to estimate the MFPT and thus to scale up the computations to large networks. Given this methodology, we study the landscape of 3 Boolean models of Arabidopsis thaliana root development and find that the latest one (a 2020 model) does not respect the biologically expected hierarchy of cell states based on their relative stabilities. Therefore we developed an iterative greedy algorithm that searches for models which satisfy the expected hierarchy of cell states. By applying our algorithm to the 2020 model, we find many Boolean models that do satisfy the expected hierarchy. Our methodology thus provides new tools that can enable reconstruction of more realistic and accurate Boolean models of DGRNs.

Boolean models are a well-established framework to model developmental gene regulatory networks (DGRNs) for acquisition of cellular identities. During the reconstruction of Boolean DGRNs, even if the network structure is given, there is generally a large number of combinations of Boolean functions that will reproduce the different cell fates (biological attractors). Here we leverage the developmental landscape to enable model selection on such ensembles using the relative stability of the attractors. First we show that previously proposed measures of relative stability are strongly correlated and we stress the usefulness of the one that captures best the cell state transitions via the mean first passage time (MFPT) as it also allows the construction of a cellular lineage tree. A property of great computational importance is the insensitivity of the different stability measures to changes in noise intensities. That allows us to use stochastic approaches to estimate the MFPT and thereby scale up the computations to large networks. Given this methodology, we revisit different Boolean models of Arabidopsis thaliana root development, showing that a most recent one does not respect the biologically expected hierarchy of cell states based on relative stabilities. We therefore developed an iterative greedy algorithm that searches for models which satisfy the expected hierarchy of cell states and found that its application to the root development model yields many models that meet this expectation. Our methodology thus provides new tools that can enable reconstruction of more realistic and accurate Boolean models of DGRNs.

Boolean network (BN) models of gene regulatory networks (GRNs) have gained widespread traction as they can easily recapitulate cellular phenotypes via their attractor states. The overall dynamics of such models are embodied in the system's state transition graph (STG) which is highly informative. Indeed, even if two BN models have the same network structure and recover the same attractors, their STGs can be drastically different depending on the type of regulatory logic rules or Boolean functions (BFs) employed. A key objective of the present work is to systematically delineate the effects of different classes of regulatory logic rules on the structural features of the STG of reconstructed Boolean GRNs, while keeping BN structure and biological attractors fixed. Furthermore, we ask how such global features might be driven by characteristics of the underlying BFs. For that, we draw from ideas and concepts proposed in cellular automata for both the structural features and their associated proxies. We use the network of 10 reconstructed Boolean GRNs to generate ensembles that differ in the type of logic used while keeping their structure fixed and recovering their biological attractors, and compute quantities associated with the structural features of the STG: 'bushiness' and 'convergence', that are based on the number of garden-of-Eden (GoE) states and transient times to reach attractor states when originating at them. We find that ensembles employing biologically meaningful BFs have higher 'bushiness' and 'convergence' than those employing random ones. Computing these 'global' measures gets expensive with larger network sizes, stressing the need for more feasible proxies. We thus adapt Wuensche's Z-parameter to BFs in BNs and provide 4 natural variants, which along with the network sensitivity, comprise our descriptors of local dynamics. One variant of the network Z-parameter as well as the network sensitivity correlate particularly very well with the bushiness, serving as a good proxy for the same. Finally, we provide an excellent proxy for the 'convergence' based on computing transient lengths originating at random states rather than GoE states.