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

All the ODE models we use in this work are downloaded from the BioModels database (Glont, 2018, Malik-Sheriff et al., 2020) in SBML format. From all models referenced on the website, we only consider the ones that are curated, that have at least three 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 Etcheverry et al., 2023b). 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 Supplementary file 1). 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 s and ΔT=0.1). We discarded simulation results containing invalid values (NaN or negative concentrations) or those that took an excessive amount of time (>1 s). 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 s). 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 (20%) exhibited invalid concentrations (NaN or negative) or unsettled behaviors (∃t≥2400suchthat|y(t)−y(T)|≥0.02×|y(T)−y(0)|) or periodic patterns (∃f>0suchthat|S(f)|≥40whereS=DFT([y(T2),⋯,y(T)])); or [filter F2] the reached space in Z was too small ((maxk=1⋯50yk(T)−mink=1⋯5yk(T))<0.1) to discard cases where ‘diversity’ could result from floating point rounding errors; or [filter F3] the number of attractors was less than four {ykT}k=1⋯50cover≤4binsa20×20binningofZ.

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

Following the methodology proposed in Biswas et al., 2023, 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 Biswas et al., 2023. To achieve this, we initially generated 300 networks based on the transcriptional gene circuit model (Reinitz and Sharp, 1995), 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,τ of these networks were randomized W∼−30,30n×n,B∼−10,10n,τ∼1,15 where model step is defined as yt+1=ΔTτ×sigmoidWy+B+1−ΔTτ×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.

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’ (Topp et al., 2000) and BIOMD0000000454 as described in the ‘Example One’ of the ‘Metabolic Control Analysis: Rereading Reder’ paper (Smallbone, 2013). 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 3—figure supplement 1 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 y1(t),...y5(t) evolve in time. We also illustrate results of our curiosity-driven exploration method and of a random search method on this more complex variant in the figure.

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 1c. 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∼N0,0.1∗y0,max−y0,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 (copy archived at Etcheverry et al., 2023a) that accompanies this paper can be used to implement this and other IMGEP variants in JAX.

We define three family of perturbations: (1) the noise perturbation Unσn,pn∨y 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 Upmp,np∨y parametrized by its magnitude (proportional to the extent of y) and number of occurrences; (3) the wall perturbation Uwlw,nw∨y 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 Figure 5—figure supplement 1.

To assess the robustness of the GRN systems in our database, we employ an evaluation procedure, as depicted in Figure 1d. 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 {yk,k=1..K} to s=18 different perturbation distributions, each representing various levels of difficulty:(σn,pn)∈{(0.001,5),(0.005,5),(0.1,5),(0.005,10),(0.005,5),(0.005,1)}, (mp,np)∈{(0.05,1),(0.1,1),(0.15,1),(1,0.1),(2,0.1),(3,0.1)}, (lw,nw)∈{(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∗r perturbations. For each perturbation in the set {up,p=1...P}, we re-run the trajectory starting from the same initial state i but with the sampled perturbation applied i,up, and observe the resulting outcome op and reached endpoint zp.

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

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 (Reinke et al., 2020; Etcheverry et al., 2020; Pugh et al., 2015) or distance-based metric from ecology (Scheiner, 2019) 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 (Benureau, 2015). Another approach, called the threshold coverage, measures diversity as the volume of the union of the set of hyperballs of radius ϵ that have for centers the observed effects {z∈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 (Benureau, 2015).

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

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 (Kitano, 2007b) or phenotypic trait (Félix and Barkoulas, 2015) 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∈Z discovered by curiosity search, (2) the set of perturbation {up} is defined in previous section and conditioned on the GRN goal-reaching trajectory i→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.

In this study, we introduce the terms ‘diversity’ and ‘versatility’ to characterize the competencies of the exploration agent (IMGEP) and the GRN agent, 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×πϵ21+2ϵ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.

This section details the additional experiments conducted on the RKIP-ERK signaling pathway (Lehman and Stanley, 2011). We refer to the accompanying notebook tutorial for reproducing these experiments: https://developmentalsystems.org/curious-exploration-of-grn-competencies/tuto1.html.

For Figure 4, clustering in behavior space was performed using the HDBSCAN algorithm (McInnes et al., 2017) 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 (Pugh et al., 2015; Gillies, 2022).

In Figure 6, we generated trajectory-based energy landscapes following the method proposed in Venkatachalapathy et al., 2021. 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, b, c 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×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=yMEKPP1,⋯,yMEKPP10 every 10 s for 100 s. 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.

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 Hiscock, 2019, where parameters of a gene circuit model are optimized to produce target oscillator patterns. The gene circuit model employed in Hiscock, 2019 is the same than the one used for the random networks database (Eq 1), with τ=1. Hence the intervention space is a n2+2n dimensional space defined as I=yt=0,min,yt=0,max⊕Wmin,Wmax⊕Bmin,Bmax, with y0,min=0,y0,max=1,Wmin=−30,Wmax=30,Bmin=−10,Bmax=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 y0 that would produce a periodic pattern y=yn=00,⋯,yn=0T with target amplitude A, frequency w and offset b. Here, the target A,ω,b are sample randomly with A∼U0.1,0.5,b∼UA,1−A,ω∼Betaα=2,β=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 Hiscock, 2019, 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 1c, 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, for a set of parameters i=W,B,y0, 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,ω,b : L=∑t(y(t)−(Acos(2πωt)+b))2. We then sample a random parameter i∼UI and use Adam optimizer with lr=10−3,b1=0.02,b3=0.001,ϵ=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.