© 2026 Michal Kobiela, Diego A. Oyarzún, Michael U.Gutmann. This is the author’s version of the work. It is posted here for personal use. Not for redistribution. This manuscript has been submitted for publication.
Sequential Design of Genetic Circuits Under Uncertainty With Reinforcement Learning
arXiv:2605.06552v1 [cs.LG] 7 May 2026
MICHAL KOBIELA, School of Informatics, The University of Edinburgh, UK DIEGO A. OYARZÚN, School of Informatics, School of Biological Sciences, The University of Edinburgh, UK MICHAEL U. GUTMANN, School of Informatics, The University of Edinburgh,, UK The design of biological systems is hindered by uncertainty arising from both intrinsic stochasticity of biomolecular reactions and variability across laboratory or experimental conditions. In this work, we present a sequential framework to optimize genetic circuits under both forms of uncertainty. By employing simulator models based on differential equations or Markov jump processes alongside a reinforcement learning (RL) policy-based approach, our method suggests experiments that adapt to unknown laboratory conditions while accounting for inherent stochasticity. While previous Bayesian methods address uncertainty through iterative experiment–inference–optimization cycles, they typically require computationally expensive inference and optimization steps after each experimental round, leading to delays. To overcome this bottleneck, we propose an amortized approach trained up-front across a distribution of possible uncertain parameters. This strategy sidesteps the need for explicit parameter inference during the design cycle, enabling immediate, observation-based adaptation. We demonstrate our framework on models for heterologous gene expression and a repressilator circuit, showing that it efficiently handles both molecular noise and cross-laboratory variability. CCS Concepts: • Applied computing → Systems biology; Computational biology; • Computing methodologies → Uncertainty quantification; Sequential decision making; Artificial intelligence; Control methods; Modeling and simulation. Additional Key Words and Phrases: Automated, Sequential, Design, Genetic circuits, Amortization, Reinforcement, Policy-based, POMDP, Epistemic, Uncertainty ACM Reference Format: Michal Kobiela, Diego A. Oyarzún, and Michael U. Gutmann. 2026. Sequential Design of Genetic Circuits Under Uncertainty With Reinforcement Learning. 1, 1 (May 2026), 23 pages. https://doi.org/10.1145/nnnnnnn.nnnnnnn
1
Introduction
The ability to design robust biological systems with specified functions is a cornerstone of engineering biology [4]. Engineered systems have applications ranging from the production of therapeutic drugs [19, 40] and industrial chemicals [21] to environmental remediation [18]. Reliable design of such functionalities is hindered by uncertainty about the Authors’ Contact Information: Michal Kobiela, [email protected], School of Informatics, The University of Edinburgh, Edinburgh, UK; Diego A. Oyarzún, [email protected], School of Informatics, and School of Biological Sciences, The University of Edinburgh, Edinburgh, UK; Michael U. Gutmann, [email protected], School of Informatics, The University of Edinburgh,, Edinburgh, UK. Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than the author(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from [email protected]. © 2026 Copyright held by the owner/author(s). Publication rights licensed to ACM. Manuscript submitted to ACM Manuscript submitted to ACM
1
Michal Kobiela, Diego A. Oyarzún, and Michael U. Gutmann
max
Action Uncertainty in θ
Action
Action
θ2
Action
Goal
Expression
Expression
Action
D
Action
Expression
Expression Expression
B
Host overworked
θ1
E
Lab 2
Observation
θ
Lab 1
x
Observation
Observation
Expression
C
Observation
Model
Expression
A
Action
2
Action
x
x
Action
Fig. 1. Overview of the approach using the example task of maximising gene expression in a biological system. A and B: A parameterized model predicts the expected gene expression level in a biological system in response to an action and the goal is to find the action that maximizes gene expression. C: Model parameters can be uncertain and vary between laboratories. Depending on the specific parameters, the optimal action may change. D: Different colored curves represent different possible values of 𝜃 , i.e., the uncertainty, and arrows indicate the corresponding optimal actions. By performing an action and observing the result, uncertainty is reduced, ultimately improving knowledge of the optimal action. E: Traditional Bayesian approaches require updating the belief each time an observation is received, which can be computationally costly (illustrated by the inner loop on the top). In contrast, a policy-based approach can immediately recommend new actions after receiving an observation, effectively bypassing repeated inference and optimization steps. The policy’s goal is to sequentially adapt to unknown parameter values 𝜃 to ultimately achieve maximal expression.
biological system and the exact experimental conditions. A design that works well in one lab may fail in another, making it challenging to predict which design will achieve the desired outcome. To guide the design of genetic systems, computational simulations based on mechanistic models are widely used [4]. These models are typically formulated as systems of differential equations or as Markov jump processes (MJPs) [13] to capture intrinsic stochasticity. They allow experimenters to explore the behavior of genetic systems in silico before performing costly wet-lab experiments. To identify promising system designs, simulations are often combined with optimization strategies [5, 8, 14, 22, 24, 25, 27, 30, 41, 44]. As an example, consider the expression of a foreign gene, where one can control the promoter strength (action) and must balance increasing gene expression against overloading the host (see Figure 1A). Optimization can then be applied to maximize the expected yield of the foreign protein (see Figure 1B). Such optimization approaches enable efficient exploration of the design spaces, but a key limitation lies in handling the uncertainty [20, 33, 36]. Uncertainty can be broadly categorized into two types: aleatoric and epistemic. Aleatoric uncertainty reflects the intrinsic stochasticity of biological systems, e.g. molecular reactions, and is irreducible; it is present in every experiment Manuscript submitted to ACM
Sequential Design of Genetic Circuits Under Uncertainty With Reinforcement Learning
3
and inevitably alters observed results. Epistemic uncertainty, by contrast, arises from incomplete knowledge of the system and the experimental conditions. Epistemic uncertainty is relevant in practice because, for instance, the true effective kinetic rates governing a biological system are rarely known precisely and can vary across experimental settings due to differences in growth conditions, reagent batches, or host strains. This uncertainty causes designs optimized under one set of parameters to fail or perform suboptimally when implemented under slightly different conditions (see Figure 1C). Epistemic uncertainty can be reduced through experimentation, which motivates sequential design, where information from earlier experiments is used to refine predictions and guide the selection of subsequent designs. In some cases, a single iteration where a prior belief on parameter values is updated using experimental data, followed by optimization, is sufficient to obtain good designs [11, 20]. But one single iteration is often not sufficient to disambiguate between several possibly optimal designs, and a sequence of experiments is needed to determine truly optimal one. This is illustrated in Fig. 1D, where a single observation (the purple cross) does not sufficiently reduce the uncertainty to disambiguate the optimal design (arrows indicate possible locations). A natural strategy is therefore to iterate this procedure: after each experiment, update the parameter belief using the newly collected data and re-optimize the design. While conceptually straightforward, this approach introduces computational overhead and delays between experimental cycles, as inference and optimization must be repeated after every observation (see Fig. 1E). This not only results in delays between experiments (as illustrated in Fig. 2), but may also require fine-tuning the inference and optimization algorithms during deployment to ensure that they work reliably for the specific problem at hand. To address this limitation, we propose a naturally sequential alternative that handles uncertainty while eliminating the need for inference and optimization between experimental rounds. We cast experimental design as a partially observable Markov decision process (POMDP), where the true system parameters constitute a latent state that is specific to each laboratory setting, while actions are the controllable components of the system. Although training the policy may require substantial upfront computation, once deployed the algorithm can immediately select the next action after each observation is received (see Fig. 2B). This is enabled by amortization: during training, the policy is exposed to many possible realizations of the model parameters sampled from the prior, allowing it to learn a direct mapping from experimental histories to actions. We evaluate the proposed approach on two representative synthetic biology tasks: maximization of heterologous gene expression and tuning an oscillator toward a desired frequency. Using synthetic datasets, we examine how the learned policy adapts its experimental interventions over sequential rounds in the presence of latent environmental variation, without performing explicit parameter inference. We also investigate the role of prior knowledge by incorporating mechanistic models during training and compare the resulting behavior with uninformed optimization strategies such as Bayesian optimization. 2
Previous work
Previous work applying reinforcement learning (RL) to synthetic biology has largely focused on bioreactor control [39] and optimal experimental design [38], rather than genetic circuit design. In this latter area, early studies used RL with a discrete action space to design a deterministic dynamical system to obtain oscillatory dynamics from random initializations [12]. RL played the role of an optimizer, with the policy operating on model-derived quantities (e.g. bifurcation properties). In our work, the model, in form of simulator, is used during training to learn a policy, but, at deployment, the policy operates directly on experimental observations without requiring access to the underlying Manuscript submitted to ACM
4
Michal Kobiela, Diego A. Oyarzún, and Michael U. Gutmann
Fig. 2. Timelines comparing inference-optimization based approach [20] (A) and amortized policy-based design (B). While the policy requires an initial time investment for training, it can be deployed immediately after each experiment produces an observation, unlike the traditional approach that must perform inference and optimization after every experiment. .
model, enabling direct application in laboratory settings. This formulation is reflected in how actions and observations are defined. Rather than using model-derived quantities (e.g., bifurcation properties), we restrict observations to experimentally measurable outputs, such as gene expression trajectories. This allows the policy trained in simulation to transfer to the deployment phase where experimental data and not model-generated data is being used. We also consider design under both aleatoric and epistemic uncertainty. Aleatoric uncertainty arises from intrinsic stochasticity in biochemical reactions and measurement noise; in our framework, it is modeled either through observational noise or through Markov jump processes that explicitly capture stochastic reaction dynamics. We consider reducible epistemic uncertainty over model parameters, whose values are implicitly inferred from the full history of observations and actions that is directly provided as input to the policy (see Fig. 3A), as opposed to the most recent observation. Furthermore, treating the design parameters as actions allows the policy to select the initial design and explore the entire continuous design space, rather than operating on a discrete action space as e.g. in Giannantoni et al. [12]. Consequently, the policy can leverage all available information to adapt its decisions in the presence of epistemic uncertainty from the very first experiments. This capability is particularly important during early stages of experimentation, where only a limited number of experiments can be performed and prior mechanistic knowledge encoded in the biological model can provide substantial guidance. While general-purpose optimization methods such as Bayesian optimization (BO; see, e.g., Garnett [10]) may become attractive for larger experimental budgets, our empirical results show that, in early-stage experimentation, approaches that exploit prior biological knowledge achieve better-performing designs than methods that rely primarily on data collected during the optimization process. Further contrasting with previous work, we note that conditioning decisions on histories of actions and observations is common in sequential Bayesian optimal experimental design (BOED) approaches [3, 9, 16, 38] and in history-based approaches for handling partial observability in RL [37]. However, the objectives of BOED differ fundamentally from ours; those methods aim to maximize information gain about uncertain system parameters, and therefore may deliberately select experiments that are highly informative about the model but not necessarily useful for achieving a desired functional behavior. As a consequence, the resulting designs may not perform well with respect to the Manuscript submitted to ACM
Experiment1 Action1 x
Policy 1
Yield (reward)
Action Action1 Observation1 Action2 Observation2
Expression
Expression
Sequential Design of Genetic Circuits Under Uncertainty With Reinforcement Learning
1
Policy
+
Experiment2 Action2 x
2
Action next Action
Policy
max Eθ E Policy
5
2
+ ... +
n
Fig. 3. Overview of the methodology. We train a neural network policy that takes all previously tried actions and corresponding observations to recommend the next action. The training objective is to maximize the total yield across 𝑛 experiments by following the policy.
underlying engineering objective. In contrast, our approach directly optimizes for functional performance of the system. From the very beginning of the experimental process, the policy selects designs that aim to perform well under the latent system parameters while simultaneously gathering information that helps refine future decisions. This focus on performance-oriented experimentation aligns more closely with practical engineering goals, where the ultimate objective is to identify designs that reliably achieve a desired behavior rather than to fully characterize the underlying model. 3
Methods
3.1
Genetic circuit design as a POMDP
We study the problem of optimizing the design of a genetic circuit when a simulator model is available that predicts system behavior as a function of a latent parameter vector 𝜃 (e.g., binding and dissociation kinetic rate constants) and controllable design parameters 𝑎 (e.g., promoter strengths or inducer concentrations). We formulate this problem as a Partially Observable Markov Decision Process (POMDP). Some background on POMDPs is provided in Appendix A.1, and Fig. 4 provides a graphical representation of our approach. The simulator defines a mapping 𝑜𝑡 ∼ Ω(𝑎𝑡 , 𝜃 ),
(1)
where 𝑜𝑡 denotes the observable output of the system, such as a gene expression trajectory (e.g., via reporter gene fluorescence). The simulator may be stochastic and non-differentiable (e.g., when implemented via a Markov Jump Process). We assume a prior distribution over latent parameters, 𝜃 ∼ 𝑝 (𝜃 ),
(2)
representing uncertainty about the biological system before experiments are performed. At the start of each episode, parameters are sampled as 𝜃 𝑗 ∼ 𝑝 (𝜃 ) and remain fixed throughout. Manuscript submitted to ACM
6
Michal Kobiela, Diego A. Oyarzún, and Michael U. Gutmann
Fig. 4. The training loop. In each episode, different values of the simulator parameters 𝜃 are sampled from the prior distribution. At each step, the policy uses all previous actions and corresponding observations ℎ𝑖 to select an action 𝑎𝑖 . The simulator predicts the resulting observation 𝑜𝑖 , and the reward function evaluates the outcome, generating a reward 𝑟𝑖 . The reward function is computed solely from observations, except for the host-aware growth case study, where 𝜃 is used directly in the reward (green dotted arrow).
The latent state in the POMDP corresponds to these unknown parameters 𝜃 . At each experimental iteration 𝑡, the agent selects design parameters 𝑎𝑡 and observes 𝑜𝑡 ∼ Ω(𝑎𝑡 , 𝜃 𝑗 ), During deployment, simulated observations are replaced by experimental measurements. We assume a task-specific reward function 𝑅(𝑎𝑡 , 𝜃 𝑗 ) that computes a scalar score 𝑟𝑡 , 𝑟𝑡 = 𝑅(𝑎𝑡 , 𝜃 𝑗 ),
(3)
based on the action taken and the latent state of the environment. While this is the general case, the reward is often only computed from the observation 𝑜𝑡 , e.g., the total expression yield as in Fig. 3 or the deviation from a desired oscillation frequency in Fig. 7. During training, 𝜃 may also be used when defining the reward in certain case studies (see 4.3), but the policy itself does not observe 𝜃 or the reward. Instead, it conditions only on the history of actions and observations, ℎ𝑡 = (𝑎 1, 𝑜 1, . . . , 𝑎𝑡 −1, 𝑜𝑡 −1 ),
(4)
𝑎𝑡 ∼ 𝜋𝜓 (𝑎𝑡 | ℎ𝑡 ),
(5)
and selects actions according to
where 𝜋𝜓 is a policy parameterised by 𝜓 . This formulation ensures that the learned policy can be deployed in real laboratory settings where 𝜃 is unknown and only experimental observations are available. We assume time, budget, or other practical constraints allow for a finite number of design iterations 𝑇 . The goal is to sequentially select actions 𝑎 1, . . . , 𝑎𝑇 that achieve the design objective while accounting for both epistemic uncertainty in 𝜃 and intrinsic stochasticity in the simulator. For small 𝑇 , biologically uninformed methods often fail to gather sufficient information, whereas our approach leverages the simulator to enable more sample-efficient design (see Section 4.1). Additional formal details are provided in Appendix A.2. Manuscript submitted to ACM
Sequential Design of Genetic Circuits Under Uncertainty With Reinforcement Learning 3.2
7
Policy optimization
We optimize the policy using Proximal Policy Optimization (PPO) [35]. The policy 𝜋𝜓 (𝑎𝑡 | ℎ𝑡 ) maps the interaction history to a distribution over design parameters, while a value network 𝑉𝜙 (ℎ𝑡 ) estimates the expected cumulative reward, when the actions of 𝜋𝜓 (𝑎𝑡 | ℎ𝑡 ) are followed. The objective is to maximize the expected return over an episode of 𝑇 experimental iterations, "𝑇 # ∑︁ 𝐽 (𝜓 ) = E𝜃 𝑗 ∼𝑝 (𝜃 ), 𝑎1:𝑇 ∼𝜋𝜓 𝑟𝑡 .
(6)
𝑡 =1
Because the simulator can be stochastic and non-differentiable, gradients cannot be computed by differentiating through it. Instead, the policy gradient estimator is used: "𝑇 # ∑︁ ∇𝜓 𝐽 (𝜓 ) = E ∇𝜓 log 𝜋𝜓 (𝑎𝑡 | ℎ𝑡 ) 𝐴ˆ𝑡 ,
(7)
𝑡 =1
where 𝐴ˆ𝑡 is an advantage estimate computed using the value network. The advantage is defined as 𝐴𝜙 (ℎ𝑡 , 𝑎𝑡 ) := 𝑄𝜙 (ℎ𝑡 , 𝑎𝑡 ) − 𝑉𝜙 (ℎ𝑡 ),
(8)
where 𝑄𝜙 (ℎ𝑡 , 𝑎𝑡 ) is the expected cumulative return starting from history ℎ𝑡 , taking action 𝑎𝑡 at time 𝑡, and thereafter following the policy, and 𝑉𝜙 (ℎ𝑡 ) is the expected cumulative return from ℎ𝑡 under the policy. Intuitively, 𝐴𝜙 (ℎ𝑡 , 𝑎𝑡 ) measures how much better taking action 𝑎𝑡 is compared to the average performance of the policy from the same history. PPO further stabilizes learning through a clipped surrogate objective that prevents large policy updates. Training proceeds over many simulated episodes with parameters 𝜃 𝑗 ∼ 𝑝 (𝜃 ), allowing the policy to implicitly infer hidden system parameters from past observations and adapt experimental designs accordingly. Additional algorithmic details are provided in Appendix A.3. 4
Results
4.1
Maximizing heterologous gene expression with a host-aware simulator
We consider the design of a gene expression system that produces a foreign protein in the Escherichia coli bacterium, which is a central tasks in many synthetic biology applications [26]. The objective is to tune the induction of the expression system to maximize protein output while accounting for limited cellular resources and variability in host physiology, captured by the latent parameters 𝜃 , see Figs. 5A and B. The system is modeled using a host-aware mechanistic framework from Weiße et al. [43] and Nikolados et al. [26], which explicitly models the dynamics of major proteome components 𝑥 ∈ {𝑟, 𝑡, 𝑚, 𝑞, gfp}. Here, 𝑟 denotes ribosomes, 𝑡 transcriptional machinery, 𝑚 metabolic proteins, 𝑞 housekeeping proteins, and gfp the heterologous protein of interest. The model incorporates metabolism, transcription, ribosome binding, translation, dilution, degradation processes: Metabolism 𝑣imp
𝑣cat
𝑠 −−−→ 𝑠 int −−→ 𝑛𝑠 𝑒
Transcription
Ribosome Binding
𝑤𝑥
𝑝𝑟 + 𝑚𝑥 ⇌ 𝑐 𝑥
∅ −−→ 𝑚𝑥
𝑘𝑏𝑥
𝑘𝑢𝑥
Degradation 𝑑𝑚,𝑥
𝑚𝑥 −−−−→ ∅,
𝑑𝑥
𝑝𝑥 −−→ ∅
Translation 𝑣𝑥
𝑛𝑥 𝑒 + 𝑐 𝑥 −−→ 𝑝𝑟 + 𝑐 𝑥 + 𝑝𝑥
(9)
Dilution 𝜆
𝑒, 𝑚𝑥 , 𝑐 𝑥 , 𝑝𝑥 → − ∅ Manuscript submitted to ACM
8
Michal Kobiela, Diego A. Oyarzún, and Michael U. Gutmann
for each protein 𝑥 ∈ {𝑟, 𝑡, 𝑚, 𝑞, gfp}. Some of the rates are functions of other parameters and species: 𝑝𝑡 𝑣 𝑡 𝑠 , 𝐾𝑡 + 𝑠 𝑤 𝑥,max 𝑒 𝑤𝑥 = , 𝜃𝑥 + 𝑒
𝑣 imp =
𝑣 cat = 𝑥 ∈ {𝑟, 𝑡, 𝑚, gfp},
𝑐 𝑥 𝛾 max𝑒 , 𝑛𝑥 𝑒 + 𝐾𝛾 ∑︁ ∑︁ 𝑀= 𝑛𝑥 𝑝 𝑥 + 𝑛𝑟 𝑐 𝑥 ,
𝑣𝑥 =
𝑥
𝑥
𝑝𝑚 𝑣𝑚 𝑠𝑖 , 𝐾𝑚 + 𝑠𝑖
1 𝑒 , 𝜃𝑞 + 𝑒 1 + (𝜌𝑞 /𝐾𝑞 )ℎ𝑞 ∑︁ 𝛾 (𝑒) 𝑐𝑥 , 𝜆= 𝑀
𝑤𝑞 = 𝑤𝑞,max
(10) (11) (12)
𝑥 ∈ {𝑟,𝑡,𝑚,𝑞,gfp}
𝛾 (𝑒) =
𝛾 max e 𝑒 + 𝐾𝛾 .
(13)
In this model, the external nutrient 𝑠 is imported into the cell at rate 𝑣 imp , which depends on the abundance of transport proteins 𝑝𝑡 and follows Michaelis–Menten kinetics with maximum rate 𝑣𝑡 and half-saturation constant 𝐾𝑡 . The internalised substrate 𝑠 int is then catabolized into energy 𝑒 at rate 𝑣 cat by metabolic enzymes 𝑝𝑚 with catalytic parameters 𝑣𝑚 and 𝐾𝑚 , producing 𝑛𝑠 units of energy per substrate molecule. Gene expression is initiated by transcription, where mRNA 𝑚𝑥 is synthesized at rate 𝑤 𝑥 , which increases with energy availability 𝑒 and saturates according to 𝜃 𝑥 ; for housekeeping genes 𝑞, transcription is further regulated by a negative feedback term involving the proteome fraction 𝜌𝑞 , threshold 𝐾𝑞 , and Hill coefficient ℎ𝑞 . Ribosomes 𝑝𝑟 reversibly bind mRNA to form translation complexes 𝑐 𝑥 with association and dissociation rates 𝑘𝑏𝑥 and 𝑘𝑢𝑥 , respectively. Protein synthesis occurs at rate 𝑣 𝑥 , which depends on the number of complexes 𝑐 𝑥 , the elongation rate 𝛾 (𝑒), and the protein length 𝑛𝑥 , with elongation limited by energy through a saturating function with parameters 𝛾 max and 𝐾𝛾 . Both mRNA and proteins are subject to degradation with rates 𝑑𝑚,𝑥 and 𝑑𝑥 , respectively. In addition, all intracellular species are diluted by cellular growth at rate 𝜆, which is computed from the translational activity across the proteome through Eq. (13). Together, these processes capture the interplay between metabolism, gene expression, and growth, and explicitly account for competition for shared resources such as ribosomes and energy. All default parameter values are given in Nikolados et al. [26], which are based on the original study by Weiße et al. [43]. The key parameters of interest in our case study are highlighted in blue above. They are 𝑤 gfp,max , which controls induction of the heterologous protein; the nutrient efficiency 𝑛𝑠 ; and the ribosome–mRNA binding and unbinding rates, 𝑘𝑏 and 𝑘𝑢 of gfp, denoted without protein index for simplicity. In the next section we will also impose growth rate 𝜆 constraints. We built a simulator model based on Ordinary Differential Equations (ODE) derived from the reaction system above to predict heterologous gene expression as a function of the design variable (action) 𝑤 gfp,max and the latent parameters 𝜃 = {𝑛𝑠 , 𝑘𝑢 , 𝑘𝑏 }. The observation 𝑜 is the steady state expression of gfp, calculated by simulating for a sufficiently long time horizon using Rosenbrock23 method from DifferentialEquations.jl library [31]. This steady state value is also the reward, as we aim to maximize the heterologous protein yield. Priors over the latent parameters are specified as 𝑛𝑠 ∼ U (0, 1), 𝑘𝑏 ∼ U (0, 2), and 𝑘𝑢 ∼ N (1, 1) truncated to the interval (0, 2). They are centered around the original values reported in Nikolados et al. [26]. The Normal prior for 𝑘𝑢 is intended to emulate a scenario in which a parameter is, in principle, controllable and partially characterized but not precisely known due to implicit regulation and/or limited knowledge. For a fixed 𝜃 , the model predicts that increasing induction beyond a critical threshold can overload the host, reducing heterologous expression due to resource competition (Fig. 5B). Sampling 𝜃 from the prior reveals substantial variability in the induction–expression relationship (Fig. 5C), demonstrating that both the location of the optimal action and the Manuscript submitted to ACM
Sequential Design of Genetic Circuits Under Uncertainty With Reinforcement Learning
A
Heterologous gene
Expression
B
C
Protein
D θ1
G
H
E θ2
Host overworked
Induction
9
1
F
2
3
4
5
θ3 θ3 Observation Mean
Fig. 5. Design of heterologous gene expression system. A: The host organism E. coli is inserted with foreign (heterologous) gene e.g. GFP. B: Example prediction of the model for fixed parameters 𝜃 , it turns out that increasing gene induction can overload the organism and result in lower expression. C: Variability of induction (action) to expression (observation and reward) curves depending on the sample of 𝜃 from prior belief. D, E and F: The plot showing the actions taken by the policy for 5-step horizon, for three different test values 𝜃 ∼ 𝑝 (𝜃 ). Aggregated results for 1000 test values 𝜃 ∼ 𝑝 (𝜃 ), y-axis shows the normalized regret, i.e. difference between actual expression for each action and optimal expression. H: Comparison with a Bayesian optimization baseline for 100 test samples 𝜃 ∼ 𝑝 (𝜃 ).
scale of the induction–expression curves can vary. To address this, we trained a policy capable of adapting the induction parameter 𝑤 gfp,max to unknown values of 𝜃 to maximize the cumulative yield of the heterologous protein. We trained the policy over a five-step experimental horizon. Fig. 5D–F shows the actions selected by the trained policy for three example latent parameter vectors drawn from the prior 𝜃 𝑖 ∼ 𝑝 (𝜃 ). The policy consistently selects the same initial action, reflecting the lack of prior data to inform the latent parameters at that point. Once the first observation becomes available, the policy begins to adapt to the true value of 𝜃 , balancing exploration and exploitation in subsequent steps and generally converging toward the optimal action by the final step. Notably, while the policy is trained to maximize the cumulative protein yield rather than identify the optimal action, it turns out that the policy does discover the optimal action. Aggregated performance over 1000 test samples of 𝜃 is quantified using the normalized regret, which quantifies the difference between the expression achieved by the policy and the optimal expression for the same system parameters, divided by the optimal expression value. Fig. 5G summarizes this metric across all test cases for the five-step experimental horizon. The low regret observed in the final three steps indicates that the recommended actions closely approach the optimal solution. Manuscript submitted to ACM
10
Michal Kobiela, Diego A. Oyarzún, and Michael U. Gutmann Finally, Fig. 5H compares our approach with a Bayesian optimization baseline (see Appendix C) evaluated on 100
test samples. We plot for each of those samples differences in protein expression obtained with both approaches, aggregated in a box plot. As the differences between our approach and BO tends to be positive for all steps. This comparison demonstrates that our method can efficiently leverage prior knowledge provided by the mechanistic model to recommend better actions compared to an uninformed approach. These results suggest that such an approach holds promise for guiding time-consuming or costly experiments, aiding the design process during early stages of experimentation. 4.2
Maximizing heterologous gene expression with growth rate constraints
We extend the previous study to account not only for heterologous protein production but also for host growth 𝜆 (Fig. 6A), using the same mechanistic model (see Equation 12 for the definition of 𝜆). In this scenario, both the steady state heterologous protein expression 𝑜 expr and the host growth rate 𝑜 growth are assumed observable and the action is the induction parameters 𝑤 gfp,max . The design objective is to maximize protein expression while ensuring that host growth does not fall below a predefined threshold, here set to 0.8. This constraint limits the feasible design space, excluding regions where growth is too low (shaded areas in Fig. 6B–C). The optimal induction action 𝑎 ∗ is determined by balancing high protein expression against the growth constraint. Variability in latent host parameters 𝜃 shifts both the location of the expression peak and the effective growth threshold, making the optimization problem more problematic than unconstrained expression maximization (Fig. 6B–C). To address this, we first trained a regressor that predicts the optimal action 𝑎 ∗ given a value of 𝜃 . For the training of the regressor, we sampled 100 values of 𝜃 from the prior and computed the corresponding induction-expression curves using the model. For each curve, 𝑎 ∗ was identified via grid search while taking into account the growth constraint, which is easily done since 𝑎 is a scalar. Next, we trained a policy using as reward function 𝑅(𝑎, 𝜃 𝑗 ) the absolute distance between 𝑎 and the optimal action 𝑎 ∗ predicted from 𝜃 𝑗 by the regressor. Note that this approach does not enforce that growth is strictly larger than the cutoff, just that it is close-by. Importantly, the history (which serves as the input to the policy) does not include rewards, and therefore 𝜃 is not available to the policy. Instead, the policy must sequentially adapt based only on observed outputs. The reward function is used solely during training via training objective (6) and is not employed during deployment (ensuring no leakege of 𝜃 ). We trained the policy over a five-step experimental horizon. Fig. 6D–F illustrate the actions selected for three example latent parameter configurations 𝜃 ∼ 𝑝 (𝜃 ), showing trajectories of both protein expression and host growth. The results demonstrate that the policy adapts induction levels to maintain growth above the threshold, even when the gene expression peak occurs in a region of low growth. Notably, the effective location of the growth constraint varies with 𝜃 , highlighting the need for sequential adaptation to host-specific dynamics. Aggregated results over 100 test samples of 𝜃 are shown in Fig. 6G–H. Fig. 6G quantifies the deviation from the optimal action, while Fig. 6H shows the corresponding host growth values. Across all cases, growth remains near or above the threshold. Overall, the method remains promising, effectively adapting induction levels to maximize heterologous protein expression even when constraints on host growth are imposed. 4.3
Designing a genetic oscillator under stochastic and parametric uncertainty
We consider the design of a genetic oscillator based on the system in Fig. 7A, in which three genes repress each other in a ring configuration. This circuit, termed the repressilator, has been shown to display oscillatory protein expression Manuscript submitted to ACM
Sequential Design of Genetic Circuits Under Uncertainty With Reinforcement Learning
11
θ1
Growth
1
H
3
2
3
4
5
Growth < 0.8
1
C
2
4
5
Growth < 0.8
Fig. 6. Host-aware heterologous gene expression with minimal impact. A: The host-aware model can predict not only protein expression but also the growth of the host. Here, we assume that both expression and growth are observable. B and C: In this work, we assume that growth should not fall below a certain threshold (specifically 0.8). Therefore, the goal is to maximize expression outside the region where growth is too low, indicated by the shaded regions. The optimal action is denoted by 𝑎 ∗ . D, E, and F: These plots show the actions taken by the policy over a 5-step horizon for three different test values 𝜃 ∼ 𝑝 (𝜃 ), displaying both expression and growth. Depending on the uncertain parameters 𝜃 , the location of the threshold can change, making the problem challenging. The optimal location can shift both due to differences in the position of the expression peak and variations in the threshold location. G and H: Aggregated results from 100 test samples are presented, showing both the distance from the optimal action and the growth value. Growth generally remains close to the threshold.
dynamics [6]. In this case study, the objective is to select design parameters that yield oscillations with a target frequency (Fig. 7B), while accounting for both intrinsic stochasticity and uncertainty in system parameters. As a preliminary step, we first verified that policy gradient methods can effectively optimize repressilator designs in the absence of parametric uncertainty. In this setting, the task reduces to optimizing circuit behavior under intrinsic stochasticity alone, similar to prior work [36]. We found that our approach reproduces comparable solutions while remaining computationally efficient, converging in competitive time using only a single CPU thread (see Appendix B for details). This confirms that policy gradient methods provide a viable alternative to traditional optimization approaches in stochastic biochemical systems. Building on this result, we extended the setting to incorporate uncertainty over system parameters. To this end, the system is modeled using a stochastic reaction network adapted from Sequeiros et al. Manuscript submitted to ACM
12
Michal Kobiela, Diego A. Oyarzún, and Michael U. Gutmann
[36], capturing transcription, translation, and degradation processes for each gene 𝑖 ∈ {1, 2, 3}: mRNA Production 𝛼1
Protein Production
mRNA Degradation
Protein Degradation
𝛾𝑚
𝑚 1 −−→ ∅
𝑝 1 −−→ ∅
𝑚 2 −−→ ∅
𝛾𝑚
𝑝 2 −−→ ∅
𝛾𝑚
𝑝 3 −−→ ∅, !
𝑘𝑋
𝑘𝑋
𝑘𝑋
𝑚 3 −−→ ∅
Gene 1
∅ −−→ 𝑚 1
𝑚 1 −−→ 𝑚 1 + 𝑝 1
Gene 2
∅ −−→ 𝑚 2
𝛼2
𝑚 2 −−→ 𝑚 2 + 𝑝 2
Gene 3
∅ −−→ 𝑚 3
𝛼3
𝑚 3 −−→ 𝑚 3 + 𝑝 3 !
with 𝛼 1 = 𝑘𝑚 𝜖 + (1 − 𝜖)
1 𝐻 𝑝 1+ 𝐾3
𝛾𝑋
𝛾𝑋
𝛾𝑋
! , 𝛼 2 = 𝑘𝑚 𝜖 + (1 − 𝜖)
1 𝐻 𝑝 1+ 𝐾1
(14)
and 𝛼 3 = 𝑘𝑚 𝜖 + (1 − 𝜖)
1 𝐻 𝑝 1+ 𝐾2
.
The system is extended to a stochastic reaction network describing the coupled dynamics of transcription, translation, and degradation for three genes 𝑖 ∈ {1, 2, 3}. Gene expression is modeled through repressing regulatory interactions, where mRNA production for each gene occurs at rate 𝛼𝑖 , with repression mediated by the protein product of another gene. Here, 𝑘𝑚 denotes the maximal transcription rate, 𝐾 is the repression threshold, 𝐻 is the Hill coefficient controlling the steepness of the regulatory response, and 𝜖 represents basal transcription that ensures a nonzero expression level 𝑘𝑋
even under full repression. Protein synthesis follows translation reactions of the form 𝑚𝑖 −−→ 𝑚𝑖 +𝑝𝑖 , where each mRNA acts catalytically to produce protein 𝑝𝑖 at rate 𝑘𝑋 . mRNA and proteins are linearly diluted by cell growth with rate constants 𝛾𝑚 and 𝛾𝑋 , respectively. The full system is simulated as a Markov jump process, thereby capturing intrinsic (aleatoric) stochasticity arising from the probabilistic nature of molecular reaction events [7]. We used Catalyst.jl [23] and JumpProcesses.jl [31] to derive and simulate the system form the above reactions using stochastic simulation algorithm (SSA) [13]. In addition to this intrinsic noise, the model incorporates parametric (epistemic) uncertainty by treating the parameter vector 𝜃 = {𝐻, 𝛾𝑋 , 𝛾𝑚 , 𝜖} as unknown but bounded, with prior distributions reflecting limited prior knowledge: 𝐻 ∼ U (3, 7), 𝛾𝑋 ∼ U (0.8, 1.1), 𝛾𝑚 ∼ U (4, 50), and 𝜖 ∼ U (0.05, 0.15). The design variables (actions) are given by 𝑎 = {𝑘𝑋 , 𝑘𝑚 , 𝐾 }, which are selected by the policy within predefined ranges: 𝑘𝑋 ∈ [100, 1000],
𝑘𝑚 ∈ [3, 120],
𝐾 ∈ [10, 200].
The simulator generates stochastic gene expression trajectories conditioned on the chosen design 𝑎 and latent parameters 𝜃 . The observation 𝑜 consists of the time series of protein concentration (we use 𝑝 1 (𝑡)). From this trajectory, we estimate the oscillation frequency 𝑓 (𝑜) by identifying the location 𝜏 of the second peak of the normalized autocorrelation function 𝐶 (𝜏). The reward function is defined as the weighted sum of two terms: 2 𝑅(𝑎, 𝜃 ) = − 𝑓 (𝑜) − 𝑓 ★ + 𝜆 𝐶 (𝜏2 ), with 𝜆 = 0.3. The first terms penalises the deviation of 𝑓 (𝑜) from the target frequency 𝑓 ★. The second term is the value of the second peak, as used by [36]; it encourages more regular and sustained oscillatory behavior of the design. To improve training stability, we sampled 500 different values of the uncertain latent parameters 𝜃 and use them to normalize (rescale) the loss. This ensures that the optimization objective is well-balanced across parameter variability and prevents certain regions of the latent space from dominating the gradients. Sampling 𝜃 from the prior reveals substantial variability in oscillatory behavior (Fig. 7C), demonstrating that both the achievable frequency and the regularity of oscillations depend strongly on the underlying parameters. Even for a fixed 𝜃 , repeated simulations exhibit variability due to intrinsic noise (Fig. 7D). Manuscript submitted to ACM
Sequential Design of Genetic Circuits Under Uncertainty With Reinforcement Learning
B
A
C
D
θ2
θ1 ω2
a1 θ1
a1 ω1
θ3
ω3
D Observation1
CNN
+ Action1
Observation2
CNN
+ Action2
Observation5
CNN
+ Action5
Step 1
Step 3
next Action
ANN Step 5
H
E
13
F
G
Fig. 7. A: Schematic representation of the oscillator circuit – repressilator. Three genes repress each other in a circular fashion. B: The goal is to design an oscillator with a specific frequency. C: The same action applied to different uncertain parameters can result in very diverse responses (epistemic uncertainty). D: Even if the uncertain parameters are fixed to specific values, the same action can still result in highly variable responses depending on the irreducible randomness due to molecular noise, represented by different random seeds 𝜔. In this case, however, the frequency and amplitude are much less variable compared to epistemic uncertainty. E: Neural network architecture of the policy. The observations are first processed with a shift-invariant convolutional layer to reduce phase dependency, which are then concatenated with corresponding actions and fed into a standard MLP, which outputs the next action. F,G,H: Three examples of the design process, showing the observations at each step. In the first step, the policy does not receive any information yet and suggests a solution in the presence of both aleatoric and epistemic uncertainty. In later steps, observations are provided to the policy, reducing epistemic uncertainty and allowing the policy to recommend better actions. The final step does not always result in the desired frequency; however, it is much closer to the target, and oscillations are often more regular than in the initial designs. I: Aggregated results showing how the reward function evolves for 100 test points of 𝜃 . In each step, both the median (middle lines) and mean (triangles) improve.
To achieve the desired design goal in presence this kind of uncertainty, we trained a sequential decision-making policy that adapts to unknown values of 𝜃 while accounting for stochastic observations. Observations were processed Manuscript submitted to ACM
14
Michal Kobiela, Diego A. Oyarzún, and Michael U. Gutmann
using a convolutional neural network with shift-invariant layers to ensure phase invariance, and the extracted features are concatenated with the corresponding actions before being passed to a fully connected network that outputs the next action (Fig. 7D). We trained the policy over a 5-step experimental horizon. Fig. 7E–G shows the actions selected by the policy for three example latent parameter configurations 𝜃 ∼ 𝑝 (𝜃 ). As in the previous case study, the policy initially selects the same action due to the absence of information about 𝜃 . Once the first observation becomes available, the policy begins to adapt, progressively reducing epistemic uncertainty and refining its recommendations. While early actions may produce irregular or off-target oscillations, later steps yield trajectories that more closely match the desired frequency and exhibit more regular behavior. In this case, due to broad priors, although the final designs improve upon the initial ones, it may not be possible to completely match the frequency. In the appendix, we present a simplified case study in which it is more feasible to achieve the target frequency across the prior range. We also provide additional results demonstrating robustness to random initializations, along with comparisons to oracle policies informed by the ground-truth parameters, which allow us to assess how well the policy adapts relative to the optimal value. Aggregated performance over 100 test samples of 𝜃 is summarized in Fig. 7H. We report the reward value at each step, showing consistent improvement in both mean and median performance. These results demonstrate that sequential, closed-loop design enables robust control of stochastic genetic circuits, effectively handling both intrinsic noise and uncertainty in system parameters, and extending the applicability of model-based design to realistic biological systems. 5
Discussion
In this paper, we present an reinforcement learning approach to sequentially optimize genetic circuits under epistemic and aleatoric uncertainty. The case studies presented here illustrate that amortized sequential design policies can provide an effective and practical approach for optimizing biological systems under uncertainty, as demonstrated using synthetic data. Compared to more traditional approaches that iterate experiment-inference-optimization blocks, our approach avoids repeated inference and optimization during the experimental phase, decoupling the computational effort from the experimental process. Once trained, the policy can propose design actions in real time, circumventing potentially costly and computationally difficult steps of Bayesian inference and optimization. One limitation to our study is the amortization gap: the policy may not perfectly match the optimal design for every possible latent parameter configuration. One way to mitigate this is through semi-amortized approaches, similar in spirit to work [17] in Bayesian optimal design, where the policy is updated sequentially using posterior information. While this can improve adaptation to the true system, it comes at the cost of additional computation, as inference must be performed and the policy potentially retrained. Exploring such semi-amortized strategies is an interesting avenue for future work. Another important consideration is that all experiments here were conducted using synthetic data. While synthetic evaluations allow controlled validation of the approach, real experimental systems may introduce additional sources of model mismatch or unmodeled dynamics. Handling potential model mismatch robustly is therefore an important future direction, particularly in biological settings where simulator fidelity may be limited. Our approach provides a natural framework for the low-data, uncertain scenarios, as they can combine prior knowledge and simulation-based inference with sparse observations. Furthermore, policy gradient methods offer an alternative to standard optimization approaches for designing genetic circuits under molecular noise such as Sequeiros et al. [36] when there is no epistemic uncertainty. This can be particularly useful in systems like the repressilator, where Markov Jump Process simulations are relatively fast but intrinsic noise is high, making traditional optimization Manuscript submitted to ACM
Sequential Design of Genetic Circuits Under Uncertainty With Reinforcement Learning
15
challenging. Because policy gradient methods are local and do not require simulator gradients, they likely can scale to high-dimensional design spaces. Overall, our results suggest that amortized sequential design provides an attractive tool, enabling rapid decisionmaking while maintaining flexibility to handle stochasticity, parametric uncertainty, and limited experimental data, while utilizing prior knowledge in form of the mechanistic model. Future work addressing model mismatch, semi-amortized updates, and validation in real experimental systems will further strengthen the applicability of this approach. Code Availability The code, data, and trained models used in this study are publicly available at: https://github.com/MichalKobiela/ GeneCircuitsRL . Acknowledgments This work was supported by the United Kingdom Research and Innovation (grant EP/S02431X/1), UKRI Centre for Doctoral Training in Biomedical AI at the University of Edinburgh, School of Informatics. For the purpose of open access, the author has applied a creative commons attribution (CC BY) licence to any author accepted manuscript version arising. During the preparation of this work the authors used ChatGPT in order to increase the readability of the text. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication. References [1] Joshua Achiam. 2018. Spinning Up in Deep Reinforcement Learning. (2018). [2] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah. 2017. Julia: A fresh approach to numerical computing. SIAM Review 59, 1 (2017), 65–98. doi:10.1137/141000671 [3] Tom Blau, Edwin V Bonilla, Iadine Chades, and Amir Dezfouli. 2022. Optimizing sequential experimental design with deep reinforcement learning. In International conference on machine learning. PMLR, 2107–2128. [4] Jennifer AN Brophy and Christopher A Voigt. 2014. Principles of genetic circuit design. Nature methods 11, 5 (2014), 508–520. [5] Madhukar S Dasika and Costas D Maranas. 2008. OptCircuit: an optimization based method for computational design of genetic circuits. BMC systems biology 2, 1 (2008), 1–19. [6] Michael B Elowitz and Stanislas Leibler. 2000. A synthetic oscillatory network of transcriptional regulators. Nature 403, 6767 (2000), 335–338. [7] Michael B Elowitz, Arnold J Levine, Eric D Siggia, and Peter S Swain. 2002. Stochastic gene expression in a single cell. Science 297, 5584 (2002), 1183–1186. [8] Maurice Filo, Nicolò Rossi, Zhou Fang, and Mustafa Khammash. 2026. GenAI-Net: A Generative AI Framework for Automated Biomolecular Network Design. arXiv preprint arXiv:2601.17582 (2026). [9] Adam Foster, Desi R Ivanova, and Ilyas Malik. [n. d.]. Deep Adaptive Design: Amortizing Bayesian Experimental Design. Variational, Monte Carlo and Policy-Based Approaches to Bayesian Experimental Design ([n. d.]), 139. [10] Roman Garnett. 2023. Bayesian optimization. Cambridge University Press. [11] Adam R Gerlach, Andrew Leonard, Jonathan Rogers, and Chris Rackauckas. 2020. The Koopman expectation: An operator theoretic method for efficient analysis and optimization of uncertain hybrid dynamical systems. arXiv preprint arXiv:2008.08737 (2020). [12] Leonardo Giannantoni, Alessandro Savino, and Stefano Di Carlo. 2023. Optimization of synthetic oscillatory biological networks through Reinforcement Learning. In 2023 IEEE International Conference on Bioinformatics and Biomedicine (BIBM). IEEE, 2624–2631. [13] Daniel T Gillespie. 1977. Exact stochastic simulation of coupled chemical reactions. The journal of physical chemistry 81, 25 (1977), 2340–2361. [14] Tom W Hiscock. 2019. Adapting machine-learning algorithms to design gene circuits. BMC bioinformatics 20 (2019), 1–13. [15] J. D. Hunter. 2007. Matplotlib: A 2D graphics environment. Computing in Science & Engineering 9, 3 (2007), 90–95. doi:10.1109/MCSE.2007.55 [16] Desi R Ivanova, Adam Foster, Steven Kleinegesse, Michael U Gutmann, and Thomas Rainforth. 2021. Implicit deep adaptive design: Policy-based experimental design without likelihoods. Advances in neural information processing systems 34 (2021), 25785–25798. [17] Desi R Ivanova, Marcel Hedman, Cong Guan, and Tom Rainforth. 2024. Step-DAD: Semi-amortized policy-based bayesian experimental design. In ICLR 2024 Workshop on Data-centric Machine Learning Research (DMLR), Vol. 2. 21. Manuscript submitted to ACM
16
Michal Kobiela, Diego A. Oyarzún, and Michael U. Gutmann
[18] Tian Jiang, Veronica A Montgomery, Karuna Jetty, Vijaydev Ganesan, Matthew R Incha, John M Gladden, Nathan J Hillson, and Di Liu. 2025. Metabolic engineering and synthetic biology for the environment: from perspectives of biodetection, bioremediation, and biomanufacturing. Biotechnology for the Environment 2, 1 (2025), 14. [19] Linda M Keefer, Marie-Agnès Piron, and Pierre De Meyts. 1981. Human insulin prepared by recombinant DNA techniques and native human insulin interact identically with insulin receptors. Proceedings of the National Academy of Sciences 78, 3 (1981), 1391–1395. [20] Michal Kobiela, Diego A Oyarzún, and Michael U Gutmann. 2026. Risk-averse optimization of genetic circuits under uncertainty. Cell Systems 17, 1 (2026). [21] Sang Yup Lee, Hyun Uk Kim, Tong Un Chae, Jae Sung Cho, Je Woong Kim, Jae Ho Shin, Dong In Kim, Yoo-Sung Ko, Woo Dae Jang, and Yu-Sin Jang. 2019. A comprehensive metabolic map for production of bio-based chemicals. Nature Catalysis 2, 1 (Jan. 2019), 18–33. doi:10.1038/s41929-018-0212-4 [22] Zhengda Li, Shixuan Liu, and Qiong Yang. 2017. Incoherent inputs enhance the robustness of biological oscillators. Cell systems 5, 1 (2017), 72–81. [23] Torkel E. Loman, Yingbo Ma, Vasily Ilin, Shashi Gowda, Niklas Korsbo, Nikhil Yewale, Chris Rackauckas, and Samuel A. Isaacson. 2023. Catalyst: Fast and flexible modeling of reaction networks. PLOS Computational Biology 19, 10 (10 2023), 1–19. doi:10.1371/journal.pcbi.1011530 [24] Wenzhe Ma, Ala Trusina, Hana El-Samad, Wendell A Lim, and Chao Tang. 2009. Defining network topologies that can achieve biochemical adaptation. Cell 138, 4 (2009), 760–773. [25] Charlotte Merzbacher, Oisin Mac Aodha, and Diego A Oyarzún. 2023. Bayesian Optimization for Design of Multiscale Biological Circuits. ACS Synthetic Biology (2023). [26] Evangelos-Marios Nikolados, Andrea Y Weiße, and Diego A Oyarzún. 2021. Prediction of Cellular Burden with Host–Circuit Models. In Synthetic Gene Circuits. Springer, 267–291. [27] Irene Otero-Muras and Julio R Banga. 2017. Automated design framework for synthetic biology exploiting pareto optimality. ACS Synthetic Biology 6, 7 (2017), 1180–1193. [28] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. 2019. Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems 32 (2019). [29] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. 2011. Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research 12 (2011), 2825–2830. [30] Lingxia Qiao, Wei Zhao, Chao Tang, Qing Nie, and Lei Zhang. 2019. Network topologies that can achieve dual function of adaptation and noise attenuation. Cell systems 9, 3 (2019), 271–285. [31] Christopher Rackauckas and Qing Nie. 2017. Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia. Journal of Open Research Software 5, 1 (2017), 15. [32] Antonin Raffin, Ashley Hill, Adam Gleave, Anssi Kanervisto, Maximilian Ernestus, and Noah Dormann. 2021. Stable-Baselines3: Reliable Reinforcement Learning Implementations. Journal of Machine Learning Research 22, 268 (2021), 1–8. http://jmlr.org/papers/v22/20-1364.html [33] Tobias Schladt, Nicolai Engelmann, Erik Kubaczka, Christian Hochberger, and Heinz Koeppl. 2021. Automated Design of Robust Genetic Circuits: Structural Variants and Parameter Uncertainty. ACS Synthetic Biology 10, 12 (Dec. 2021), 3316–3329. doi:10.1021/acssynbio.1c00193 [34] John Schulman, Philipp Moritz, Sergey Levine, Michael Jordan, and Pieter Abbeel. 2015. High-dimensional continuous control using generalized advantage estimation. arXiv preprint arXiv:1506.02438 (2015). [35] John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov. 2017. Proximal policy optimization algorithms. arXiv preprint arXiv:1707.06347 (2017). [36] Carlos Sequeiros, Carlos Vázquez, Julio R Banga, and Irene Otero-Muras. 2023. Automated Design of Synthetic Gene Circuits in the Presence of Molecular Noise. ACS Synthetic Biology 12, 10 (2023), 2865–2876. [37] David Silver and Joel Veness. 2010. Monte-Carlo Planning in Large POMDPs. In Advances in Neural Information Processing Systems, J. Lafferty, C. Williams, J. Shawe-Taylor, R. Zemel, and A. Culotta (Eds.), Vol. 23. Curran Associates, Inc. https://proceedings.neurips.cc/paper_files/paper/2010/ file/edfbe1afcf9246bb0d40eb4d8027d90f-Paper.pdf [38] Neythen J Treloar, Nathan Braniff, Brian Ingalls, and Chris P Barnes. 2022. Deep reinforcement learning for optimal experimental design in biology. PLOS Computational Biology 18, 11 (2022), e1010695. [39] Neythen J Treloar, Alex JH Fedorec, Brian Ingalls, and Chris P Barnes. 2020. Deep reinforcement learning for the control of microbial co-cultures in bioreactors. PLoS computational biology 16, 4 (2020), e1007783. [40] Jean-Yves Trosset and Pablo Carbonell. 2015. Synthetic biology for pharmaceutical drug discovery. Drug design, development and therapy (2015), 6285–6302. [41] Babita K Verma, Ahmad A Mannan, Fuzhong Zhang, and Diego A Oyarzún. 2021. Trade-offs in biosensor optimization for dynamic pathway engineering. ACS synthetic biology 11, 1 (2021), 228–240. [42] Michael L. Waskom. 2021. seaborn: statistical data visualization. Journal of Open Source Software 6, 60 (2021), 3021. doi:10.21105/joss.03021 [43] Andrea Y Weiße, Diego A Oyarzún, Vincent Danos, and Peter S Swain. 2015. Mechanistic links between cellular trade-offs, gene expression, and growth. Proceedings of the National Academy of Sciences 112, 9 (2015), E1038–E1047. [44] Mae L Woods, Miriam Leon, Ruben Perez-Carrasco, and Chris P Barnes. 2016. A statistical approach reveals designs for the most robust stochastic gene oscillators. ACS synthetic biology 5, 6 (2016), 459–470. Manuscript submitted to ACM
Sequential Design of Genetic Circuits Under Uncertainty With Reinforcement Learning A
17
Methods details
A.1
Formal POMDP background
A Partially Observable Markov Decision Process (POMDP) [37] is defined by the tuple (S, A, O,𝑇 , Ω, 𝑅, 𝜌 0, 𝛾), where • S is the set of latent states 𝑠 ∈ S, • A is the set of actions 𝑎 ∈ A, • O is the set of observations 𝑜 ∈ O, • 𝜌 0 (𝑠) is the initial state distribution, • 𝑇 (𝑠𝑡 +1 | 𝑠𝑡 , 𝑎𝑡 ) is the state transition model, • Ω(𝑜𝑡 +1 | 𝑠𝑡 +1, 𝑎𝑡 ) is the observation model, • 𝑅(𝑠𝑡 , 𝑎𝑡 , 𝑠𝑡 +1 ) is the reward function, • 𝛾 ∈ [0, 1] is the discount factor. At the start of an episode the latent state is sampled as 𝑠 ∼ 𝜌 0 (·).
(15)
At each timestep 𝑡, the agent selects an action 𝑎𝑡 , after which the environment evolves according to 𝑠𝑡 +1 ∼ 𝑇 (· | 𝑠𝑡 , 𝑎𝑡 ),
(16)
𝑜𝑡 +1 ∼ Ω(· | 𝑠𝑡 +1, 𝑎𝑡 ),
(17)
𝑟𝑡 +1 = 𝑅(𝑠𝑡 , 𝑎𝑡 , 𝑠𝑡 +1 ).
(18)
Because the latent state is not directly observable, decisions are conditioned on the action–observation history ℎ𝑡 = (𝑎 0, 𝑜 1, 𝑎 1, 𝑜 2, . . . , 𝑎𝑡 −1, 𝑜𝑡 ),
(19)
𝑎𝑡 ∼ 𝜋 (𝑎𝑡 | ℎ𝑡 ).
(20)
and actions are sampled from a policy
A.2
Genetic circuit design as a POMDP
In our setting the latent state 𝑎 corresponds to the vector of model parameters 𝜃 . At the beginning of each episode 𝑗, parameters are sampled from a prior 𝜃 𝑗 ∼ 𝑝 (𝜃 ).
(21)
These parameters remain fixed during the entire episode, reflecting that the underlying biological system does not change between experimental iterations. Consequently, 𝑠𝑡 = 𝜃 𝑗
∀𝑡
(22)
and the transition model becomes deterministic, 𝑇 (𝑠𝑡 +1 | 𝑠𝑡 , 𝑎𝑡 ) = 𝛿 (𝑠𝑡 +1 − 𝑠𝑡 ).
(23)
In the standard POMDP formulation the observation distribution is written as 𝑜𝑡 +1 ∼ Ω(· | 𝑠𝑡 +1, 𝑎𝑡 ).
(24) Manuscript submitted to ACM
18
Michal Kobiela, Diego A. Oyarzún, and Michael U. Gutmann However, since 𝑠𝑡 +1 = 𝑠𝑡 = 𝜃 𝑗 in our problem, the observation depends only on the current action and the fixed latent
parameters. We therefore adopt the equivalent but simpler indexing convention 𝑜𝑡 ∼ Ω(· | 𝜃 𝑗 , 𝑎𝑡 ),
(25)
where 𝑜𝑡 denotes the experimental outcome produced by action 𝑎𝑡 . Under this convention (and by starting indexing from 1 rather than 0), the interaction history becomes ℎ𝑡 = (𝑎 1, 𝑜 1, . . . , 𝑎𝑡 −1, 𝑜𝑡 −1 ).
(26)
Similarly, the reward is defined with a shifted index as 𝑟𝑡 = 𝑅(𝜃 𝑗 , 𝑎𝑡 ),
(27)
which evaluates how well the observed circuit behavior satisfies the design objective. Shifting the reward index is commonly done in some frameworks, e.g. Achiam [1].
A.3
Policy optimization with PPO
We optimize the policy using Proximal Policy Optimization (PPO) [35]. The policy 𝜋𝜓 (𝑎𝑡 | ℎ𝑡 ) is parameterized by 𝜓 and maps the interaction history to a distribution over actions, while a value network 𝑉𝜙 (ℎ𝑡 ) estimates the expected cumulative reward from history ℎ𝑡 . The expected return of a policy is 𝐽 (𝜓 ) = E𝜃 𝑗 ∼𝑝 (𝜃 ), 𝑎1:𝑇 ∼𝜋𝜓
"𝑇 ∑︁
# 𝑟𝑡 (𝜃 𝑗 , 𝑎𝑡 ) .
(28)
𝑡 =1
Because the simulator may be stochastic and non-differentiable, gradients of the reward with respect to actions cannot be computed directly. Instead, PPO uses the policy gradient estimator "𝑇 # ∑︁ ∇𝜓 𝐽 (𝜓 ) = E ∇𝜓 log 𝜋𝜓 (𝑎𝑡 | ℎ𝑡 )𝐴ˆ𝑡 ,
(29)
𝑡 =1
where 𝐴ˆ𝑡 denotes an estimate of the advantage function, typically computed using generalized advantage estimation [34]. To stabilize training, PPO maximizes the clipped surrogate objective h i 𝐿 CLIP (𝜓 ) = E min 𝑟𝑡 (𝜓 )𝐴ˆ𝑡 , clip(𝑟𝑡 (𝜓 ), 1 − 𝜖, 1 + 𝜖)𝐴ˆ𝑡 ,
(30)
where 𝑟𝑡 (𝜓 ) =
𝜋𝜓 (𝑎𝑡 | ℎ𝑡 ) 𝜋𝜓old (𝑎𝑡 | ℎ𝑡 )
(31)
is the probability ratio between the updated and previous policies. Training proceeds by generating trajectories from the simulator using parameters sampled from 𝑝 (𝜃 ) and updating the policy and value networks using stochastic gradient ascent on the surrogate objective 𝐿 CLIP (𝜓 ).
Manuscript submitted to ACM
Sequential Design of Genetic Circuits Under Uncertainty With Reinforcement Learning A.4
19
Neural network architecture for the genetic oscillator
The history is encoded as a stacked sequence of past actions and observations: 𝑥𝑡 = [𝑜 1, 𝑎 1, 𝑜 2, 𝑎 2, . . . , 𝑜 𝑁 , 𝑎 𝑁 ].
(32)
Each 𝑜𝑏𝑠𝑖 is a protein expression trajectory simulated with the MJP, and each 𝑎𝑖 is a vector of design parameters. This structure enables the policy to capture dependencies between past designs and outcomes and infer latent 𝜃 𝑗 . The policy and value networks share a CNN-based feature extractor. Each observation 𝑜𝑖 is processed with a one-dimensional CNN and flattened: 𝑓𝑖 = CNN(𝑜𝑖 ),
(33)
then concatenated with the corresponding action 𝑎𝑖 to form 𝑧𝑡 = [𝑓1, . . . , 𝑓𝑁 , 𝑎 1, . . . , 𝑎 𝑁 ].
(34)
This embedding is fed to fully connected layers to produce the final feature for PPO. A.5
Implementation and training
Training uses Stable Baselines3 [32] PPO with Gaussian policies. The CNN feature extractor has latent dimension 128. Training runs for millions of steps, corresponding to thousands of simulated design–experiment iterations, allowing the policy to progressively infer hidden 𝜃 𝑗 from history ℎ𝑡 and adapt its design strategy (Fig. 3D). Stable Baselines is designed for policies that are conditioned on a single observation. In our setting, however, the policy needs access to a sequence of past observations and actions. To handle this, at each timestep we construct an “effective observation” by stacking the history into a fixed-size representation and padding it with zeros where necessary. This zero-padded history is then passed to the policy. For visualization, we employed Seaborn [42] and Matplotlib [15] within the Python programming language. For the simulations, we utilized the Julia programming language [2], together with the packages DifferentialEquations.jl (with JumpProcesses.jl) [31], and Catalyst.jl [23]. For regressor training, we employed PyTorch [28]. B
Genetic oscillator with biomolecular noise only—no epistemic uncertainty
In this experiment, we replicate the repressilator case study from Sequeiros et al. [36], but instead of relying on a mixed-integer optimization framework and partial integro-differential equations, we employ reinforcement learning (RL) for design optimization and Markov jump processes (MJPs) for simulation. The objective is to identify parameter values that maximize the location of the second peak of the normalized autocorrelation function, used as a proxy for oscillation robustness. Using RL, we obtained a policy that includes the optimal design reported by Sequeiros et al. [36], as shown in Figure 8. Moreover, the distribution of the reward under our learned policy closely matches the distribution obtained using their optimal design. Since our policy is Gaussian, it is also possible to use its mean as a deterministic design instead of sampling. We found that this mean-based design achieves performance comparable to the one reported by Sequeiros et al. [36]. Furthermore, our approach required only 46 minutes on a single CPU thread, while [36] reported 5.5 hours using a GPU. This demonstrates that approaches based on the score derivative trick, such as Proximal Policy Optimization (PPO), can effectively optimize genetic oscillators under biomolecular noise simulated via MJPs. While Manuscript submitted to ACM
20
Michal Kobiela, Diego A. Oyarzún, and Michael U. Gutmann 1.0
1.0
0.8 0.6
0.8
Reward
reward
0.4 0.2 0.0
0.6
0.4
0.2 0.4
0.2
H
0.6 0
20
5.5
6.0
40
60
80
thousand steps
6.5
7.0
100
120
140
Ours (stochastic)
Ours (deterministic)
Sequiros et. al
7.00 6.75 6.50 6.25 6.00 5.75 5.50 27.5
m
25.0 22.5 20.0 17.5 15.0 200
kX
180 160 140 120 100 120
km
110 100 90 200
K
180 160 140 120
H
15
20
m
25 Ours (Stochastic)
100 125 150 175 200 kX Sequiros et. al
90
100 110 km
120 120
140
160 K
180
200
Ours (Deterministic)
Fig. 8. RL-based optimization of the repressilator under biomolecular noise. Top left: Training curve showing improvement of the reward function over time. Top right: Comparison between our stochastic policy, its deterministic mean, and the design reported by Sequeiros et al., showing comparable distributions of the reward. Bottom: Heatmap of sampled designs from the trained policy, its mean and Sequiros et. al design, which is contained within the policy.
such simulations are non-differentiable and thus incompatible with standard gradient-based optimization methods, reinforcement learning provides a viable alternative by enabling efficient optimization in this setting. Manuscript submitted to ACM
Sequential Design of Genetic Circuits Under Uncertainty With Reinforcement Learning C
21
Bayesian optimization baseline
We compare our approach against a Bayesian Optimization (BO) baseline implemented with a Gaussian Process (GP) surrogate model. To make the comparison fair we use the mechanistic model to inform the hyperparameters of the Gaussian process. Namely, we generated 100 input-output curves by sampling latent parameters 𝜃 ∼ 𝑝 (𝜃 ) and evaluating the simulator at 20 equally spaced action-values for each sample 𝜃 . The empirical mean of all observations was computed and subtracted from the data prior to fitting, resulting in a zero-mean dataset; this value was retained as the GP prior mean. The GP was equipped with a radial basis function kernel. The signal variance, length scale, and observation noise variance were estimated by maximizing the marginal likelihood on the generated dataset. BO was run over a five-step horizon. The first action was selected uniformly at random. The next three actions were obtained by maximizing an acquisition function using the default gp_hedge strategy in scikit-learn [29], which randomly alternates between Probability of Improvement, Expected Improvement, and Lower Confidence Bound. The final action was selected as the maximizer of the GP posterior mean. Evaluation was performed over 100 independently sampled latent parameter settings 𝜃 ∼ 𝑝 (𝜃 ) using different random seeds. For each instance, we recorded the objective values obtained by both methods and computed their differences. D
Additional results for the genetic oscillator
Figure 9 illustrates the effect of different random seeds used to initialize the policy and environment. It also compares the multi-step policy to a single-step policy that optimizes the reward marginalized over the prior, highlighting their close agreement in the first step and supporting the interpretation that the policy adapts to the observations rather than only improving performance in a marginal (prior-averaged) sense across steps. Figure 10 shows a comparison between the adaptive policy and four “oracle” baselines trained with access to groundtruth parameter values 𝜃 . Across representative settings of 𝜃 , the adaptive policy performs comparably to these oracle baselines, with some degradation in extreme regions of the parameter space. The figure further highlights that even the oracle policies do not always achieve fully functional designs, reflecting the inherent difficulty of the task.
seed 1
seed 2
single step
Fig. 9. We present the differences between two random seeds used to initialize the policy and environment. Although the results are not identical, their performance is very similar. We also compare these to a single-step policy that maximizes the reward marginalized over the prior. This comparison shows that the first step of the multi-step policy achieves performance comparable to the solution obtained when the prior is marginalized out. This suggests that the policy is not merely finding progressively better solutions in a marginal sense (i.e., averaged over the prior) across steps, but is instead adapting to the unknown values of the uncertain parameters, see also Figure 10.
Manuscript submitted to ACM
22
Michal Kobiela, Diego A. Oyarzún, and Michael U. Gutmann Ground truth
Adaptive
Oracle
θ Prior range
θ
θ
θ
Fig. 10. We compare the performance of the adaptive policy to four “oracle” policies, i.e., policies trained using the ground-truth parameter values rather than a prior distribution. We evaluate these across four different ground-truth settings of 𝜃 : one at the center of the normalized parameter space, (0, 0, 0, 0); one on the boundary but not at a corner, (1, −1, 0, 0); and two at the corners, (−1, −1, −1, −1) and (1, 1, 1, 1). Note that 𝜃 is four dimensional; the figures show a two-dimensional representation. Overall, the adaptive policy discovers solutions that are close in performance to the corresponding oracle policies, with the exception of the corner case (−1, −1, −1, −1). This discrepancy is likely due to the amortization gap, as the volume of the surrounding neighborhood intersected with the parameter space is smallest at the corners (grey shaded areas); in the discussion we suggest possible remedies. Importantly, the oracle policies do not always achieve fully functional designs (i.e., rewards close to 1), indicating that for some values of 𝜃 , perfectly optimal designs may not be attainable.
E
Genetic oscillator: Simplified case study
In this case study, we consider a simplified version of the repressilator case study with more control over the parameters and less epistemic uncertainty. Specifically, here we take the parameter 𝛾𝑋 to be controllable within the interval Manuscript submitted to ACM
Sequential Design of Genetic Circuits Under Uncertainty With Reinforcement Learning
23
[100, 1000], the leakage is fixed to 𝜖 = 0.05 and protein degradation rate is fixed to 𝛾𝑥 = 1.0. The setup is otherwise the same as in the main case study. Figure 11 summarizes the results. The top-left panel shows the learning curve, indicating that the final reward improves steadily as the policy is trained over multiple episodes. The top-right panel plots the reward for each design step in the sequential setting. Unlike in the main case study, the final reward is concentrated around 1, where the desired number of oscillations is generally reasonably achieved. The middle and bottom rows present heatmaps of gene expression levels under specific values of the uncertain parameters (which are unknown to the policy). These figures demonstrate how the policy progressively adapts its strategy: for instance, in the middle row, the first design produces too many oscillations, but the policy adjusts its suggestions in subsequent steps and converges to a design yielding the desired target of seven oscillations.
Fig. 11. Simplified repressilator study. The policy adapts to unknown components of the system to achieve seven oscillations. Top left: Final reward increases as the policy is trained. Top right: Reward for each design in the sequential setting (design 1: before any observations; design 2: after observing first outcome, etc.). Middle: Heatmaps showing gene expression for the first, third, and fifth designs under one realization of the uncertain parameters. Bottom: Same as middle, but under a different realization of 𝜃 .
Received n/a; revised n/a; accepted n/a Manuscript submitted to ACM