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

Figure 9: Comparison of four alternative strategies for the design of oscillator circuits: curiosity search (blue), random search (pink), gradient descent (orange) and an evolutionary algorithm (green). (a-c) Given a budget of 5000 experiments, curiosity search is able to find 1167 oscillator circuits (ones showing sustained oscillations), whereas random search only finds 42 oscillators, and optimization-driven search fail to discover them (only one discovered by CMA-ES and none for gradient descent when starting from a random initialization). (a) 3D scatter plot of the 42 random search discoveries (pink) and 1167 curiosity search ones (blue) in the (amplitude, main frequency, offset) analytic behavior space. (b) Box plots projecting points from the 3D scatter plot into the respective (amplitude, main frequency, offset) axes. (c) Diversity discovered throughout exploration, where diversity is measured with a binning-based space coverage metric (20 bins per dimension). (d-e-f-g) Best discoveries (for which $L$ is minimal) made by the four exploration strategies. (h) Evolution of the optimization loss $L$ for the four algorithm variants. (i) Evolution of local training loss when finetuning the best IMGEP (blue) and RS (pink) discoveries with gradient descent, with the finetuned results displayed in (j-k).

The final type of reuse we explore is not a direct reuse of the constructed behavioral catalogs, but rather a reuse of the proposed automated tools to reveal different kinds of behaviors in a bioengineering context. A common problem in synthetic biology is to optimize the configuration and parameters of a gene model network to optimally perform a desired functionality, also known as gene circuit engineering [119]. Recent approaches rely on optimization-driven machine learning strategies, such as evolutionary algorithms and stochastic gradient descent. However, choosing the right loss function and parameter initialization for these optimization methods is a well-known problem in machine learning. These issues can lead to optimization algorithms getting trapped in local minima within the complex landscape of possibilities. In response to these challenges, we propose to investigate whether the curiosity-driven exploration strategy can be employed as an alternative (diversity-driven) strategy. Whereas traditionally-employed for exploratory purposes, these exploration strategies were also shown to facilitate the resolution of external, pre-defined tasks characterized by sparse or deceptive rewards [74], by effectively exploring solution space.

Here, we consider the target application of oscillator circuit engineering followed in [133], where parameters of a gene circuit model are optimized to produce oscillation patterns with target amplitude $A$, frequency $w$ and offset $b$. This time, the intervention space includes both genetic interventions (setting kinematic parameters of regulatory rules) and environmental interventions (setting the initial state $y_{0}$). We then compare several exploration strategies: a random search, two optimization-driven strategies, one using gradient descent as proposed in [133] and one using an evolutionary algorithm, and finally a diversity-driven strategy using the proposed curiosity search. All algorithms are given the same experimental budget $(N = 5000)$. For curiosity search, the behavior space $Z$ is defined as the image space of the discrete Fourier transform of the observation. We then use the exact same IMGEP algorithm as before, but operating within the novel problem spaces $(I,Z)$. For gradient descent, we follow the procedure proposed in [133]. We define a loss function which measures the mean square error between the observed node activation levels $y$ and the target oscillation (represented as a cosine wave). We then randomly initialize the parameters $i \sim U(I)$ and use Adam optimizer for N=5000 optimization steps. For the evolutionary algorithm, we use the CMA-ES algorithm [148], with the negative of the loss as fitness function. In addition, we also use gradient descent for local refinement of the best discoveries made by the other exploration strategies (curiosity search and random search), this time with a limited budget of $N = 100$ optimization steps.

In Figure 9, we show that curiosity search is again significantly more efficient than random search in revealing a diversity of possible oscillator behaviors. Out of 5000 trials, random search was able to find only 42 configurations leading to sustained oscillations whereas curiosity search was able to find 1167 (and gradient descent did not find any). Without focusing on the target objective, curiosity search is able to efficiently cover the analytic $(A,\omega,b)$ space (Figure 9, a-c), thus discovering oscillators close to the target one (Figure 9-d). Instead, when starting from a random initial condition, gradient descent is very likely to get trapped in a local minimum where it converges to the target offset b but fails to produce oscillations (Figure 9-f and 9-h). While CMA-ES explores more the solution space at the beginning of optimization than SGD does, it also ultimately converges into a similar local minima. However, whereas the global optimization strategies are unsuccessful in this example, they seem to be useful to locally refine close-enough solutions, as can be seen here when refining the best discoveries made by curiosity search and random search with gradient descent (Figure 9-i, 9-j). These results suggest that a diversity-driven exploration strategy, eventually combined with a more advanced local optimization strategy, can offer promising and cost-effective alternatives for the design of synthetic gene networks. More generally, as our framework only relies on empirical investigation for inferring the mapping between interventions and behaviors (treating them as abstract variables in observable problem spaces), we believe it offers an exciting perspective to be deployed across various problem spaces and classes of organisms.

Discussion

This paper presents a novel framework aimed at uncovering the navigation competencies of gene regulatory networks (GRNs). The framework conceptualizes GRNs as agents actively navigating the transcriptional space and provides a set of tools, leveraging computational models of curiosity-driven learning and exploration, with a battery of empirical tests inspired from behaviorist tradition, for automated experimentation and behavioral characterization. The proposed framework is novel in two central ways. First, it introduces a novel AI-based toolbox to the field of biological network analysis. We show how this toolbox, leveraging the successful ingredients of recent intrinsically motivated learning algorithms - originally developed to enable robotic AI agents to explore and learn diverse skills in novel and unstructured environments [136][65] - can be transposed to assist efficient discovery of behavioral abilities within biological pathway models like GRNs. Secondly, rather than merely mapping the attractor states [9][129][103] or analyzing their sensitivity to model parameter changes [145][25] as extensively proposed in conventional GRN analysis methods, our framework investigates the dynamic adaptability of these networks' navigation competencies in response to various changing environmental conditions. With this approach, our aim is to uncover whether diverse competencies, analogous to the ones exhibited by living agents, can be found within physiological network dynamics. Notably, these competencies are discovered without necessitating structural alterations to network properties or connectivity. Importantly, our framework and its associated tools do not make any assumptions about the structure or origin of the biological network, making it in theory adaptable to the study of diverse unconventional intelligences across various domains.

By applying this framework to a curated database of GRN models, we discovered a diverse range of behavioral responses that GRN can exhibit under different initial conditions and characterized their robustness to various perturbations. Notably, our analysis revealed a number of interesting aspects of navigation of the state space which can be leveraged in several contexts. These automated tools form the first step towards cost-effective in silico simulation and interrogation platforms; as the "behavioral catalogs" produced by this process can be a first stepping stone for better understanding the GRN functionalities as well as for designing drug-driven interventions in a biomedical or bioengineering context.

There are several limitations and avenues for future work to this study. First, these networks are studied as a model in isolation and it is possible that some of the ODE models (or solvers) provide spurious behaviors within certain parameter ranges which might not map to observable phenotypes in vitro. Interestingly, this limitation also suggests an interesting further direction to this work: using the automated discovery toolbox to assist model inference, allowing to efficiently identify the rare or unexpected behaviors of the ODE model and suggest whether further refinement is needed or not. Another interesting direction for future work, as our framework considers the GRN model as a black-box and works with limited experimental budget, would be to directly apply it to in vitro GRN models at the bench. One could for instance integrate experimental constraints to the search by defining families of empirically-testable interventions and perturbations, as well as specify clinically-relevant goal spaces and perturbations. Even if in a biological setting versatility and robustness phenomena may be harder to detect, or harder to alter, these results can be used to (1) design synthetic biology circuits with advanced capabilities [138], and (2) conduct studies of subcellular proto-cognitive phylogenetics, to help understand the evolutionary pressures for and against reprogrammability in cell regulatory machinery. Another limitation of our work is that we consider predefined problem spaces, here the space of GRN steady states (or Fourier descriptors of the dynamics in the bioengineering example). The dynamics of gene regulatory networks are relatively simple (usually converge to stable points or periodic orbits) allowing such hand-defined descriptors. To scale the framework to higher-dimensional and more complex problem spaces, recent works from the IMGEP literature suggest using unsupervised learning of goal space representations [66][72]. Whereas these works were applied to abstract models of multicellular patterning, similar works could be envisaged in more realistic systems, such as sophisticated model of multicellular morphogen and/or bioelectrical patterning which were used to suggest in-vitro experimental manipulations [127][80].

The tools presented here, and the behavioral repertoire we identified, are just the beginning, and much work remains. Future efforts must test additional competencies across the spectrum of cognition (memory, creative problem-solving, valence, etc.) and extend the tools we presented here to explore them . The predictions made by our computational tools can now be tested in real cells, using emerging tools for physiological profiling in the living state and a diverse set of biochemical, biomechanical, and bioelectrical perturbations. We anticipate a tight and productive feedback loop between computational theory that suggests new experiments, and results in living cells that greatly extend our computational perspective on what they can do [115][112][114][108][113]. Such interdisciplinary work, pulling together concepts and techniques across fields, is likely to have major implications for fundamental understanding of evolution, intelligence, and dynamical control, as well as drive novel kinds of therapeutics that leverage the innate behavioral competencies of living matter[97][59].

Methods

Read more

GRN models and numerical simulation

This study employs ordinary differential equation (ODE) models to represent molecular pathways, with nodes representing pathway components and edges capturing their interactions. The continuous node states, encompassing variables like gene expression levels and protein concentrations, are interconnected through a system of ODEs, enabling the modeling of complex regulatory dynamics. ODE models are often available in the Systems Biology Markup Language (SBML), a standardized format that contains essential information about variables, parameters, equations, and model metadata in XML files.

To perform numerical simulations of ODE SBML models, we rely on the SBMLtoODEjax python library, a recent development that automates the parsing and conversion of SBML models into python models written entirely in JAX [27]. Taking advantage of JAX computing capabilities, SBMLtoODEjax enables efficient and parallel numerical solutions for gene expression levels and other node states by recursively invoking the generated python models to integrate the ODE equations with current gene expression levels. Additionally, we have developed a python library (https://github.com/flowersteam/autodiscjax) comprising additional modules and pipelines that facilitate interventions on the GRN models such as genome or drug interventions, as well as other perturbations such as noise, pushes, and walls that can be applied to the states and kinematic parameters of gene regulatory networks.

Given the model species initial state $y(t = 0)$, the desired rollout length $T(secs)$ and step size $\Delta T$, as well as the chosen intervention $i$ and/or perturbation $u$, the model rollout iteratively 1) integrates the system of ODE-governed equations that specifies the rate of species changes $\frac{dy}{dt}$ using JAX odeint solver to update model species $y(t) \rightarrow y(t + \Delta T)$, 2) calls the model assignment rules to update kinematic parameters if needed, and 3) apply the intervention and/or perturbation function to update $\left( y(t + \Delta T),w(t + \Delta T),c \right)$ accordingly. In this paper we use $T = 2500secs$ and $\Delta T = 0.1$ (25 001 time points per rollout including $t_{0}$). The ODE solver uses an absolute tolerance of $1e^{- 6}$ and relative tolerance of $1e^{- 12}$, with maximum number of solver steps of $1000$. For a step-by-step guide on utilizing these libraries within the proposed framework, we refer interested readers to our tutorial (https://developmentalsystems.org/curious-exploration-of-grn-competencies/tuto1.html), which offers practical examples and detailed instructions.

Experimental setup

In our computational models, we are able to record the activities of all nodes during a model rollout. The observation space $O \subset {\mathbb{R}}_{+}^{n \times \frac{T}{\Delta T}}$ is such that $o = \left( y(0),\ldots,y(T) \right)$ where $y(t)$ represents the n-dimensional vector of node states at each time step, with T being the total reaction time. The boundaries of the observation space are not known.

Regarding the exploration of problem spaces, namely the intervention space I and behavior space Z, we specify them as follows.

For the main experiments on biological networks, the intervention space $I \subset {\mathbb{R}}_{+}^{n}$ consists of initial node states sampled from the hyper-rectangle $\left\lbrack y_{0,\min},y_{0,\max} \right\rbrack$ where $y_{0,\min} = \frac{1}{r} \times y_{d,\min}$ and $y_{0,\max} = r \times y_{d,\max}$ with $r = 20$ and $\left( y_{d,\min},y_{d,\max} \right)$ the minimum and maximum of each node of the model over the default time course simulation (with initial conditions provided in the SBML file and T=25000). On the other hand, the behavior space $Z \subset {\mathbb{R}}_{+}^{2}$ endpoint states $z = \left( y_{i}(T),y_{j}(T) \right)$ where $(i,j)$ corresponds to the target phenotype nodes. We ensure that most trajectories have reached stable states at T = 2500 (as elaborated in the next section) such that Z can be viewed as the space of reachable endpoints, whose boundaries are not known.

Database creation

Biological networks database

All the ODE models we use in this work are downloaded from the BioModels database [120][121] in SBML format. From all models referenced on the website, we only consider the ones that are curated, that have at least 3 nodes, and that are handled by the SBMLtoODEjax simulator (as SBMLtoODEjax does not handle models with discrete events, custom functions or other specific cases as detailed in [27]. Note that the original models 37, 262, 263, 284, 454, 455, 459, 461 and 624 have been previously analyzed while clamping certain nodes at fixed values (as detailed in table S1). Here, we relaxed this condition for a more realistic simulation in which all of the nodes' concentration are free to vary.

To ensure the inclusion of models suitable for our analyses, we then applied specific filters to the collected models.

First, we simulated the default model rollout for each model to obtain the concentration profiles of the pathway components over a short time span (T=10 secs and $\Delta T = 0.1$). We discarded simulation results containing invalid values (NaN or negative concentrations) or those that took an excessive amount of time (>1sec). While it is acceptable that a rollout sometimes returns NaN values (when there are no solutions given ODE tolerance options for specific initial conditions), we consider the model invalid if this occurs for the default initial conditions provided in the SBML file.

For the remaining models, we conducted further simulations with an extended time span (T=2500) and 50 random initial conditions uniformly sampled within the model's intervention space $I$ (as defined before). Once again we discarded models whose batch simulations took an excessive amount of time ($>15$ secs). From the remaining models, we derived the resulting 50 trajectories for each node pair (i, j) and subjected them to additional filters to refine the database. We removed node pairs where either 1) [filter F1] a substantial proportion of trajectories $\left( \ge 20\% \right)$ exhibited invalid concentrations (NaN or negative) or unsettled behaviors $(\exists t \ge 2400 \text{ such that } |y(t)-y(T)| \ge0.02 \times |y(T)-y(0)|)$ or periodic patterns $(\exists f>0 \text{ such that } |S(f)| \ge 40 \text{ where } S=DFT(y(\frac{T}{2}, \dots, y(T)))$; or [filter F2] the reached space in $Z$ was too small$\left( \left( \max_{k = 1\cdots 50}y^{k}(T) - \min_{k = 1\cdots 5}y^{k}(T) \right) < 0.1 \right)$ to discard cases where "diversity" could result from floating point rounding errors; or [filter F3] the number of attractors was less than four $($ {$ y^k(T)$}$_{k = 1 \dots 50}\text{ cover } \leq 4 \text{ bins over a }20 \times 20 \text{ binning of } Z)$.

Upon completion of the filtering process, our final database comprised 30 models, consisting of a total of 432 systems, as detailed in Supplementary Table S1. These curated models and systems served as the foundation for our subsequent analyses and investigations into the navigation competencies of the molecular pathways.

Random networks database

Following the methodology proposed in [63], we aimed to create a database of synthetic networks with topologies similar to those of the biological networks, but with random regulatory rules instead of evolved ones. The objective was to compare the versatility and robustness competencies between biological and random networks, akin to the approach used for memory competencies in [63]. To achieve this, we initially generated 300 networks based on the transcriptional gene circuit model [54], ensuring that they had the same distribution of network size (number of nodes) and connectivity (nodes in-degree) as the biological network database (using fitted gaussian distributions). The kinematic parameters $W,b,\tau$ of these networks were randomized $( W \sim [-30,~30]^{n\times n}, ~ B\sim [-10,~10]^n, ~ \tau \sim [1,~15])$ where model step is defined as $y(t + 1) = \frac{\Delta T}{\tau} \times sigmoid(Wy + B) + \left( 1 - \frac{\Delta T}{\tau} \right) \times y$ and in-degree connectivity is enforced by setting some weights of $W$ to zero. However, during the creation process, we observed that none of the generated networks met the criterion for exhibiting a sufficient number of steady states (criterion F3). This limitation arose from the inherent constraints imposed by the gene circuit model's shape of ODE equations, limiting the diversity of possible dynamical behaviors. As our focus was on networks with a possible spectrum of steady states, akin to the biological network database, we decided not to pursue further analyses on these networks.

Instead, we selected the systems (models and pairs of nodes) that demonstrated the highest versatility (metric detailed below) from among all the generated systems that passed the filters F1 and F2. The selected networks' versatility is presented in Figure 7, but for future research, it would be interesting to explore broader and more complex classes of equations to assess their potential for achieving higher behavioral diversity.

Sanity Check

We also tested our exploration pipeline on two simple models, namely BIOMD0000000341 as described in the paper "A model of beta-cell mass, insulin, and glucose kinetics: pathways to diabetes" [149] and BIOMD0000000454 as described in the "Example One" of the "Metabolic Control Analysis: Rereading Reder" paper [150]. For both these models, the ground truth goal states (attractors of the models) are already known, because easy to find analytically or numerically as described in the corresponding papers. As a sanity check, we validate that our exploration pipeline is able to re-discover those goal states and added those results in Figure S4 of the supplementary for completeness. Note that a more complex variant of BIOMD0000000454 is also included in our main biological networks database where we let the metabolite concentrations $y_1(t), \dots, y_5(t)$ evolve in time as described in the caption of table S1. We illustrate results of our curiosity-driven exploration method and of a random search method on this more complex variant in Figure S4-c, S4-d.

Curiosity-driven exploration

This section provides additional information about the internal models and hyperparameters of the intrinsically-motivated goal exploration process. The overall IMGEP pipeline is illustrated in Figure 1-c. To sample a goal, the IMGEP uses a uniform sampling strategy within the bounding hyper-rectangle of currently reached goals (scaled by a factor 1.3). Hence sampling bounds adapt to the discoveries and do not need to be predefined via expert knowledge. The volume of the hyper-rectangle is larger compared to the cloud of currently-reached goals, which incentivizes targeting unexplored areas outside of the cloud and promotes diversity in the exploration process. Then, to generate an intervention for achieving the sampled goal, the IMGEP selects the nearest previously reached goal in $Z$, identifies its associated intervention, and performs a local random step from that point $(stepsize \sim \mathcal{N}(0, 0. 1 \cdot [y_{0,\max} - y_{0,\min}])$ in the intervention space.

While our implementation choices for the IMGEP goal representation, goal generation, and goal-conditioned optimization are relatively straightforward, it is worth noting that alternative strategies could be considered for each of these components for more complex problems. The python library AutoDiscJax (https://github.com/flowersteam/autodiscjax) that accompanies this paper can be used to implement this and other IMGEP variants in JAX.

Robustness tests

We define 3 family of perturbations: 1) the noise perturbation $U_{n}\left( \sigma_{n},p_{n} | y \right)$ which is parametrized by its standard-deviation (scaled proportionally to the extent of the observed trajectory $y$ prior perturbation) and period (secs); 2) the push perturbation $U_{p}\left( m_{p},n_{p} | y \right)$ parametrized by its magnitude (proportional to the extent of $y$) and number of occurrences; 3) the wall perturbation $U_{w}\left( l_{w},n_{w} | y \right)$ parametrized by its length (proportional to the extent of $y$) and number, and where walls are generated in locations of the space that the GRN would "naturally" visit without the induced perturbation. Details about the implementation of walls are provided in Supplementary Figure S3.

To assess the robustness of the GRN systems in our database, we employ an evaluation procedure, as depicted in Figure 1-d. For each system $(I,Z)$ in the database with its corresponding behavioral catalog $H$ discovered using the curiosity-search algorithm, we perform the following steps. We first retrieve $K$ representative trajectories out of the $N$ discoveries, i.e. ones that cover well the reachable space. To do so, we randomly sample tuples of K discoveries (among N) 500 times, and select the one with the maximum diversity. One could test all trajectories with K=N but here we use K=N/10 mainly for compute reasons, as we run the experimental campaign on all 432 systems. Next, we subject each of these K trajectories {${ y}_{k},k = 1\ldots K$} to s=18 different perturbation distributions, each representing various levels of difficulty:
$\left( \sigma_{n},p_{n} \right) \in$ {$(0.001,5),(0.005,5),(0.1,5),(0.005,10),(0.005,5),(0.005,1)$},
$\left( m_{p},n_{p} \right) \in $ {$(0.05,1),(0.1,1),(0.15,1),(1,0.1),(2,0.1),(3,0.1)$},
$\left( l_{w},n_{w} \right) \in $ {$(0.05,1),(0.1,1),(0.15,1),(1,0.1),(2,0.1),(3,0.1)$}.
In each perturbation distribution, we sample r=3 random perturbations, resulting in $P = s \ast r$ perturbations. For each perturbation in the set {${ u}_{p},p = 1\ldots P$}, we re-run the trajectory starting from the same initial state $i$ but with the sampled perturbation applied $\left( i,u_{p} \right)$, and observe the resulting outcome $\left( o_{p} \right)$ and reached endpoint $\left( z_{p} \right)$.

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

Evaluation Metrics

Diversity measure

Diversity is measured by the area that explored observations cover in behavior space Z. Each single exploration results in a new point in this space, such that diversity measures how much area the algorithms explored in those spaces.

In general, existing approaches in the NS, QD and IMGEP literature use binning-based metrics [66][72][102] or distance-based metric from ecology [147] to quantify the diversity of a set of explored instances. However, those metrics are sensitive to the binning strategy, or fail to discriminate between qualitatively significantly different explorations [26]. Another approach, called the threshold coverage, measures diversity as the volume of the union of the set of hyperballs of radius $\epsilon$ that have for centers the observed effects
{$z \in Z$}. This diversity measure, while difficult to compute in high-dimensional spaces, avoids the pitfalls of bin-based and distance-based metrics and is easily computable in 2-dimensional spaces [26].

Threshold coverage quantifies the area of the space that has been reached at a given precision $\epsilon$(the threshold), and is what we used in Figure 3 to compare random search and curiosity-driven exploration strategies.

Sensitivity measure

In general, existing approaches in systems biology and evolutionary genetics measure sensitivity (opposite of robustness) in a relative manner with respect to 1) a functionality [3] or phenotypic trait [38] of interest, 2) specific perturbations (environmental or genetic changes), and 3) a measure of the degree of variation. Here, we adopt a similar metric where 1) the phenotypic trait of interest is defined as a goal state $z \in Z$ discovered by curiosity search, 2) the set of perturbation {$ u_{p}$} is defined in previous section and conditioned on the GRN goal-reaching trajectory $i \rightarrow z$, and 3) variation is measured as the Euclidean distance in behavior space, normalized by the extent of the trajectory prior perturbation in Z.

This distance-based sensitivity measure proves straightforward as we explicitly use "spaces" to observe and analyze behaviors. The results of this sensitivity analysis are presented in Figure 5.

Versatility-Robustness measure

In this study, we introduce the terms "diversity" and "versatility" to characterize the competencies of the exploration agent (IMGEP) and the gene regulatory network agent (GRN), respectively. Diversity refers to the ability of the IMGEP agent to reveal a wide range of behaviors in the GRN, while versatility refers to the capability of the GRN agent to reach diverse goal states. The GRN versatility is unknown, and can only be approximated via proxy metric. Here, we consider that the diversity of the IMGEP (measured with the threshold coverage metric) is a good approximation of the versatility of a given GRN, as the IMGEP was shown to efficiently drive the GRN into diverse possible goal states.
In Figure 7a, we employ this diversity metric to categorize the versatility of surveyed networks based on the class of organism they belong to. For the random networks, as they all have less or equal than 4 attractors, the versatility remains below $0.026 = 4 \times\frac{\pi\epsilon^{2}}{(1+2\epsilon)^{2}}.$

Figure 7b, we introduce the versatility-robustness metric, which conditions the diversity metric on a sensitivity threshold. Only goal states with sensitivity to perturbations below this threshold are considered when computing the reached area of the space. A high versatility-robustness score indicates that diverse goal states are achieved with a high level of precision.

Experiments on the RKIP-ERK signaling pathway

This section details the additional experiments conducted on the RKIP-ERK signaling pathway [56]. We refer to the accompanying notebook tutorial for reproducing these experiments: [https://developmentalsystems.org/curious-exploration-of-grn-competencies/tuto1.html]{.underline}.

For Figure 4, clustering in behavior space was performed using the HDBSCAN algorithm [73] with hyperparameters set as min_cluster_size=10 and cluster_selection_epsilon=0.1. Points in the 10-dimensional intervention space are visualized by applying a TSNE 2-dimensional reduction. To visualize the clusters in behavior space (and corresponding clusters in intervention space), we fitted polygons on the cluster points using shapely library unary_union, dilatation, and erosion operations [24].

In Figure 6, we generated trajectory-based energy landscapes following the method proposed in [1]. Energy landscapes provide an intuitive way to understand how a system with multiple steady states behave, by picturing it as a ball rolling downhill towards low-energy valleys (steady states). Given a set of trajectories in behavior space Z, we constructed a probability distribution (P) of system states and converted it into a pseudopotential energy surface (U = −ln(P)). This energy surface was smoothed using cubic spline interpolation and visualized using Plotly 3D surface plots. Figure 6a, 6b, and 6c differed by the input set of trajectories used for generating the landscape: a) employed the set of trajectories discovered by random search, b) used the set of trajectories discovered by curiosity search, and c) utilized the set of trajectories generated by robustness tests.

In Figure 8, the "healthy" and "disease" clusters were the same as in Figure 4 and visualized similarly. We displayed trajectories with the lowest sensitivity (averaged over all $P = 3 \times 18$ perturbations). The stimuli-based intervention shown in Figure 8b was found using a simple random search procedure. First, we defined an arbitrary target node and a stepwise node-activation function, clamping MKEPP values to desired values $x = \left\lbrack {y_{MEKPP}}^{(1)},\cdots,{y_{MEKPP}}^{(10)} \right\rbrack$ every 10 seconds for 100 seconds. Then, we randomly sampled x within a range of values near the MKEPP current steady states (endpoints from the 6 "disease" trajectories, assuming that the drug intervention cannot drastically remodel those values). For each candidate x, we ran new trajectories starting from the disease states and applying the intervention x under a distribution of noise, push, and wall perturbations. Finally, we selected the intervention x that most successfully brought ERK-RKIP levels back to the target setpoint (centroid of the healthy region). The resulting intervention (shown in Figure 8b) succeeds to robustly reset all 6 disease state points despite perturbations, as shown in Figure 8c. We refer to the notebook for reproducing the experiments.

Experiments on synthetic gene networks

This section details the additional experiments conducted on the synthetic gene networks (Figure 9). We refer to the second accompanying tutorial for the full codebase: https://developmentalsystems.org/curious-exploration-of-grn-competencies/tuto2.html.

In these experiments, we consider the target application of gene circuit engineering followed in [133], where parameters of a gene circuit model are optimized to produce target oscillator patterns. The gene circuit model employed in [133] is the same than the one used for the random networks database (Eq 1), with $\tau = 1$. Hence the intervention space is a $n^{2} + 2n$ dimensional space defined as $I = \left\lbrack y_{t = 0,\min},y_{t = 0,\max} \right\rbrack \oplus \left\lbrack W_{\min},W_{\max} \right\rbrack \oplus \left\lbrack B_{\min},B_{\max} \right\rbrack$, with $y_{0,\min} = 0,$ $y_{0,\max} = 1,$ $W_{\min} = - 30,$ $W_{\max} = 30,$ $B_{\min} = - 10,$ $B_{\max} = 10$. Here we consider networks of n=3 nodes, with the first node being the target phenotype node. Thus, what we seek here is kinematic parameters $(W,B)$ and initial concentrations $y_{0}$ that would produce a periodic pattern $y =[y_{n=0}(0), ~\cdots~y_{n=0}(T)]$ with target amplitude $A$, frequency $w$ and offset $b$. Here, the target $(A,\omega,b)$ are sample randomly with $A \sim U\left( \lbrack 0.1,0.5\rbrack \right),$ $b \sim U\left( \lbrack A,1 - A\rbrack \right),$ $\omega \sim Beta(\alpha = 2,\beta = 8)$.

We then compare three alternative exploration strategies: 1) curiosity search, 2) random search and 3) gradient descent, i.e. pure optimization-driven search as proposed in [133], all given the same experimental budget $N = 5000$.

For curiosity search, the behavior space $Z$ is defined as the image space of the discrete Fourier transform of the 1d-signal $y$, where distance in the space measures average difference in spectral amplitude. The IMGEP algorithm is then the same that the one previously used, as detailed in Figure 1-c, but operating within the novel problem spaces $(I,Z)$.

For random search, interventions are sample uniformly $\left( i_1,~ \cdots ~,i_{N} \right) \sim U(I)$.

For gradient descent, we follow the procedure proposed in [133]. We define a loss function which, for a set of parameters $i = \left( W,B,y_{0} \right)$, measures the mean square error between the phenotype node activation levels $y$ and the target oscillation represented as a cosine wave with the desired $(A,\omega,b)$ : $L = \sum_{t}(y(t) - (A\cos(2\pi\omega \times t) + b))^2$. We then sample a random parameter $i \sim U(I)$ and use Adam optimizer with $l_{r} = 10^{- 3},$ $b1 = 0.02,$ $b3 = 0.001,$ $\epsilon = 10^{- 8}$ for $N = 5000$ optimization steps (same number of model rollouts allowed than for curiosity search and random search).

In addition, we use gradient descent for local refinement of the best discoveries made by the other exploration strategies (curiosity search and random search), this time with a limited budget of $N = 100$ optimization steps.

Visualizations in Figure 9 show: (a-b) the oscillators discovered by random search and curiosity search (gradient descent did not find any oscillator in this example) in the $(A,\omega,b)$ space, (c) the corresponding diversity (using this time a binning-based space coverage measure with $20^{3}$bins as the space is 3-dimensional), (d) the evolution of the training loss $L$ throughout the N=5000 trials for the three exploration strategies, (e-f-g) the corresponding best discoveries (for which $L$ is minimal) for the three exploration strategies, and (h-i) the local training loss and resulting finetuning of the best discoveries with gradient descent.

Statistics

Statistical analyses were performed using custom code in python 3.9 (relevant libraries include jax 0.4.8 and scipy 1.10 and plotly 5.16.1). Welch’s two-sample t-tests reported in Figure 3 are two-tailed, and additional details on statistical analyses are provided in the figure legend.

Data Availability

Source code is available on GitHub at https://github.com/flowersteam/curious-exploration-of-grn-competencies. It contains experimental data and an executable notebook version of the paper to reproduce all paper figures, as well as additional step-by-step tutorials to reproduce results from scratch for Figures 4, 6 and 8 (tutorial 1) and Figure 9 (tutorial 2), as well as the codebase to reproduce the whole experimental campaign. All our codebase is open-source under MIT License.

The ODE models we use in this work are downloaded from the BioModels database, and ID references of these models are detailed in Supplementary Table S1.

Acknowledgments

We thank Patrick Erickson and Randall Jordan Ellis for review and discussion, as well as Tom Cirrito, Wesley Clawson and Santosh Manicka for useful discussions. We also thank Alexander Mordvinstev for providing the executable paper template, as well as Julia Poirier for assistance with the manuscript.

The authors acknowledge support from the biotechnology company Poietis and the French National Association of Research and Technology (ANRT), as well as from the French National Research Agency (ANR, DeepCuriosity AI chair project). M.L. gratefully acknowledges funding support from Astonishing Labs, and from the Templeton World Charity Foundation via grant TWCF0606.

This work also benefited from the use of the Jean Zay supercomputer associated with the Genci grant A0151011996.

Author Contributions

M.E, M.L, PY.O and C.MF conceived the project together. M.E. wrote the code and performed the computational experiments and data analysis. M.L, PY.O and C.MF supervised the project. M.E, M.L, PY.O and C.MF wrote the paper.

Corresponding Author

Correspondence to Michael Levin: michael.levin@tufts.edu

Competing Interests Statement

M.L.’s lab receives funding in the form of a sponsored research agreement from, and consults for, Astonishing Labs, which has interest in therapeutic applications of GRN learning behaviors.

Supplementary Material

Figure S1: Examples of interventions that can be implemented within the accompanying AutodiscJax software. All those examples can be reproduced in the accompanying tutorial 1. (a) Numerical simulations with interventions can be performed in parallel by vectorizing simulations over different intervention parameters, simply using the jax vmap operator. This offers a convenient (and fast) way to test several interventions in the biological network, as shown here for testing the network under various initial conditions in batch mode. 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. In this example, despite the numerous interventions, the GRN trajectories still converge to the same endpoint B. (b) Example intervention where species amounts are clamped to specific values. In this example the node MEKPP is clamped to $2.5\mu M$ for 10 seconds at t=0 and then to $1\mu M$ for 10 additional seconds at t=400. In this example, after the first clamping the GRN trajectory still follows a similar S-shape curve and arrives close to the original endpoint B but after the second clamping, ERK expression levels are shifted to a higher steady state B'. (c) Example intervention where the numerical value of one kinematic parameter of the model (k5) is changed from 0.0315 to 0.1. In this example we can see that changing the parameter k5 shifts the trajectory end point quite significantly, but qualitatively the trajectory seems to preserve a similar S-shape.

BioModel ID reference	number of nodes	organism class	systems (observed nodes pairs)
BIOMD0000000524	16	Homo sapiens	('PrER', 'mGFP'), ('DISC', 'tBid'), ('p18inactive', 'PrER'), ('p18inactive', 'tBid'), ('p18inactive', 'mGFP'), ('p18inactive', 'PrNES'), ('DISC', 'mCherry'), ('PrNES', 'PrER'), ('Bid', 'mCherry'), ('tBid', 'PrNES'), ('p18inactive', 'PrER_mGFP'), ('PrNES', 'mGFP'), ('mCherry', 'PrER_mGFP'), ('PrER_mGFP', 'PrER'), ('FADD', 'PrNES'), ('tBid', 'mGFP'), ('DISC', 'PrER'), ('FADD', 'mGFP'), ('mCherry', 'mGFP'), ('FADD', 'mCherry'), ('Bid', 'PrER'), ('FADD', 'tBid'), ('DISC', 'p18inactive'), ('Bid', 'PrNES'), ('p18inactive', 'Bid'), ('tBid', 'PrER_mGFP'), ('FADD', 'PrER'), ('tBid', 'mCherry'), ('PrNES', 'PrER_mGFP'), ('FADD', 'p18inactive'), ('FADD', 'Bid'), ('Bid', 'mGFP'), ('mCherry', 'PrER'), ('PrNES', 'mCherry'), ('DISC', 'Bid'), ('DISC', 'PrNES'), ('Bid', 'tBid'), ('DISC', 'PrER_mGFP'), ('FADD', 'DISC'), ('DISC', 'mGFP'), ('p18inactive', 'mCherry'), ('PrER_mGFP', 'mGFP'), ('FADD', 'PrER_mGFP'), ('tBid', 'PrER')
BIOMD0000000050	14	n/a	('Gly', 'Fru'), ('FA', 'Mel'), ('FA', 'MG'), ('FA', 'Fru'), ('Mel', 'Fru'), ('Gly', 'Mel'), ('FA', 'Glu'), ('Gly', 'AA'), ('AA', 'Fru'), ('Glu', 'Mel'), ('MG', 'Fru'), ('AA', 'Glu'), ('FA', 'AA'), ('Glu', 'Fru'), ('Gly', 'MG'), ('Glu', 'MG'), ('AA', 'MG'), ('AA', 'Mel'), ('Mel', 'MG'), ('Gly', 'FA')
BIOMD0000000647	11	n/a	('RKIPP', 'MEKPP_ERK'), ('RKIPP', 'RKIPP_RP'), ('Raf1_RKIP', 'ERK'), ('ERK', 'RP'), ('ERK', 'MEKPP'), ('Raf1', 'MEKPP'), ('MEKPP', 'RKIPP_RP'), ('MEKPP', 'RP'), ('RKIP', 'RP'), ('ERKPP', 'RKIPP'), ('Raf1', 'RKIPP'), ('ERKPP', 'RP'), ('Raf1_RKIP_ERKPP', 'ERK'), ('Raf1_RKIP', 'MEKPP'), ('Raf1', 'ERKPP'), ('ERKPP', 'MEKPP_ERK'), ('ERK', 'MEKPP_ERK'), ('RKIP', 'RKIPP_RP'), ('Raf1', 'ERK'), ('ERKPP', 'Raf1_RKIP_ERKPP'), ('RKIP', 'ERKPP'), ('Raf1_RKIP_ERKPP', 'RKIPP_RP'), ('Raf1_RKIP_ERKPP', 'MEKPP'), ('Raf1_RKIP', 'ERKPP'), ('RP', 'RKIPP_RP'), ('MEKPP_ERK', 'RKIPP_RP'), ('ERK', 'RKIPP_RP'), ('ERKPP', 'ERK'), ('Raf1_RKIP_ERKPP', 'RP'), ('Raf1', 'Raf1_RKIP_ERKPP'), ('Raf1', 'Raf1_RKIP'), ('MEKPP', 'MEKPP_ERK'), ('Raf1_RKIP', 'RKIPP'), ('Raf1_RKIP_ERKPP', 'RKIPP'), ('Raf1_RKIP_ERKPP', 'MEKPP_ERK'), ('MEKPP_ERK', 'RP'), ('RKIP', 'MEKPP'), ('Raf1', 'RKIP'), ('RKIPP', 'RP'), ('RKIP', 'MEKPP_ERK'), ('Raf1_RKIP', 'Raf1_RKIP_ERKPP'), ('Raf1', 'RP'), ('Raf1_RKIP', 'RKIPP_RP'), ('Raf1', 'RKIPP_RP'), ('Raf1_RKIP', 'MEKPP_ERK'), ('ERKPP', 'MEKPP'), ('RKIP', 'RKIPP'), ('Raf1', 'MEKPP_ERK'), ('RKIP', 'ERK'), ('RKIPP', 'MEKPP'), ('ERKPP', 'RKIPP_RP'), ('RKIP', 'Raf1_RKIP'), ('ERK', 'RKIPP'), ('Raf1_RKIP', 'RP'), ('RKIP', 'Raf1_RKIP_ERKPP')
BIOMD0000000520	3	Rodents	('N0', 'N1'), ('N1', 'N2'), ('N0', 'N2')
BIOMD0000000523	16	Homo sapiens	('tBid', 'mCherry'), ('PrER', 'mGFP'), ('PrNES', 'mGFP'), ('DISC', 'PrNES'), ('DISC', 'tBid'), ('FADD', 'tBid'), ('p18inactive', 'tBid'), ('FADD', 'mCherry'), ('mCherry', 'mGFP'), ('p18inactive', 'mCherry'), ('tBid', 'PrNES'), ('tBid', 'mGFP'), ('FADD', 'DISC'), ('p18inactive', 'mGFP'), ('mCherry', 'PrER'), ('p18inactive', 'PrER'), ('FADD', 'p18inactive'), ('DISC', 'p18inactive'), ('DISC', 'PrER'), ('FADD', 'PrER'), ('DISC', 'mGFP'), ('PrNES', 'mCherry'), ('p18inactive', 'PrNES'), ('FADD', 'PrNES'), ('tBid', 'PrER'), ('FADD', 'mGFP'), ('DISC', 'mCherry'), ('PrNES', 'PrER')
BIOMD0000000454	8	Generic	('x3', 'y3'), ('y4', 'y2'), ('y1', 'y5'), ('x3', 'y2'), ('y4', 'y5'), ('x2', 'y4'), ('x1', 'y3'), ('x1', 'y5'), ('y2', 'y3'), ('y1', 'y3'), ('x1', 'y4'), ('x1', 'x3'), ('y4', 'y3'), ('y1', 'y2'), ('y1', 'y4'), ('y1', 'x3'), ('y5', 'y2'), ('y1', 'x1'), ('x2', 'x1'), ('x1', 'y2'), ('x2', 'x3'), ('x3', 'y4'), ('x3', 'y5'), ('y1', 'x2'), ('x2', 'y2'), ('x2', 'y5'), ('x2', 'y3'), ('y5', 'y3')
BIOMD0000000069	10	Homo sapiens	('srco', 'srca'), ('srci', 'srca'), ('srco', 'srcc'), ('srca', 'Cbp_P'), ('srcc', 'Cbp_P_CSK'), ('srci', 'Cbp_P'), ('srci', 'PTP'), ('Cbp_P', 'PTP_pY789'), ('srcc', 'PTP_pY789'), ('CSK_cytoplasm', 'PTP'), ('srci', 'PTP_pY789'), ('srca', 'PTP_pY789'), ('srcc', 'CSK_cytoplasm'), ('CSK_cytoplasm', 'Cbp_P_CSK'), ('srca', 'CSK_cytoplasm'), ('srca', 'PTP'), ('srcc', 'Cbp_P'), ('srco', 'Cbp_P_CSK'), ('CSK_cytoplasm', 'PTP_pY789'), ('srcc', 'PTP'), ('srca', 'Cbp_P_CSK'), ('srci', 'CSK_cytoplasm'), ('Cbp_P', 'Cbp_P_CSK'), ('Cbp_P', 'PTP'), ('srci', 'srco'), ('srco', 'PTP'), ('CSK_cytoplasm', 'Cbp_P'), ('srca', 'srcc'), ('srci', 'Cbp_P_CSK'), ('Cbp_P_CSK', 'PTP_pY789'), ('PTP', 'PTP_pY789'), ('Cbp_P_CSK', 'PTP'), ('srci', 'srcc'), ('srco', 'Cbp_P'), ('srco', 'PTP_pY789'), ('srco', 'CSK_cytoplasm')
BIOMD0000000455	9	Generic	('x2', 'y5'), ('y5', 'y2'), ('x2', 'y2'), ('y1', 'y3'), ('x3', 'y3'), ('y1', 'x2'), ('x2', 'x1'), ('x2', 'y4'), ('x3', 'y2'), ('x2', 'y3'), ('y1', 'y4'), ('y5', 'y3'), ('y1', 'y5'), ('x1', 'x3'), ('x1', 'y4'), ('y4', 'y5'), ('x1', 'y3'), ('y1', 'x3'), ('y4', 'y3'), ('x3', 'y4'), ('x3', 'y5'), ('x1', 'y5'), ('y1', 'x1'), ('y1', 'y2'), ('y2', 'y3'), ('y4', 'y2'), ('x2', 'x3'), ('x1', 'y2')
BIOMD0000000526	16	Homo sapiens	('FADD', 'PrER_mGFP'), ('PrNES', 'mCherry'), ('p18inactive', 'PrER_mGFP'), ('DISC', 'p18inactive'), ('FADD', 'p18inactive'), ('p18inactive', 'mCherry'), ('DISC', 'mGFP'), ('FADD', 'mGFP'), ('p18inactive', 'tBid'), ('p18inactive', 'mGFP'), ('FADD', 'DISC'), ('tBid', 'PrER'), ('DISC', 'tBid'), ('p18inactive', 'PrER'), ('p18inactive', 'PrNES'), ('PrNES', 'PrER_mGFP'), ('PrNES', 'PrER'), ('FADD', 'PrNES'), ('DISC', 'PrER'), ('DISC', 'mCherry'), ('FADD', 'tBid'), ('PrER_mGFP', 'PrER'), ('PrNES', 'mGFP'), ('FADD', 'PrER'), ('mCherry', 'mGFP'), ('DISC', 'PrNES'), ('tBid', 'PrNES'), ('mCherry', 'PrER'), ('PrER_mGFP', 'mGFP'), ('tBid', 'mCherry'), ('FADD', 'mCherry'), ('tBid', 'mGFP'), ('mCherry', 'PrER_mGFP'), ('DISC', 'PrER_mGFP'), ('PrER', 'mGFP'), ('tBid', 'PrER_mGFP')
BIOMD0000000284	9	Bacteria	('D', 'E'), ('A', 'D'), ('X', 'C'), ('B', 'E'), ('C', 'Z'), ('X', 'E'), ('A', 'Z'), ('X', 'Y'), ('C', 'D'), ('Y', 'Z'), ('E', 'F'), ('F', 'Z'), ('D', 'F'), ('B', 'F'), ('C', 'F'), ('X', 'Z'), ('A', 'Y'), ('X', 'D'), ('A', 'F'), ('B', 'Z'), ('X', 'F'), ('A', 'B'), ('A', 'C'), ('B', 'D'), ('Y', 'E'), ('B', 'C'), ('B', 'Y'), ('C', 'Y'), ('Y', 'F'), ('E', 'Z'), ('X', 'B'), ('D', 'Z'), ('D', 'Y'), ('A', 'E'), ('C', 'E'), ('X', 'A')
BIOMD0000000084	8	n/a	('Rin', 'x2'), ('x1', 'x3'), ('x2', 'x3'), ('x1', 'x2'), ('Rin', 'x3'), ('Rin', 'x1')
BIOMD0000000052	11	n/a	('Triose', 'lys_R'), ('Cn', 'lys_R'), ('Acetic_acid', 'lys_R'), ('lys_R', 'Melanoidin'), ('Triose', 'Melanoidin'), ('Cn', 'Melanoidin'), ('Formic_acid', 'lys_R'), ('C5', 'lys_R'), ('C5', 'Formic_acid'), ('Triose', 'Acetic_acid'), ('C5', 'Acetic_acid'), ('Acetic_acid', 'Melanoidin'), ('C5', 'Cn'), ('C5', 'Melanoidin'), ('Cn', 'Acetic_acid'), ('C5', 'Triose'), ('Formic_acid', 'Triose'), ('Formic_acid', 'Cn'), ('Triose', 'Cn'), ('Formic_acid', 'Melanoidin'), ('Formic_acid', 'Acetic_acid')
BIOMD0000000271	6	Rodents	('EpoR', 'dEpoe'), ('Epo_EpoRi', 'dEpoe'), ('Epo_EpoR', 'dEpoe'), ('Epo_EpoRi', 'dEpoi'), ('EpoR', 'dEpoi'), ('Epo', 'dEpoe'), ('Epo', 'dEpoi'), ('Epo_EpoR', 'dEpoi'), ('dEpoi', 'dEpoe')
BIOMD0000000461	4	Bacteria	('lacz', 'x'), ('IPTG', 'sigb'), ('sigb', 'x'), ('sigb', 'lacz'), ('IPTG', 'x'), ('IPTG', 'lacz')
BIOMD0000000525	16	Homo sapiens	('p18inactive', 'PrNES'), ('p18inactive', 'PrER'), ('tBid', 'PrNES'), ('DISC', 'PrER'), ('mCherry', 'PrER'), ('PrNES', 'mCherry'), ('DISC', 'tBid'), ('FADD', 'mCherry'), ('p18inactive', 'tBid'), ('FADD', 'PrER'), ('FADD', 'p18inactive'), ('tBid', 'mGFP'), ('DISC', 'p18inactive'), ('mCherry', 'mGFP'), ('DISC', 'PrNES'), ('PrER', 'mGFP'), ('DISC', 'mGFP'), ('FADD', 'DISC'), ('tBid', 'PrER'), ('p18inactive', 'mCherry'), ('tBid', 'mCherry'), ('PrNES', 'mGFP'), ('PrNES', 'PrER'), ('FADD', 'PrNES'), ('DISC', 'mCherry'), ('p18inactive', 'mGFP'), ('FADD', 'tBid'), ('FADD', 'mGFP')
BIOMD0000000521	4	Homo sapiens	('P', 'Q')
BIOMD0000000010	8	Amphibians	('MKK_P', 'MAPK'), ('MKK', 'MAPK_PP'), ('MKK', 'MKK_PP'), ('MKK_P', 'MKK_PP'), ('MKK', 'MAPK'), ('MKK_PP', 'MAPK_P'), ('MAPK', 'MAPK_PP'), ('MKK_PP', 'MAPK'), ('MKK', 'MAPK_P'), ('MKK_PP', 'MAPK_PP'), ('MAPK_P', 'MAPK_PP'), ('MKK_P', 'MAPK_PP'), ('MKK', 'MKK_P'), ('MKK_P', 'MAPK_P')
BIOMD0000000029	4	Amphibians	('M', 'MpT'), ('M', 'MpY'), ('MpY', 'MpT')
BIOMD0000000197	5	n/a	('x1', 'x2'), ('x1', 'x4'), ('x1', 'x3'), ('x1', 'x5'), ('x5', 'x4'), ('x3', 'x5'), ('x5', 'x2')
BIOMD0000000272	6	Rodents	('SAv_EpoR', 'SAv_EpoRi'), ('EpoR', 'SAv_EpoR'), ('EpoR', 'SAv_EpoRi')
BIOMD0000000167	7	Generic	('Pstat_nuc', 'stat_nuc'), ('stat_nuc', 'species_test'), ('Pstat_nuc', 'species_test')
BIOMD0000000262	11	Rodents	('Akt', 'S6')
BIOMD0000000240	6	Bacteria	('DegU', 'mDegU'), ('Dim', 'mDegU'), ('DegUP', 'mDegU')
BIOMD0000000037	12	Slime Mold	('Xi', 'Ya'), ('Xi', 'Pi'), ('Ya', 'Pi')
BIOMD0000000263	11	Rodents	('Akt', 'pro_TrkA'), ('S6', 'pro_TrkA'), ('Akt', 'S6')
BIOMD0000000641	5	Homo sapiens	('CellCact', 'Effectoract')
BIOMD0000000413	5	Plants	('TIR1', 'auxinTIR1'), ('auxinTIR1', 'VENUS'), ('TIR1', 'VENUS')
BIOMD0000000624	7	Homo sapiens	('APAP', 'APAPconj_Glu')
BIOMD0000000945	5	Homo sapiens	('L_m', 'H_m')
BIOMD0000000459	4	Bacteria	('IPTG', 'sigb')

Table S1: List of biological networks from Biomodels used in this study. The resulting database includes 30 biological networks (one row per network) and a total of 432 systems, which is defined as a (GRN model, behavior space (Z)) tuple and where the pairs of observed nodes (used as behavior spaces) per network are given in the last column. Please note that we let the following chemical species evolve in time, despite being labeled as "constant" or "boundary" species in the original SBML files: model #262: ['EGF', 'pro_EGFR'], model #284: ['X', 'Y', 'Z'], model #37: ['Gluc'], model #263: ['NGF', 'pro_TrkA'], model #454: ['y1', 'y2', 'y3', 'y4', 'y5'], model #455: ['y1', 'y2', 'y3', 'y4', 'y5', 'y6'], model #459: ['IPTG'], model #461: ['IPTG'], model #624: ['X1'].

Figure S2: Additional results complementing Figure 8 of the main paper. This figure shows the resulting trajectories after applying the discovered stimuli-based intervention (shown in Figure 8-b) to the example RKIP-ERK signaling pathway [56] for the 6 "disease" trajectories originally discovered in the behavioral catalog (shown in Figure 8-a). (a) For each trajectory (one per row), we see that the intervention successfully re-sets the disease setpoint (startpoint of the trajectory shown in red in the orange region) to a healthy set-point (endpoint of the trajectory shown in cyan in the green region). (b-c) Similar results are achieved despite adding push perturbations (b) or wall perturbations (c) in addition to the stimuli-based intervention.

Figure S3: Wall implementation. Walls are implemented within the 2D space spanned by the 2 observed nodes. Within that space, we can interpret the node activation levels $y(0),\cdots,y(t)$ as the trajectory of a particle moving. In order to simulate the interaction with "walls" in that space, several implementations could be envisaged. Within the accompanying software AutoDiscJax we propose two possible variants: perfectly elastic collision (equivalent to a discontinuous force field) and some continuous force field variant. The second variant (continuous force field) is employed for the main results of this paper. (a) For the first variant, we consider a perfectly elastic collision without loss of kinetic energy. In that case, when the trajectory is touching the wall at position $p$ with speed $v = v_{\bot} + v_{\parallel}$ we deviate the trajectory in such a way that is "bouncing" against the wall such that $v_{\parallel}$ is unchanged and $v_{\bot} \leftarrow - v_{\bot}$. To implement it, we simply check whether the segment $\left\lbrack y(t),y(t + \Delta t) \right\rbrack$ intersect the wall at each time step. It it does, we compute the intersection point $p$ and time $t_{1}$, and set the activation level $y(t + \Delta t)$ to $p + (\Delta t - t1) \cdot \left\lbrack - v_{\bot} + v_{\parallel} \right\rbrack$. (b) For the second variant, we implement walls as energy barriers acting as a new force field in the environment, constraining the GRN traversal of the space. This time, instead of having a discontinuous effect on the perpendicular speed $v_{\bot}$ we define a wall force $f_{\bot} = \pm \alpha v_{\bot}$ (+ if $v_{\bot}$ is going toward wall, - otherwise ) and use it to update the perpendicular component of the trajectory speed as $v_{\bot} \leftarrow v_{\bot} + f_{\bot} \cdot \Delta T$. Here $\alpha \in \lbrack 0, - 2\rbrack$ and is calculated as a function of the distance between the point and the wall. As illustrated in the figure, this basically results in a stadium-shaped force field around the wall.

Figure S4: Sanity Check: Discoveries made by our exploration pipeline on simple systems for which the ground truth goal states are known (a-b) as well as for more complex variants of these systems (c-d). Time-course evolution of the discovered trajectories are shown for the nodes given in x and y axis. Trajectories start from the initial conditions (red points) sampled by the exploration method and end in the discovered steady states (cyan points). (a) Discoveries made by our curiosity search method for BIOMD0000000341 which models the dynamics of beta-cell mass, insulin and glucose dynamics as described in Topp et al.'s paper [149]. This simple model is known to have two attractors, representing physiological (B=300, I=10, G=100) and pathological (B=0, I=0, G=600) steady states, as described in the paper. We show here that our exploration pipeline is able to re-discover those two steady states. (b) Discoveries made by our curiosity search method for BIOMD0000000454 with the same assumptions as described in the "Example One" from the "Metabolic Control Analysis: Rereading Reder" paper [150], which is a small metabolic pathway whose (single) steady state can be derived analytically as shown in the paper. As expected our exploration method is again able to re-discover this steady state (x1=0.056, x2=0.769, x3=4.231). (c-d) Discoveries made by our curiosity search method and a random search method for BIOMD0000000454 but this time relaxing some of the assumptions made in "Example One" of [150]. We notably relax the assumption that metabolite concentrations $y_1(t), ... y_5(t)$ are fixed to certain values, and instead let these quantities evolve according to the provided chemical reactions. In that more complex case, the system exhibits a spectrum of possible steady states that are not analytically easy to find and unknown a priori. Here we see that the curiosity search method (c) is able to uncover the distribution of possible steady states more efficiently than a random search method (d).

References