Introduction

Understanding how biochemical networks produce biological behavior is one of the central goals of systems and synthetic biology. Studying the design principles of these networks, or circuits, provides fundamental insights into their computational capabilities and enables the application of synthetic biology to problems such as cell-based therapeutics (Lim et al. (2013); Williams et al. (2020)) and engineered microbial communities (Jones et al. (2024); Purnick and Weiss (2009)). The topology, or qualitative architecture, of a synthetic gene circuit is a critical determinant of its behavior, and choosing a topology is the first step in engineering circuits with novel functions such as switching in response to an input (Gardner et al. (2000)), spontaneous oscillations (Elowitz and Leibler (2000)), or multistability (Zhu et al. (2022)). Circuit topologies with a particular desired phenotype have traditionally been discovered de novo using a method we call computational enumeration. Enumeration involves listing every topology in a predefined space of possibilities and simulating its behavior with an exhaustive set of random perturbations of the modeling parameters (Ma et al. (2009); Cotterell and Sharpe (2010); Chau et al. (2012); Schaerli et al. (2014); Gerardin et al. (2019)). This method identifies the topologies that are the most robust to parameter choice, as well as the network motifs of the phenotype, which are topological patterns overrepresented among highly robust circuits (Alon (2007, 2019)).

Currently, computational complexity has limited the role of automation in circuit design. Assuming 3 possible types of interactions (activation, inhibition, or no interaction), the number of distinct topologies that can be made with N circuit components scales as , not accounting for symmetries. Because the enumeration method requires simulating each topology in this potential space with many parameter perturbations (typically ≥ 10⁴), designing a five-component circuit with this approach would require on the order of 10¹⁵-10¹⁶ simulations, a computationally intractable number. This “curse of dimensionality” constitutes the main barrier to understanding the structure-function relationship for large circuits. The search space can be narrowed by only considering combinations of smaller motifs (Qiao et al. (2019)), but this approach risks excluding novel, non-modular designs that appear at higher complexity (Jiménez et al. (2017)). For large circuits, a functional topology may be found with evolutionary optimization (François et al. (2007); François and Siggia (2010)) or recurrent neural network-guided inference (Shen et al. (2021)), but these methods can be sensitive to the optimization parameters and the nature of the fitness land-scape, a fact that complicates inference of network motifs and other general design principles. The ideal design framework for circuit topologies would be: (i) Efficient: Using relatively few samples, it returns a set of reasonably robust solutions. (ii) Generalizable: It generates solutions that are representative of the ground truth, enabling inference of motifs and design principles. (iii) Scalable: It maintains features (i) and (ii) when searching large spaces of potential designs.

Our inspiration for a circuit design approach with these features comes from algorithms developed in modern artificial intelligence to solve decision problems, which share core computational similarities with genetic circuit design. In a decision problem, or game, one starts from an initial position (the root state) and makes a sequence of decisions ending in success or failure. Algorithms for optimal game playing sample the branching tree of possible decisions in search of the best decision sequences, or paths. Large game trees, characteristic of long, complex games like chess (∼ 10⁴³ states) and Go (10¹⁷⁰ states), can overwhelm classical exhaustive search methods (e.g., flat search, A∗, and minimax) and branch-and-bound methods (e.g., α-β pruning) (Shannon (1950); Tromp and Farnebäck (2007)). In contrast, most modern probabilistic search algorithms use reinforcement learning (RL) to sample decision paths; for example, Monte Carlo tree search (MCTS) uses an RL decision rule that optimistically estimates the success of each possible decision (Kocsis and Szepesvári (2006)). MCTS has proven to be instrumental in solving large games (Silver et al. (2016, 2018)) due to its low memory footprint, ability to be parallelized (Enzenberger and Müller (2010)), balance of exploitation and exploration (Auer et al. (2002)), and myriad variations (Świechowski. (2022); Browne et al. (2012)).

In this work, we present a circuit design platform based on tree search. This platform, which we call CircuiTree, approaches circuit design as a sequence of assembly decisions and optimizes the circuit topology for a given phenotype by using MCTS to search the tree of possible circuit assemblies. We first used CircuiTree to search a trial space – all 3,325 connected circuits with three transcription factors – for topologies that exhibit spontaneous, sustained oscillations. As it dynamically explores the search space, CircuiTree discovers sparse oscillators before complex ones and reliably identifies the most robust topologies. Next, we demonstrate that CircuiTree accurately identifies assembly motifs, which are overrepresented search states analogous to network motifs found by enumeration. The best assembly motifs for 3-component oscillators also contain the activator-inhibitor (AI) and/or repressor (Rep) motifs, a result consistent with previous studies of 2- and 3-component oscillators. Finally, we use a parallelized implementation of CircuiTree to study mutational tolerance in large circuits. We discover that highly fault-tolerant five-component oscillator circuits are “multiplexed” in that they contain multiple AI and Rep motifs that are syner-gistically interleaved such that the oscillations are robust to complete and partial knockdowns of single genes. Our results demonstrate that CircuiTree’s RL-based framework enables the efficient, generalizable, and scalable design of biochemical circuits, which could be used to engineer and understand complex intracellular and multicellular biological systems.

Results

Circuit design as an assembly game

We define a circuit as a system of k dynamically interacting biochemical species x = [x₀, …, x_k] characterized by its qualitative topology s, a member of the set of possible topologies 𝒮, and its parameters θ ∈ ℝ^n(s), which describe the rates of biochemical reactions. We also assume the existence of a phenotype function f (x(t)) = q ∈ {0, 1} that returns a yes-or-no verdict as to whether a given simulation x(t) exhibits the desired phenotype. This function, which we conceptualize quite generally, may represent the presence of a dynamical phenotype such as adaptation or oscillation or a more general property such as multistability. Like many prior computational studies of design principles (Ma et al. (2009); Chau et al. (2012); Cotterell and Sharpe (2010); Schaerli et al. (2014)), each simulation of the system is performed with a set of randomly chosen parameters, and we define robustness Q(s) as the mean phenotype score from these samples Q(s) = 𝔼_θ[q(s)] = 𝔼_θ [f (x(t; s, θ))]. The most robust topology s∗ is therefore the solution to a combinatorial optimization problem. (Note that Q in this context should not be confused with the value function used in RL.)

In this work, for convenience, we consider models where all topologies are governed by the same reaction parameters (i.e. n(s_i) = n(s_j)∀s_i, s_j ∈ 𝒮), but CircuiTree does not require this assumption.

To solve Equation 1, we approach circuit topology design as a multi-step decision problem, or game (Figure 1A). Beginning with an “empty” circuit s₀ with a given set of components and without interactions, each step of the game can add an activating or inhibitory (auto-)regulatory interaction. The game ends when the builder chooses to “terminate” the assembly at some topology s_i, at which point the outcome is decided by a simulation with a win probability Q(s_i). In the language of Markov decision processes (MDPs), each topology s represents a state of the game, and the addition of an interaction or the termination of assembly represents an action a that converts a state s_i to a different state s_j. Each playout of the game traces a path on the tree of possible decisions T rooted at s₀. The nodes in the tree are circuit topologies 𝒮, and each directed edge (s_i, s_j) ∈ ε represents the assembly of s_j from s_i by an action. Note that because multiple decision paths can converge on the same s_j, T is technically a directed acyclic graph (DAG), but we will call it a tree for simplicity.

Figure 1.

For large design spaces, it is computationally infeasible to estimate Q(s) for every s ∈ T. Instead, we iteratively sample from T using Monte Carlo tree search (MCTS), an RL algorithm commonly used in MDPs and game-playing artificial intelligences (Silver et al. (2016, 2017, 2018)) to bound the complexity of searching large decision trees. In these applications, tree search is used during game play (and sometimes also during training) to find the best move by simulating many hypothetical scenarios forward in time.

Each MCTS iteration starts at the root state s₀ and undergoes four phases, shown in Figure 1B and described in Algorithm 1. The key innovation of MCTS occurs during the first phase (selection), when it uses the UCT criterion from RL to choose moves while balancing (i) exploitation of moves that have yielded high rewards in the past and (ii) exploration of new moves that may yield higher rewards. At each branch (Figure 1C, left), from a state s_i, multiple actions a_j are available that have yielded cumulative rewards r_ij over υ_ij previous visits. MCTS decides the optimal action using the decision policy UCT(i, j) (Figure 1C, right), originally developed to solve the classic multi-armed bandit problem Auer et al. (2002).

The first term of the UCT criterion estimates the mean reward of the move s_i → s_j, while the second term, derived from Hoeffding’s inequality (Kocsis and Szepesvári (2006)), is an estimated upper bound on the sampling error. The second term ensures that a move with a historically low reward but relatively few samples (υ_j << ∑_j υ_ij) will still be tried occasionally to account for insufficient sample size. Thus, UCT is an optimistic predictor of a move’s true mean reward. Note that the two terms encourage exploitation and exploration, respectively. (By definition, UCT(i, j) = +∞ if s_i → s_j has not been visited before).

Using the UCT criterion, actions are selected recursively, breaking any ties randomly, until the new state s_j is terminal (i.e. the assembly is finished) or has not been visited before during sampling. An unvisited s_j is added to T. Then, during the simulation phase, the outcome q of the trial is determined by taking random actions until the “terminate” action is chosen and simulating the resulting topology with randomly chosen reaction parameters and initial conditions. Finally, the history of rewards and visits is updated for each selected edge in T (backpropagation).

In addition to finding individual circuit topologies, an important goal of circuit design is to identify specific structural features, or network motifs, that lead to successful designs (Milo et al. (2002); Alon (2007); Ma et al. (2009); Cotterell and Sharpe (2010); Shah and Sarkar (2011); Chau et al. (2012); Lim et al. (2013); Schaerli et al. (2014)). Similarly, in our game-playing paradigm, a motif of the assembly game is a series of moves (an assembly strategy) that leads to topologies with a high rate of phenotypic success, on average. Thus, when a motif is present in the space of possible designs, it creates an enriched region in T, which can be exploited during the tree search to rapidly discover successful designs (Figure 1D).

Establishing a ground-truth for three-node stochastic oscillators with enumeration

Oscillations appear in diverse cellular contexts such as metabolism (Chandra et al. (2011)), DNA damage (Geva-Zatorsky et al. (2006, 2010)), the cell cycle (Ferrell et al. (2011)), and circadian rhythms (Tyson et al. (1999)), and consequently, many of their design principles have been elucidated (Novák and Tyson (2008)). We first benchmark CircuiTree on the well-studied problem of designing an oscillator circuit with three nodes (Figure 2A, part i). Specifically, we consider a system of symmetric transcription factors (TFs) modeled as a stochastic birth-and-death process of individual mRNAs, TFs, and TF-response element (TF-RE) complexes, with elementary reactions of transcription, translation, degradation, binding, and unbinding (Figure S1A; reactions and their rate parameters are summarized in Table 1). TFs regulate transcription by binding to one of two REs in the regulatory region of each promoter, and cooperativity is introduced by assuming that TF-RE binding is stronger when both sites are occupied by the same TF. This assumption reduces the computational cost associated with modeling every homo- and heterodimer and their dimer-RE complexes. For each stochastic simulation, the system was randomly initialized and stochastically simulated using the Gillespie method for t_max = 4 × 10⁴ sec = 11.1 hrs (see Methods for details of initialization). Oscillations were quantified by computing the normalized autocorrelation function (ACF) and identifying its lowest minimum value across all TFs ACF_min, excluding the bounds, and a Boolean reward value was assigned based on a cutoff value (Figure 2A, part ii).

Table 1.

Figure 2.

To compare with enumeration, all 3,325 unique, fully connected topologies were simulated with 10⁴ randomly sampled parameter sets and initial conditions. The 10 rate parameters of the model were reduced to 8 dimensionless variables, which were sampled uniformly from a range of values determined based on experimentally measured rates (Milo et al. (2010)) and known requirements for oscillations (Elowitz and Leibler (2000); Novák and Tyson (2008)). (See Table 2 for variable definitions and sampling ranges and Methods for details of parameter sampling). The robustness to parameter perturbation Q was calculated as the fraction of parameter sets that produced oscillation, and topologies with Q > 0.01 were considered oscillators. This procedure uncovered 221 oscillator topologies (6.65% of the search space) (see Figures S1B and S1C for a summary of the best topologies). These results generally agree with prior studies of two- and three-node oscillators (Qiao et al. (2022); Elowitz and Leibler (2000); Novák and Tyson (2008); Woods et al. (2016); Stricker et al. (2008)). For instance, almost all oscillators contain a repressilator (Rep) loop, activator-inhibitor (A-I) loop, or a combination thereof. Positive autoregulation (PAR), which has been shown to buffer extrinsic and intrinsic noise (Qiao et al. (2022)), is also ubiquitous among these oscillators and is required for robust A-I loop oscillations (Novák and Tyson (2008); Stricker et al. (2008)). Further discussion of these motifs can be found in the section below. Interestingly, stochastic switching between stable states was occasionally mistaken for low-frequency oscillations (Figure S2), causing some toggle switch circuits to be classified as oscillators (highest Q = 0.024, ranked #139).

Table 2.

Algorithm 1 Circuit topology search using MCTS Algorithm 1 Circuit topology search using MCTS

Across all oscillators, oscillation was favored by tight TF-RE binding (K_D,1 < 10³), low basal transcription (k_m,unbound < 0.3 sec⁻¹), strong repression (k_m,rep/k_m,unbound < 10⁻¹), and a high activated transcription rate (k_m,act, monotonic effect) (Figure S3A). Protein and mRNA degradation rates had a non-monotonic effect with a peak at γ_m ≈ γ_p ≈ 0.2 molec⁻¹ sec⁻¹ (Figure S3A), and oscillation period depended strongly on these rates and their ratio (Figures S3B and S3C). However, the best oscillators were exceptionally robust to the choice of parameters. For example, the most robust oscillator (the repressilator with PAR) has a Q of 0.767, meaning that for each of the 8 sampled variables, on average, 96.7% of the possible values would permit oscillations (0.967⁸ ≈ 0.767).

CircuiTree efficiently and systematically discovers three-node oscillators

MCTS masters a game with a limited number of trials by balancing deep sampling of promising regions with broad sampling of under-explored regions. To understand how this strategy performs for circuit topology design, we use CircuiTree to search for 3-node oscillators given N = 10⁵ iterations of MCTS (n = 50 replicates; see Methods for MCTS implementation specifics). In the average run, the first putative result is discovered after just 0.69 samples/topology (2,280 iterations), and by 19.6 samples/topology (65,360 iterations), a top-5 oscillator has been sampled >100 times. A topology is discovered as a “successful” oscillator if its estimated robustness exceeds the threshold value Q_thresh = 0.01, which was chosen heuristically prior to the search. Figure 2B shows a representative example of how, over sampling time, CircuiTree incrementally builds the search tree from the root state (top to bottom), encounters more oscillators (shown in orange), and improves its prediction for the most robust design (shown as a black circuit diagram). By the end of sampling time, the average run saturates the tree, sampling 99.97% of topologies at least once.

CircuiTree balances efficiency and comprehensiveness by first finding sparse solutions before exploring deeper areas of the tree. In Figure 2C, each row of the heatmap is one of the 221 oscillators, and the rows are sorted first by the circuit’s complexity, or the number of interactions, then by decreasing robustness (Q). The color scale indicates the likelihood of discovery, or recall, of each oscillator measured as the proportion of replicates that discovered it. Because sparse topologies require fewer assembly steps, MCTS encounters the sparsest oscillators (such as the repressilator) first, and it discovers increasingly complex solutions over time until oscillators with 9 interactions (the maximum) are found at N ≈ 5 × 10⁴. As shown by the line plot in Figure 2D, CircuiTree has a very high recall of 94.7% for the top 10% of oscillators (95% CI: 86.4% − 100.0%) and has a recall of 35.9% (95% CI: 30.8% − 41.5%) for oscillators in general. This is fairly high, considering that it gets an average of samples/topology, and a majority of oscillators have a Q of less than 1/30. Additionally, CircuiTree has a precision of 81.5% (95% CI: 74.1% − 88.3%), indicating a low rate of false positives.

The competing goals of breadth and depth manifest during tree search as distinct temporal phases in which MCTS first explores a broad set of topologies before focusing on a narrow, enriched subset. Among topologies containing a combination of AI and Rep feedback loops, 27.2% (47/173) are oscillators, a very high proportion compared to 6.7% (221/3,325) for the entire design space and 9.3% (146/1,575) and 8.0% (6/75) for the AI or Rep loops alone (Figure S4A). Figure S4B shows how MCTS allocates samples to each of these categories over a total sampling time of 10⁵ and 10⁶ iterations (mean of n = 50 and n = 12 replicates, respectively). During an initial exploratory phase, samples are taken broadly across categories; however, at N ≈ 6 × 10⁴, the AI-Rep combination becomes heavily favored, followed by Rep alone, for the rest of sampling time. This transition can be observed directly from the search history by measuring regret. Defined as

regret is the difference between the expected reward from sampling the best topology (with robustness Q∗) and the actual reward over N iterations. As MCTS moves from exploration into exploitation, the reward rate becomes higher and regret accumulates more slowly (Figures S4C and S4D). Thus, the discovery of an enriched region during search triggers a shift to more focused sampling that can be identified in real-time without a priori knowledge about design principles or motifs.

CircuiTree infers oscillator motifs from search results

In addition to finding individual circuit topologies, an important goal of circuit design is to identify specific structural features, or network motifs, that lead to successful designs (Milo et al. (2002); Alon (2007); Ma et al. (2009); Cotterell and Sharpe (2010); Shah and Sarkar (2011); Chau et al. (2012); Lim et al. (2013); Schaerli et al. (2014)). Similarly, in our game-playing paradigm, a motif of the assembly game is a series of moves (an assembly strategy) that leads to a subset of T highly enriched with successful design results. Specifically, we define an assembly motif as any successful topology (Q > Q_thresh) that is also overrepresented as a sub-pattern within other successful topologies, compared to a set of randomly assembled designs. Using a statistical test for overrepresentation (described in detail in the Methods section and shown schematically in Figure S5A), CircuiTree can mine the results of a tree search for putative assembly motifs.

Looking for 3-node oscillator motifs with ≤ 3 interactions, CircuiTree identifies the same four motifs as enumeration, shown at the top of Figure 3: the repressilator motif (Rep) and the activator-inhibitor (AI) loop with either PAR of the activator (AIPAR), constitutive inhibition of the inhibitor (AICI), or constitutive activation of the activator (AICA). CircuiTree finds these minimal motifs in 100% (Rep), 86% (AIPAR), 94% (AICI), and 62% (AICA) of replicates. To see how these motifs are situated in the overall design space, we plot all 221 oscillators in Figure 3 as a “complexity atlas” (Cotterell and Sharpe (2010); Schaerli et al. (2014)), a graph-of-circuits where every oscillator is a node and nodes are organized in layers according to their complexity. Edges connect nodes in adjacent layers if they differ by the addition of a single circuit interaction, analogous to a move in the assembly game. Notably, 216/221 oscillators constitute a large connected component of the graph originating from the four minimal motifs, indicating that almost all 3-node oscillators live in a subset of design space defined by Rep and/or AI motifs. The size of each node in Figure 3 denotes its motif robustness Q_motif, the overall oscillation probability for all circuits containing this motif (in other words, the mean reward once reaching this state in the game). The color of each node reflects its motif discovery rate, measured as the percentage of n = 50 replicates that labeled it as a motif.

Figure 3.

Note that these features correlate visually (large circles tend to be red and vice versa) and quan-titatively, as shown by the scatterplot in Figure S5B, indicating that higher quality motifs are more likely to be found. Additionally, the motifs that are discovered most reliably (discovery rate <80%, indicated by bold edges) form an optimal subset corresponding to the best motifs (the largest nodes). Surprisingly, the most robust oscillator, the repressilator with PAR on all three TFs (Rep+3xPAR, Q = 0.767), is discovered at a rate of 98%. Because Rep+2xPAR is a poor oscillator (ranked #187) and Rep+1xPAR does not oscillate at all, Rep+3xPAR cannot be assembled from Rep (or any other intermediate) without breaking oscillations. Nonetheless, its high motif discovery rate suggests that CircuiTree is capable of identifying special design strategies that require a specific combination of moves. Overall, we find that CircuiTree reliably infers the most robust motifs for a given phenotype, even finding obscure but optimal solutions, and the resulting motifs for 3-node oscillation form a single cluster of optimal designs. Note that, in contrast to classic network motifs which compare the frequency of a pattern in successful topologies against the entire set of topologies (or a comparable null model), we compare frequencies between circuit assemblies that were successful during search and random circuit assemblies. Therefore our method returns motifs that are overrepresented among search results and thereby accounts for the bias inherent to sampling the space of topologies by assembly rather than by flat enumeration.

Motif multiplexing allows five-node oscillators to compensate for deletions and single-gene knockdowns

While synthetic circuits generally rely on a single, minimal circuit module, naturally occurring circuits often use many functional modules. For instance, circadian oscillators across divergent taxa contain two or more oscillatory feedback loops (Lee et al. (2000); Cheng et al. (2001); Bell-Pedersen et al. (2005); Pokhilko et al. (2012)). Could these additional modules make circadian oscillators more resistant to deletions during evolution (Wagner (2005))? To explore this possibility, we implemented a parallelized version of CircuiTree with pruning (described in details in the Methods and illustrated in Figure S6) and used it to design 5-node circuit topologies that oscillate despite a 50% chance of deletion of a random TF (Figure 4A). Please see Methods for a description of the parallel implementation and Figure S6 for a schematic of the pruning.

Figure 4.

After 5 million search iterations (less than four days of real time using 1,000 parallel search threads and 300 CPUs), CircuiTree finds 1,386 putative oscillators with and > 100 samples, the first result arriving in just 1.3 hours. Due to the computational cost of exact simulations of large stochastic systems, only topologies with up to 15 interactions were included. As shown by the heatmap in Figure S7A, when clustered based on their structural similarity (measured as the pairwise graph edit distance between topologies), these topologies separate into sparse and dense clusters (Figure S7B). In Figure 4B, each of these topologies is plotted as a circle based on their observed robustness with or without a random deletion ( and , respectively). The overall robustness is indicated by the color gradient, and the circle size represents the number of samples. The circuit’s fault-tolerance (FT), the proportion of components that can be deleted without losing oscillation, is estimated as , and contour lines for , and are shown as dashed lines. The putative oscillators included 42 topologies with 2 or 3 nodes, all of which had low robustness to deletions . For instance, the second-best oscillator in the 3-node search, the AI+Rep circuit (labeled (i) in Figure 4B and Figure 4C), was found to have a fault tolerance of , consistent with oscillations that persist after deletion of the two unused TFs (D and E) but attenuating with deletion of any of the active components (A, B, or C).

In contrast, fault-tolerant oscillators contain many interleaved oscillatory motifs that activate under different deletion scenarios, a design feature we call motif multiplexing. As shown in Figure 4D, the number of repressilator (Rep) and activator-inhibitor (AI) motifs is a strong predictor of robustness in the presence of deletion, measured as , and higher complexity in general is associated with higher (Figure S7C). To see this design principle, consider the AI+3Rep circuit (labeled (ii) in Figures 4B and 4C), which is an extension of the AI+Rep topology with two additional backup repressilator loops (A-C-D and A-C-E, highlighted in yellow). This circuit has a higher fault tolerance of because, while deletion of gene B is fatal for oscillations in the AI+Rep circuit, the A-C-D motif takes over to rescue oscillations in the AI+3Rep circuit (representative simulated trajectories shown in Figure 4E, upper and middle row; see Table 1 for parameter values). The pattern extends to the 3AI+3PAR circuit (labeled (iii) in Figures 4B and 4C, right diagram) which is similar to the AI+3Rep except for the addition of two activator-inhibitor motifs (B-D and D-C, highlighted in yellow). This circuit similarly activates a repressilator motif (A-C-E) upon deletion of gene B. Now, however, oscillations are rescued after deletion of gene C by activating the B-D motif. Consequently, the 3AI+3PAR circuit has a high fault tolerance of (Figure 4E, bottom row).

During evolution, genomic mutation may lead to a partial reduction in transcription rate rather than a complete knockout. Do multiplexed oscillators maintain their mutational robustness in these conditions? To explore this question, we simulated the 3AI+3Rep circuit under conditions where a single gene is partially knocked down by multiplying its transcription rate by a factor (100 − KD)/100, KD being the percent of knockdown. In Figure 5A, we visualize a representative stochastic simulation of the wild-type (WT) 3AI+3Rep circuit by reducing the system of five TFs (a five-dimensional phase space) to two composite dimensions using principal components analysis and plotting a simulated trajectory on the first and second principal component axes. At each time point, the dominantly expressed TF is indicated by a transparent marker of the same color, and the inset diagram shows the order of TFs activated in the limit cycle (A-B-D-C). As gene A is knocked down from KD=0% to KD=100% (Figure 5B, upper middle panel), the limit cycle gradual drifts in phase space until the trajectory eventually flattens, consistent with the lack of oscillations after deletion of gene A (Figure 4E). During knockdown of gene B (Figure 5B, top right panel), in contrast, the limit cycle drifts before discontinuously jumping to a new limit cycle (the A-E-C repressilator motif) between 75% and 100% KD, as shown in the inset diagram. Similarly, knockdown of gene C induces a gradual drift followed by a jump between 75% and 100% KD to a limit cycle driven by the B-D activator-inhibitor motif. Unlike genes B and C, knockdown of gene D produces no obvious discontinuities as the limit cycle gradually adapts from A-B-D-C to the A-B-C repressilator. Knockdown of gene E produces no appreciable changes in the limit cycle. Overall, the multiplexed oscillator 3AI+3Rep maintains oscillatory behavior during partial knockdown of genes B, C, D, or E and transitions between disparate limit cycles in qualitatively unique ways.

Figure 5.

To quantitatively understand how the system responds to knockdowns, single-gene knock-downs were simulated with a range of values of KD with n = 50 replicates with different random seeds. The quality and frequency of oscillations were then assessed by computing the ACF_min and, if oscillations were detected, the frequency of oscillation. In Figure 5C (upper panel), the ACF_min (mean ± 95% confidence interval) is plotted as a function of KD, and the dashed line indicates the value ACF_min = −0.4 used as a threshold between oscillatory and non-oscillatory dynamics (below and above the line, respectively. Knockdown of gene A leads to a rise in ACF_min and bomes lethal for oscillations above KD = 60%. For gene B, ACF_min increases before suddenly crossing the threshold to indicate loss of oscillations between KD = 85% and 95%, and for gene C, oscillations disappear in the range KD = 65%to90%. For genes D and E, no detrimental effects are observed during knockdown. Thus, oscillations persist in almost all cases (91.9% of samples) when partially knocking down genes B, C, D, or E (Figure 5C, middle panel). Partial knockdowns of A, B, C, and D each have distinct effects on oscillation frequency, a phenomenon called frequency pulling (Heltberg et al. (2021)). The bottom panel of Figure 5C shows how, for KD < 60%, the WT resonant frequency of 6.0hr⁻¹ is pulled higher when knocking down A or B and lower when knocking down C or D (no change for E). At higher values of KD, knockdowns of B and C discontinuously jump from the drifted frequency to a new resonant frequency (2.6hr⁻¹ and 6.6hr⁻¹, respectively), while the D knockdown transitions smoothly to its new frequency of 3.8hr⁻¹. Thus, the 3AI+3Rep oscillator accommodates most partial knockdowns of genes B, C, D, or E by modulating, yet retaining, oscillatory dynamics.

Up to this point, our modeling has assumed that the rates of reactions such as binding, transcription, translation, and degradation are identical between circuit components. To investigate whether the oscillators we discover are fine-tuned for this symmetry, we perturbed the rate parameters for each TF individually. A Gaussian perturbation kernel with standard deviation σ_param was applied to each parameter, scaled to the pre-defined ranges of those parameters (see Table 2 for parameter ranges and Methods for details of the perturbation kernel). Stochastic simulations were performed for a range of values of σ_param between 1 · 10⁻³ and 3 · 10⁻². To assess the effects of parameter perturbation on oscillation, the power spectral density (PSD) of the dynamics of TF A (mean of n = 50 replicates with different random seeds) was calculated for each value of σ_param and plotted as a heatmap (Figure 5D). For σ_param < 10⁻², the PSD shows peaks at the oscillator’s fundamental frequency of 6hr⁻¹ and first and second harmonics (12hr⁻¹ and 18hr⁻¹, respectively). As σ_param increases further to 5 · 10⁻², the fundamental frequency peak gradually diffuses, and the oscillation rate drops from 1.0 to ∼ 0.5, as shown in the line plot in Figure 5D (envelope indicates 95% confidence interval). above σ_param = 5 · 10⁻², the oscillation rate drops to near-zero, and the mean PSD shows no visible peaks. Therefore, this oscillator can tolerate a minimal amount of heterogeneity between TFs, above which it appears somewhat sensitive to asymmetry in at least one parameter of the model.

Discussion

Natural biological networks contain a number of biochemical components that exceeds our current capabilities of engineering by orders of magnitude, underscoring the importance of scalable computational methods for synthetic circuit design and analysis of biological design principles. Inspired by state-of-the-art artificial intelligence, we approach circuit design as a game of step-by-step topology assembly where success is determined by the achievement of a target phenotype in simulation. Similar to game-playing platforms that have achieved superhuman mastery of complex games (Silver et al. (2017)), CircuiTree searches the space of possible circuit topologies using Monte Carlo tree search (MCTS), which balances exploitation of promising circuit assembly moves with exploration of other possibilities (Figure 1). This search strategy is comprehensive enough to infer motifs for a given phenotype (Figure 3) and efficient enough to search large spaces of designs fruitfully with limited samples (Figures 4A and 4B). Finally, we demonstrate CircuiTree’s scalability by characterizing a novel class of five-gene oscillators that use a strategy we call motif multiplexing to resist deletion and single-gene knockdown mutations by densely interleaving multiple sub-oscillators (Figures 4C-4E and 5A-5C). This design principle may grant evolutionary robustness to eukaryotic circadian clocks, which have been observed across many phyla to contain multiple oscillatory loops (Cheng et al. (2001); Bell-Pedersen et al. (2005); Pokhilko et al. (2012)), and it leaves open the question of how multiple sub-oscillators can be coupled without triggering chaos (Heltberg et al. (2021)).

CircuiTree is distinguished by its general framework. Developed for general game-playing and planning problems, MCTS can query very large spaces in search of robust topologies and assembly motifs for any measurable phenotype, without restrictions on the modeling or simulation frame-work, and without needing to enumerate all possible topologies. This property bridges a gap in computational circuit design where generally, methods for finding single topologies (such as evolution (François and Hakim (2004); François and Siggia (2008)), mixed-integer optimization (Otero-Muras and Banga (2016)), and recurrent neural networks (Shen et al. (2021))) do not generalize easily to design principles and/or require a specific mathematical form, while more comprehensive methods (such as enumeration (Chau et al. (2012); Schaerli et al. (2014)) and Bayesian sampling (Woods et al. (2016))) require an explicit list of topologies. Notably, while CircuiTree currently implements a “vanilla” version of the algorithm, MCTS can be modified in many ways to suit different problems (Browne et al. (2012); (Świechowski. (2022).

This work also highlights the possibilities of using reinforcement learning to address difficult combinatorial design and inference problems in biology. In future work, CircuiTree could be applied to combinatorial design of therapeutics, both in the sense of optimizing chimeric receptor design (Daniels et al. (2022)) or combinatorial antigen recognition (Dannenfelser et al. (2020); Williams et al. (2020)), and in the sense of generating novel multicellular therapies consisting of multiple engineered cell types that collaborate in situ, analogous to a natural immune system. Our work could also be used to generate synthetic morphogenesis circuits, which are difficult to design due to the many possible combinations of chemical and physical components (Davies (2017); Toda et al. (2018); Fleischer and Barr (1997)) and a high computational cost per simulation. Given a reward function that measures goodness-of-fit, CircuiTree could also be extended as an inference tool — for instance, for inferring transcriptional regulatory networks from data. More generally, this computational platform extends the study of design principles (Lim et al. (2013)), which has been limited to small motifs of 2 or 3 components, to a larger space of biological networks, presenting an opportunity to dissect the algorithms and strategies used by biology to assemble complex networks at scale.

Methods

Modeling and simulation

Transcription factor (TF) cooperativity was modeled as a slower unbinding rate when both response elements (REs) for the a promoter are bound by the same TF. For K genes (each with mRNA and protein species) with cumulative A activation reactions and R repression reactions, the elimination of explicit dimerization reduces the number of species from 2K + K² + K³ (mRNAs, TFs, TF-TF complexes, and TF-TF-RE complexes) to 2K +A+R (mRNAs, TFs, and TF-RE complexes) and the number of reactions accordingly from 4K + 2K² + 3K³ to 4K + 3A + 3R.

Stochastic trajectories were simulated using Gillespie’s exact method (Gillespie (1977)) and saved at n_t = 2000 time points at intervals of dt = 20 sec. The separation of time-scales between fast binding-unbinding kinetics and the other reactions creates a long simulation time (>2 mins) that cannot be alleviated with, for example, traditional τ-leaping (Cao et al. (2004)). Thus, instead of using realistic but impractical binding and unbinding rates, these rates were set to virtual values determined from the equilibrium constant for first-binding K_D,1 by the solving the equations

where the value Z = 100 was chosen heuristically to be large enough to maintain a separation of time-scales at low quantities. For each protein species, high-frequency noise was filtered from the stochastic signal using a 9-point binomial filter (Aubury and Luk (1996)), and the quality of oscillations overall was determined by computing the normalized autocorrelation function for each TF and finding the lowest minimum among TFs (ACF_min), excluding the bounds. Oscillations are considered present if this quantity, related to the dissipation constant for oscillations (Otero-Muras and Banga (2016)), is below a cutoff ACF_thresh = −0.4.

For simulations with partial knockdown of a gene, all transcriptional rate parameters for that gene were multiplied by a coefficient on [0, 1] (for example, KD=80% was achieved using a coefficient of 0.2.), and all species were initialized with zero quantity. For all other simulations, the initial quantity of each TF was selected from a Poisson distribution with a mean of 10 proteins, and all other species were initialized with zero quantity.

For perturbation studies, the eight sampled dimensionless variables were converted to values on [0, 1] by normalization to the upper and lower values in Table 2. These values were stored in a K × 8 matrix, where each row represented one of the K = 5 TFs, and each value was perturbed independently with a Gaussian kernel of standard deviation σ_param, truncated to the range [0, 1] to prevent values outside the reasonable range. The perturbed values were then converted to rate parameters for each TF as outlined above.

Random sampling

Random sets of the 10 rate parameters listed in Table 1 were generated by drawing samples of the 8 dimensionless variables described in Table 2. Latin Hypercube sampling was used to draw 10⁴ samples from a multivariate uniform distribution, with bounds shown in Table 2. For the 3-node and 5-node cases, every possible initialization (with a unique set of random generator seed, initial protein quantity, and parameters) was stored in a table of 10⁴ rows, and each simulation was initialized with a random row of this table.

Monte Carlo tree search

For all MCTS runs, the hyperparameter in Equation 2 was set to c = 2.00 to encourage exploration (the default value is ). Replicate runs for an experiment were performed using different random seeds, which are used to break ties during the selection phase and draw a random parameter set during the simulation step. For the 3-node search, each reward value q was drawn as a Bernoulli trial with a success probability equal to the Q value found using enumeration rather than running fresh simulations.

MCTS was parallelized using the lock-free method Enzenberger and Müller (2010) and implemented with the Python utility ‘celery’ (https://docs.celeryq.dev/en/stable/index.html). Tree search was parallelized over 1,000 green threads which dispatched simulation jobs to run in parallel on multiple cloud computing instances totaling 300 CPUs. Reward values for each topology-parameter set pair were stored in an in-memory cache to speed up subsequent training runs (epochs). Back-propagation was executed asynchronously with virtual loss, assuming a reward of 0 until the actual reward is returned. To prevent excessive sampling of local optima during the 5-node search, we introduced a form of decision tree pruning we term “node exhaustion.” Once a terminal node (a completed topology) is visited > 10⁴ times, it is considered exhaustively sampled and pruned from the search graph, and a non-terminal node is pruned once all its successors have been pruned. As illustrated in Figure S5, the visits and reward for each in-edge to an exhausted node are subtracted from the parent node. Thus, once a node is marked exhausted, its history is “forgotten” by its predecessors.

Please see the documentation at https://pranav-bhamidipati.github.io/circuitree/index.html for all additional details of implementation, as well as code tutorials and descriptions of the API.

Motif identification

Before testing for motifs based on a tree search, we first determine whether each terminated topology s_j discovered during the search is a successful oscillator by comparing its empirical robustness to the predefined threshold Q_thresh = 0.01. This segregates all topologies into disjoint successful and unsuccessful sets (X and Y, respectively). We then take samples from the null distribution, drawing n_sample = 10⁵ samples from the tree of topologies using random assembly, each time starting at the root of the design tree and choosing random actions until a terminal topology is reached (Figure S5A). The procedure is then repeated, this time rejecting the result unless it is a member of the successful set X (note that this can be computationally expensive if solutions are sparse). Next, a contingency table is constructed for each successful oscillator x_i ∈ X to compare the frequency of observing x_i within samples of the successful subspace against the frequency of observation in the overall space. (see Figure S5A for an illustration of this process). Finally, statistically significant overrepresentation is determined by the χ² independence test with a significance threshold of p < α = 0.05 after Bonferroni correction. To identify motifs based on the results of enumeration, we conduct the same hypothesis test on each oscillator, except using ground-truth frequencies found by enumeration.

Code availability

CircuiTree, written in Python 3.10, is available on GitHub (https://doi.org/10.5281/zenodo.11285522). The code used to run computational experiments, perform analyses, and plot results is available separately at https://doi.org/10.5281/zenodo.11285550.

Acknowledgements

We thank M. Elowitz, J. Bois, A. Barr, L. Morsut, J. Courte, and J. Paulsson for scientific discussions. We thank J. Bois and G. Manela for advice and assistance on scientific software development. We thank J. Gornet for assistance with distributed computing. We are grateful to I.-M. Strazhnik for assistance with figures and illustrations.

Additional files

Supplementary Figures