Tuto2 • Tuto1 • Paper

AI-driven Automated Discovery Tools Reveal Diverse Behavioral Competencies of Biological Networks

Abstract

Many applications in biomedicine and synthetic bioengineering rely on understanding, mapping, predicting, and controlling the complex behavior of chemical and genetic networks. The emerging field of diverse intelligence investigates the problem-solving capacities of unconventional agents. However, few quantitative tools exist for exploring the competencies of non-conventional systems. Here, we view gene regulatory networks (GRNs) as agents navigating a problem space and develop automated tools to map the robust goal states GRNs can reach despite perturbations. Our contributions include: (1) Adapting curiosity-driven exploration algorithms from AI to discover the range of reachable goal states of GRNs, and (2) Proposing empirical tests inspired by behaviorist approaches to assess their navigation competencies. Our data shows that models inferred from biological data can reach a wide spectrum of steady states, exhibiting various competencies in physiological network dynamics without requiring structural changes ofin network properties or connectivity. We also explore the applicability of these ‘"behavioral catalogs’" for comparing evolved competencies across biological networks, for designing drug interventions in biomedical contexts and synthetic gene networks for bioengineering. These tools and the emphasis on behavior-shaping open new paths for efficiently exploring the complex behavior of biological networks. For the interactive version of this paper, please visit https://developmentalsystems.org/curious-exploration-of-grn-competencies.

Introduction

Developing methods to recognize, map, predict, and control the complex, context-sensitive behavior of chemical and genetic networks is an essential frontier of research in science and engineering. These systems, such as gene regulatory networks and protein pathways, are known to be instructive drivers of embryogenesis, cell behavior, and complex physiology [91][77][116]. Understanding the control properties of these systems is critical not only for the study of evolutionary developmental biology [88][85][84][71][122], but also for comprehending and intervening in various disease states, including cancer [95][11][144], and for the construction of novel synthetic biologicals in bioengineering contexts[17][35][86][45][18].

Thus, much work has gone into mathematical modeling and computational inference of both protein pathways and gene regulatory network models [50][98][47][104], which has resulted in the development of large collections of publicly-available models such as the Biomodels database [120][121]. Yet, despite the wealth of available models, scientists still largely lack an effective understanding of the range of possible behaviors that these models can exhibit under different initial conditions and environmental stimuli, and are in search of systematic methods to reveal and optimize those behaviors via external interventions. The full extent of the computational and control properties of such networks are not yet well-understood; while dynamical systems theory has been extensively used to characterize their behavior [9][129], it is not known what other sets of tools might reveal and exploit interesting properties of this ubiquitous biological substrate. The field of diverse intelligence (also known as basal cognition) has suggested that strong functional symmetries between pathway networks and neural networks could imply the existence of learning and other kinds of behavior in this unconventional substrate [124][43][109][2][89]. Specifically, it has been hypothesized that gene regulatory networks (GRNs) and other molecular networks could be endowed with surprising navigation competencies allowing them to robustly reach diverse homeostatic or allostatic states despite a wide range of perturbations [62][134][105][130], and that exploiting these innate competencies could provide a promising roadmap for the design of interventions in regenerative medicine and bioengineering contexts [111][78].

However, significant challenges remain in practice for the exploration and behavior-shaping of these innate competencies, which presents a barrier to the use of these ideas in regenerative medicine and bioengineering. Because of the non-linearity and redundancy in pathway dynamics, passive exploration strategies such as random screening are likely to either fail in uncovering the full range of potential behaviors or require time and energy beyond the available resources. Here, we formalize and investigate a view of gene regulatory networks as agents navigating a problem space. We propose a framework and automated tools, leveraging (1) curiosity-driven goal-directed exploration algorithms coming from recent advances in machine learning and (2) a battery of empirical tests inspired from behaviorist approaches, for mapping the repertoire of robust goal states that GRNs can reach within this problem space despite various perturbations. A key novelty of this work is the use of AI-based exploration tools to map the space of possible behaviors in biological networks, which opens interesting avenues for efficient mapping of unfamiliar system behaviors, yielding transferable insights for diverse problem-solving once such a map is discovered.

The challenge of exploring and mapping spaces of complex and self-organized behaviors appears in many fields such as diverse intelligence in biological systems, minimal active matter or robotics: many systems in these areas provide a rich space of evolved, engineered, and hybrid systems that offer many of the same fundamental problems of behavior and control regardless of specific composition or provenance [87]. These span many orders of spatio-temporal scale, from molecular assemblies to swarms of complex organisms [2][139][94][49]. One set of approaches seeks to develop tools to identify the optimal level of control, ranging from physical rewiring to various methods from cybernetics and behavioral sciences, to reveal and exploit the native competencies and computational capacities of these systems [17]. Specifically, it is increasingly realized that the level of competency (and thus the appropriate level of control) often cannot be guessed by inspection of a system's components, and that its position on a spectrum ranging from passive matter to complex metacognition must be determined empirically [87][126][58][16]. This is critical not only for fundamental understanding of evolution of bodies and minds [43][14][42][52][135][97], but also for the design of interventions in biomedicine and synthetic morphology contexts [32][5]. Yet, a common property in many of these systems is that it is expensive in time and energy to conduct experiments: empirical exploration needs to be made under limited resources. Thus, methods for automating efficient exploration and discovery of a diversity of behaviors in these spaces may be widely useful. As explained below, we will here leverage methods from developmental artificial intelligence initially designed for the specific purpose of exploring a diversity of behaviors using a limited budget of experiments.

One especially fascinating set of systems concerns cellular molecular pathways, or gene regulatory networks (GRNs). In the lab or clinic, these pathways are usually treated as simple machines, with intervention strategies focusing on rewiring their structure to achieve a desired outcome: adding or removing nodes (gene therapy), or changing connection weights (by targeting promoter sequences or protein structures) [31][100][75][68]. However, the emergent, generative nature of development and physiology ensure that it is often very hard to know which genes/proteins to modify, and how, in order to reach a complex desired system-level outcome [142]. Moreover, the responses of cells and tissues to drugs changes over time, making it even more difficult to infer specific interventions (e.g. drugs) that will induce a stable improvement in pathway state in vivo. Indeed, with the exception of antibiotics and surgery, most available treatment modalities do not solve the underlying problem -- they seek to mitigate symptoms, which recur (or expand) once the drug is withdrawn. This is because current therapeutics function bottom-up, attempting to force specific molecular states, as it has been challenging to develop methods for shifting complex tissues and organs towards a stable health profile. Next-generation solutions, which would offer true healing (stable correction), require an understanding of the homeostatic and allostatic properties of networks with respect to how they traverse the space of transcriptional, physiological and anatomical states. An understanding of the behavior policies of networks as they dynamically navigate these problem spaces is essential for predicting what stimuli can be used to re-set their setpoints and guide them to autonomously maintain a healthy state. In the language of behavioral neuroscience, this strategy corresponds to exploiting their native robustness, decision-making, and navigational competencies to induce predictable, long-lasting changes in functionality.

Significant challenges remain in revealing and controlling the range of behaviors that can self-organize in these cellular and molecular pathways . To characterize steady-state concentrations and responses to small perturbations, conventional methods rely on piecewise-linear approximation of the system behavior [123][145][25][40][103], but struggle with higher-dimensional systems or wider parameter ranges which limits their applicability [29]. Other works have proposed the porting of tools from network control theory to identify sets of control nodes allowing to drive the network behavior toward target steady states [60]. These methods typically exploit the network topology [60][106][117][19][101] or regulatory structure [70][8][51] to identify control strategies based either on permanent knockout/activation of genes or on temporary perturbations, the latter being preferable in biomedical context.

However, these approaches often require prior knowledge of target attractor states or are limited to Boolean network models. Other works have explored the use of machine learning tools, such as evolutionary search [69][83][82] and gradient-descent optimization [133][79], for controlling continuous ODE biomolecular networks with high-dimensional parameter spaces, mainly in the context of synthetic circuit engineering [46][119]. While providing powerful optimization tools, these approaches tend to focus on rewiring network structure and connectivity. Moreover, the choice of a predefined fitness function and parameter range initialization is not only critical to the success of optimization [83] but largely restricts exploration of the behavior space [79].

In contrast, an alternative line of research proposes exploring and leveraging the inherent molecular mechanisms of adaptivity and robustness in cellular pathways as a promising approach for drug interventions that do not rely on genomic editing or gene therapy [62][140]. Recently, a broad, substrate-independent behavior science perspective suggests novel properties of gene regulatory networks (GRNs) and other biological networks [12][124]. This perspective views GRNs as agents that convert activation levels of specific genes (inputs) to those of effector genes (outputs), with intermediate nodes in between, leading to strategies for controlling network behavior based on a specific history of inputs (experience) rather than through network rewiring. Notably, the concept of training a chemical pathway using pulsed input stimuli (node activation or suppression drugs) has been formalized, and several networks have been analyzed to establish a taxonomy of memory types found in biological GRNs and pathways [76][63].

Here, building upon recent research [105][76][63], we take the next step and investigate a view of gene regulatory networks as agents navigating a problem space toward target goal states with varying degrees of competency (Figure 1-a). We seek to implement a definition of goal that abstracts it from conventional associations with human or other advanced brains and facilitates the use of tools from cybernetics, behavior science, and control theory to understand broader aspects of biological regulation. Here we use the term “goal” state to refer to a system’s steady state, which it expends effort to reach despite interventions or barriers - a definition appropriate to the study of basal (or minimal) proto-cognitive regulatory systems.. Our definition of goal does not imply “purpose” (high-level goals where an agent has the meta-cognition to think about having goals and what they might be), and we do not attribute high-level competencies (such as re-setting one’s own goals) to GRNs.

Our particular focus lies in investigating two types of navigation competencies: versatility, which refers to the capacity to reach diverse goal states under different interventions, and robustness, which refers to the ability to reach a goal state despite various perturbations. The primary scientific question we aim to address is: What is the repertoire of robust goal states that a GRN can actively reach through minimal and non-genetic interventions within a navigation task context, and can we develop systematic methods and automated tools to aid scientists in discovering this repertoire?

Figure 1: Overview of the proposed framework. (a) MOTIVATION: We often focus on studying the navigation and behavior of organisms in conventional three-dimensional environments, neglecting the intelligence underlying competencies at sub-organismal scales [105]. To better understand navigation competencies in unconventional organisms solving problems in unconventional spaces (e.g., embryos in morphological space), it is essential to construct comprehensive "behavioral catalogs" for these novel entities, which in turn requires sophisticated exploration methods to discover the extent of possible behaviors. Images are taken and adapted with permissions from Figure 2-C in [16], Figure 6-A in [152], graphical abstract of [53], Figure 3-D of [151] and Figure 4-A of [6].(b) EXPERIMENTAL DESIGNS: We formalize GRN behavior as a navigation task and propose to investigate it by defining abstract and observer-dependent "problem spaces" that we use to organize the observed biological behaviors and their exploration in practice. (c) AUTOMATED EXPERIMENTATION: Pseudo-code of the curiosity-driven goal exploration process we use to automate the discovery of behavioral abilities that the GRN can exhibit in behavior space. (d) EMPIRICAL TESTS: We use a battery of empirical tests to identify the robust goal states of the systems, i.e. the one that can be attained under a wide variety of perturbation (including noise in gene expression, and pushes or walls during traversal of transcription space). (e) PERSPECTIVES: We explore several potential reuses of the discovered "behavioral catalog" and proposed framework across evolutionary biology, biomedicine and bioengineering contexts.

To address this question in practice, our experimental framework revolves around the definition of "problem spaces", which we use as tractable components of the GRN's overall state space (Figure1-b), and on a set of methodological contributions which we organize around three sub-questions:

1. Automated discovery of diverse behavioral abilities with autotelic curiosity search (Figure 1-c): What is the range of possible goal states that GRNs can exhibit and how can we devise efficient exploration strategies to automatically identify these goal states? Defining goal states as attractor states of the underlying gene regulatory network, we show that traditional screening methods can be very inefficient in discovering the range of possible goal states. To address this, we propose to use intrinsically-motivated goal exploration processes (IMGEP) [136][65], a recent family of diversity-driven machine learning approaches also known as autotelic curiosity search which was recently shown to form a useful discovery assistant for revealing the behavioral diversity of unfamiliar systems such as chemical oil-droplet systems[146], physical non-equilibrium systems [99] and models of continuous cellular automata [66][72][61].

2. Evaluation of the navigation competencies (Figure 1-d): How competent is the GRN, in terms of robustness to perturbations, in attaining the diverse previously-identified goal states? Prior studies have offered definitions of robustness in biological networks, characterized as the degree of variation in functionality [3] or phenotypic trait [38] under specific environmental or genetic changes. However, these studies often consider a predefined functionality and random perturbations in network parameters [4][28][82] or specific gene knockouts [48]. Environmental perturbations on the other hand are often limited to random variations in initial conditions within a predefined range [29][7]. Here, inspired from behaviorist approaches, we test hypotheses about non-genetic resistance with respect to various navigation competencies that living agents often exhibit, and that do not require structural changes of network properties or connectivity. Those tests assess the system's ability to maintain robustness despite various perturbations encountered during traversal, including developmental noise in gene expression levels, sudden "pushes" within transcriptional space, and the presence of energy barriers or "walls" acting as force fields in the environment.

3. Potential reuses of the discovered "behavioral catalog" and framework (Figure 1-e): Can the constructed behavioral catalogs be useful for fundamental research and practical therapeutic applications, and can the framework be easily applied to other systems and problem spaces? We propose that the discovered competencies may provide valuable insights for understanding evolvability and developmental robustness, and provide a fertile source for the design of interventions in biomedicine and synthetic morphology contexts. We also suggest that the framework and automated tools, which are observer-focused and substrate-independent, could be transposed to other systems and problem spaces.

The overall framework is summarized in Figure 1. Applying it on a database of 30 continuous (ODE) models from the Biomodels website, consisting of a total of 432 systems defined as GRN model-behavior space tuples, revealed several interesting insights. First, results suggested that most of the surveyed systems are capable of reaching a surprisingly wide spectrum of steady states depending on their initial state. Interestingly, random screening strategies were not able to reveal this diversity of reachable states (or at least not in a sample efficient way), confirming the need for more advanced exploration strategies like curiosity search. Secondly, among the discovered steady states, we were able to identify several robust goal states i.e. ones that the system consistently reaches despite various perturbations during traversal of transcriptional space. Altogether, these findings seem to suggest that cell phenotype and functionality could be the result of a multi-step program [106] that could be flexibly and robustly reprogrammed by appropriate stimuli [16]. Finally, we demonstrate possible reuses of this "behavioral catalog" for comparing the network's competencies across different classes of organisms, as well as for the design of non-genetic drug interventions. We also demonstrate an alternative reuse of the framework to reveal new kinds of reachable "goals" in synthetic gene networks, suggesting alternative strategies for the design of gene networks in a bioengineering context.

An interactive executable version of the paper, as well as step-by-step tutorials and notebooks can be found online at [https://developmentalsystems.org/curious-exploration-of-grn-competencies]{.underline}. The full codebase of the proposed automated experimentation pipeline is written end-to-end in JAX, a high-performance numerical computing library that we leverage for parallel experimentation and computational speedups of the ODE models time-course simulations.

Results

Generalizing GRN behavior as a navigation task

Table 1: Glossary of terms used in this paper, with the proposed isomorphism which investigates concepts from dynamical complex systems and behavioral sciences under a common navigation task perspective.

The GRNs analyzed in this study are biological pathway networks taken from the BioModels repository [120][121]. The term "GRN" is used broadly to include protein interaction, gene regulatory, and metabolic networks. In these mathematical models, the dynamic interactions between nodes of the network (molecular species) are modeled with a system of ordinary differential equations, enabling to quantitatively simulate time-course behavior (model rollouts) and observe the dynamics of node activities over time (Figure 2a). Here, following a terminology which aims to integrate concepts from dynamical complex systems with concepts from behavioral sciences, we propose to conceptualize GRN behavior as a navigation task (Table 1). Model rollouts are viewed as "trajectories" in transcriptional space where network steady states are "goal states" (endpoints) that the "agent" (GRN) can reach with varying levels of competencies. As for living agents, these competencies may range from unstable locomotion patterns to more advanced forms of goal-directed behavior like path following, obstacle avoidance, or even forms of spatial memory and foresight. In this paper, we are particularly interested in investigating two forms of navigation competencies that we refer to as versatility, the capacity to reach diverse goal states under various interventions, and robustness, the capacity to reach a goal state despite various perturbations. Note that versatility and robustness are studied with respect to different sources of incoming environmental variation, respectively interventions and perturbations.

Table 2: Problem spaces used in this study

To investigate these competencies in practice, our experimental framework is based on the definition of "problem spaces", which include the observation space (O), behavior space (Z), intervention space (I) and perturbation space (U) as defined in Table 2. To be consistent with our navigation task terminology introduced in Table 1, we refer to a behavior z as the reached "goal state" of a GRN trajectory. However these "goals" may lie on a continuum between complete robustness and high sensitivity, and our primary interest lies in identifying robust goals of the system. Whereas several choices could be made for the intervention space I and perturbation space U, we intentionally consider minimal and non-genetic interventions to investigate the "native" goal states of the GRN, and environmental obstacles to investigate for navigation competencies classically observed in other living agents. Examples of simulations, interventions, and perturbations are illustrated in Figure 2.

Then, a typical analysis using our framework relies on a 2-step procedure, detailed in the subsequent sections. First, to assess the versatility of the GRN, we define an exploration strategy which organizes the sequence of interventions $i_{1},\ldots,i_{N}$ used to drive the system toward a maximally diverse set of reachable endpoints { $z_{k} \in Z$ }$_{k = 1,N}$ , while being given a limited budget of experiments N. Secondly, to assess the robustness of the discovered goal states {${ z}_{k} \in Z$ }, we conduct a battery of empirical tests to characterize their degree of sensitivity to novel perturbations, with a fixed experimental budget of P perturbations per selected behavior z. At the end of this 2-step procedure, we obtain the "behavioral catalog" (H) of the studied GRN, which includes the history of experiments $H = $ {$( i_{k},o_{k},z_{k},~$ { $( u_{p},o_{p},z_{p} ),~p = 1...P$ } $ ),~~k = 1\ldots N$ }.

Following this framework, the behavioral catalog is constructed for a database of 30 biological networks consisting of a total of 432 systems, where a system is defined as a (GRN model, intervention space (I), behavior space (Z)) tuple, as described in Materials and Methods and Table S1. These catalogs provide valuable empirical observations and insights into the navigation competencies of the studied GRNs, particularly in their ability to consistently achieve diverse goal states under various tested perturbations. Statistical analyses of the results are presented in Figures 3, 5, and 7, and specific results for the RKIP-ERK signaling pathway [56] are shown in Figures 2, 4, 6, and 8.

Figure 2: Illustration of the experimental setup and chosen problem spaces on an example GRN model which has 10 nodes and models the influence of RKIP on the ERK Signaling Pathway [56]. (a) Time-course evolution of the different nodes y1, ..., y10 (one color per node) when starting from the default initial conditions (as provided in [56]). The observation captures the states taken through time $o=[y(t=0), ..., y(t=T)]$ where $y=[y_1, ..., y_{10}]$. (b) Corresponding trajectory in transcriptional space (phase space), for two target nodes (ERK, RKIPP_RP), from t=0 (A, in red) to T=1000 seconds (B, in cyan). We can see that the trajectory converges to endpoint B in less than 100 seconds, and then stay there. The behavior (or reached goal state) is the endpoint $B = \left\lbrack y_{ERK}(T),y_{RKIPRP}(T) \right\rbrack$, where T is chosen big enough to ensure convergence. (c) The intervention is setting the initial state of the system trajectory (for all nodes): $i = [y_1(t=0), ..., y_{10}(t=0)]$. (d-e) Example of perturbations used in this paper. (d) Noise perturbation, here applied to all 10 nodes every 5 secs until t=80 secs. (e) Push perturbation, here applied to the two target nodes (ERK, RKIPP_RP) at t=3 seconds. (f) Wall perturbation, also applied to the two target nodes (ERK, RKIPP_RP), here at 10% and 90% of the total distance traveled. Supplementary Figure S1 shows examples of other possible "drug" or "genome" interventions that can be implemented in the accompanying software, as well as the possibility to perform interventions (or perturbations) in parallel using batched computations.

Curiosity Search Uncovers a Diversity of Reachable Goal States

Figure 3: Curiosity search uncovers a wide spectrum of reachable states in behavior space Z. (a) Diversity of endpoints discovered by random search (pink) and curiosity search (blue) for the 432 systems. Diversity is measured as the volume of the union of the set of hyperballs of radius $\epsilon$ that have for centers the discovered endpoints {$z \in Z$} as depicted by the shaded area in (b-c) with $\epsilon = 0.05$. (a-left) Mean and standard deviation curves of the diversity of behaviors discovered throughout exploration (with random search having twice more experiments n=900). Dots indicate significance (p<0.05) when testing curiosity search (n) against random search (n) in brown, and against random search (n=900) in black, with a Welch's t-test. Standard deviation is divided by 4 for visibility. (a-right) Detail of the diversity obtained in the left plot for all 432 systems at n=450 and n=900, where *** indicate significance (p<0.001). (b-c) Discovered endpoints at the end of exploration (n=450) by random search (pink) and curiosity search (blue) for 6 example systems of our database. (b) Examples of systems for which curiosity search is much more sample-efficient than random search in finding a diversity of reachable states in behavior space Z. (c) Examples of systems with low-redundancy mapping $I \rightarrow Z$ such that random search in $I$ is already quite efficient in covering behavior space Z, and curiosity search performs equivalently.

One advantage of modeling GRN behavior within a tractable behavior space Z is that we can then deploy strategies to efficiently discover and map that space. Notably, recent diversity-driven machine learning techniques such as Novelty Search [81][137], Quality Diversity [30][34] and Intrinsically-Motivated Goal Exploration Processes (IMGEP) [136][65] are explicitly designed to efficiently explore a so-called "behavior space" or "goal space" which is basically a (predefined or learned) model of the overall state space. In particular IMGEPs, which were originally developed for the learning of inverse models of highly-redundant mapping in robotics context [136], were recently shown to successfully assist discovery in complex self-organizing systems [146][99][66][72].

Here, we propose to use an IMGEP to control GRN initial states and maximize the diversity of discovered endpoints {$z \in Z$ } within a limited budget of $N$ experiments. The IMGEP operates in two phases: initially, $N_{int}$ interventions are sampled randomly from $I$ to populate history $H$, then remaining interventions are generated through a goal-directed process which relies on several key internal models. Those including a goal-embedding module $(R)$ that encodes observations $(o)$ into the IMGEP goal space $(Ƭ)$, a goal generator module $(G)$ that samples goals from the goal space based on intrinsic motivation incentives (e.g. to promote novelty or learning progress), and a goal-conditioned optimization policy $(\Pi)$ that generates candidate intervention parameters to achieve the current goal. Given those internal models, the goal-directed phase iterates through 1) sample a target goal $g \sim G(H)$, 2) infer intervention parameters to achieve that goal $i \sim \Pi(g,H)$, 3) conduct an experiment with the intervention i, observe the outcome o, and compute the reached goal $z = R(o)$, and 4) store the tuple $(i,o,z)$ in history $H$. Here, we use the GRN behavior space Z as the IMGEP goal space $Ƭ = Z$. Hence "target goal" refers to a goal sampled by IMGEP while "reached goal" refers to an actual endpoint of the GRN trajectory, discovered by IMGEP while targeting a potentially different point in Z. Throughout exploration, the IMGEP dynamically refines its Z-traversal strategy based on the knowledge acquired by its discoveries. Here we opt for a simple IMGEP variant such that the exploration process can be seen as performing novelty search in behavior space Z [44]. The pseudocode of our IMGEP pipeline is shown in Figure 1-c and details about the internal models are provided in Materials and Methods. The final outcome is a "behavioral catalog" of the GRN, containing the diverse goal states discovered by IMGEP: $H =$ { $( i_{k},o_{k},z_{k} ),~~k = 1\ldots N$}.

We deploy the IMGEP, also known as "curiosity search," on all 432 systems in the biological network database. Our evaluation focuses on two related competencies: the IMGEP agent's ability to empirically reveal a diversity of reachable goal states in the (GRN, I, Z) system, referred to as "discovered diversity," and the GRN agent's competency to naturally reach diverse goal states, referred to as "versatility." The true versatility of the GRN is unknown and can only be inferred through empirical exploration and proxy metrics.

For evaluating diversity, we measure the area covered in Z by the IMGEP discoveries using the threshold-coverage metric [26] and compare it with the area covered by the diversity of a naive random screening strategy (which uniformly samples interventions in $I$). In Figure 3, the diversity discovered by the two exploration variants is shown for the 432 $(GRN,I,Z)$ systems, where random search is given a budget of experiments $N$ which is twice bigger (N=900) as the one given to the curiosity-search algorithm (N=450). Interestingly we see that, on average, at n=290 the curiosity search already significantly outperforms the final diversity achieved by random search, while only utilizing one-third of its experimental budget (N=900). Whereas we are dealing with numerical systems and our codebase allow for efficient and parallel execution, each experiment still consists of $\frac{T}{\Delta T} = 25000$ model steps, where each step integrates the ODE system. Repeating that N times for each of the 432 systems starts to be very costly, which is why having efficient exploration strategies is very valuable (and would be even more valuable when scaling the framework to larger databases). Even more critical, as illustrated in Figure 3-b, it seems that, for some systems, random search is not able to discover the "latent" regions revealed by the IMGEP in Z, or it would need an extremely large budget of experiments. On the other hand, as illustrated in Figure 3-c, there are some systems for which random search is already quite efficient in revealing diverse behaviors in Z, and for which IMGEP performs equivalently.

In fact, the goal-directed strategy of the IMGEP is particularly beneficial for $(GRN,I,Z)$ systems with high nonlinearity or redundancy in their $I \rightarrow Z$ mapping, as seen in Figure 4 and studied in robotics contexts [26]. Redundancy implies that many interventions in $I$ lead to similar effects in $Z$, as illustrated in Figure 2 where various interventions and perturbations converge to the same endpoint. In these systems, random search will preferentially discover points in areas of high redundancy in Z whereas the IMGEP, whose exploration is directed uniformly in goal space, should cover different levels of redundancy equally. In general, when dealing with large intervention spaces and limited experimental budgets, curiosity search can be particularly useful for efficiently navigating Z-space.

Figure 4: Illustration of the non linearity and redundancy of the $I \rightarrow Z$ mapping, and of the interest of using goal-directed exploration strategies. Plot shows the reachable points discovered by curiosity search (a) and by random search (b) in the behavior space Z and their corresponding starting points in the intervention space I, for the RKIP-ERK signaling pathway system [56]. The intervention space is 10-dimensional, and here we show the TSNE reduction in 2D. We apply HDBSCAN clustering [73] on the points discovered in Z, which produced the 4 clusters for curiosity search (displayed in gray, green, purple and orange; non-assigned points are displayed in light blue) and 2 clusters for random search (displayed in light and dark orange). We then visualize where those regions in behavior space mapped back in the intervention space, by applying the same coloring. (a) Looking at the curiosity search discoveries, we can see the non-linearity of the $I \rightarrow Z$ mapping, where small regions of intervention space can map to large regions of the behavior space (like the orange area) and reversely (gray area). We can also see the redundancy of the behavior space which is clearly concentrated in the left border of the space (ERK close to zero) which can seemingly be reached from very large portions of the intervention space (gray area). (b) Looking at random search discoveries, we can understand that it is very inefficient as it spends most of its exploration budget in the region of intervention space that converges to the left border in Z, and fails to explore the orange, purple and green regions discovered by curiosity search which seemingly lead to the more novelty in Z.

Finally, as the IMGEP efficiently drives the GRN into diverse goal states with minimal interventions, we propose that the diversity achieved by the IMGEP can serve as a good proxy metric of the GRN versatility. Notably, analysis of example systems in Figure 3 reveals that many GRNs can reach a broad spectrum of steady states. Whereas our database is limited to certain systems (see Materials and Methods) and might not be representative of all biological pathways, this observation underlines the existence of various phenotypes that can be realized. It also highlights the critical importance of identifying salient interventions that can effectively control cellular states within this spectrum of possibilities, notably as many cancer types are due to epigenetically non-identical cells [36].

Empirical Tests Reveal Robust Navigation Competencies

Figure 5: Identification of robust traversal strategies in transcriptional space. (a) Violin plots show, for each of the 432 systems (one point per system), the median sensitivity (over the K representative goal states) to the noise (green), push (gray) and wall (yellow) perturbation families. Violin plots on the right detail the median sensitivity for the 18 sub-families. (b-g) Each row provides examples of perturbed trajectories of either extremely-robust or extremely-sensitive example (GRN, Z) system (on average over the K goal states) for the three families of perturbations, as shown by annotations in (a). For instance, the first row (b) shows perturbed trajectories of the (model #10, nodes (3,7)) system which has the highest sensitivity to noise whereas the last row (g) shows trajectories of the (model #272, nodes (2,3)) system which has a nearly perfect robustness to walls. Each image contains an example trajectory for a given $(i,u)$, and one $u$ per sub-family is shown per column. For instance, in the first row (b), the trajectories are perturbed with the different sub-families of noise $\left( \sigma_{n} \in \lbrack 0.001,0.005,0.1\rbrack,p_{n} \in \lbrack 10,5,1\rbrack \right)$ which can be seen as various levels of difficulty. For each trajectory we annotate the starting position (A), endpoint prior perturbation (B), and endpoint after perturbation (B'), and show the original trajectory in black. The perturbed trajectory is shown in colorscale (from red at t=0 to cyan at t=3000 secs). (b) Except for few cases (trajectory #43), the system (model #10, nodes (3,7)) system is not robust to noise as its trajectories are easily deviated from the original endpoint. (c) The (model #52, nodes (4,7)) system however, except for rare cases (trajectory #35), consistently reaches its original target despite encountering various amounts of noise. Interestingly, trajectories #36 and #40 consistently follows a complex up->right-down->left path, despite the induced noise. (d) The (model #647, nodes (2,10)) system, except for few cases (trajectory #0), is typically deviated from its original trajectory when being pushed away. Interestingly though, it seems to follow similar (parallel) trajectories. (e) The (model #284, nodes (4,6)) system, is an example of an extremely robust system which, despite many push configurations (in magnitude and number), consistently returns to its original trajectory. Interestingly, the trajectories of this system are relatively complex with several loops and detours. (f) The (model #84, nodes (4,6)) system is not very robust to walls, and typically deviates or blocked when it encounters a wall. (g) The (model #272, nodes (2,3)) system is another example of an extremely robust system which, despite many wall configurations (in length and number), consistently returns to its original path. Once again interestingly, the trajectories of this system are relatively complex with several loops and detours.

We are then interested in characterizing the degree of robustness of the previously-discovered "goal states" in order to identify the ones that can consistently be reached by the GRN despite encountering various perturbations. Whereas many studies have proposed rigorous analysis of the "robustness" of biological networks [3][38], the generated perturbations often target variations in the regulatory rules (i.e. variations at the hardware level) and variations are often sampled independently (and prior) to observations of the GRN dynamical behaviors [29][82][4][28][7][143]. Here instead, we propose to conduct a battery of empirical tests that draw inspiration from classical "displacement experiments" [37][15] and "barrier experiments" [90] commonly used in behavioral sciences to assess the navigation competencies of various animals. As illustrated in Figure 2, we consider environmental perturbations that perturb the GRN trajectory with 1) various degree of noise in the gene expression levels, 2) sudden "pushes" during the GRN traversal of transcriptional space, and 3) energy barriers or "walls" acting as new force fields that constrain the GRN traversal. Importantly, those perturbations are conditioned on the observed behavior of the GRN. The magnitude of the noise and of the pushes is scaled proportionally to the extent of the observed trajectories, and the walls are generated in locations of the space that the GRN would "naturally" visit without the induced perturbation. While intuitive from a behaviorist point of view, where one would adapt experimentation when testing animals in different contexts (e.g. to study homing behavior of an ant and of a sea turtle, or of an ant in food deprivation and in reproduction phase) [141], robustness studies in systems biology often neglect those aspects. We propose that a behaviorist lens on robustness can help understanding forms of non-genetic resistance in transcriptional space, which is crucial for the development of therapeutic strategies [36].

To assess the degree of robustness of the discovered goal states, our evaluation procedure is the following. For each (GRN, I,Z) system of the database, we retrieve a representative set of trajectories previously discovered using the curiosity-search algorithm and subject these trajectories to $P = s \times r$ perturbations conditioned on the GRN goal-reaching trajectory $i \rightarrow z$ prior perturbation. Here, s represents the different perturbation distributions which correspond to various "tests" and "levels of difficulty" (e.g. noise magnitude and frequency, number of walls, etc.) and $r$ is the number of (stochastic) perturbations sampled per family. The pseudocode is illustrated in Figure 1-c and details about the different family of perturbations are provided in Materials and Methods. At the end of this process, the behavioral catalog is augmented with the perturbed trajectories
$H =$ {$( i_{k},o_{k},z_{k},~${$( u_{p},o_{p},z_{p}),~p = 1\ldots P$}$),~~k = 1\ldots K$}.

As the use of "spaces" comes with the notion of similarity and distance, we can then easily evaluate the sensitivity of a goal state $z$ with respect to a set of perturbation {$u_{p},p = 1\ldots P$} as the average distance in behavior space Z between the original trajectory endpoint $z$ and the perturbed trajectories endpoints {${ z}_{p}$}. Here our distance is simply the Euclidean distance, normalized by the extent of the trajectory prior perturbation in Z. We can then identify the so-called "robust goals" of the systems as the ones that have the lower sensitivity to perturbations. These sensitivity analyses can be useful in two important ways. On the one hand, they allow us to quickly identify the "extreme" examples of robustness, both at the system-level and at the goal-level, providing several insights into the degree of "competencies" that some biological networks might exhibit in their relative space (Figure 5). On the other hand, these analyses also allow us to map the heterogeneity of cellular responses and to better understand how non-genetic perturbations might modulate the landscape of reachable cell phenotypes (Figure 6).

Figure 5 shows the median sensitivity, over the representative goal states, for the 432 systems of our database and for the noise, push and wall perturbations families (as well as for the s=18 sub-families which correspond to varying degrees of perturbations). Overall, even though we observe varying degrees of sensitivity between systems (and between magnitudes of perturbations, which is expected), one first and interesting observation is that the median sensitivity remains relatively low, suggesting that GRNs could not only exhibit versatility (with respect to the considered interventions) but also robustness (with respect to the considered perturbations). In fact, looking at the "extreme" examples, we can identify quite impressive examples of complex and yet highly-robust space traversal strategies, with non-linear trajectories exhibiting many "detours" and "loops" but yet consistently reaching the same endpoint despite several pushes (Figure 5-e) or walls (Figure 5-g) on the way.

Figure 6: Energy landscape visualization based on the trajectory-based landscape generation method here, and constructed from different set of GRN trajectories, respectively trajectories generated (a) by the random search exploration, (b) by the curiosity-driven exploration, and (c) by the robustness tests experiments.

Figure 6 shows how the constructed catalog $H$ can be used to generate the energy landscape of the studied system. In biology, landscape formalisms have been used to comprehend the underlying dynamics of several systems, such as cell cycles and cell differentiation [64][33]. It is believed that such system-level visualizations could be particularly useful to apprehend non-genetic heterogeneity in the context of cancer treatment and stem cell differentiation [36][1]. A recent landscape-generation method only proposes to approximate the pseudopotential energy through simulation trajectories obtained throughout exploration of the system [1], making it a widely applicable method which we can directly apply here. However, the paper relied on Monte Carlo simulation to generate the trajectories. Due to the previously mentioned non-linearity and redundancy of the $I \rightarrow Z$ mapping, this can lead to poor estimation of the overall energy landscape (Figure 6-a). Instead, when generating the landscape from the trajectories discovered by our curiosity search exploration, we are able to reveal a new and wide "valley" of reachable states (Figure 6-b). Interestingly, the landscape-generation method can also be used to better comprehend the effect of external cues on the gene regulatory network, by visualizing how much they deform the energy landscape for instance leading to new shaped valleys (Figure 6-c). For the example system RKIP-ERK pathway [56], results highlighted a specific region of behavior space (with low RKIP and high ERK activation levels) that seems to be particularly robust, i.e. consistently reached by the GRN from certain initial conditions, and that might be associated with tumor development [57].

Possible reuses of the behavioral catalog and framework

Our framework generated a catalog of stimuli, responses, and navigation tests for the different GRN models contained in our database. Creating and sharing such a "behavioral catalog" with the scientific community is possibly one of the more exciting aspects of the work with new organisms. Furnished with such an empirically based data-set and detailed observations, one can 1) conduct statistical analysis across the population of studied organisms to inform fundamental research questions and 2) reuse the acquired knowledge to design specific behavior-shaping experiments in organisms of interest. As our framework focuses on observable behavior and is agnostic about the internal construction of the organism, another exciting perspective is to deploy it to different problem spaces and other classes of natural, chimeric or synthetic organisms. This section illustrates preliminary experiments along those three types of reuse.

To develop insights on the degree of sophistication of the different GRNs

Figure 7: Analysis and comparison of the degree of sophistication, in terms of versatility and robustness, between different classes of GRN. We categorize the GRNs by class of organism they belong to: plant, bacteria, slime mold, amphibian, rodent, homo sapiens, or generic. "n/a" refers to network models for which this information is not available. (a) Violin plots show the versatility of the 432 systems (one point per system) for each class. Versatility of one system is measured as the area covered by all the goal states discovered by curiosity search (equivalent to what we call diversity in Figure 3). (b) Trade-off (aka Pareto) mean and standard deviation curves that represent the trade-off among versatility and wall robustness performances as taken by the different classes of GRNs (standard deviation is divided by 4 for visibility). For each system, versatility (y-value) is measured as the area covered by the set of robustly achieved goal states, where the criterion of goal-achievement is a binary which tests whether the goal-reaching sensitivity (on average overall wall perturbations) is below a certain threshold (x-values). Violin plots in (a) are ordered in ascending order according to the class mean y-value at x=0.4 in (b).

A first use-case we explore is to conduct statistical analysis to categorize versatility and robustness in the surveyed networks on the basis of species in evolutionary strata. We consider seven categories, namely, plant, bacteria, slime mold, amphibian, rodent, homo sapiens, or generic. Here, generic corresponds to the networks not associated with any species but related to generalized biological processes. Please note that the surveyed database is relatively small with respect to the wealth of available models and biological pathways, so we can hardly claim that these results represent the true distribution of competencies across these organism categories. Still, as shown in Figure 7, results suggested interesting patterns.

First, on average, generic and Homo sapiens GRNs exhibit higher versatility (mean 0.228 and 0.238) compared to rodent and amphibian GRNs (mean 0.163 and 0.169), which in turn show higher versatility than bacteria and plant GRNs (mean 0.136 and 0.117). These findings are particularly intriguing in the context of the recently-formulated hypothesis of multi-scale competency architecture [16]: could the observed variation in versatility among different classes of GRNs contribute to the degree of versatility observed at higher-level scales? Collecting such experimental data for broader classes of organisms and GRNs will be crucial to understand how competencies at the molecular scale can impact the overall functionality and adaptability of organisms at higher scales, and to understand how evolution might have exploited this modular architecture for shaping the observed adaptivity and reprogrammability of biological systems.

Secondly, when comparing with the versatility of random networks (in black), generated to follow the same distributions of network size and connectivity as biological networks (as proposed in [63], see Materials and Methods), we observe that random network versatility is much lower (<0.026) than the versatility observed in biological networks. Once again, it is difficult to draw strong conclusions as the gene circuit model used for the random networks is relatively limited, whilst generic and studied across a range of biological contexts [54][93][131][39], and it will be interesting to scale the comparison to a broader and more complex range of ODE-based random models. Still, these findings hint that versatility prevalence might be a strong invariant of biological intelligence shaped by evolutionary processes.

Finally, we categorize the versatility-robustness tradeoff in the different categories of organisms. The idea is to compare the GRN competencies to robustly achieve diverse goal states, for different robustness thresholds. In Figure 7-b, we plot the mean and standard deviation pareto curves for the different categories of organisms and observe that, in average, the pareto-optimal solutions are mostly achieved by generic cell GRNs, even though bacteria GRNs can robustly reach more goal states for exigent robustness criteria (high x-values). The slime mold GRN can reach highly diverse goal states but the tradeoff quickly drops with wall perturbations, and there is only one system in our database belonging to this category so results might be not representative. Once again, those results are very interesting as generic cells GRNs are a building block that has been extensively reused by evolution across several organisms and contexts, bacteria have evolved to be very resistant (e.g. to antibiotics), and slime molds are a unicellular organism well known for its diverse capabilities, especially navigational ones [22][118][132][67].

For the development of therapeutic interventions

Figure 8: Identification of stimuli-based stepwise intervention triggering robust re-set of disease states into healthy physiological states. (a) The 10 most robust identified goal states (average sensitivity <0.05) and the corresponding reaching trajectories are displayed for the example RKIP-ERK signaling pathway [56]. We can see that most of them converge toward attractors in the "disease" region (orange). (b) Discovered stepwise stimuli intervention on MEKPP which we apply on states stuck in the "disease" region for 100 seconds. (c) The discovered intervention successfully brings back all points from the "disease" region closer to the target endpoint in the "healthy" region, and this under various tested perturbations (as shown in Supplementary Figure S2). The optimization procedure that led to the discovery of this intervention is described in the main text.

Understanding forms of non-genetic resistance and non-genetic heterogeneity is crucial across a wide range of cancer and treatment contexts [36]. Here, we illustrate how the constructed behavioral catalog can provide a fertile source for the design of therapeutic strategies, notably in the context of network control, using again the example of the RKIP-ERK signaling pathway [56]. In Figure 4, we saw that curiosity search revealed four clusters of reachable steady states for this system. From a clinical perspective, one might denote the green cluster as "healthy" region of behavior space and the orange cluster as "disease" region of the behavior space, as high levels of ERK and low-levels of RKIP are often linked to tumor development [57]. In Figure 8-a, we plot those two clusters as well as the 10 more robust goal-reaching behaviors from the behavioral catalog of this system, i.e. the goal states with the lower average sensitivity to the induced perturbations. We see that 6 out of the 10 more robust trajectories end up in the "disease" region, suggesting that certain configurations of initial state are very likely to reach that region despite chemical blockers (here pushes, walls, and noise), which was also visible on the system's energy landscape in Figure 6-c. Looking at the six trajectories, it seems that they all follow similar patterns where RKIP activation level increases past a certain threshold, and only then converge toward the disease region. This might already provide an interesting biomarker for prediction of tumor development, but what we are really interested here is to build upon that knowledge to develop stimuli-based interventions allowing to re-set the GRN setpoints from the identified "disease" steady states back to steady states within the identified "healthy" region . To do so, we define a parameterized stimuli-based intervention and a performance function, and search for parameters that optimize this performance. For the intervention function, we use a piecewise constant function that determines which nodes to intervene on (here MEKPP), when to apply the intervention (here every 10 seconds for 100 seconds), and with what amplitude (which are the parameters that we are seeking to optimize). The choice of the intervention function, which is arbitrary in this example, would typically depend on the experimental constraints, e.g. which nodes can be targeted with drugs and at which precision. For the performance function, we define the centroid of the "healthy" region as the target setpoint and compute performance of the stepwise intervention as the average distance of the novel setpoints (after intervention when starting from the 6 disease setpoints) to the target setpoint, and under a distribution of stochastic walls, pushes and noise perturbations. Hence a successful intervention should re-set the disease setpoints to healthy setpoints for all discovered disease states and robustly across the various tested perturbations. For optimization, we simply perform random search as this was sufficient here to discover one intervention (as shown in Figure 8-b) that successfully reset the setpoints (as shown in Figure 8-c) under various tested perturbations (as shown in Supplementary Figure S2). Here random search was sufficient to find a successful intervention, but more advanced optimization strategies like evolutionary algorithms or stochastic gradient descent could be envisaged for harder problems. Overall, mapping the "latent" behavioral abilities of GRNs in healthy physiology and disease states may have important implications for the identification of robust stimuli-based interventions that focus on behavior shaping instead of micromanaging all molecular states, and that can be exploited in therapeutic contexts.

As an alternative strategy to gene circuit engineering