Systems and methods for using ai-driven automated discovery tools to reveal diverse behavioral competencies of biological networks

Abstract

Disclosed are apparatus and methods for identifying a behavioral catalog of a gene regulatory network. The disclosed apparatus and methods use a curiosity driven machine learning algorithm to identify goal behavioral states and determine how controlled sources of variation connect to goal behavioral states.

Description

Tufts Ref. T002751 Quarles Ref.166118.01450 SYSTEMS AND METHODS FOR USING AI-DRIVEN AUTOMATED DISCOVERY TOOLS TO REVEAL DIVERSE BEHAVIORAL COMPETENCIES OF BIOLOGICAL NETWORKS CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] The present application is based on and claims priority from U.S. Patent Application Serial No.63 / 582,442, filed on September 13, 2023, and U.S. Patent Application Serial No.63 / 563,222, filed on March 8, 2024. The entire disclosures of each are incorporated herein by reference. STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT

[0002] N / A. FIELD OF THE INVENTION

[0003] The field of the invention relates to machine learning methods to investigate biological networks. More specifically, the patent disclosure relates to apparatus and methods for using a curiosity-driven algorithm to uncover details of gene regulatory networks, and approaches to manipulate gene regulatory networks. The disclosed methods include constructing behavioral catalogs of gene regulatory networks and using behavioral catalogs for therapeutic development and synthetic biology. BACKGROUND

[0004] Developing methods to recognize, map, predict, and control the complex, context- sensitive behavior of chemical and genetic networks is an important frontier of research in science and engineering. These systems, such as gene regulatory networks (GRNs) and protein pathways, are known to be instructive drivers of embryogenesis, cell behavior, and complex physiology. Understanding the control properties of these systems is important not only for the study of evolutionary developmental biology, but also for comprehending and intervening inTufts Ref. T002751 Quarles Ref.166118.01450 various disease states, including cancer, and for the construction of novel synthetic biologicals in bioengineering contexts.

[0005] Nevertheless, significant challenges remain 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. 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. These span many orders of spatio- temporal scale, from molecular assemblies to swarms of complex organisms. SUMMARY

[0006] Accordingly, disclosed herein are apparatus, systems, and methods for identifying a behavioral catalog of a gene regulatory network (GRN). The disclosed apparatus, systems, and methods include or utilize curiosity-driven machine learning algorithms to identify goal behavioral states and determine the versatility and robustness of said states.

[0007] In one aspect, a method of identifying a behavioral catalogue of a gene regulatory network (GRN) is provided, where the behavioral catalog includes at least one goal behavioral state and a plurality of trajectories connecting controlled sources of variation in the GRN and each of the at least one goal behavioral state. The method includes: obtaining a data set including a mathematical model of the GRN; using a curiosity driven machine learning algorithm to identify at least one goal behavioral state and a plurality of corresponding trajectories connecting controlled sources of variation in the GRN and each of the at least one goal behavioral state; testing a robustness of the at least one goal behavioral state by perturbing each of the plurality of corresponding trajectories and evaluating, for each of the plurality of corresponding trajectories, whether the respective goal behavioral state is reached; selecting a subset of the goal behavioralTufts Ref. T002751 Quarles Ref.166118.01450 states which are robust to perturbations; and saving the subset of goal behavioral states and the corresponding plurality of trajectories.

[0008] In another aspect, a system for identifying a behavioral catalogue of a gene regulatory network (GRN) is provided, where the behavioral catalog includes at least one goal behavioral state and a plurality of trajectories connecting controlled sources of variation in the GRN and each of the at least one goal behavioral state. The system includes at least one hardware processor that is programmed to perform one or more of the following steps: obtaining a data set including a mathematical model of the GRN; executing a curiosity driven machine learning algorithm to identify at least one goal behavioral state and a plurality of corresponding trajectories connecting controlled sources of variation in the GRN and each of the at least one goal behavioral state; testing a robustness of the at least one goal behavioral state by perturbing each of the plurality of corresponding trajectories and evaluating, for each of the plurality of corresponding trajectories, whether the respective goal behavioral state is reached; selecting a subset of the goal behavioral states which are robust to perturbations; saving the subset of goal behavioral states and the corresponding plurality of trajectories; and outputting the subset of goal behavioral states and the corresponding plurality of trajectories to a user.

[0009] In yet another aspect, a method of identifying a therapeutic intervention is provided. The method includes: obtaining a behavioral catalog of a gene regulatory network (GRN) associated with a disease, wherein the behavioral catalog includes a plurality of behavioral end states and a plurality of trajectories connecting controlled sources of variation in the GRN and each behavioral end state in the plurality of behavioral end states; selecting a disease associated with the GRN; classifying each behavioral end state in the plurality of behavioral end states as a first setpoint, wherein the first setpoint can be a disease setpoint or the first setpoint can be a healthy setpoint, resulting in a plurality of first setpoints, wherein a disease setpoint is associated with a phenotype or gene expression pattern associated with the disease, and wherein the healthy setpoint is not associated with a phenotype or gene expression pattern associated with the disease; defining a plurality of stimulus-based interventions to apply to each disease setpoint; introducing each stimulus-based intervention in the plurality of stimulus-based interventions to each first setpoint in the plurality of first setpoints and identifying a plurality of trajectories, wherein each trajectory in the plurality of trajectories corresponds to a response of aTufts Ref. T002751 Quarles Ref.166118.01450 first setpoint in the plurality of first setpoints to a stimulus-based intervention in the plurality of stimulus-based interventions, and determining a second setpoint of each trajectory in the plurality of trajectories; selecting a subset of trajectories in which the first setpoint is a disease state and the second setpoint after introducing the stimulus-based interventions is a healthy state; determining the subset of stimulus-based interventions that corresponds with the subset of trajectories in which the first setpoint is a disease state, and the second setpoint is a healthy state; and outputting the subset of stimulus-based interventions.

[0010] In still another aspect, a system for identifying a therapeutic intervention is provided. The system includes at least one hardware processor that is programmed to perform one or more of the following steps: obtaining a behavioral catalog of a gene regulatory network (GRN) associated with a disease, wherein the behavioral catalog includes a plurality of behavioral end states and a plurality of trajectories connecting controlled sources of variation in the GRN and each behavioral end state in the plurality of behavioral end states; selecting a disease associated with the GRN; classifying each behavioral end state in the plurality of behavioral end states as a first setpoint, wherein the first setpoint can be a disease setpoint or the first setpoint can be a healthy setpoint, resulting in a plurality of first setpoints, wherein a disease setpoint is associated with a phenotype or gene expression pattern associated with the disease, and wherein the healthy setpoint is not associated with a phenotype or gene expression pattern associated with the disease; defining a plurality of stimulus-based interventions to apply to each disease setpoint; introducing each stimulus-based intervention in the plurality of stimulus-based interventions to each first setpoint in the plurality of first setpoints and identifying a plurality of trajectories, wherein each trajectory in the plurality of trajectories corresponds to a response of a first setpoint in the plurality of first setpoints to a stimulus-based intervention in the plurality of stimulus-based interventions, and determining a second setpoint of each trajectory in the plurality of trajectories; selecting a subset of trajectories in which the first setpoint is a disease state and the second setpoint after introducing the stimulus-based interventions is a healthy state; determining the subset of stimulus-based interventions that corresponds with the subset of trajectories in which the first setpoint is a disease state, and the second setpoint is a healthy state; and outputting the subset of stimulus-based interventions.Tufts Ref. T002751 Quarles Ref.166118.01450

[0011] In yet another aspect, a method of identifying at least one of a configuration or a parameter of a gene regulatory network (GRN) to perform a desired function is provided. The method includes: obtaining a data set including a mathematical model of a GRN; using a curiosity driven machine learning algorithm to identify a plurality of goal behavioral states and a plurality of corresponding trajectories connecting controlled sources of variation in the GRN and each of the plurality of goal behavioral states; and optionally using a local optimization strategy to identify an additional goal behavioral state and add the additional goal behavioral state to the plurality of goal behavioral states.

[0012] In still another aspect, a system for identifying at least one of a configuration or a parameter of a gene regulatory network (GRN) to perform a desired function is provided, where the system includes at least one hardware processor that is programmed to perform one or more of the following steps: obtaining a data set including a mathematical model of a GRN; using a curiosity driven machine learning algorithm to identify a plurality of goal behavioral states and a plurality of corresponding trajectories connecting controlled sources of variation in the GRN and each of the plurality of goal behavioral states; and optionally using a local optimization strategy to identify an additional goal behavioral state and add the additional goal behavioral state to the plurality of goal behavioral states.

[0013] In yet another aspect, a method of priming a biological system for therapeutic treatment is provided. The method includes administering a compound to a biological system to move the biological system to a desired state. The desired state may be a state the biological system has not yet achieved (e.g., a state that is more differentiated than the biological system’s original state).

[0014] In still another aspect, a method of administering a therapeutic treatment to a biological system is provided. The method includes administering a compound to the biological system to move the biological system to a desired state. The desired state may be a more differentiated state of the biological system (e.g., more differentiated state than the biological system’s original state).Tufts Ref. T002751 Quarles Ref.166118.01450 BRIEF DESCRIPTION OF THE DRAWINGS

[0015] The disclosure will hereafter be described with reference to the accompanying drawings, wherein like reference numerals denote like elements.

[0016] The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee.

[0017] FIGS.1A-1E show an overview of the proposed framework. (FIG.1A) MOTIVATION: Functionalist approach. 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. 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 from. (FIG.1B) EXPERIMENTAL DESIGNS: Investigating GRN navigation in problem spaces. 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. (FIG.1C) AUTOMATED EXPERIMENTATION: Curious exploration of behavioral abilities. 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. (FIG.1D) EMPIRICAL TESTS: Assess navigation competencies. 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). (FIG.1E) PERSPECTIVES: Behavioral catalog reuse. We explore several potential reuses of the discovered “behavioral catalog” and proposed framework across evolutionary biology, biomedicine and bioengineering contexts.

[0018] FIGS.2A-2F provide an illustration of the experimental setup and chosen problem spaces on an example GRN model which has 10 nodes and which models the influence of RKIP on the ERK Signaling Pathway. (FIG.2A) Time-course evolution of the different nodesTufts Ref. T002751 Quarles Ref.166118.01450 y1, …, y10 (one color per node) when starting from the default initial conditions. The observation captures the states taken through time o=[y(t=0), …, y(t=T)] where y=[y1, …, y10]. (FIG.2B) 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 ^^ = [ ^^ERK( ^^), ^^RKIPRP( ^^)], where T is chosen big enough to ensure convergence. (FIG.2C) The intervention is setting the initial state of the system trajectory (for all nodes): i = [y1(t=0), …, y10(t=0)]. (FIG.2D-2E) Example of perturbations used in this paper. (FIG.2D) Noise perturbation, here applied to all 10 nodes every 5 secs until t=80 secs. (FIG.2E) Push perturbation, here applied to the two target nodes (ERK, RKIPP_RP) at t=3 seconds. (FIG.2F) Wall perturbation, also applied to the two target nodes (ERK, RKIPP_RP), here at 10% and 90% of the total distance traveled.

[0019] FIGS.3A-3D demonstrate how curiosity search uncovers a wide spectrum of reachable states in behavior space Z. (FIG.3A) 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 ^^ that have for centers the discovered endpoints { ^^ ∈ ^^} as depicted by the shaded area in (b-c) with ^^ = 0.05. (a-left) Mean and std- deviation curves of the diversity of behaviors discovered throughout exploration (with random search run twice longer until 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. (FIG.3B) 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). (FIGS.3C-3D) Discovered endpoints at the end of exploration (n=450) by random search (pink) and curiosity search (blue) for 6 example systems of our database. (FIG.3C) 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. (FIG.3D) Examples of systems with low-redundancy mapping I -> Z such that random search in ^^ is already quite efficient in covering behavior space Z, and curiosity search performs equivalently.

[0020] FIGS.4A-4B provide an illustration of the non-linearity and redundancy of the I- >Z mapping, and of the interest of using goal directed exploration strategies. Plot shows theTufts Ref. T002751 Quarles Ref.166118.01450 reachable points discovered by curiosity search (FIG.4A) and by random search (FIG.4B) in the behavior space Z and their corresponding starting points in the intervention space I, for the RKIP-ERK signaling pathway system. The intervention space is 10-dimensional, so what we show here is a 2D TSNE reduction. We apply HDBSCAN clustering on the points discovered in Z, which led to 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. (FIG.4A) Looking at the curiosity search discoveries, we can see the non-linearity of the I->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). (FIG.4B) 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 led to more novelty in Z.

[0021] FIGS.5A-5H show identification of robust traversal strategies in transcriptional space. (FIG.5A-5B) Robustness test statistics. 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. (FIG.5B) Violin plots on the right detail the median sensitivity for the 18 sub-families. (FIGS.5C-5H) 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 (FIG.5A). (FIG.5C) shows perturbed trajectories of the (model #10, nodes (3,7)) system which has the highest sensitivity to noise, whereas (FIG.5H) 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 ( ^^, ^^), and one ^^ per sub-family is shown per column. (FIG.5C) the trajectories are perturbed with the different sub-families of noise ( ^^n ∈ [0.001,0.005,0.1], ^^n ∈ [10,5,1]) which can be seen as various levels of difficulty. For each trajectory we annotate the starting position (A), endpoint prior perturbation (B), and endpointTufts Ref. T002751 Quarles Ref.166118.01450 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). (FIG.5C) 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. (FIG.5D) The (model #52, nodes (4,7)) system however, except for rare cases (trajectory #35), can consistently reach its original target despite encountering various amounts of noise. Trajectories #36 and #40 consistently follows a complex up->right-down->left path, despite the induced noise. (FIG.5E) The (model #647, nodes (2,10)) system, except for few cases (trajectory #0), is typically deviated from its original trajectory when being pushed away, and it seems to follow similar (parallel) trajectories. (FIG.5F) 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 gets back to its original trajectory. The trajectories of this system are relatively complex with several loops and detours. (FIG.5G) The (model #84, nodes (4,6)) system is not very robust to walls, and typically deviates or blocked when it encounters a wall. (FIG.5H) 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 gets back to its original path. The trajectories of this system are relatively complex with several loops and detours.

[0022] FIGS.6A-6C provide an energy landscape visualization based on the trajectory- based landscape generation method, and constructed from different set of GRN trajectories, respectively trajectories generated (FIG.6A) by the random search exploration, (FIG.6B) by the curiosity-driven exploration, and (FIG.6C) by the robustness tests experiments.

[0023] FIGS.7A-7B show 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. (FIG.7A) 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 (FIG.7B) Trade-off (aka Pareto) mean and std-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,Tufts Ref. T002751 Quarles Ref.166118.01450 versatility (y-value) is measured as the area covered by the set of robustly achieved goal states, where the criteria of goal-achievement is a binary and which tests whether the goal-reaching sensitivity (on average overall wall perturbations) is below a certain threshold (x-values). Violin plots in (FIG.7A) are ordered in ascending order according to the class mean y-value at x=0.4 in (FIG.7B).

[0024] FIGS.8A-8C show identification of stimuli-based stepwise intervention triggering robust re-set of disease states into healthy physiological states. (FIG.8A) 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. Most trajectories converge toward attractors in the “disease” region (orange). (FIG.8B) Discovered stepwise stimuli intervention on MEKPP which we apply on states stuck in the “disease” region during a 100 second period. (FIG.8C) 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.

[0025] FIGS.9A-9I provide a comparison of three alternative strategies for the design of oscillator circuits: curiosity search (blue), random search (pink), and gradient descent (orange). (FIGS.9A-9C) 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 gradient descent does not discover any (starting from a single random initialization). (FIG.9A) 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. (FIG.9B) Box plots projecting points from the 3D scatter plot into the respective (amplitude, main frequency, offset) axes. (FIG.9C) Diversity discovered throughout exploration, where diversity is measured with a binning-based space coverage metric (20 bins per dimension). (FIG. 9D) Evolution of the training loss ^^ for the three exploration strategies. Corresponding best discoveries (for which ^^ is minimal) for (FIG.9E) curiosity search, (FIG.9F) random search, and (FIG.9G) gradient descent. Local training loss and resulting finetuning of the best discoveries of (FIG.9H) curiosity search and (FIG.9I) random search with gradient descent.

[0026] FIGS.10A-10C provide examples of interventions that can be implemented within the accompanying AutodiscJax software. (FIG.10A) Simulation in batch mode.Tufts Ref. T002751 Quarles Ref.166118.01450 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. (FIG.10B) Example drug intervention. Intervention where species amounts are clamped to specific values. In this example the node MEKPP is clamped to 2.5 ^^ ^^ for 10 seconds at t=0 and then to 1 ^^ ^^ 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’. (FIG.10C) Example genome intervention. 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.

[0027] FIGS.11A-11F show the resulting trajectories after applying the discovered stimuli-based intervention to the example RKIP-ERK signaling pathway for the 6 “disease” trajectories originally discovered in the behavioral catalog. Each figure (FIGS.11A-11F) shows 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). Left: without perturbation; center: push perturbation; right: wall perturbation.

[0028] FIGS.12A-12B show 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 ^^(0),⋯, ^^( ^^) 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. (FIG.12A) For the first variant, we consider a perfectly elastic collision without loss of kinetic energy. In that case, whenTufts Ref. T002751 Quarles Ref.166118.01450 the trajectory is touching the wall at position ^^ with speed ^^ = ^^+ ^^∥ we deviate the trajectory in such a way that is “bouncing” against the wall such that ^^∥ is unchanged and ^^ / ← − ^^ / . To implement it, we simply check whether the segment [ ^^( ^^), ^^( ^^ + ^^ ^^)] intersect the wall at each time step. If it does, we compute the intersection point ^^ and time ^^", and set the activation level ^^( ^^ + ^^ ^^) to ^^ + ( ^^ ^^ − ^^1) ⋅ [− ^^⊥ + ^^∥]. (FIG.12B) 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 ^^ / we define a wall force ^^⊥ = ± ^^ ^^⊥ (+ if ^^⊥ is going toward wall, - otherwise) and use it to update the perpendicular component of the trajectory speed as ^^⊥ ← ^^⊥ + ^^⊥ ⋅ ^^ ^^. Here ^^ ∈ [0, −2] and is calculated as a function of the distance between the point and the wall. This basically results in a stadium-shaped force field around the wall.

[0029] FIGS.13A-13D provide a list of biological networks from Biomodels used in this study. The resulting database included 32 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.

[0030] FIG.14 shows an example of a process for identifying a behavioral catalog in accordance with some embodiments of the disclosed subject matter.

[0031] FIG.15 shows an example of a process for identifying a therapeutic intervention in accordance with some embodiments of the disclosed subject matter.

[0032] FIG.16 shows an example of a process for classifying healthy setpoints and disease setpoints in accordance with some embodiments of the disclosed subject matter.

[0033] FIG.17 shows an example of a process for identifying a configuration of a gene regulatory network to perform a desired function in accordance with some embodiments of the disclosed subject matter.

[0034] FIG.18 shows an example of a system for identifying a behavioral catalog, identifying therapeutic interventions, or identifying a configuration of gene regulatory network to perform a desired function in accordance with some embodiments of the disclosed subject matter.Tufts Ref. T002751 Quarles Ref.166118.01450 DETAILED DESCRIPTION

[0035] The disclosures of these patents, patent applications, and publications in their entireties are hereby incorporated by reference into this application in order to more fully describe the state of the art as known to those skilled therein as of the date of the invention described and claimed herein. The instant disclosure will govern in the instance that there is any inconsistency between the patents, patent applications, and publications and this disclosure.

[0036] Before any embodiments of the disclosure are explained in detail, it is to be understood that the disclosure is not limited in its application to the details of construction and the arrangement of components set forth in the following description or illustrated in the following drawings. The disclosure is capable of other embodiments and of being practiced or of being carried out in various ways.

[0037] It is to be understood that the phraseology and terminology used herein is for the purpose of description and should not be regarded as limiting. The use of “including,” “comprising,” or “having” and variations thereof herein is meant to encompass the items listed thereafter and equivalents thereof as well as additional items. Unless specified or limited otherwise, the terms “mounted,” “connected,” “supported,” and “coupled” and variations thereof are used broadly and encompass both direct and indirect mountings, connections, supports, and couplings. Further, “connected” and “coupled” are not restricted to physical or mechanical connections or couplings.

[0038] In one aspect, the disclosure provides methods and systems of identifying a behavioral catalog of a gene regulatory network (GRN). The behavioral catalog can include a library of all accessible states a GRN can achieve, characterization of how sensitive each state is, and / or a record of the history of experiments completed involving the GRN. In some embodiments, a data set of the GRN includes a mathematical model (e.g., a mathematical model of a biological and / or biomedical system such as an ODE model from the Biomodels website). In some embodiments, an algorithm is used to identify at least one goal behavioral state (e.g., a steady-state or attractor, or a stable state towards which a system tends to evolve over time) and one or more trajectories (e.g., the set or sequence of states exhibited by a system when starting from an initial condition) that connect controlled sources of variation in the GRN and each of theTufts Ref. T002751 Quarles Ref.166118.01450 goal behavioral states. In some embodiments, the algorithm uses a protocol based on machine learning, artificial intelligence, autonomous vehicle navigation, behavioral science, and / or artificial life. In some embodiments, a curiosity driven machine learning algorithm (e.g., an intrinsically-motivated goal exploration process, or IMGEP, algorithm) is used. The protocol the algorithm follows can be based on active interference framework, perceptual control theory, training for habituation, training for sensitization, training for anticipation, training for instrumental learning, training for Pavlovian conditioning, training for planning, predictive coding, principle of empowerment, or casual information theory. Throughout the specification, “an algorithm” can be used in place of “a curiosity driven algorithm.”

[0039] In some embodiments, the robustness (e.g., the ability to reach a goal state despite various perturbations) of the goal behavioral states is tested by perturbing the trajectories that connect controlled sources of variation and a corresponding goal behavioral state and determining whether the corresponding goal behavioral state is reached. In some embodiments, the versatility (e.g., the ability to reach diverse goal states under various interventions) of the goal behavioral state is tested by intervening in the system and determining the diversity of goal states that can be reached. In general, interventions refer to different initial states of functionalities of a model rollout before an experiment. In contrast, perturbations refer to stochastic events that can occur during the experiment. In some embodiments, a subset of the goal behavioral states and their corresponding trajectories which are robust to perturbation are selected and saved.

[0040] In some embodiments, the mathematical model includes a plurality of nodes representing molecular species (e.g., proteins and nucleotides) in the GRN. In some embodiments, one or more differential equations are used to model dynamic interactions (e.g., complex interactions between nodes which can change in response to different stimuli or over time) between the nodes. In some embodiments, the dynamic interactions are observed over time (e.g., over seconds, minutes, or other durations corresponding to the length of an experiment). In some embodiments, the goal behavioral states are provided as steady states of the dynamic interactions observed over time.

[0041] In some embodiments, using the algorithm further includes randomly sampling more than one source of variation (e.g., different functionality of the network, perturbations toTufts Ref. T002751 Quarles Ref.166118.01450 network parameters, initial conditions, or external stimuli) and determining a corresponding goal behavioral state and corresponding trajectory to populate an initial mapping of the system, including sources of variation, trajectories, and goal behavioral states. In some embodiments, the initial mapping is further populated by sampling the goal behavioral states of the initial mapping and optimizing additional goal states and additional trajectories. In some embodiments, the goal behavioral states and corresponding trajectories are saved.

[0042] In some embodiments, randomly sampling further includes using a uniform sampling strategy (e.g., a naïve random screening strategy) within a bounding hyper-rectangle of currently reached goal behavioral states, in which a volume of the hyper-rectangle is larger than a distribution of currently reached goal behavioral states. The bounding hyper-rectangle can be described as an intervention space, or a space where interventions represent controlled sources of incoming variation. Practically, this sets the initial state of a model rollout.

[0043] In some embodiments, testing the robustness includes perturbing the corresponding trajectories with at least one of noise perturbations, push perturbations, or wall perturbations. In some embodiments, the noise perturbations are parameterized by their standard deviation and period. In some embodiments, the push perturbations are parameterized by their magnitude and number of occurrences. In some embodiments, the wall perturbations are parameterized by their length and number. In some embodiments, the walls are generated in locations of a space that the GRN would visit without induced perturbations.

[0044] In some embodiments, each goal behavioral state has more than one corresponding trajectory and each of the corresponding trajectories reaches the goal behavioral state when perturbed. In some embodiments, a goal behavioral state is a cell phenotype, such as cell metabolism, cell proliferation, cell differentiation, cell migration, cell shape, cell signaling, programmed cell death, gene expression, or senescence. In other embodiments, the goal behavioral state is a multi-cellular phenotype, such as morphogenesis, remodeling, wound healing, regeneration, or embryonic development. In still other embodiment, a goal behavioral state is potentiation of long-term memory.

[0045] In another aspect, a method or system of identifying therapeutic interventions (e.g., suggested genetic therapies, drugs or medicines, or the timing with which to administer a specific treatment) is disclosed. In some embodiments, a behavioral catalog of a GRN associatedTufts Ref. T002751 Quarles Ref.166118.01450 with a disease is obtained, the behavioral catalog including behavioral end states and corresponding trajectories connecting controlled sources of variation and each behavioral end state. In some embodiments, a disease associated with the GRN is selected. In some embodiments, each behavioral end state is classified as a first setpoint, and the first setpoint can be a disease setpoint or a healthy setpoint. In some embodiments, a disease setpoint is associated with a phenotype or gene expression pattern associated with the disease, and a healthy setpoint is not associated with a phenotype or gene expression pattern associated with the disease. In some embodiments, more than one stimulus-based interventions (e.g., perturbations to one or nodes of a system, with pre-determined timing and amplitude) to apply to each disease setpoint are defined. In some embodiments, each stimulus-based intervention is introduced to each first setpoint and trajectories are identified, each trajectory corresponding to a response of a first setpoint to a stimulus-based intervention and a second setpoint of each trajectory is determined. In some embodiments, a subset of trajectories in which the first setpoint is a disease state and the second setpoint after introducing the stimulus-based intervention is a healthy state. In some embodiments, the subset of stimulus-based interventions that correspond with the subset of trajectories in which the first setpoint is a disease state and the second setpoint is a healthy state, and the subset of stimulus-based interventions are output. Stimulus-based interventions are optimized based on which nodes to intervene on, when to apply the intervention, and with what amplitude to apply the intervention.

[0046] In some embodiments, classifying healthy setpoints and disease setpoints includes constructing an energy landscape (e.g., a graph of the energy function across the configuration space of a system) of the behavioral catalog of the GRN and determining which regions of the energy landscape correspond with a disease region, in which the phenotype or gene expression pattern of the region is associated with a disease, and determining which regions of the energy landscape correspond with a healthy region, in which the phenotype or gene expression pattern of the region is not associated with the disease. In some embodiments, the method or system further includes identifying a subset of endpoints in the behavioral catalog in which each endpoint falls within the disease region and classifying the endpoint as a disease setpoint, and identifying a subset of endpoints in the behavioral catalog in which endpoint falls within the healthy region and classifying the endpoint as a healthy setpoint. In some embodiments, the method or system further includes selecting a subset of trajectories associated with a diseaseTufts Ref. T002751 Quarles Ref.166118.01450 setpoint and a subset of trajectories associated with a healthy setpoint that exhibit a low level of sensitivity to perturbations (e.g., high robustness and versatility, or an ability to reach the healthy setpoint under various perturbations or interventions) and outputting the subset of trajectories and associated end setpoints.

[0047] In some embodiments, each stimulus-based intervention is a piecewise constant function that determines, which nodes in the GRN to intervene on, when to apply the intervention, and with what amplitude the stimulus-based intervention should be applied. In some embodiments, the nodes in the GRN are selected based on which nodes can be targeted with drugs that treat the disease, and with what precision the drugs can be administered. In particular embodiments, introducing the stimulus-based interventions further includes applying a plurality of perturbations, in which the perturbations include at least one of noise perturbations, push perturbations, or wall perturbations. In some embodiments, selecting the subset of trajectories further includes selecting trajectories which reach the healthy state under all perturbations applied to the trajectory.

[0048] In some embodiments, the GRN is the RKIP-ERK signaling pathway. In some embodiments, the phenotype or gene expression pattern of the disease region is associated with tumor development. In some embodiments, the disease region has a phenotype or gene expression pattern associated with tumor development. In some embodiments, the healthy end state is a state in which a cell has an elevated susceptibility to therapies, such as cancer therapies.

[0049] Also disclosed herein is a method or system of identifying a configuration (e.g., an initial condition and / or external stimuli) or parameter of a gene regulatory network to perform a desired function. In some embodiments, the method includes obtaining a data set including a mathematical model of a GRN. The method may include using an algorithm to identify one or more goal behavioral states and corresponding trajectories connecting controlled sources of variation in the GRN and each of the goal behavioral states. The method optionally may include using a local optimization strategy to identify additional goal behavioral states. In some embodiments, the GRN is an oscillator circuit (e.g., a biological oscillator with an inhibitory feedback loop which includes one or more oscillating variables, and a source of delay in the feedback loop). In some embodiments, the controlled sources of variation can include one or more of a genetic intervention or an environmental intervention. In another aspect, a method ofTufts Ref. T002751 Quarles Ref.166118.01450 priming a biological system for therapeutic treatment is disclosed. The method may include administering a compound to the system to move the system to a desired state, such as a state that the biological system has not yet achieved (e.g., a more differentiated state than the biological system initially exhibits). In some embodiment, the biological system is in a disease state. In some embodiments, a desired state may be a biological process with a phenotype associated with acquired sensitivity to at least one therapeutic treatment. A desired state may also be restored sensitivity, acquired sensitivity, reduced resistance to, or prevented resistance to at least one therapeutic agent. A desired state may also be a cell phenotype that is susceptible to a therapeutic agent that the biological system was previously not susceptible to. In some embodiments, an administered therapeutic treatment includes administering at least one drug. The drug can be selected from a group of chemical classes such as propanolamine, C-amidated tetrapeptide, aryl sulfide, alkyl-phenylketone, benzazepine, secondary amine, benzenoid, monomethoxylflavone, monohydroxypyrrolidine, benzodiazepine, quinazoline, quinazoline derivative, phenyltriazine, phenol, propanolamine, non-imidazole heterocyclic, and / or hydroxamate. In some cases, more than one drug is administered. In some embodiments, it is then determined whether the biological system has reached the desired state.

[0050] Also provided herein are methods of administering a therapeutic treatment to a biological system by administering a compound to a biological system to move the biological system to a desired state. The desired state may be a more differentiated state of the biological system. It can also be determined whether the biological system has achieved the desired state. The method may include administering multiple compounds to the biological system. In some embodiments, a new therapeutic treatment is selected based on the achieved desired state. In some embodiments, the biological system is a patient, such as a patient with cancer. The cancer may be a hematologic cancer (e.g., leukemia or lymphoma). In some embodiments, the compound may be selected from a group of chemical classes including propanolamine, C- amidated tetrapeptide, aryl sulfide, alkyl-phenylketone, benzazepine, secondary amine, benzenoid, monomethoxylflavone, monohydroxypyrrolidine, benzodiazepine, quinazoline, quinazoline derivative, phenyltriazine, phenol, propanolamine, non-imidazole heterocyclic, and / or hydroxamate.Tufts Ref. T002751 Quarles Ref.166118.01450

[0051] The ability to understand, map, predict, and control chemical and genetic networks is critical for a deeper understanding of evolutionary biology, comprehending and intervening in various disease states, and constructing novel synthetic biologicals in bioengineering. A significant amount of work into mathematical modeling of gene regulatory networks (GRNs) has been done previously. However, there is still a lack of understanding of how different initial states and environmental stimuli can affect the behaviors GRNs can exhibit. It has been hypothesized that GRNs and other molecular networks could be endowed with surprising navigation competencies which allows them to robustly reach diverse homeostatic of allostatic states despite a wide range of perturbations.

[0052] 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. 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. This is critical not only for fundamental understanding of evolution of bodies and minds, but also for the design of interventions in biomedicine and synthetic morphology contexts.

[0053] 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. Recently, a broad, substrate-independent behavior science perspective has been suggested to reveal novel properties of GRNs and other biological networks. 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.Tufts Ref. T002751 Quarles Ref.166118.01450

[0054] Investigating GRNs in the lab has many practical challenges. In the lab or clinic, GRNs 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). However, it is difficult to identify candidate genes and proteins to modify, and how to modify them in order to reach a desired system-level state due to the emergent, generative nature of development and physiology. Experimental techniques, particularly experiments that need to be carried out in vivo, inherently require a large amount of time and resources, and have a large degree of variation. 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 even expand, once the drug is withdrawn. Current therapeutics function bottom- up, attempting to force specific molecular states. This has made it challenging to develop methods for shifting complex tissues and organs towards a stable health profile.

[0055] Next-generation solutions, which would offer stable correction of a GRN and thus true healing, 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 stimuli is essential for predicting what stimuli can be used to re-set their setpoints and guide them to autonomously maintain a healthy state. Therefore, much work has gone into mathematical modeling and computational inference of both protein pathways and gene regulatory network models, which has resulted in the development of large collections of publicly available models. 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 data 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, it is not known what other sets of tools might reveal and exploit interesting properties of this ubiquitous biological substrate.Tufts Ref. T002751 Quarles Ref.166118.01450

[0056] Mathematical modeling has the potential to characterize biomolecular pathways at steady state and to determine controlled sources of variation and determine the dynamic responses of the pathways. Prior research has been limited due to many challenges. To characterize steady-state concentrations and responses to small perturbations, conventional methods rely on piecewise-linear approximation of the system behavior but struggle with higher- dimensional systems or wider parameter ranges which limits their applicability. 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. These methods typically exploit the network topology or regulatory structure to identify control strategies based either on permanent knockout / activation of genes or on temporary perturbations, the latter being preferable in biomedical context. Identifying how external stimuli may affect a GRN is similarly difficult. Pathway dynamics are non-linear and redundant, which makes passive exploration strategies likely to either fail in uncovering the full range of potential behaviors or require a prohibitive amount of time and energy.

[0057] Existing 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 and gradient-descent optimization, for controlling continuous ordinary differential equation (ODE) models of biomolecular networks with high-dimensional parameter spaces, mainly in the context of synthetic circuit engineering. 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 but largely restricts exploration of the behavior space. Prior studies use random perturbations in network parameters or specific gene knockouts. Environmental perturbations are often limited to random variations in initial conditions within a predefined range. Overall, fitness functions and methods of introducing perturbations to the system are currently greatly limited.

[0058] 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, there remains a need for methods for automating efficient exploration and discovery of a diversity of behaviors in these spaces.Tufts Ref. T002751 Quarles Ref.166118.01450

[0059] Here, we formalize and investigate a view of gene-regulatory networks as agents navigating a problem space. We develop automated tools to efficiently map the repertoire of robust goal states that GRNs can reach despite perturbations, using a limited budget of experiments. These tools rely on two main advancements: (1) Using curiosity driven exploration algorithms, originating from the AI community to explore the range of behavioral abilities of a given system, that we adapt and leverage to automatically discover the range of reachable goal states of GRNs and (2) Proposing a battery of empirical tests inspired from behaviorist tradition to assess their navigation competencies.

[0060] We investigate a view of gene-regulatory networks as agents navigating a problem space toward target goal states and with varying degrees of competency (FIG.1A). 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. This disclosure focuses on two initial goals: identifying the repertoire of robust goal states that a GRN can actively reach through minimal and non-genetic interventions within a navigation task context, and developing systematic methods and automated tools to aid scientists in discovering this repertoire.

[0061] To address these goals, the experimental framework revolves around the definition of ‘problem spaces,’ which we use as tractable components of the GRN’s overall state space (FIG.1B), and on a set of methodological contributions. This can be divided further into three subcategories. First, we develop methods of automated discovery of behavioral abilities within autotelic curiosity search. 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 use intrinsically-motivated goal exploration processes (IMGEP), a recent family of diversity-driven machine learning approaches also known as autotelic curiosity search. Second, we evaluate the navigation competencies. Specifically, we examine how robust and versatile GRNs are in attaining the diverse previously identified goal states. 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 maintainTufts Ref. T002751 Quarles Ref.166118.01450 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. Third, we find reuses of the discovered behavioral catalog and framework. We demonstrate 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".

[0062] Applying IMGEP on a database of 32 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. 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 that could be flexibly and robustly reprogrammed by appropriate stimuli. 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.

[0063] 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 – can be transposed to assist efficient discovery of behavioral abilitiesTufts Ref. T002751 Quarles Ref.166118.01450 within biological pathway models like GRNs. Secondly, rather than merely mapping the attractor states or analyzing their sensitivity to model parameter changes 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. Notably, these competencies are discovered without necessitating structural alterations to network properties or connectivity. An advantage of our framework and its associated tools is that this method does not make any assumptions about the structure or origin of the biological network, making it adaptable to the study of diverse unconventional intelligences across various domains.

[0064] The disclosed procedures have many practical applications that can further advance scientific discovery. 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 may be an advancement for understanding the GRN functionalities as well as for designing drug-driven interventions in a biomedical or bioengineering context.

[0065] The methods presented herein can be used in using the automated discovery toolbox to assist model inference, allowing to efficiently identify the rare or unexpected behaviors of ordinary differential equation models and suggest whether further refinement is needed or not. The disclosed procedures can be directly applied to in vitro GRN models at the bench. One could 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, and (2) conduct studies of subcellular proto- cognitive phylogenetics, to help understand the evolutionary pressures for and against reprogrammability in cell regulatory machinery.Tufts Ref. T002751 Quarles Ref.166118.01450

[0066] The predictions made by our computational tools can be tested in vivo, 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. 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.

[0067] This method can be used in higher-dimensional and more complex problem spaces by using unsupervised learning of goal space representations. Whereas this disclosure is applied here 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.

[0068] This disclosure improves on existing methods in several ways. Previous methods utilized a predefined fitness function and strict parameter range initialization which greatly restricts exploration of the behavior space. We show that traditional screening methods can be very inefficient in discovering the range of possible goal states. In contrast, we can deploy strategies to discover and map behavior space because of our approach of modeling GRN behavior within a tractable behavior space Z. Notably, recent diversity-driven machine learning techniques such as Novelty Search, Quality Diversity, and IMGEP 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. IMGEPs, which were originally developed for the learning of inverse models of highly redundant mapping in robotics context, were recently shown to successfully assist discovery in complex self-organizing systems.

[0069] Having efficient exploration strategies like curiosity search can be highly valuable, especially when scaling to larger databases. The goal-directed strategy of the IMGEP is particularly beneficial for ( ^^ ^^ ^^, ^^, ^^) systems with high nonlinearity or redundancy in their ^^ → ^^ mapping. Redundancy implies that many interventions in ^^ lead to similar effects in ^^. 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 coverTufts Ref. T002751 Quarles Ref.166118.01450 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. Overall, the disclosed procedures efficiently drive a GRN into diverse goal states with minimal interventions.

[0070] In existing methods, perturbations of GRNs 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. Here instead, we conduct a battery of empirical tests that draw inspiration from classical “displacement experiments” and “barrier experiments” commonly used in behavioral sciences to assess the navigation competencies of various animals. This offers a better characterization of the versatility and robustness of pathways.

[0071] Additionally, this method offers improvements as an alternative approach in synthetic biology, particularly as an alternative to gene circuit engineering. Our 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, our framework only relies on empirical investigation for inferring the mapping between interventions and behaviors (treating them as abstract variables in observable problem spaces). Therefore, it offers an exciting perspective to be deployed across various problem spaces and classes of organisms.

[0072] The following examples set forth, in detail, ways in which the present disclosure may be used or implemented and will enable one of ordinary skill in the art to more readily understand the principles thereof. The following examples are presented by way of illustration and are not meant to be limiting in any way. EXAMPLES EXAMPLE 1: Materials and Methods

[0073] 1.1 GRN models and numerical simulation

[0074] This study employs ordinary differential equation (ODE) models to represent molecular pathways, with nodes representing pathway components and edges capturing theirTufts Ref. T002751 Quarles Ref.166118.01450 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.

[0075] 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

[0135] . 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.

[0076] Additionally, we have developed a python library (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.

[0077] Given the model species initial state ^^( ^^ = 0), the desired rollout length ^^ ( ^^ ^^ ^^ ^^) and step size ^^ ^^, as well as the chosen intervention ^^ and / or perturbation ^^, the model rollout iteratively 1) integrates the system of ODE-governed equations that specifies the rate of species changes using JAX odeint solver to update model species ^^( ^^) → ^^( ^^ + ^^ ^^), 2) calls the model assignment rules to update kinematic parameters if needed, and 3) applies the intervention and / or perturbation function to update ( ^^( ^^ + ^^ ^^), ^^( ^^ + ^^ ^^), ^^) accordingly. In this paper we use ^^ = 2500 ^^ ^^ ^^ ^^ and ^^ ^^ = 0.1 (25001 time points per rollout including ^^,). The ODE solver uses an absolute tolerance of 1 ^^-6and relative tolerance of 1 ^^-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 (developmentalsystems.org / curious-exploration-of- grncompetencies / tuto1.html), which offers practical examples and detailed instructions.

[0078] 1.2 Experimental set upTufts Ref. T002751 Quarles Ref.166118.01450

[0079] In our computational models, we are able to record the activities of all nodes a model rollout. The observation space ^^ ⊂ ℝ^ൈ ^during∆^ା is such that ^^ = ( ^^(0),…, ^^( ^^)) 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.

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

[0081] For the main experiments on biological networks, the intervention space ^^ ⊂ ℝ^ାconsists of initial node states sampled from the hyper-rectangle [ ^^0,min, ^^0,max] where^^ ൈ ^^ௗ,^^^ and ^^^,^^௫ ൌ ^^ ൈ ^^ௗ,^^௫ with r – 20 and (yd,min, yd,max) the minimum and maximumof 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 ^^ ⊂ ℝ⊥2endpoint states ^^ = (yi(T), yj(T)) where ( ^^, ^^) 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.

[0082] 1.3 Database creation

[0083] Biological networks database

[0084] All the ODE models we use in this work are downloaded from the BioModels database

[0021] ,

[0022] 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

[0135] ). To ensure the inclusion of models suitable for our analyses, we applied specific filters to the collected models.

[0085] 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 ^^ ^^ = 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 forTufts Ref. T002751 Quarles Ref.166118.01450 specific initial conditions), we consider the model invalid if this occurs for the default initial conditions provided in the SBML file.

[0086] 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 ^^ (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 (> 20%) exhibited invalid concentrations (NaN or negative) or unsettled behaviors (∃ ^^ ≥ 2400 ^^ ^^ ^^ℎ ^^ℎ ^^ ^^ | ^^( ^^) − ^^( ^^)| ≥ 0.02 × | ^^( ^^) − ^^(0)| ) or periodic patterns (∃ ^^ > 0 ^^ ^^ ^^ℎ ^^ℎ ^^ ^^ | ^^( ^^)| ≥ 40 ^^ℎ ^^ ^^ ^^ ^^ = ^^ ^^ ^^([ ^^(T / 2),⋯, ^^( ^^)]); or [filter F2] the reached space in ^^ was too small ^൬^^௫ ^^ୀ^∙∙∙ହ^^^ ^ ^^^ െ^^^ ^^ୀ^∙∙∙ହ^^^ ^ ^^^^ ^ 0.1^ to discard cases where “diversity” could result [filter F3] the number ofattractors was less than four ^^ ^^^^ ^^^^^ୀ^∙∙∙ହ^^^ ^^ ^^ ^^ ^^ ^ 4 ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ 20 ൈ 20 binning of Z^.

[0087] Upon completion of the filtering process, our final database comprised 32 models, consisting of a total of 432 systems, as detailed in FIG.13. These curated models and systems served as the foundation for our subsequent analyses and investigations into the navigation competencies of the molecular pathways.

[0088] Random networks database

[0089] Following the methodology proposed in

[0079] , 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

[0079] ). To achieve this, we initially generated 320 networks based on the transcriptional gene circuit model

[0117] , 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 ^^, ^^, ^^ of these networks wererandomized ( ^^ ∼ [−30,30]n×n, ^^ ∼ [−10,10]n, ^^ ∼ [1,15]) where model step is defined as ^^^ ^^ ^1^ ൌ ∆் ൈ ∆்ఛ ^^ ^^ ^^ ^^ ^^ ^^ ^^^ ^^ ^^ ^ ^^^ ^ ^1 െ ఛ ^ ൈ ^^ and in-degree connectivity is enforced byTufts Ref. T002751 Quarles Ref.166118.01450 setting some weights of ^^ to zero. However, during the creation process, we observed that none of the generated networks met the criteria 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.

[0090] 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 FIG.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.

[0091] 1.4 Curiosity-driven exploration

[0092] 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 FIG.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 ^^, identifies its associated intervention, and performs a local random step from that point ( ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ∼\ ^^ ^^ ^^ℎ ^^ ^^ ^^{ ^^}(0, 0.1 ∗ [y0,max – y0,min]) in the intervention space.

[0093] 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 (github.com / flowersteam / autodiscjax) that accompanies this paper can be used to implement this and other IMGEP variants in JAX.

[0094] 1.5 Robustness testsTufts Ref. T002751 Quarles Ref.166118.01450

[0095] We define 3 family of perturbations: 1) the noise perturbation ^^n( ^^n, ^^n| ^^) which is parametrized by its standard-deviation (scaled proportionally to the extent of the observed trajectory ^^ prior perturbation) and period (secs); 2) the push perturbation ^^p(mp, np| ^^) parametrized by its magnitude (proportional to the extent of ^^) and number of occurrences; 3) the wall perturbation ^^w( ^^w, ^^w= | ^^) parametrized by its length (proportional to the extent of ^^) 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 FIG.12.

[0096] To assess the robustness of the GRN systems in our database, we employ an evaluation procedure, as depicted in FIG.1D. For each system ( ^^, ^^) in the database with its corresponding behavioral catalog ^^ discovered using the curiosity-search algorithm, we perform the following steps. We first retrieve ^^ representative trajectories out of the ^^ 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 to reduce the computation burden, as we run the experimental campaign on all 432 systems. Next, we subject each of these K trajectories { ^^k, ^^ = 1.. ^^} to s=18 different perturbation distributions, each representing various levels of difficulty:( ^^n, ^^n) ∈ {(0.001,5), (0.005,5), (0.1,5), (0.005,10), (0.005,5), (0.005,1)}, ( ^^p, np^ ∈ {(0.05,1), (0.1,1), (0.15,1), (1,0.1), (2,0.1), (3,0.1)}, ( ^^w, ^^w) ∈ {(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 ^^ = ^^ ∗ ^^ perturbations. For each perturbation in the set { ^^p, ^^ = 1... ^^}, we re-run the trajectory starting from the same initial state ^^ but with the sampled perturbation applied ( ^^, ^^p), and observe the resulting outcome ( ^^p) and reached endpoint ( ^^p).

[0097] At the end of this process, the behavioral catalog is augmented with the perturbed trajectories ^^ = {( ^^k, ^^k, ^^k, {( ^^p, ^^p, ^^p), ^^ = 1... ^^}), ^^ = 1… ^^}.

[0098] 1.6 Evaluation Metrics

[0099] Diversity measureTufts Ref. T002751 Quarles Ref.166118.01450

[0100] 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.

[0101] In general, existing approaches in the NS, QD and IMGEP literature use binning- based metrics

[0090] ,

[0091] ,

[0136] or distance-based metric from ecology

[0137] 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

[0105] . 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 { ^^ ∈ ^^}. 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

[0105] .

[0102] Threshold coverage quantifies the area of the space that has been reached at a given precision ^^ (the threshold). Threshold coverage was used in FIG.3 to compare random search and curiosity-driven exploration strategies.

[0103] Sensitivity measure

[0104] 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

[0093] or phenotypic trait

[0094] 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 ^^ ∈ ^^ discovered by curiosity search, 2) the set of perturbation { ^^p} is defined in previous section and conditioned on the GRN goal-reaching trajectory ^^ → ^^, and 3) variation is measured as the Euclidean distance in behavior space, normalized by the extent of the trajectory prior perturbation in Z.

[0105] 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 FIG.4.

[0106] Versatility-Robustness measure

[0107] In this study, we introduce the terms "diversity" and "versatility" to characterize the competencies of the exploration agent (IMGEP) and the gene regulatory network agentTufts Ref. T002751 Quarles Ref.166118.01450 (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.

[0108] In FIG.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 ൈగఢ^^ାଶఢ^మ.

[0109] In FIG.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.

[0110] 1.7 Experiments on the RKIP-ERK signaling pathway

[0111] This section details the additional experiments conducted on the RKIP-ERK signaling pathway

[0099] . We refer to the accompanying notebook tutorial for reproducing these experiments: developmentalsystems.org / curious-exploration-of-grncompetencies / tuto1.html.

[0112] For FIG.4, clustering in behavior space was performed using the HDBSCAN algorithm

[0106] 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

[0115] .

[0113] In FIG.6, we generated trajectory-based energy landscapes following the method proposed in

[0113] . 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 energyTufts Ref. T002751 Quarles Ref.166118.01450 surface (U = −ln(P)). This energy surface was smoothed using cubic spline interpolation and visualized using Plotly 3D surface plots. FIGS.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.

[0114] In FIG.8, the "healthy" and "disease" clusters were the same as in FIG.4 and visualized similarly. We displayed trajectories with the lowest sensitivity (averaged over all ^^ =3 × 18 perturbations). The stimuli-based intervention shown in FIG.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 ^^ = [ ^^MEKPP^1^,⋯, ^^MEKPP^10^] every 10 seconds during 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 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 FIG.8B) succeeds to robustly reset all 6 disease state points despite perturbations, as shown in FIG.8C. We refer to the notebook for reproducing the experiments.

[0115] 1.8 Experiments on synthetic gene networks

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

[0117] In these experiments, we consider the target application of gene circuit engineering followed in

[0072] , where parameters of a gene circuit model are optimized to produce target oscillator patterns. The gene circuit model employed in

[0072] is the same than the one used for the random networks database (Eq 1), with ^^ = 1. Hence the intervention space is a ^^1 + 2 ^^ dimensional space defined as ^^ = [ ^^t=0,min, yt=0,max]⊕[ ^^min, ^^max]⊕ [ ^^min, ^^max^, with ^^0,min= 0, ^^0max= 1, ^^min= −30, ^^max= 30, ^^min= −10, ^^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 kinematicTufts Ref. T002751 Quarles Ref.166118.01450 parameters ( ^^, ^^) and initial concentrations ^^, that would produce a periodic pattern ^^ = [ ^^nൌ0(0),⋯, ^^n=0(T)] with target amplitude ^^, frequency ^^ and offset ^^. Here, the target ( ^^, ^^, ^^) are sample randomly with ^^ ∼ ^^([0.1,0.5]), ^^ ∼ ^^([ ^^, 1 − ^^]), ^^ ∼ ^^ ^^ ^^ ^^( ^^ = 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

[0072] , all given the same experimental budget ^^ = 5000.

[0118] For curiosity search, the behavior space ^^ is defined as the image space of the discrete Fourier transform of the 1d-signal ^^, 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 FIG.1C, but operating within the novel problem spaces ( ^^, ^^).

[0119] For random search, interventions are sample uniformly ( ^^1,⋯, ^^N) ∼ ^^( ^^).

[0120] For gradient descent, we follow the procedure proposed in

[0072] . We define a loss function which, for a set of parameters ^^ = ( ^^, ^^, ^^0), measures the mean square error between the phenotype node activation levels ^^ and the target oscillation represented as a cosine wave with the desired ( ^^, ^^, ^^) : ^^ = Σt( ^^( ^^) − ( ^^ ^^ ^^ ^^ (2 ^^ ^^ × ^^) + ^^))2. We then sample a random parameter ^^ ∼ ^^( ^^) and use Adam optimizer with ^^r = 10-3, ^^1 = 0.02, ^^3 = 0.001, ^^ = 10-8for ^^ = 5000 optimization steps (same number of model rollouts allowed than for curiosity search and random search).

[0121] 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 ^^ = 100 optimization steps.

[0122] Visualizations in FIG.9 show: (FIG. A, FIG. B) the oscillators discovered by random search and curiosity search (gradient descent did not find any oscillator in this example) in the ( ^^, ^^, ^^) space, (FIG. C) the corresponding diversity (using this time a binning-based space coverage measure with 203bins as the space is 3-dimensional), (FIG. D) the evolution of the training loss ^^ throughout the N=5000 trials for the three exploration strategies, (FIG. E, FIG. F, FIG. G) the corresponding best discoveries (for which ^^ is minimal) for the three exploration strategies, and (FIG. H, FIG. I) the local training loss and resulting finetuning of the best discoveries with gradient descent.Tufts Ref. T002751 Quarles Ref.166118.01450 EXAMPLE 2: Generalizing GRN behavior as a navigation task Dynamical Systems Behavioral Science Proposed Navigation Task Terminology Terminology Isomorphism Terminology System: a set of Organism: a living Both are collections Agent or GRN interconnected being that responds to of lower-level elements that interact stimuli and adapts to elements that interact to produce emergent its environment to produce emergent behavior behavior and can adapt at the system level Phase-space Behavioral Both represent the trajectory trajectory: set of trajectory: the sequence of states or states taken by the sequence of states behaviors that a system when starting that an organism system or individual from one particular exhibits in response experiences over time initial condition to stimuli Initial condition: Stimuli: an event hat Both represent Intervention or initial state of a might (or might not) incoming variations perturbation system’s variables trigger a response in that set a system or and parameters that an organism organism in motion condition its dynamics Critical parameter: Salient stimuli: Both represent the Effective a parameter or stimuli that are incoming variations intervention condition that, if particularly relevant that have a significant changed, can cause a or meaningful to an impact on a system’s system to undergo a organism, either steady-state or qualitative change or because they are organism’s response phase transition associated with reward or punishment or because they are novel or unexpected Steady-state (or Observed response: Both represent the Reached endpoint or attractor): a stable outcome or endpoint endpoint that a goal state (or set of states), of a behavioral system or organism is towards which the trajectory towards moving towardsTufts Ref. T002751 Quarles Ref.166118.01450 system tends to which an organism evolve over time converges Robust attractor: Target goal: it is Both represent a Robust goal stable attractor assumed that an stable endpoint or toward which the organism engages in goal that the system system tends to a goal-directed successfully attains evolve under various manner when it under various initial conditions and exhibits new ways or perturbations perturbations actions to achieve a similar outcome when faced with novel circumstances Controllability: Trainability: degree Both represent the Versatility degree to which the to which an capacity of a system system’s dynamics organism’s behavior or individual to be (and resulting steady can be modified or influenced or states) can be shaped by experience changed by controlled or or conditioning controlled manipulated) interventions Table 1: Glossary of terms used in this application, with the proposed isomorphism which generalizes concepts from dynamical complex systems and behavioral sciences under a common navigation task perspective.

[0123] The GRNs analyzed in this study are biological pathway networks taken from the BioModels repository

[0021] ,

[0022] . 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 (FIG.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,Tufts Ref. T002751 Quarles Ref.166118.01450 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. Problem Space Generic Definition Specific definition in this application Observation Space (O) Space of raw observations Records node activities over made during the GRN model time so as o – (y(0), …, rollout to measure its state or y(T)), where y(t) is an n- behavior. dimensional vector (n = number of nodes) and T is the measured reaction time. Behavior Space (Z) A projection of the Encodes the trajectory observation space used by the endpoint of a model rollout. experimenter to encode the Represents a cell phenotype ‘goal states’ of a model defined by the state values of rollout into a tractable (lower- some nods (relevant dimensional) space. biological markers), such that z = (yi1(T), … yim(T)) (we use m=2 in this study for simplicity and visualization). Intervention Space (I) A space where interventions Sets the initial state i = (y1(0), represent controlled sources …, yn(0)) of a model rollout. of incoming variation that the Defined as an hyper-rectangle experimenter can exert on the ^^ ⊆ ℝ^where the boundaries GRN model rollout to drive it are 20 times larger than the toward novel or targeted min and max values taken by states. the respective nodes from default initial conditions. Perturbation Space (U) A space where perturbations Includes three classes of represent external sources of (stochastic) perturbations incoming variation, used by including noise perturbation the experimenter to Un, push perturbation Up, and characterize the robustness of wall perturbation Uw. a given goal state. Table 2: Problem spaces used in this study.

[0124] 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 beTufts Ref. T002751 Quarles Ref.166118.01450 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 nongenetic 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 FIG.2.

[0125] 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 is in charge of organizing the sequence of interventions ^^1,…, ^^N used to drive the system toward a maximally diverse set of reachable endpoints { ^^k ∈ ^^}k=1..N, while being given a limited budget of experiments N. Secondly, to assess the robustness of the discovered goal states { ^^k ∈ ^^}, 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 ^^ = {( ^^k, ^^k, ^^k, {( ^^p, ^^p, ^^p), ^^ = 1... ^^}5, ^^ = 1… ^^}.

[0126] Following this framework, the behavioral catalog is constructed for a database of 32 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 FIG.13. 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 FIG.3, FIG.5, and FIG.7, and specific results for the RKIP-ERK signaling pathway

[0099] are shown in FIG.2, FIG.4, FIG.6, and FIG.8. EXAMPLE 3: Curiosity search uncovers a diversity of reachable goal states

[0127] 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, recentTufts Ref. T002751 Quarles Ref.166118.01450 diversity driven machine learning techniques such as Novelty Search

[0100] ,

[0101] , Quality Diversity

[0102] ,

[0103] and Intrinsically-Motivated Goal Exploration Processes (IMGEP)

[0086] ,

[0087] 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

[0086] , were recently shown to successfully assist discovery in complex self-organizing systems

[0088] -

[0091] .

[0128] Here, we propose to use an IMGEP to control GRN initial states and maximize the diversity of discovered endpoints { ^^ ∈ ^^} within a limited budget of ^^ experiments. The IMGEP operates in two phases: initially, ^^'() interventions are sampled randomly from ^^ to populate history ^^, then remaining interventions are generated through a goal-directed process which relies on several key internal models. Those including a goal-embedding module ( ^^) that encodes observations ( ^^) into the IMGEP goal space (Ƭ), a goal generator module ( ^^) 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 ( ^^) 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 ^^~ ^^( ^^), 2) infer intervention parameters to achieve that goal ^^ ∼ ^^( ^^, ^^), 3) conduct an experiment with the intervention i, observe the outcome o, and compute the reached goal ^^ = ^^( ^^), and 4) store the tuple ( ^^, ^^, ^^) in history ^^. Here, we use the GRN behavior space Z as the IMGEP goal space Ƭ = ^^. 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

[0104] . The pseudocode of our IMGEP pipeline is shown in FIG.1C 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: ^^ = {( ^^k, ^^k, ^^k), ^^ = 1… ^^}.

[0129] 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: theTufts Ref. T002751 Quarles Ref.166118.01450 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.

[0130] For evaluating diversity, we measure the area covered in Z by the IMGEP discoveries using the threshold-coverage metric

[0105] and compare it with the area covered by the diversity of a naive random screening strategy (which uniformly samples interventions in ^^). In FIG.3, the diversity discovered by the two exploration variants is shown for the 432 ( ^^ ^^ ^^, ^^, ^^) systems, where random search is given a budget of experiments ^^ which is twice bigger (N=900) as the one given to the curiosity-search algorithm (N=450).

[0131] 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 our dealing with numerical systems and our codebase allow for efficient and parallel execution, each experiment still consists of்∆்= 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 FIG.3B, 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 FIG.3C, there are some systems for which random search is already quite efficient in revealing diverse behaviors in Z, and for which IMGEP performs equivalently.

[0132] We evaluate diversity using the threshold-coverage metric

[0105] . FIG.3 compares the diversity discovered by the curiosity-search algorithm (N=450) and random search (N=900) for the 432 (GRN, I, Z) systems. Remarkably, at n=290, curiosity search already significantly outperforms the final diversity achieved by random search, using only one third of its experimental budget (N=900). Although our codebase allows for efficient and parallel execution, each experiment consists of்= 25000 model steps, making the repetition for all N experiments and all systems in ourcostly. Having efficient exploration strategies like curiosityTufts Ref. T002751 Quarles Ref.166118.01450 search can be highly valuable, especially when scaling to larger databases. In fact, the goal- directed strategy of the IMGEP is particularly beneficial for ( ^^ ^^ ^^, ^^, ^^) systems with high nonlinearity or redundancy in their ^^ → ^^ mapping, as seen in FIG.4 and studied in robotics contexts

[0105] . Redundancy implies that many interventions in ^^ lead to similar effects in ^^, as illustrated in FIG.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.

[0133] 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 FIG.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 diseases are due to epigenetically non-identical cells

[0107] . EXAMPLE 4: Curiosity search uncovers a diversity of reachable goal states

[0134] 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

[0093] ,

[0094] , 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

[0098] ,

[0095] ,

[0096] ,

[0071] ,

[0108] ,

[0060] . Here instead, we propose to conduct a battery of empirical tests that draw inspiration from classical “displacement experiments”

[0109] ,

[0110] and “barrier experiments”

[0111] commonly used in behavioral sciences to assess the navigation competencies of various animals. As illustrated in FIG.2, we consider environmental perturbations that perturbTufts Ref. T002751 Quarles Ref.166118.01450 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)

[0112] , 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

[0107] .

[0135] 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) the trajectories previously discovered using the curiosity-search algorithm and subject these trajectories to ^^ = ^^ × ^^ perturbations conditioned on the GRN goal-reaching trajectory ^^ → ^^ prior perturbation. Here, s represents the different perturbation which correspond to various “tests” and “levels of difficulty” (e.g., noise magnitude and frequency, number of walls, etc.) and ^^ is the number of (stochastic) perturbations sampled per family. The pseudocode is illustrated in FIG.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 ^^ = {( ^^k, ^^k, ^^k, {( ^^p, ^^p, ^^p), ^^ = 1... ^^}), ^^ = 1… ^^}.

[0136] As the use of “spaces” comes with the notion of similarity and distance, we can then evaluate the sensitivity of a goal state ^^ with respect to a set of perturbation { ^^p, ^^ = 1... ^^} as the average distance in behavior space Z between the original trajectory endpoint ^^ and the perturbed trajectories endpoints { ^^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 theTufts Ref. T002751 Quarles Ref.166118.01450 goal-level, providing several insights into the degree of “competencies” that some biological networks might exhibit in their relative space (FIG.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 (FIG.6).

[0137] FIG.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 exhibit 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” yet consistently reaching the same endpoint despite several pushes (FIG.5E) or walls (FIG.5G) on the way.

[0138] FIG.6 shows how the constructed catalog ^^ 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

[0114] ,

[0115] . It is believed that such system-level visualizations could be particularly useful to apprehend epigenetic heterogeneity in the context of cancer treatment and stem cell differentiation

[0107] ,

[0113] . A recent landscape-generation method only proposes to approximate the pseudopotential energy through simulation trajectories obtained throughout exploration of the system

[0113] , 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 (FIG.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 (FIG.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 instanceTufts Ref. T002751 Quarles Ref.166118.01450 leading to new shaped valleys (FIG.6C). For the example system RKIP-ERK pathway

[0099] , 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

[0116] .

[0139] Our framework generated a catalog of stimuli, responses, and navigation tests situations 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 dataset 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.

[0140] FIG.14 shows an example 1400 of a process for identifying a behavioral catalog in accordance with some embodiments of the disclosed subject matter. As shown in FIG.12, at 1402 process 1400 can obtain a data set of a GRN. The data set 1402 can be a mathematical model of the GRN. At 1404, process 1400 can identify at least one goal behavioral state and corresponding trajectory. The goal behavioral state can be a state at which the system reaches steady state. The corresponding trajectory can connect an initial state and a goal behavioral state. At 1406, process 1400 can test the robustness of one of more goal behavioral states. Testing the robustness can involve introducing perturbations during the experiment. The perturbations can be noise perturbations, push perturbations, or wall perturbations. At 1408, the process 1400 can select a subset of goal behavioral states. The subset can be goal behavioral states with a higher robustness than the average of all goal behavioral states. At 1410, the process 1400 can save the goal behavioral states. EXAMPLE 5: Develop insights on the degree of sophistication of the different GRNs

[0141] 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, homoTufts Ref. T002751 Quarles Ref.166118.01450 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 FIG.7, results suggested interesting patterns.

[0142] 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

[0042] : 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.

[0143] 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

[0079] , 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

[0117] -

[0120] , 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.

[0144] 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 FIG.7B, we plot the mean and std- 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 GRNsTufts Ref. T002751 Quarles Ref.166118.01450 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

[0121] -

[0124] . EXAMPLE 6: Development of therapeutic interventions

[0145] Understanding forms of non-genetic resistance and epigenetic heterogeneity is crucial across a wide range of cancer and treatment contexts

[0107] . 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

[0099] . In FIG.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

[0116] . In FIG.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 FIG.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 GRN setpoints from the identified “disease” steady states back to steady states within the identified “healthy” region.Tufts Ref. T002751 Quarles Ref.166118.01450

[0146] 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 during 100 seconds), and with what amplitude (which are the parameters that we are seeking to optimize). The choice of the intervention function, whereas 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 FIG.8B) that successfully reset the setpoints (as shown in FIG.8C) under various tested perturbations (as shown in FIG.11). 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.

[0147] FIG.15 shows an example 1500 of a process for identifying therapeutic interventions. As shown in FIG.15, at 1502, process 1500 can obtain a behavioral catalog of a gene regulatory network. The behavioral catalog can include behavioral end states, at which the system reaches steady state. The behavioral catalog can also include a history of what stimuli can produce a trajectory that ends in a behavioral end state. At 1504, process 1500 can select a disease associated with the GRN. At 1506, process 1500 can classify behavioral end states in the behavioral catalog, classifying each as a disease first setpoint or a healthy first setpoint. At 1508, process 1500 can define stimulus-based interventions. Stimulus based interventions can correspond with how a drug or genetic therapy would act on a GRN. At 1510, process 1500 can introduce stimulus-based interventions to at least one first setpoint in the behavioral catalog. InTufts Ref. T002751 Quarles Ref.166118.01450 some cases, the stimulus-based interventions are introduced to all disease first setpoints. At 1512, process 1500 can identify trajectories after stimulus-based interventions are introduced, and at 1514 the second setpoint of the trajectories can be determined. The second setpoint can be a state at which the system reaches steady state. At 1516, the process 1500 can select a subset of trajectories. In some cases, the subset of trajectories may be those in which the first setpoint is a disease setpoint, and the second setpoint is a healthy setpoint. At 1518, process 1500 can determine the stimulus-based interventions associated with the subset of trajectories, and at 1520 process 1500 can output the subset of stimulus-based interventions.

[0148] FIG.16 shows an example of a process for classifying healthy setpoints and disease setpoints. At 1602, process 1600 can construct an energy landscape of a behavioral catalog. Process 1600 can determine which regions of the energy landscape correspond with a disease region at 1604a, and can determine which regions of the energy landscape correspond with a healthy region at 1604b. Process 1600 can identify endpoints of the behavioral catalog that fall within a disease region and classify them as a disease setpoint 1606a and identify endpoints of the behavioral catalog that fall within a healthy region and classify them as a healthy setpoint. At 1608, process 1600 can select a subset of trajectories associated with a disease setpoint. In some cases, the subset of trajectories are trajectories with a higher than average robustness. At 1610, process 1600 can output a subset of trajectories and associated setpoints. EXAMPLE 7: Strategy for gene circuit engineering

[0149] 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

[0075] . 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-Tufts Ref. T002751 Quarles Ref.166118.01450 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

[0125] , by effectively exploring solution space.

[0150] Here, we consider the target application of oscillator circuit engineering followed in

[0072] , where parameters of a gene circuit model are optimized to produce oscillation patterns with target amplitude ^^, frequency ^^ and offset ^^. This time, the intervention space includes both genetic interventions (setting kinematic parameters of regulatory rules) and environmental interventions (setting the initial state ^^,). We then compare three alternative exploration strategies: curiosity search, random search and a global optimization strategy using gradient descent as proposed in

[0072] , all given the same experimental budget ( ^^ = 5000). For curiosity search, the behavior space ^^ is defined as the image space of the discrete Fourier transform of the observation. We then use the exact same IMGEP algorithm than before, but operating within the novel problem spaces ( ^^, ^^). For gradient descent, we follow the procedure proposed in

[0072] . We define a loss function which measures the mean square error between the observed node activation levels ^^ and the target oscillation (represented as a cosine wave). We then randomly initialize the parameters ^^ ∼ ^^( ^^) and use Adam optimizer for N=5000 optimization steps. 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 ^^ = 100 optimization steps.

[0151] In FIG.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 ( ^^, ^^, ^^) space (FIG.9A-9C), thus discovering oscillators close to the target one (FIG.9-E). Instead, when starting from a random initial condition, gradient descent is very likely to get trapped in a local minima where it converges to the target offset b but fails to produce oscillations (FIG.9D and FIG.9G). However, whereas the global optimization is unsuccessful in this example, gradient descent seems to be useful to locally refine close-enough solutions, as can be seen here whenTufts Ref. T002751 Quarles Ref.166118.01450 refining the best discoveries made by curiosity search and random search (FIG.9H, FIG.9I). 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.

[0152] FIG.17 shows an example process for identifying a configuration of parameter of a gene regulatory network to perform a desired function. At 1702, process 1700 can obtain a data set of a GRN. At 1704, process 1700 can identify goal behavioral states and corresponding trajectories. In some cases, a curiosity driven machine learning algorithm can be used to identify the goal behavioral states. In some cases, a goal behavioral state is a state in which the system reaches steady state. In some cases, a trajectory connects an initial state with a goal behavioral state. At 1706, process 1700 can use a local optimization strategy to identify additional behavioral states.

[0153] In FIG.18, an example 1800 of a system (e.g., a data processing system) for identifying and using behavioral catalogs is shown in accordance with some embodiments of the disclosed subject matter.

[0154] In some embodiments, computing device 1804 and / or server 1816 can be any suitable computing device or combination of devices, such as a desktop computer, a laptop computer, a smartphone, a tablet computer, a wearable computer, a server computer, a virtual machine being executed by a physical computing device, etc. As described herein, system 1800 can present information about the behavioral catalog or its applications to a user (e.g., a researcher and / or a physician).

[0001] In some embodiments, communication network 1802 can be any suitable communication network or combination of communication networks. In some embodiments, communication network 1006 can be any suitable communication network or combination of communication networks. For example, communication network 1006 can include a Wi-Fi network (which can include one or more wireless routers, one or more switches, etc.), a peer-to- peer network (e.g., a Bluetooth network), a cellular network (e.g., a 4G network, a 5G network,Tufts Ref. T002751 Quarles Ref.166118.01450 etc., complying with any suitable standard, such as CDMA, GSM, LTE, LTE Advanced, WiMAX, etc.), a wired network, etc. In some embodiments, communication network 1006 can be a local area network, a wide area network, a public network (e.g., the Internet), a private or semi-private network (e.g., a corporate or university intranet), any other suitable type of network, or any suitable combination of networks. Communications links shown in FIG.18 can each be any suitable communications link or combination of communications links, such as wired links, fiber optic links, Wi-Fi links, Bluetooth links, cellular links, etc.

[0155] FIG.18 additionally shows an example of hardware that can be used to implement computing device 1804 and server 1816 in accordance with some embodiments of the disclosed subject matter. In some embodiments, computing device 1802 can be used to execute one or more set of instructions to identify a behavioral catalog. In other embodiments, computing device 1802 can be used to identify therapeutic interventions. In still other embodiments, computing device 1802 can be used to identify a configuration of parameter of a gene regulatory network to perform a desired function.

[0156] As shown in FIG.18, computing device 1804 can include one or more hardware processor 1806, one or more displays 1808, one or more inputs 1810, one or more communications 1812, and / or memory 1814. In some embodiments, processor 1806 can be any suitable hardware processor or combination of processors, such as central processing unit, a graphics processing unit, etc. In some embodiments, display 1808 can include any suitable display devices, such as a computer monitor, a touchscreen, a television, etc. In some embodiments, inputs 18010 can include any suitable input device and / or sensors that can be used to receive user input, such as a keyboard, a mouse, a touchscreen, a microphone, etc.

[0157] In some embodiments, communications systems 1812 can include any suitable hardware, firmware, and / or software for communicating information over communication network 1802 and / or any other suitable communication networks. For example, communications systems 1812 can include one or more transceivers, one or more communication chips and / or chip sets, etc. In a more particular example, communications systems 1812 can include hardware, firmware and / or software that can be used to establish a Wi-Fi connection, a Bluetooth connection, a cellular connection, an Ethernet connection, etc.Tufts Ref. T002751 Quarles Ref.166118.01450

[0158] In some embodiments, memory 1814 can include any suitable storage device or devices that can be used to store instructions, values, etc., that can be used, for example, by processor 1806 to present content using display 1808, to communicate with server 1816 via communications system(s) 1812, etc.

[0159] Memory 1814 can include any suitable volatile memory, non-volatile memory, storage, or any suitable combination thereof. For example, memory 1814 can include RAM, ROM, EEPROM, one or more flash drives, one or more hard disks, one or more solid state drives, one or more optical drives, etc. In some embodiments, memory 1814 can have encoded thereon a computer program for controlling operation of computing device 1814. In such embodiments, processor 1804 can execute at least a portion of the computer program to present content (e.g., images, user interfaces, graphics, tables, etc.), receive content from server 1816, transmit information to server 1816, etc.

[0160] In some embodiments, server 1816 can include a processor 1818, a display 1820, one or more inputs 1822, one or more communications systems 1824, and / or memory 1826. In some embodiments, processor 1818 can be any suitable hardware processor or combination of processors, such as a central processing unit, a graphics processing unit, etc. In some embodiments, display 1820 can include any suitable display devices, such as a computer monitor, a touchscreen, a television, etc. In some embodiments, inputs 1822 can include any suitable input devices and / or sensors that can be used to receive user input, such as a keyboard, a mouse, a touchscreen, a microphone, etc.

[0161] In some embodiments, communications systems 1824 can include any suitable hardware, firmware, and / or software for communicating information over communication network 1802 and / or any other suitable communication networks. For example, communications systems 1824 can include one or more transceivers, one or more communication chips and / or chip sets, etc. In a more particular example, communications systems 1824 can include hardware, firmware and / or software that can be used to establish a Wi-Fi connection, a Bluetooth connection, a cellular connection, an Ethernet connection, etc.

[0162] In some embodiments, memory 1826 can include any suitable storage device or devices that can be used to store instructions, values, etc., that can be used, for example, by processor 1818 to present content using display 1820, to communicate with one or moreTufts Ref. T002751 Quarles Ref.166118.01450 computing devices 1804, etc. Memory 1826 can include any suitable volatile memory, non- volatile memory, storage, or any suitable combination thereof. For example, memory 1826 can include RAM, ROM, EEPROM, one or more flash drives, one or more hard disks, one or more solid state drives, one or more optical drives, etc. In some embodiments, memory 1826 can have encoded thereon a server program for controlling operation of server 1816. In such embodiments, processor 1818 can execute at least a portion of the server program to transmit information and / or content (e.g., results of a tissue identification and / or classification, a user interface, etc.) to one or more computing devices 1804, receive information and / or content from one or more computing devices 1804, receive instructions from one or more devices (e.g., a personal computer, a laptop computer, a tablet computer, a smartphone, etc.), etc.

[0163] In some embodiments, any suitable computer readable media can be used for storing instructions for performing the functions and / or processes described herein. For example, in some embodiments, computer readable media can be transitory or non-transitory. For example, non-transitory computer readable media can include media such as magnetic media (such as hard disks, floppy disks, etc.), optical media (such as compact discs, digital video discs, Blu-ray discs, etc.), semiconductor media (such as RAM, Flash memory, electrically programmable read only memory (EPROM), electrically erasable programmable read only memory (EEPROM), etc.), any suitable media that is not fleeting or devoid of any semblance of permanence during transmission, and / or any suitable tangible media. As another example, transitory computer readable media can include signals on networks, in wires, conductors, optical fibers, circuits, or any suitable media that is fleeting and devoid of any semblance of permanence during transmission, and / or any suitable intangible media.

[0164] Citations to a number of patent and non-patent references may be made herein. The cited references are incorporated by reference herein in their entireties. In the event that there is an inconsistency between a definition of a term in the specification as compared to a definition of the term in a cited reference, the term should be interpreted based on the definition in the specification.Tufts Ref. T002751 Quarles Ref.166118.01450

[0165] References

[0166] [1] J. J. Sanz-Ezquerro, A. E. Münsterberg, and S. Stricker, “Editorial: Signaling Pathways in Embryonic Development,” Front. Cell Dev. Biol., vol.5, 2017, Accessed: Jun.12, 2023. [Online]. Available: frontiersin.org / articles / 10.3389 / fcell.2017.00076

[0167] [2] E. R. A.-B. Padilla-Longoria Enrique Balleza, Mariana Benítez, Carlos Espinosa-Soto, Pablo, “Gene regulatory network models: A dynamic and integrative approach to development,” in Practical Systems Biology, Taylor & Francis, 2008.

[0168] [3] S. Huang, G. Eichler, Y. Bar-Yam, and D. E. Ingber, “Cell Fates as High- Dimensional Attractor States of a Complex Gene Regulatory Network,” Phys. Rev. Lett., vol.94, no.12, p.128701, Apr.2005, doi: 10.1103 / PhysRevLett.94.128701.

[0169] [4] E. H. Davidson, “Emerging properties of animal gene regulatory networks,” Nature, vol.468, no.7326, pp.911–920, Dec.2010, doi: 10.1038 / nature09645.

[0170] [5] I. S. Peter and E. H. Davidson, “Evolution of Gene Regulatory Networks Controlling Body Plan Development,” Cell, vol.144, no.6, pp.970–985, Mar.2011, doi: 10.1016 / j.cell.2011.02.017.

[0171] [6] K. H. ten Tusscher and P. Hogeweg, “Evolution of Networks for Body Plan Patterning; Interplay of Modularity, Robustness and Evolvability,” PLOS Comput. Biol., vol.7, no.10, p. e1002208, Oct.2011, doi: 10.1371 / journal.pcbi.1002208.

[0172] [7] H. Kim and H. Sayama, “How Criticality of Gene Regulatory Networks Affects the Resulting Morphogenesis under Genetic Perturbations,” Artif. Life, vol.24, no.02, pp.85–105, May 2018, doi: 10.1162 / ARTL_a_00262.

[0173] [8] M. Srivastava, “Beyond Casual Resemblance: Rigorous Frameworks for Comparing Regeneration Across Species,” Annu. Rev. Cell Dev. Biol., vol.37, no.1, pp.415– 440, 2021, doi: 10.1146 / annurev-cellbio-120319-114716.

[0174] [9] A. J. Singh, S. A. Ramsey, T. M. Filtz, and C. Kioussi, “Differential gene regulatory networks in development and disease,” Cell. Mol. Life Sci., vol.75, no.6, pp.1013– 1025, Mar.2018, doi: 10.1007 / s00018-017-2679-6.Tufts Ref. T002751 Quarles Ref.166118.01450

[0175]

[0010] G. Qin, L. Yang, Y. Ma, J. Liu, and Q. Huo, “The exploration of disease- specific gene regulatory networks in esophageal carcinoma and stomach adenocarcinoma,” BMC Bioinformatics, vol.20, no.22, p.717, Dec.2019, doi: 10.1186 / s12859-019-3230-6.

[0176]

[0011] H. Fazilaty et al., “A gene regulatory network to control EMT programs in development and disease,” Nat. Commun., vol.10, no.1, Art. no.1, Nov.2019, doi: 0.1038 / s41467-019-13091-8.

[0177]

[0012] J. Davies and M. Levin, “Synthetic morphology via active and agential matter.” OSF Preprints, Jun.10, 2022. doi: 10.31219 / osf.io / xrv8h.

[0178]

[0013] S. Toda, L. R. Blauch, S. K. Y. Tang, L. Morsut, and W. A. Lim, “Programming self-organizing multicellular structures with synthetic cell-cell signaling,” Science, vol.361, no.6398, pp.156–162, Jul.2018, doi: 10.1126 / science.aat0271.

[0179]

[0014] S. Toda, W. L. McKeithan, T. J. Hakkinen, P. Lopez, O. D. Klein, and W. A. Lim, “Engineering synthetic morphogen systems that can program multicellular patterning,” Science, vol.370, no.6514, pp.327–331, Oct.2020, doi: 10.1126 / science.abc0033.

[0180]

[0015] C. Ho and L. Morsut, “Novel synthetic biology approaches for developmental systems,” Stem Cell Rep., vol.16, no.5, pp.1051–1064, May 2021, doi: 10.1016 / j.stemcr.2021.04.007.

[0181]

[0016] M. Santorelli, C. Lam, and L. Morsut, “Synthetic development: building mammalian multicellular structures with artificial genetic programs,” Curr. Opin. Biotechnol., vol.59, pp.130–140, Oct.2019, doi: 10.1016 / j.copbio.2019.03.016.

[0182]

[0017] H. de Jong, “Modeling and Simulation of Genetic Regulatory Systems: A Literature Review,” J. Comput. Biol., vol.9, no.1, pp.67–103, Jan.2002, doi: 10.1089 / 10665270252833208.

[0183]

[0018] T. Schlitt and A. Brazma, “Current approaches to gene regulatory network modelling,” BMC Bioinformatics, vol.8, no.6, p. S9, Sep.2007, doi: 10.1186 / 1471-2105-8-S6- S9.Tufts Ref. T002751 Quarles Ref.166118.01450

[0184]

[0019] J. S. Fetrow and P. C. Babbitt, “New computational approaches to understanding molecular protein function,” PLOS Comput. Biol., vol.14, no.4, p. e1005756, Apr.2018, doi: 10.1371 / journal.pcbi.1005756.

[0185]

[0020] F. M. Delgado and F. Gómez-Vela, “Computational methods for Gene Regulatory Networks reconstruction and analysis: A review,” Artif. Intell. Med., vol.95, pp. 133–145, Apr.2019, doi: 10.1016 / j.artmed.2018.10.006.

[0186]

[0021] M. Glont et al., “BioModels: expanding horizons to include more modelling approaches and formats,” Nucleic Acids Res., vol.46, no. D1, pp. D1248–D1253, Jan.2018, doi: 10.1093 / nar / gkx1023.

[0187]

[0022] R. S. Malik-Sheriff et al., “BioModels—15 years of sharing computational models in life science,” Nucleic Acids Res., vol.48, no. D1, pp. D407–D415, Jan.2020, doi: 10.1093 / nar / gkz1055.

[0188]

[0023] S. A. Kauffman, The origins of order: Self-organization and selection in evolution. Oxford University Press, USA, 1993.

[0189]

[0024] S. A. Kauffman, At home in the universe: The search for laws of self- organization and complexity. Oxford University Press, USA, 1995.

[0190]

[0025] C. I. Abramson and M. Levin, “Behaviorist approaches to investigating memory and learning: A primer for synthetic biology and bioengineering,” Commun. Integr. Biol., vol.4, no.11, pp.230–247, Jan.2021, doi: 10.1080 / 19420889.2021.2005863.

[0191]

[0026] F. Baluška and M. Levin, “On Having No Head: Cognition throughout Biological Systems,” Front. Psychol., vol.7, 2016, Accessed: Jun.12, 2023. [Online]. Available: frontiersin.org / articles / 10.3389 / fpsyg.2016.00902

[0192]

[0027] G. Dodig-Crnkovic, “Cognition as Morphological / Morphogenetic Embodied Computation In Vivo,” Entropy Basel Switz., vol.24, no.11, p.1576, Oct.2022, doi: 10.3390 / e24111576.

[0193]

[0028] Y. Timsit and S.-P. Grégoire, “Towards the Idea of Molecular Brains,” Int. J. Mol. Sci., vol.22, no.21, Art. no.21, Jan.2021, doi: 10.3390 / ijms222111868.Tufts Ref. T002751 Quarles Ref.166118.01450

[0194]

[0029] Y. Katz, M. Springer, and W. Fontana, “Embodying probabilistic inference in biochemical circuits.” arXiv, Jun.26, 2018. doi: 10.48550 / arXiv.1806.10161.

[0195]

[0030] P. Csermely et al., “Learning of Signaling Networks: Molecular Mechanisms,” Trends Biochem. Sci., vol.45, no.4, pp.284–294, Apr.2020, doi: 10.1016 / j.tibs.2019.12.005.

[0196]

[0031] D. M. Gyurkó, D. V. Veres, D. Módos, K. Lenti, T. Korcsmáros, and P. Csermely, “Adaptation and learning of molecular networks as a description of cancer development at the systems-level: Potential use in anti-cancer therapies,” Semin. Cancer Biol., vol.23, no.4, pp.262–269, Aug.2013, doi: 10.1016 / j.semcancer.2013.06.005.

[0197]

[0032] C. Fields and M. Levin, “Competency in Navigating Arbitrary Spaces: Intelligence as an Invariant for Analyzing Cognition in Diverse Embodiments.” PsyArXiv, Apr. 24, 2022. doi: 10.31234 / osf.io / 87nzu.

[0198]

[0033] R. Watson, C. L. Buckley, R. Mills, and A. Davies, “Associative memory in gene regulation networks,” H. Fellerman, M. Dörr, M. M. Hanczyc, L. Ladegaard Laursen, S. Maurer, D. Merkle, P.-A. Monnard, K. Stoy, and S. Rasmussen, Eds., MIT Press, 2010, pp.659– 666. Accessed: Sep.05, 2023. [Online]. Available: eprints.soton.ac.uk / 339763 /

[0199]

[0034] J. Mathews, A. (Jaelyn) Chang, L. Devlin, and M. Levin, “Cellular signaling pathways as plastic, proto-cognitive systems: Implications for biomedicine,” Patterns, vol.4, no. 5, p.100737, May 2023, doi: 10.1016 / j.patter.2023.100737.

[0200]

[0035] E. Lagasse and M. Levin, “Future medicine: from molecular pathways to the collective intelligence of the body,” Trends Mol. Med., vol.29, no.9, pp.687–710, Sep.2023, doi: 10.1016 / j.molmed.2023.06.007.

[0201]

[0036] W. P. Clawson and M. Levin, “Endless forms most beautiful 2.0: teleonomy and the bioengineering of chimaeric and synthetic organisms,” Biol. J. Linn. Soc., p. blac073, Jul.2022, doi: 10.1093 / biolinnean / blac073.

[0202]

[0037] K. T. Krist, A. Sen, and W. G. Noid, “A simple theory for molecular chemotaxis driven by specific binding interactions,” J. Chem. Phys., vol.155, no.16, p.164902, Oct.2021, doi: 10.1063 / 5.0061376.Tufts Ref. T002751 Quarles Ref.166118.01450

[0203]

[0038] J. Čejková, T. Banno, M. M. Hanczyc, and F. Štěpánek, “Droplets As Liquid Robots,” Artif. Life, vol.23, no.4, pp.528–549, Nov.2017, doi: 10.1162 / ARTL_a_00243.

[0204]

[0039] M. M. Hanczyc, F. Caschera, and S. Rasmussen, “Models of Minimal Physical Intelligence,” Procedia Comput. Sci., vol.7, pp.275–277, Jan.2011, doi: 10.1016 / j.procs.2011.09.058.

[0205]

[0040] A. Rosenblueth, N. Wiener, and J. Bigelow, “Behavior, Purpose and Teleology,” Philos. Sci., vol.10, no.1, pp.18–24, Jan.1943, doi: 10.1086 / 286788.

[0206]

[0041] J. Bongard and M. Levin, “Living Things Are Not (20th Century) Machines: Updating Mechanism Metaphors in Light of the Modern Science of Machine Behavior,” Front. Ecol. Evol., vol.9, 2021, Accessed: Jun.12, 2023. [Online]. Available: www.frontiersin.org / articles / 10.3389 / fevo.2021.650726

[0207]

[0042] M. Levin, “Technological Approach to Mind Everywhere: An Experimentally-Grounded Framework for Understanding Diverse Bodies and Minds,” Front. Syst. Neurosci., vol.16, 2022, Accessed: Jun.12, 2023. [Online]. Available: frontiersin.org / articles / 10.3389 / fnsys.2022.768201

[0208]

[0043] P. Lyon, “The biogenic approach to cognition,” Cogn. Process., vol.7, no.1, pp.11–29, Mar.2006, doi: 10.1007 / s10339-005-0016-8.

[0209]

[0044] X. Barandiaran and A. Moreno, “On What Makes Certain Dynamical Systems Cognitive: A Minimally Cognitive Organization Program,” Adapt. Behav., vol.14, no. 2, pp.171–185, Jun.2006, doi: 10.1177 / 105971230601400208.

[0210]

[0045] F. di Primio, B. S. Müller, and J. W. Lengeler, “Minimal cognition in unicellular organisms,” Anim. Animats, pp.3–12, 2000.

[0211]

[0046] P. McGivern, “Active materials: minimal models of cognition?,” Adapt. Behav., vol.28, no.6, pp.441–451, Dec.2020, doi: 10.1177 / 1059712319891742.

[0212]

[0047] M. Levin, “Darwin’s agential materials: evolutionary implications of multiscale competency in developmental biology,” Cell. Mol. Life Sci., vol.80, no.6, p.142, May 2023, doi: 10.1007 / s00018-023-04790-z.Tufts Ref. T002751 Quarles Ref.166118.01450

[0213]

[0048] G. Pezzulo and M. Levin, “Re-membering the body: applications of computational neuroscience to the top-down control of regeneration of limbs and other complex organs,” Integr. Biol. Quant. Biosci. Nano Macro, vol.7, no.12, pp.1487–1517, Dec.2015, doi: 10.1039 / c5ib00221d.

[0214]

[0049] G. Pezzulo and M. Levin, “Top-down models in biology: explanation and control of complex living systems above the molecular level,” J. R. Soc. Interface, vol.13, no. 124, p.20160555, Nov.2016, doi: 10.1098 / rsif.2016.0555.

[0215]

[0050] D. J. Wong et al., “Revealing Targeted Therapy for Human Cancer by Gene Module Maps,” Cancer Res., vol.68, no.2, pp.369–378, Jan.2008, doi: 10.1158 / 0008- 5472.CAN-07-0382.

[0216]

[0051] T. J. Samuel, R. P. Rosenberry, S. Lee, and Z. Pan, “Correcting Calcium Dysregulation in Chronic Heart Failure Using SERCA2a Gene Therapy,” Int. J. Mol. Sci., vol. 19, no.4, Art. no.4, Apr.2018, doi: 10.3390 / ijms19041086.

[0217]

[0052] R. Krzysztoń, Y. Wan, J. Petreczky, and G. Balázsi, “Gene-circuit therapy on the horizon: Synthetic biology tools for engineered therapeutics,” Acta Biochim. Pol., vol.68, no.3, pp.377–383, Aug.2021, doi: 10.18388 / abp.2020_5744.

[0218]

[0053] C. Baum, “Insertional mutagenesis in gene therapy and stem cell biology,” Curr. Opin. Hematol., vol.14, no.4, p.337, Jul.2007, doi: 10.1097 / MOH.0b013e3281900f01.

[0219]

[0054] D. Lobo, M. Solano, G. A. Bubenik, and M. Levin, “A linear-encoding model explains the variability of the target morphology in regeneration,” J. R. Soc. Interface, vol. 11, no.92, p.20130918, Mar.2014, doi: 10.1098 / rsif.2013.0918.

[0220]

[0055] J. W. Stucki, “Stability analysis of biochemical systems— A practical guide,” Prog. Biophys. Mol. Biol., vol.33, pp.99–187, Jan.1979, doi: 10.1016 / 0079- 6107(79)90027-0.

[0221]

[0056] B. P. Ingalls, “A Frequency Domain Approach to Sensitivity Analysis of Biochemical Networks,” J. Phys. Chem. B, vol.108, no.3, pp.1143–1152, Jan.2004, doi: 10.1021 / jp036567u.Tufts Ref. T002751 Quarles Ref.166118.01450

[0222]

[0057] B. Ingalls, “Sensitivity analysis: from model parameters to system behaviour,” Essays Biochem., vol.45, pp.177–194, Sep.2008, doi: 10.1042 / bse0450177.

[0223]

[0058] A. Donzé, G. Clermont, and C. J. Langmead, “Parameter synthesis in nonlinear dynamical systems: application to systems biology,” J. Comput. Biol. J. Comput. Mol. Cell Biol., vol.17, no.3, pp.325–336, Mar.2010, doi: 10.1089 / cmb.2009.0172.

[0224]

[0059] T. Dang, C. Le Guernic, and O. Maler, “Computing reachable states for nonlinear biological models,” Theor. Comput. Sci., vol.412, no.21, pp.2095–2107, May 2011, doi:10.1016 / j.tcs.2011.01.014.

[0225]

[0060] A. Donzé, E. Fanchon, L. M. Gattepaille, O. Maler, and P. Tracqui, “Robustness Analysis and Behavior Discrimination in Enzymatic Reaction Networks,” PLOS ONE, vol.6, no.9, p. e24246, Sep.2011, doi: 10.1371 / journal.pone.0024246.

[0226]

[0061] J. Rozum and R. Albert, “Leveraging network structure in nonlinear control,” Npj Syst. Biol. Appl., vol.8, no.1, Art. no.1, Oct.2022, doi: 10.1038 / s41540-022-00249-2.

[0227]

[0062] S. N. Steinway, J. G. T. Zañudo, P. J. Michel, D. J. Feith, T. P. Loughran, and R. Albert, “Combinatorial interventions inhibit TGFβ-driven epithelial-to-mesenchymal transition and support hybrid cellular phenotypes,” Npj Syst. Biol. Appl., vol.1, no.1, Art. no.1, Nov.2015, doi: 10.1038 / npjsba.2015.14.

[0228]

[0063] J. G. T. Zañudo and R. Albert, “Cell Fate Reprogramming by Control of Intracellular Network Dynamics,” PLOS Comput. Biol., vol.11, no.4, p. e1004193, Apr.2015, doi: 10.1371 / journal.pcbi.1004193.

[0229]

[0064] J. G. T. Zañudo, G. Yang, and R. Albert, “Structure-based control of complex networks with nonlinear dynamics,” Proc. Natl. Acad. Sci., vol.114, no.28, pp.7234– 7239, Jul.2017, doi: 10.1073 / pnas.1617387114.

[0230]

[0065] L. Cifuentes Fontanals, E. Tonello, and H. Siebert, “Control Strategy Identification via Trap Spaces in Boolean Networks,” in Computational Methods in Systems Biology, A. Abate, T. Petrov, and V. Wolf, Eds., in Lecture Notes in Computer Science. Cham: Springer International Publishing, 2020, pp.159–175. doi: 10.1007 / 978-3-030-60327-4_9.Tufts Ref. T002751 Quarles Ref.166118.01450

[0231]

[0066] D. Murrugarra, A. Veliz-Cuba, B. Aguilar, and R. Laubenbacher, “Identification of control targets in Boolean molecular network models via computational algebra,” BMC Syst. Biol., vol.10, no.1, p.94, Sep.2016, doi: 10.1186 / s12918-016-0332-x.

[0232]

[0067] S.-M. Choo, B. Ban, J. I. Joo, and K.-H. Cho, “The phenotype control kernel of a biomolecular regulatory network,” BMC Syst. Biol., vol.12, no.1, p.49, Apr.2018, doi: 10.1186 / s12918-018-0576-8.

[0233]

[0068] S.-M. Choo, S.-M. Park, and K.-H. Cho, “Minimal intervening control of biomolecular networks leading to a desired cellular state,” Sci. Rep., vol.9, no.1, Art. no.1, Sep.2019, doi: 10.1038 / s41598-019-49571-6.

[0234]

[0069] S. R. Paladugu, V. Chickarmane, A. Deckard, J. P. Frumkin, M. McCormack, and H. M.Sauro, “In silico evolution of functional modules in biochemical networks,” IEE Proc. -Syst. Biol., vol.153, no.4, pp.223–235, Jul.2006, doi: 10.1049 / ip- syb:20050096.

[0235]

[0070] P. François, “Evolving phenotypic networks in silico,” Semin. Cell Dev. Biol., vol.35, pp.90–97, Nov.2014, doi: 10.1016 / j.semcdb.2014.06.012.

[0236]

[0071] N. Noman, T. Monjo, P. Moscato, and H. Iba, “Evolving Robust Gene Regulatory Networks,” PLOS ONE, vol.10, no.1, p. e0116258, Jan.2015, doi: 10.1371 / journal.pone.0116258.

[0237]

[0072] T. W. Hiscock, “Adapting machine-learning algorithms to design gene circuits,” BMC Bioinformatics, vol.20, no.1, p.214, Apr.2019, doi: 10.1186 / s12859-019-2788- 3.

[0238]

[0073] J. Shen, F. Liu, Y. Tu, and C. Tang, “Finding gene network topologies for given biological function with recurrent neural network,” Nat. Commun., vol.12, no.1, Art. no. 1, May 2021, doi: 10.1038 / s41467-021-23420-5.

[0239]

[0074] D. M. Camacho, K. M. Collins, R. K. Powers, J. C. Costello, and J. J. Collins, “Next-Generation Machine Learning for Biological Networks,” Cell, vol.173, no.7, pp. 1581–1592, Jun.2018, doi: 10.1016 / j.cell.2018.05.015.Tufts Ref. T002751 Quarles Ref.166118.01450

[0240]

[0075] M. J. Volk, I. Lourentzou, S. Mishra, L. T. Vo, C. Zhai, and H. Zhao, “Biosystems Design by Machine Learning,” ACS Synth. Biol., vol.9, no.7, pp.1514–1533, Jul. 2020, doi: 10.1021 / acssynbio.0c00129.

[0241]

[0076] H. Kitano, “A robustness-based approach to systems-oriented drug design,” Nat. Rev. Drug Discov., vol.6, no.3, Art. no.3, Mar.2007, doi: 10.1038 / nrd2195.

[0242]

[0077] S. Manicka and M. Levin, “The Cognitive Lens: a primer on conceptual tools for analysing information processing in developmental and regenerative morphogenesis,” Philos. Trans. R. Soc. B Biol. Sci., vol.374, no.1774, p.20180369, Apr.2019, doi: 10.1098 / rstb.2018.0369.

[0243]

[0078] S. Biswas, S. Manicka, E. Hoel, and M. Levin, “Gene regulatory networks exhibit several kinds of memory: Quantification of memory in biological and random transcriptional networks,” iScience, vol.24, no.3, p.102131, Mar.2021, doi: 10.1016 / j.isci.2021.102131.

[0244]

[0079] S. Biswas, W. Clawson, and M. Levin, “Learning in Transcriptional Network Models: Computational Discovery of Pathway-Level Memory and Effective Interventions,” Int. J. Mol. Sci., vol.24, no.1, Art. no.1, Jan.2023, doi: 10.3390 / ijms24010285.

[0245]

[0080] V. G. Laties, “Society for the Experimental Analysis of Behavior: The first thirty years (1957–1987),” J. Exp. Anal. Behav., vol.48, no.3, pp.495–512, Nov.1987, doi: 10.1901 / jeab.1987.48-495.

[0246]

[0081] G. Amdam and A. Hovland, “Measuring Animal Preferences and Choice Behavior,” Nat. Educ. Knowl., vol.3, p.74, Jan.2012.

[0247]

[0082] S. McLeold, “Behavioral Perspective in Psychology [Behaviorism Theory],” Nov.03, 2022. simplypsychology.org / behaviorism.html (accessed Jun.16, 2023).

[0248]

[0083] N. J. Murugan et al., “Mechanosensation Mediates Long-Range Spatial Decision-Making in an Aneural Organism,” Adv. Mater. Deerfield Beach Fla, vol.33, no.34, p.e2008161, Aug.2021, doi: 10.1002 / adma.202008161.

[0249]

[0084] A. Mikhaltsov, Paramecium bursaria.2013. [Online]. Available: commons.wikimedia.org / wiki / File:Paramecium_bursaria.jpgTufts Ref. T002751 Quarles Ref.166118.01450

[0250]

[0085] J. Bongard and M. Levin, “There’s Plenty of Room Right Here: Biological Systems as Evolved, Overloaded, Multi-Scale Machines,” Biomimetics, vol.8, no.1, Art. no.1, Mar.2023, doi: 10.3390 / biomimetics8010110.

[0251]

[0086] A. Baranes and P.-Y. Oudeyer, “Active learning of inverse models with intrinsically motivated goal exploration in robots,” Robot. Auton. Syst., vol.61, no.1, pp.49–73, Jan.2013, doi: 10.1016 / j.robot.2012.05.008.

[0252]

[0087] S. Forestier, R. Portelas, Y. Mollard, and P.-Y. Oudeyer, “Intrinsically Motivated Goal Exploration Processes with Automatic Curriculum Learning.” arXiv, May 05, 2022. Accessed: Apr.21, 2023. [Online]. Available: arxiv.org / abs / 1708.02190

[0253]

[0088] J. Grizou, L. J. Points, A. Sharma, and L. Cronin, “A curious formulation robot enables the discovery of a novel protocell behavior,” Sci. Adv., vol.6, no.5, p. eaay4237, Jan.2020, doi: 10.1126 / sciadv.aay4237.

[0254]

[0089] M. J. Falk, F. D. Roach, W. Gilpin, and A. Murugan, “Curiosity-driven search for novel non-equilibrium behaviors.” arXiv, Mar.17, 2023. Accessed: Apr.28, 2023. [Online]. available: arxiv.org / abs / 2211.02589

[0255]

[0090] C. Reinke, M. Etcheverry, and P.-Y. Oudeyer, “Intrinsically Motivated Discovery of Diverse Patterns in Self-Organizing Systems,” presented at the Eighth International Conference on Learning Representations, Apr.2020. Accessed: Apr.21, 2023. [Online]. Available: iclr.cc / virtual_2020 / poster_rkg6sJHYDr.html

[0256]

[0091] M. Etcheverry, C. Moulin-Frier, and P.-Y. Oudeyer, “Hierarchically Organized Latent Modules for Exploratory Search in Morphogenetic Systems,” in Advances in Neural Information Processing Systems, Curran Associates, Inc., 2020, pp.4846–4859. Accessed: Apr.21, 2023. [Online]. Available: proceedings.neurips.cc / paper / 2020 / hash / 33a5435d4f945aa6154b31a73bab3b73-Abstract.html

[0257]

[0092] G. Hamon, M. Etcheverry, B. W.-C. Chan, C. Moulin-Frier, and P.-Y. Oudeyer, “Learning Sensorimotor Agency in Cellular Automata,” 2022, Accessed: Apr.23, 2023. [Online]. Available: inria.hal.science / hal-03519319Tufts Ref. T002751 Quarles Ref.166118.01450

[0258]

[0093] H. Kitano, “Towards a theory of biological robustness,” Mol. Syst. Biol., vol. 3, no.1, p.137, Jan.2007, doi: 10.1038 / msb4100179.

[0259]

[0094] M.-A. Félix and M. Barkoulas, “Pervasive robustness in biological systems,” Nat. Rev. Genet., vol.16, no.8, Art. no.8, Aug.2015, doi: 10.1038 / nrg3949.

[0260]

[0095] N. T. Ingolia, “Topology and Robustness in the Drosophila Segment Polarity Network,” PLOS Biol., vol.2, no.6, p. e123, Jun.2004, doi: 10.1371 / journal.pbio.0020123.

[0261]

[0096] W. Ma, L. Lai, Q. Ouyang, and C. Tang, “Robustness and modular design of the Drosophila segment polarity network,” Mol. Syst. Biol., vol.2, no.1, p.70, Jan.2006, doi: 10.1038 / msb4100111.

[0262]

[0097] D. Deutscher, I. Meilijson, M. Kupiec, and E. Ruppin, “Multiple knockout analysis of genetic robustness in the yeast metabolic network,” Nat. Genet., vol.38, no.9, Art. no.9, Sep.2006, doi: 10.1038 / ng1856.

[0263]

[0098] G. von Dassow, E. Meir, E. M. Munro, and G. M. Odell, “The segment polarity network is a robust developmental module,” Nature, vol.406, no.6792, Art. no.6792, Jul.2000, doi: 10.1038 / 35018085.

[0264]

[0099] C. Kwang-Hyun, S. Sung-Young, K. Hyun-Woo, O. Wolkenhauer, B. McFerran, and W. Kolch, “Mathematical Modeling of the Influence of RKIP on the ERK Signaling Pathway,” in Computational Methods in Systems Biology, C. Priami, Ed., in Lecture Notes in Computer Science, vol.2602. Berlin, Heidelberg: Springer Berlin Heidelberg, 2003, pp. 127–141. doi: 10.1007 / 3-540-36481-1_11.

[0265]

[0100] J. Lehman and K. O. Stanley, “Exploiting Open-Endedness to Solve Problems Through the Search for Novelty,” presented at the IEEE Symposium on Artificial Life, 2008. Accessed: Jul.31, 2023. [Online]. Available: semanticscholar.org / paper / Exploiting-Open- Endedness-to-Solve- Problems-Through-Lehma Stanley / fb144a1d31aec3b2bece6a59bd11a876a9fafb34

[0266]

[0101] J. Lehman and K. O. Stanley, “Abandoning Objectives: Evolution Through the Search for Novelty Alone,” Evol. Comput., vol.19, no.2, pp.189–223, Jun.2011, doi: 10.1162 / EVCO_a_00025.Tufts Ref. T002751 Quarles Ref.166118.01450

[0267]

[0102] A. Cully, J. Clune, D. Tarapore, and J.-B. Mouret, “Robots that can adapt like animals,” Nature, vol.521, no.7553, Art. no.7553, May 2015, doi: 10.1038 / nature14422.

[0268]

[0103] J. K. Pugh, L. B. Soros, and K. O. Stanley, “Quality Diversity: A New Frontier for Evolutionary Computation,” Front. Robot. AI, vol.3, 2016, Accessed: Jun.13, 2023. [Online]. Available frontiersin.org / articles / 10.3389 / frobt.2016.00040

[0269]

[0104] S. Doncieux, A. Laflaquière, and A. Coninx, “Novelty search: a theoretical perspective,” in Proceedings of the Genetic and Evolutionary Computation Conference, Prague Czech Republic: ACM, Jul.2019, pp.99–106. Doi 10.1145 / 3321707.3321752.

[0270]

[0105] F. Benureau, “Self Exploration of Sensorimotor Spaces in Robots.”.

[0271]

[0106] L. McInnes, J. Healy, and S. Astels, “hdbscan: Hierarchical density based clustering,” J. Open Source Softw., vol.2, no.11, p.205, 2017.

[0272]

[0107] C. C. Bell and O. Gilan, “Principles and mechanisms of non-genetic resistance in cancer,” Br. J. Cancer, vol.122, no.4, Art. no.4, Feb.2020, doi: 10.1038 / s41416- 019-0648-6.

[0273]

[0108] A. Rizk, G. Batt, F. Fages, and S. Soliman, “A general computational method for robustness analysis with applications to synthetic gene networks,” Bioinformatics, vol.25, no.12, pp. i169–i178, Jun.2009, doi: 10.1093 / bioinformatics / btp200.

[0274]

[0109] C. Walcott, “Pigeon Homing: Observations, Experiments and Confusions,” J. Exp. Biol., vol.199, no.1, pp.21–27, Jan.1996, doi: 10.1242 / jeb.199.1.21.

[0275]

[0110] P. Luschi et al., “Testing the Navigational Abilities of Ocean Migrants: Displacement Experiments on Green Sea Turtles (Chelonia mydas),” Behav. Ecol. Sociobiol., vol.50, no.6, pp.528–534, 2001.

[0276]

[0111] S. Bisch-Knaden and R. Wehner, “Egocentric information helps desert ants to navigate around familiar obstacles,” J. Exp. Biol., vol.204, no.24, pp.4177–4184, Dec.2001, doi: 10.1242 / jeb.204.24.4177.

[0277]

[0112] C. I. Abramson and others, A primer of invertebrate learning: the behavioral perspective. American Psychological Association, 1994.Tufts Ref. T002751 Quarles Ref.166118.01450

[0278]

[0113] H. Venkatachalapathy, S. M. Azarin, and C. A. Sarkar, “Trajectory-based energy landscapes of gene regulatory networks,” Biophys. J., vol.120, no.4, pp.687–698, Feb. 2021, doi: 10.1016 / j.bpj.2020.11.2279.

[0279]

[0114] C. Li and J. Wang, “Landscape and flux reveal a new global view and physical quantification of mammalian cell cycle,” Proc. Natl. Acad. Sci., vol.111, no.39, pp. 14130–14135, Sep.2014, doi: 10.1073 / pnas.1408628111.

[0280]

[0115] C. Li and J. Wang, “Quantifying Cell Fate Decisions for Differentiation and Reprogramming of a Human Stem Cell Network: Landscape and Biological Paths,” PLOS Comput. Biol., vol.9, no.8, p. e1003165, Aug.2013, doi 10.1371 / journal.pcbi.1003165.

[0281]

[0116] H. C. Lee, B. Tian, J. M. Sedivy, J. R. Wands, and M. Kim, “Loss of Raf Kinase Inhibitor Protein Promotes Cell Proliferation and Migration of Human Hepatoma Cells,” Gastroenterology, vol.131, no.4, pp.1208–1217, Oct.2006, doi: 10.1053 / j.gastro.2006.07.012.

[0282]

[0117] J. Reinitz and D. H. Sharp, “Mechanism of eve stripe formation,” Mech. Dev., vol.49, no.1, pp.133–158, Jan.1995, doi: 10.1016 / 0925-4773(94)00310-J.

[0283]

[0118] J. Jaeger et al., “Dynamical Analysis of Regulatory Interactions in the Gap Gene System of Drosophila melanogaster,” Genetics, vol.167, no.4, pp.1721–1737, Aug.2004, doi: 10.1534 / genetics.104.027334.

[0284]

[0119] J. Cotterell and J. Sharpe, “An atlas of gene regulatory networks reveals multiple three-gene mechanisms for interpreting morphogen gradients,” Mol. Syst. Biol., vol.6, no.1, p.425, Jan.2010, doi: 10.1038 / msb.2010.74.

[0285]

[0120] E. J. Molinelli et al., “Perturbation Biology: Inferring Signaling Networks in Cellular Systems,” PLOS Comput. Biol., vol.9, no.12, p. e1003290, Dec.2013, doi: 10.1371 / journal.pcbi.1003290.

[0286]

[0121] J. Vallverdú et al., “Slime mould: The fundamental mechanisms of biological cognition,” Biosystems, vol.165, pp.57–70, Mar.2018, doi: 10.1016 / j.biosystems.2017.12.011.Tufts Ref. T002751 Quarles Ref.166118.01450

[0287]

[0122] M. Beekman and T. Latty, “Brainless but Multi-Headed: Decision Making by then Acellular Slime Mould Physarum polycephalum,” J. Mol. Biol., vol.427, no.23, pp. 3734– 3743, Nov.2015, doi: 10.1016 / j.jmb.2015.07.007.

[0288]

[0123] T. Saigusa, A. Tero, T. Nakagaki, and Y. Kuramoto, “Amoebae Anticipate Periodic Events,” Phys. Rev. Lett., vol.100, no.1, p.018101, Jan.2008, doi: 10.1103 / PhysRevLett.100.018101.

[0289]

[0124] T. Nakagaki and R. D. Guy, “Intelligent behaviors of amoeboid movement based on complex dynamics of soft matter,” Soft Matter, vol.4, no.1, pp.57–67, 2008, doi: 10.1039 / B706317M.

[0290]

[0125] C. Colas, O. Sigaud, and P.-Y. Oudeyer, “GEP-PG: Decoupling Exploration and Exploitation in Deep Reinforcement Learning Algorithms,” in Proceedings of the 35thInternational Conference on Machine Learning, PMLR, Jul.2018, pp.1039–1048. Accessed: Jul. 15, 2023. [Online]. Available: proceedings.mlr.press / v80 / colas18a.html

[0291]

[0126] A. Pandi et al., “A versatile active learning workflow for optimization of genetic and metabolic networks,” Nat. Commun., vol.13, no.1, Art. no.1, Jul.2022, doi: 10.1038 / s41467-022-31245-z.

[0292]

[0127] A. R. G. Libby et al., “Automated Design of Pluripotent Stem Cell Self- Organization,” Cell Syst., vol.9, no.5, pp.483-495.e10, Nov.2019, doi 10.1016 / j.cels.2019.10.008.

[0293]

[0128] A. Pietak and M. Levin, “Exploring Instructive Physiological Signaling with the Bioelectric Tissue Simulation Engine,” Front. Bioeng. Biotechnol., vol.4, 2016, Accessed: Jun.24, 2022. [Online]. Available: frontiersin.org / article / 10.3389 / fbioe.2016.00055

[0294]

[0129] A. Koseska and P. I. Bastiaens, “Cell signaling as a cognitive process,” EMBO J., vol.36, no.5, pp.568–582, Mar.2017, doi: 10.15252 / embj.201695383.

[0295]

[0130] F. Baluška, A. S. Reber, and W. B. Miller, “Cellular sentience as the primary source of biological order and evolution,” Biosystems, vol.218, p.104694, Aug.2022, doi: 10.1016 / j.biosystems.2022.104694.Tufts Ref. T002751 Quarles Ref.166118.01450

[0296]

[0131] F. Baluška, W. B. Miller, and A. S. Reber, “Cellular and evolutionary perspectives on organismal cognition: from unicellular to multicellular organisms,” Biol. J. Linn. Soc., vol.139, no.4, pp.503–513, Aug.2023, doi: 10.1093 / biolinnean / blac005.

[0297]

[0132] A. S. Reber and F. Baluška, “Cognition in some surprising places,” Biochem. Biophys. Res. Commun., vol.564, pp.150–157, Jul.2021, doi: 10.1016 / j.bbrc.2020.08.115.

[0298]

[0133] F. Baluška and A. S. Reber, “Cellular and organismal agency - Not based on genes: A comment on Baverstock,” Prog. Biophys. Mol. Biol., vol.167, pp.161–162, Dec. 2021, doi: 10.1016 / j.pbiomolbio.2021.11.001.

[0299]

[0134] A. Bernheim-Groswasser, N. S. Gov, S. A. Safran, and S. Tzlil, “Living Matter: Mesoscopic Active Materials,” Adv. Mater., vol.30, no.41, p.1707028, 2018, doi: 10.1002 / adma.201707028.

[0300]

[0135] M. Etcheverry, M. Levin, C. Moulin-Frier, and P.-Y. Oudeyer, “SBMLtoODEjax: efficient simulation and optimization of ODE SBML models in JAX.” arXiv, Jul.17, 2023. doi: 10.48550 / arXiv.2307.08452.

[0301]

[0136] J. K. Pugh, L. B. Soros, P. A. Szerlip, and K. O. Stanley, “Confronting the challenge of quality diversity,” in Proceedings of the 2015 annual conference on genetic and evolutionary computation, 2015, pp.967–974.

[0302]

[0137] S. M. Scheiner, “A compilation of and typology for abundance-, phylogenetic-and functional-based diversity metrics,” BioRxiv Prepr. Serv. Biol., p.530782, 2019

[0303] Other features, objects, and advantages of the present disclosure are apparent in the detailed description. It should be understood, however, that the detailed description, while indicating embodiments of the present disclosure, is given by way of illustration, not limitation. Various changes and modifications within the scope of the disclosure will become apparent to those skilled in the art from the detailed description.

[0304] Although inventive concepts have been described with reference to the embodiments illustrated in the attached drawing figures, equivalents may be employed andTufts Ref. T002751 Quarles Ref.166118.01450 substitutions made herein without departing from the scope of the claims. Components illustrated and described herein are merely examples of a system / device and components that may be used to implement embodiments of the inventive concepts and may be replaced with other devices and components without departing from the scope of the claims.

Claims

Tufts Ref. T002751 Quarles Ref.166118.01450 CLAIMS What is claimed is:

1. A method of identifying a behavioral catalogue of a gene regulatory network (GRN), the behavioral catalog comprising at least one goal behavioral state and a plurality of trajectories connecting controlled sources of variation in the GRN and each of the at least one goal behavioral state, the method comprising: a. obtaining a data set comprising a mathematical model of the GRN; b. using an algorithm to identify at least one goal behavioral state and a plurality of corresponding trajectories connecting a plurality of controlled sources of variation in the GRN and each of the at least one goal behavioral state; c. testing a robustness of the at least one goal behavioral state by perturbing each of the plurality of corresponding trajectories and evaluating, for each of the plurality of corresponding trajectories, whether the respective goal behavioral state is reached; d. selecting a subset of the goal behavioral states which are robust to perturbations; and e. saving the subset of goal behavioral states, the corresponding plurality of controlled sources of variation, and the corresponding plurality of trajectories.

2. The method of claim 1, wherein the mathematical model comprises: a. a plurality of nodes representing molecular species in the GRN; b. a plurality of differential equations to model dynamic interactions between the plurality of nodes; c. a plurality of observations of the dynamic interactions between the plurality of nodes over time; and d. a plurality of goal behavioral states provided as steady states of the dynamic interactions observed over time.Tufts Ref. T002751 Quarles Ref.166118.01450 3. The method of claim 1, wherein the algorithm comprises a protocol based on at least one of machine learning, artificial intelligence, robotics, autonomous vehicle navigation, behavioral science, cognitive neuroscience, and artificial life.

4. The method of claim 3, wherein the protocol is based on active interference framework, perceptual control theory, training for habituation, training for sensitization, training for anticipation, training for instrumental learning, training for Pavlovian condition, training for planning, predictive coding, principle of empowerment, and casual information theory.

5. The method of claim 1, wherein the algorithm comprises a curiosity driven algorithm.

6. The method of claim 1, wherein using the algorithm further comprises: a. randomly sampling a plurality of sources of variation and determining a corresponding goal behavioral state and a corresponding trajectory to populate an initial mapping comprising a plurality of sources of variation, a plurality of trajectories, and a plurality of goal behavioral states; b. further populating the initial mapping by sampling the goal behavioral states of the initial mapping and optimizing additional goal states and additional trajectories; and c. storing the goal behavioral states and corresponding trajectories.

7. The method of claim 6, wherein randomly sampling further comprises using a uniform sampling strategy within a bounding hyper-rectangle of currently reached goal behavioral states, wherein a volume of the hyper-rectangle is larger than a distribution of currently reached goal behavioral states.

8. The method of claim 1, wherein testing the robustness comprises perturbing the plurality of corresponding trajectories with at least one of noise perturbations, push perturbations, or wall perturbations.

9. The method of claim 8, wherein the noise perturbations are parameterized by their standard deviation and period.Tufts Ref. T002751 Quarles Ref.166118.01450 10. The method of claim 8, wherein the push perturbations are parameterized by their magnitude and number of occurrences.

11. The method of claim 8, wherein the wall perturbations are parameterized by their length and number, and wherein the walls are generated in locations of a space that the GRN would visit without induced perturbations.

12. The method of claim 1, wherein each controlled sources of variation in the plurality of controlled sources of variation comprises a specific node of the GRN and an action to change the behavior of the specific node of the GRN.

13. The method of claim 12, wherein the node comprises at least one of a gene, a protein, or other GRN component.

14. The method of claim 12, wherein an action to change the behavior of the node of the GRN includes at least one of perturbing parameters of the node, altering initial conditions of the node, or applying external stimuli to the node.

15. The method of claim 1, wherein the goal behavioral state is a cell phenotype of at least one of cell metabolism, cell proliferation, cell differentiation, cell migration, cell shape, cell signaling, programmed cell death, gene expression, and senescence.

16. The method of claim 1, wherein the goal behavior state is a multi-cellular phenotype of at least one of morphogenesis, remodeling, wound healing, regeneration, and embryonic development.

17. The method of claim 1, wherein the goal behavioral state is potentiation of long-term memory.Tufts Ref. T002751 Quarles Ref.166118.01450 18. The method of claim 1, wherein each goal behavioral state of the subset of goal behavioral states has more than one corresponding trajectory and each of the corresponding trajectories reaches the goal behavioral state when perturbed.

19. A system for identifying a behavioral catalogue of a gene regulatory network (GRN), the behavioral catalog comprising at least one goal behavioral state and a plurality of trajectories connecting controlled sources of variation in the GRN and each of the at least one goal behavioral state, the system comprising at least one hardware processor that is programmed to perform one or more of the following steps: a. obtaining a data set comprising a mathematical model of the GRN; b. executing an algorithm to identify at least one goal behavioral state and a plurality of corresponding trajectories connecting controlled sources of variation in the GRN and each of the at least one goal behavioral state; c. testing a robustness of the at least one goal behavioral state by perturbing each of the plurality of corresponding trajectories and evaluating, for each of the plurality of corresponding trajectories, whether the respective goal behavioral state is reached; d. selecting a subset of the goal behavioral states which are robust to perturbations; e. saving the subset of goal behavioral states, the corresponding plurality of controlled sources of variation, and the corresponding plurality of trajectories; and f. outputting the subset of goal behavioral states and the corresponding plurality of trajectories to a user.

20. The system of claim 15, wherein the mathematical model comprises: a. a plurality of nodes representing molecular species in the GRN; b. a plurality of differential equations to model dynamic interactions between the plurality of nodes; c. a plurality of observations of the dynamic interactions between the plurality of nodes over time; and d. a plurality of goal behavioral states provided as steady states of the dynamic interactions observed over time.Tufts Ref. T002751 Quarles Ref.166118.01450 21. The system of claim 19, wherein the algorithm comprises a protocol based on at least one of machine learning, artificial intelligence, robotics, autonomous vehicle navigation, behavioral science, cognitive neuroscience, and artificial life.

22. The system of claim 21, wherein the protocol is based on active interference framework, perceptual control theory, training for habituation, training for sensitization, training for anticipation, training for instrumental learning, training for Pavlovian condition, training for planning, predictive coding, principle of empowerment, and casual information theory.

23. The system of claim 19, wherein the algorithm comprises a curiosity driven algorithm.

24. The system of claim 19, wherein executing the algorithm further comprises: a. randomly sampling a plurality of sources of variation and determining a corresponding goal behavioral state and a corresponding trajectory to populate an initial mapping comprising a plurality of sources of variation, a plurality of trajectories, and a plurality of goal behavioral states; b. further populating the initial mapping by sampling the goal behavioral states of the initial mapping and optimizing additional goal states and additional trajectories; and c. storing the goal behavioral states and corresponding trajectories.

25. The system of claim 24, wherein randomly sampling further comprises using a uniform sampling strategy within a bounding hyper-rectangle of currently reached goal behavioral states, wherein a volume of the hyper-rectangle is larger than a distribution of currently reached goal behavioral states.

26. The system of claim 19, wherein testing the robustness comprises perturbing the plurality of corresponding trajectories with at least one of noise perturbations, push perturbations, or wall perturbations.

27. The system of claim 26, wherein the noise perturbations are parameterized by their standard deviation and period.Tufts Ref. T002751 Quarles Ref.166118.01450 28. The system of claim 26, wherein the push perturbations are parameterized by their magnitude and number of occurrences.

29. The system of claim 26, wherein the wall perturbations are parameterized by their length and number, and wherein the walls are generated in locations of a space that the GRN would visit without induced perturbations.

30. The system of claim 15, wherein each controlled sources of variation in the plurality of controlled sources of variation comprises a specific node of the GRN and an action to change the behavior of the specific node of the GRN.

31. The system of claim 30, wherein the node comprises at least one of a gene, a protein, or other GRN component.

32. The system of claim 30, wherein an action to change the behavior of a node of the plurality of nodes of the GRN includes at least one of perturbing parameters of the node, altering initial conditions of the node, or applying external stimuli to the node.

33. The system of claim 19, wherein the goal behavioral state is a cell phenotype of at least one of cell proliferation, cell differentiation, cell migration, cell shape, cell signaling, programmed cell death, gene expression, and senescence.

34. The system of claim 19, wherein the goal behavior state is a multi-cellular phenotype of at least one of morphogenesis, remodeling, wound healing, regeneration, and embryonic development.

35. The system of claim 19, wherein the goal behavioral state is potentiation of long-term memory.Tufts Ref. T002751 Quarles Ref.166118.01450 36. The system of claim 19, wherein each goal behavioral state of the subset of goal behavioral states has more than one corresponding trajectory and each of the corresponding trajectories reaches the goal behavioral state when perturbed.

37. A method of identifying a therapeutic intervention, the method comprising: a. obtaining a behavioral catalog of a gene regulatory network (GRN) associated with a disease, wherein the behavioral catalog comprises a plurality of behavioral end states and a plurality of trajectories connecting each controlled source of variation in the plurality of controlled sources of variation and the corresponding behavioral end state in the plurality of behavioral end states; b. selecting a disease associated with the GRN; c. classifying each behavioral end state in the plurality of behavioral end states as a first setpoint, wherein the first setpoint can be a disease setpoint or the first setpoint can be a healthy setpoint, resulting in a plurality of first setpoints, wherein a disease setpoint is associated with a phenotype or gene expression pattern associated with the disease, and wherein the healthy setpoint is not associated with a phenotype or gene expression pattern associated with the disease; d. defining a plurality of stimulus-based interventions to apply to each disease setpoint; e. introducing each stimulus-based intervention in the plurality of stimulus-based interventions to each first setpoint in the plurality of first setpoints and identifying a plurality of trajectories, wherein each trajectory in the plurality of trajectories corresponds to a response of a first setpoint in the plurality of first setpoints to a stimulus-based intervention in the plurality of stimulus-based interventions, and determining a second setpoint of each trajectory in the plurality of trajectories; f. selecting a subset of trajectories in which the first setpoint is a disease state and the second setpoint after introducing the stimulus-based interventions is a healthy state; g. determining the subset of stimulus-based interventions that corresponds with the subset of trajectories in which the first setpoint is a disease state, and the second setpoint is a healthy state; and h. outputting the subset of stimulus-based interventions.Tufts Ref. T002751 Quarles Ref.166118.01450 38. The method of claim 37, wherein classifying healthy setpoints and disease setpoints comprises: a. constructing an energy landscape of the behavioral catalog of the GRN; b. determining which regions of the energy landscape correspond with a disease region, wherein the phenotype or gene expression pattern of the region is associated with a disease, and determining which regions of the energy landscape correspond with a healthy region, wherein the phenotype or gene expression pattern of the region is not associated with the disease; c. identifying a subset of endpoints in the behavioral catalog wherein each endpoint in the subset of endpoints falls within the disease region and classifying the endpoint as a disease setpoint, and identifying a subset of endpoints in the behavioral catalog wherein each endpoint in the subset of endpoints falls within the healthy region and classifying the endpoint as a healthy setpoint; d. selecting a subset of trajectories associated with a disease setpoint and a subset of trajectories associated with a healthy setpoint that exhibit a low level of sensitivity to perturbations; and e. outputting the subset of trajectories and associated end setpoints.

39. The method of claim 37, wherein the plurality of stimulus-based interventions is a piecewise constant function that determines, which nodes in the GRN to intervene on, when to apply the intervention, and with what amplitude the stimulus-based intervention should be applied.

40. The method of claim 37, wherein the nodes in the GRN are selected based on which nodes can be targeted with drugs that treat the disease and with what precision the drugs can be administered.

41. The method of claim 37, wherein introducing the stimulus-based interventions further comprises applying a plurality of perturbations, wherein the perturbations comprise at least one of noise perturbations, push perturbations, or wall perturbations.Tufts Ref. T002751 Quarles Ref.166118.01450 42. The method of claim 37, wherein selecting the subset of trajectories further comprises selecting trajectories which reach the healthy state under all perturbations applied to the trajectory.

43. The method of claim 37 wherein the GRN is the RKIP-ERK signaling pathway.

44. The method of claim 37, wherein the phenotype or gene expression pattern of the disease region is associated with tumor development.

45. The method of claim 37, wherein the healthy end state comprises a state in which a cell has an elevated susceptibility to cancer therapies.

46. A system for identifying a therapeutic intervention, the system comprising at least one hardware processor that is programmed to perform one or more of the following steps: a. obtaining a behavioral catalog of a gene regulatory network (GRN) associated with a disease, wherein the behavioral catalog comprises a plurality of behavioral end states and a plurality of trajectories connecting controlled sources of variation in the GRN and each behavioral end state in the plurality of behavioral end states; b. selecting a disease associated with the GRN; c. classifying each behavioral end state in the plurality of behavioral end states as a first setpoint, wherein the first setpoint can be a disease setpoint or the first setpoint can be a healthy setpoint, resulting in a plurality of first setpoints, wherein a disease setpoint is associated with a phenotype or gene expression pattern associated with the disease, and wherein the healthy setpoint is not associated with a phenotype or gene expression pattern associated with the disease; d. defining a plurality of stimulus-based interventions to apply to each disease setpoint; e. introducing each stimulus-based intervention in the plurality of stimulus-based interventions to each first setpoint in the plurality of first setpoints and identifying a plurality of trajectories, wherein each trajectory in the plurality of trajectories corresponds to a response of a first setpoint in the plurality of first setpoints to a stimulus-based intervention in the plurality ofTufts Ref. T002751 Quarles Ref.166118.01450 stimulus-based interventions, and determining a second setpoint of each trajectory in the plurality of trajectories; f. selecting a subset of trajectories in which the first setpoint is a disease state and the second setpoint after introducing the stimulus-based interventions is a healthy state; g. determining the subset of stimulus-based interventions that corresponds with the subset of trajectories in which the first setpoint is a disease state, and the second setpoint is a healthy state; and h. outputting the subset of stimulus-based interventions.

47. The system of claim 46, wherein the processor, when classifying healthy setpoints and disease setpoints, is further programmed to perform the steps of: a. constructing an energy landscape of the behavioral catalog of the GRN; b. determining which regions of the energy landscape correspond with a disease region, wherein the phenotype or gene expression pattern of the region is associated with a disease, and determining which regions of the energy landscape correspond with a healthy region, wherein the phenotype or gene expression pattern of the region is not associated with the disease; c. identifying a subset of endpoints in the behavioral catalog wherein each endpoint in the subset of endpoints falls within the disease region and classifying the endpoint as a disease setpoint, and identifying a subset of endpoints in the behavioral catalog wherein each endpoint in the subset of endpoints falls within the healthy region and classifying the endpoint as a healthy setpoint; d. selecting a subset of trajectories associated with a disease setpoint and a subset of trajectories associated with a healthy setpoint that exhibit a low level of sensitivity to perturbations; and e. outputting the subset of trajectories and associated end setpoints.

48. The system of claim 46, wherein the plurality of stimulus-based interventions is a piecewise constant function that determines, which nodes in the GRN to intervene on, when to apply the intervention, and with what amplitude the stimulus-based intervention should be applied.Tufts Ref. T002751 Quarles Ref.166118.01450 49. The system of claim 46, wherein the nodes in the GRN are selected based on which nodes can be targeted with drugs that treat the disease and with what precision the drugs can be administered.

50. The system of claim 46, wherein introducing the stimulus-based interventions further comprises applying a plurality of perturbations, wherein the perturbations comprise at least one of noise perturbations, push perturbations, or wall perturbations.

51. The system of claim 46, wherein selecting the subset of trajectories further comprises selecting trajectories which reach the healthy state under all perturbations applied to the trajectory.

52. The system of claim 46, wherein the GRN is the RKIP-ERK signaling pathway.

53. The system of claim 46, wherein the phenotype or gene expression pattern of the disease region is associated with tumor development.

54. A method of identifying at least one of a configuration or a parameter of a gene regulatory network (GRN) to perform a desired function, the method comprising: a. obtaining a data set comprising a mathematical model of a GRN; b. using an algorithm to identify a plurality of goal behavioral states and a plurality of corresponding trajectories connecting a plurality of controlled sources of variation in the GRN and each of the plurality of goal behavioral states; and c. optionally using a local optimization strategy to identify an additional goal behavioral state and add the additional goal behavioral state to the plurality of goal behavioral states.

55. The method of claim 54, wherein the GRN is an oscillator circuit.

56. The method of claim 54, wherein the controlled sources of variation comprise at least one of a genetic intervention or an environmental intervention.Tufts Ref. T002751 Quarles Ref.166118.01450 57. A system for identifying at least one of a configuration or a parameter of a gene regulatory network (GRN) to perform a desired function, the system comprising at least one hardware processor that is programmed to perform one or more of the following steps: a. obtaining a data set comprising a mathematical model of a GRN; b. using an algorithm to identify a plurality of goal behavioral states and a plurality of corresponding trajectories connecting a plurality of controlled sources of variation in the GRN and each of the plurality of goal behavioral states; and c. optionally using a local optimization strategy to identify an additional goal behavioral state and add the additional goal behavioral state to the plurality of goal behavioral states.

58. The system of claim 57, wherein the GRN is an oscillator circuit.

59. The system of claim 57, wherein the controlled sources of variation comprise at least one of a genetic intervention or an environmental intervention.

60. A method of priming a biological system for therapeutic treatment, comprising administering a compound to the system to move the biological system to a desired state, wherein the desired state is a state the biological system has not yet achieved.

61. The method of claim 60, wherein the biological system is in a disease state.

62. The method of claim 60, wherein the desired state comprises a state in which the biological process has a phenotype associated with acquired sensitivity to at least one therapeutic treatment.

63. The method of claim 60, wherein the desired state is a cell phenotype of at least one of restored sensitivity to at least one therapeutic agent, reduced resistance to the at least one therapeutic agent, and prevented resistance to the at least one therapeutic agent.Tufts Ref. T002751 Quarles Ref.166118.01450 64. The method of claim 60, wherein the desired state is a cell phenotype that is susceptible to a therapeutic agent that the biological system was previously not susceptible to.

65. The method of claim 60, wherein the administered therapeutic treatment comprises administering at least one drug.

66. The method of claim 65, wherein the at least one drug is selected from a group of chemical classes comprising propanolamine, C-amidated tetrapeptide, aryl sulfide, alkyl- phenylketone, benzazepine, secondary amine, benzenoid, monomethoxylflavone, monohydroxypyrrolidine.

67. The method of claim 60, wherein the method further comprises determining that the biological system has reached the desired state.

68. The method of claim 60, wherein the method further comprises administering a plurality of compounds to the system.

69. A method of administering a therapeutic treatment to a biological system, comprising administering a compound to the biological system to move the biological system to desired state, wherein the desired state is a more differentiated state of the biological system.

70. The method of claim 69, wherein the method further comprises determining that the biological system has reached the more differentiated state.

71. The method of claim 70, wherein an appropriate new therapeutic treatment is chosen based on the determined more differentiated state.

72. The method of claim 69, wherein the biological system is a patient in need of the therapeutic treatment.

73. The method of claim 72, wherein the patient has cancer.Tufts Ref. T002751 Quarles Ref.166118.01450 74. The method of claim 73, wherein the cancer is a hematologic cancer.

75. The method of claim 69, wherein the compound is selected from a group of chemical classes comprising propanolamine, C-amidated tetrapeptide, aryl sulfide, alkyl-phenylketone, benzazepine, secondary amine, benzenoid, monomethoxylflavone, monohydroxypyrrolidine.

76. The method of claim 69, wherein the method further comprises administering a plurality of compounds to the biological system.

Images