AI-driven automated discovery tools reveal diverse behavioral competencies of biological networks

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 (Bell and Gilan, 2020).

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 (Kitano, 2007b; Félix and Barkoulas, 2015), 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 (von Dassow et al., 2000; Ingolia, 2004; Ma et al., 2006; Noman et al., 2015; Rizk et al., 2009; Donzé et al., 2011). Here instead, we propose to conduct a battery of empirical tests that draw inspiration from classical ‘displacement experiments’ (Walcott, 1996; Luschi et al., 2001) and ‘barrier experiments’ (Bisch-Knaden and Wehner, 2001) 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) (Abramson, 1994), robustness studies in systems biology often neglect those aspects. We propose that a behaviorist lens on robustness can help in understanding forms of non-genetic resistance in transcriptional space, which is crucial for the development of therapeutic strategies (Bell and Gilan, 2020).

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×r perturbations conditioned on the GRN goal-reaching trajectory i→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 1c 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={ik,ok,zk,{up,op,zp,p=1...P},k=1…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 {up,p=1...P} as the average distance in behavior space Z between the original trajectory endpoint z and the perturbed trajectories endpoints {zp}. 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 5e) or walls (Figure 5g) on the way.

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 (Li and Wang, 2014; Li and Wang, 2013). 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 (Bell and Gilan, 2020; Venkatachalapathy et al., 2021). A recent landscape-generation method only proposes to approximate the pseudopotential energy through simulation trajectories obtained throughout exploration of the system (Venkatachalapathy et al., 2021), 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->Z mapping, this can lead to poor estimation of the overall energy landscape (Figure 6a). 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 6b). 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 6c). For the example system RKIP-ERK pathway (Kwang-Hyun et al., 2003), 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 (Lee et al., 2006).

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.

The 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 (Levin, 2022): 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.

Second, 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 Biswas et al., 2023, 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 (Reinitz and Sharp, 1995; Jaeger et al., 2004; Cotterell and Sharpe, 2010; Molinelli et al., 2013), 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 7b, 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 an unicellular organism well known for its diverse capabilities, especially navigational ones (Vallverdú et al., 2018; Beekman and Latty, 2015; Saigusa et al., 2008; Nakagaki and Guy, 2008).

Understanding forms of non-genetic resistance and non-genetic heterogeneity is crucial across a wide range of cancer and treatment contexts (Bell and Gilan, 2020). 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 (Kwang-Hyun et al., 2003). 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 (Lee et al., 2006). In Figure 8a, 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 6c. 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 gene regulatory network (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 s for 100 s), 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 8b) that successfully reset the setpoints (as shown in Figure 8c) under various tested perturbations (as shown in Figure 8—figure supplement 1). 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.

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 (Volk et al., 2020). 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 (Colas et al., 2018), by effectively exploring solution space.

Here, we consider the target application of oscillator circuit engineering followed in Hiscock, 2019, 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 y0). We then compare several exploration strategies: a random search, two optimization-driven strategies, one using gradient descent as proposed in Hiscock, 2019 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 Hiscock, 2019. 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∼UI and use Adam optimizer for n=5000 optimization steps. For the evolutionary algorithm, we use the CMA-ES algorithm (Hansen et al., 2003), 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,ω,b space (Figure 9a-c), thus discovering oscillators close to the target one (Figure 9d). 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 9h-f). 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 9i-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.

This paper presents a novel framework aimed at uncovering the navigation competencies of 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 (Baranes and Oudeyer, 2013; Forestier et al., 2022) - can be transposed to assist efficient discovery of behavioral abilities within biological pathway models like GRNs. Second, rather than merely mapping the attractor states (Kauffman, 1993; Kauffman, 1995; Dang et al., 2011) or analyzing their sensitivity to model parameter changes (Ingalls, 2004; Ingalls, 2008) 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 on 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 (Pandi et al., 2022), 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 (Reinke et al., 2020; Etcheverry et al., 2020). 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 (Libby et al., 2019; Pietak and Levin, 2016; Pietak and Levin, 2016).

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 (Koseska and Bastiaens, 2017; Baluška et al., 2022; Baluška et al., 2023; Reber and Baluška, 2021; Baluška and Reber, 2021). 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 (Levin, 2023b; Bernheim-Groswasser et al., 2018).

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 (Scheiner, 2019). 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 (copy archived at Etcheverry et al., 2023a) 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 yt=0, the desired rollout length Tsecs and step size Δ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 dydt using JAX odeint solver to update model species yt→yt+ΔT, (2) calls the model assignment rules to update kinematic parameters if needed, and (3) apply the intervention and/or perturbation function to update yt+ΔT,wt+ΔT,c accordingly. In this paper, we use T=2500s and ΔT=0.1 (25,001 time points per rollout including t0). 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.

In our computational models, we are able to record the activities of all nodes during a model rollout. The observation space O⊂R+n×TΔT is such that o=(y(0),…,y(T)) 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⊂R+n consists of initial node states sampled from the hyper-rectangle y0,min,y0,max where y0,min=1r×yd,min and y0,max=r×yd,max with r=20 and yd,min,yd,max 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⊂R+2 endpoint states z=yiT,yjT 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.