Strain design optimization using reinforcement learning.

Engineered microbial cells present a sustainable alternative to fossil-based synthesis of chemicals and fuels. Cellular synthesis routes are readily assembled and introduced into microbial strains using state-of-the-art synthetic biology tools. However, the optimization of the strains required to reach industrially feasible production levels is far less efficient. It typically relies on trial-and-error leading into high uncertainty in total duration and cost. New techniques that can cope with the complexity and limited mechanistic knowledge of the cellular regulation are called for guiding the strain optimization. In this paper, we put forward a multi-agent reinforcement learning (MARL) approach that learns from experiments to tune the metabolic enzyme levels so that the production is improved. Our method is model-free and does not assume prior knowledge of the microbe's metabolic network or its regulation. The multi-agent approach is well-suited to make use of parallel experiments such as multi-well plates commonly used for screening microbial strains. We demonstrate the method's capabilities using the genome-scale kinetic model of Escherichia coli, k-ecoli457, as a surrogate for an in vivo cell behaviour in cultivation experiments. We investigate the method's performance relevant for practical applicability in strain engineering i.e. the speed of convergence towards the optimum response, noise tolerance, and the statistical stability of the solutions found. We further evaluate the proposed MARL approach in improving L-tryptophan production by yeast Saccharomyces cerevisiae, using publicly available experimental data on the performance of a combinatorial strain library. Overall, our results show that multi-agent reinforcement learning is a promising approach for guiding the strain optimization beyond mechanistic knowledge, with the goal of faster and more reliably obtaining industrially attractive production levels.

This is a PLOS Computational Biology Methods paper.

Introduction

Engineered microbial cells present a sustainable alternative to fossil-based chemical and fuel production [1]. Assembling and introducing the production routes into microbial strains is enabled by the state-of-the-art synthetic biology tools. However, to improve the initial laboratory demonstrated production to industrially feasible level requires substantial optimization of the cells. This strain optimization is typically challenging since the complex and insufficiently known cellular regulation has to be overcome to divert resources to production. Predictive models have been sought for guidance of the strain optimization process. Genome-scale metabolic models have already shown useful. Several strain design approaches using genome-scale metabolic models have been derived [2-4] but they are inherently limited to accounting for mass conservation and thermodynamic constraints to metabolism. They cannot account for instance kinetic constraints, allosteric regulation and regulatory interactions beyond metabolism, all playing important roles in metabolic flux distribution in cells.

Due to the limitations in the predictive models of cellular regulation the optimization of microbial cells commonly follows the Design, Build, Test, Learn cycle (DBTL). In this cycle, the strain is designed (D), built in a laboratory (B), measured, and tested (T) to learn (L) a model on the current strain to be exploited for the next cycle design phase (D) again. While efficient engineering solutions already exist for the build and testing phases, the design and learn phases of the DBTL cycle relies still strongly on manual evaluation by domain experts. This hinders the development of new industrially relevant productions strains [5].

In recent years, artificial intelligence and machine learning have emerged as approaches to facilitate the microbial strain design by being able to account for cellular regulation beyond established mechanistic knowledge [6-12]. For example, Hamerirad et al. (2019) [6] propose a Bayesian optimization approach using Gaussian Processes to tune the gene expression levels of three promoter of the lycopene production pathway of E. coli to maximize lycopene production. In another recent study a method called ART (Automated Recommendation Tool) [7], was proposed which leverages on a Bayesian ensemble approach to learn a mapping from the omics profiles of the strains to the outputs of interest.

In this study we propose a Reinforcement Learning (RL) [13] approach for automating the design and learn phases of the DBTL cycle for strain optimization. Our method facilitates the strain optimization by analysing the response (e.g. product yield) from the previous rounds and by suggesting most promising modifications to the strain for the next round. Our method is a model-free RL approach [13-15], that is, it does not rely on prior knowledge on the structure and dynamics of the metabolic system. Instead, the method learns to predict the responses of given strain modifications and based on the predictions, recommends the most promising strain modifications to try out next. We further propose a Multi-Agent RL (MARL) approach [16], which is a good match of the parallel experimentation in the laboratory, such as using multi-well plates or multiple bioreactors running simultaneously [17-19].

In our case, we cannot assume a stationary model which allows to estimate the value function in a reliable way, therefore, in the implementation of the model-free RL our learning approach closely follows the scheme of the multi-armed bandit problem with associative tasks, also called as contextual bandit [13, 20]. In the contextual bandit framework only the current value of the reward is used with the combination of the history of the context, the earlier observed states and the modifications executed until the current observation. To support these kind of learning approaches several theoretical models [20, 21] and practical realizations [22, 23], have already been developed.

The main contributions of this paper are summarized in the following:

We present a model-free reinforcement learning method that proposes modifications to enzyme levels based on the previous experiments.

A Multi-agent extension is presented that is able to take advantage of parallel experiments, such as screening using multi-well plates or other parallelized cultivation systems.

A comprehensive empirical evaluation is presented on the genome-scale kinetic model of E. coli (k-ecoli457), evaluating the sample complexity, noise tolerance, and stability of the designs.

The proposed MARL method is assessed in improving L-tryptophan production with yeast Saccharomyces cerevisiae, using publicly available experimental data on performance of a library of yeast strains.

An open source implementation of the algorithm is provided at https://github.com/maryamsabzevari-ai/strain-optimization

Materials and methods

Reinforcement learning for optimising strain designs

Reinforcement learning (RL) refers to a goal-oriented learning approach via interacting with an environment. In RL, an agent x aims to maximize a reward r by making iterative decisions, or actions a, based on the observations of the states s of the environment. The decision are guided by a policy π that gives suitable actions in a given state.

Our overall DBTL cycle using multi-agent reinforcement learning (MARL) is shown in Fig 1. We assume that each round of the cycle corresponds to a set of cultivations run in parallel, and for each cycle a new set of strains is generated, based on the recommendations of the agents. In each round the system is observed in a pseudo steady-state (e.g. corresponding to cells in exponential growth phase), which is estimated by the observable variables. The observable variables correspond here to metabolite concentration and enzyme expression level variables, whose values are assessed at appropriate state of the cultivations. Actions correspond to genetic engineering steps that increase or decrease the levels of the metabolic enzymes that have been chosen as the controllable variables. Rewards correspond to improvement of a target variable such as the product yield or a product of the specific production and growth rates. The policy corresponds to mapping the state of the system (estimated by the observations) to a set of changes of the enzyme levels.

Fig 1

Strain design optimization loop using reinforcement learning.

Enzyme levels corresponding to the strain i are denoted as e, and y and s, correspond to the response (used in reward) and output concentrations (used as state), respectively. The action (a), corresponding to the difference of the enzyme levels in the two consecutive iterations, is given by the policy learned with MMR.

The components of the framework are summarized below:

Actions , real valued vectors containing the changes of the enzyme levels (dimension n)

States , vectors of the steady state concentrations of metabolites and enzymes (dimension n).

Rewards : the change of the target variable y (e.g. product yield) between two consequtive rounds r = y − y

Policy : mapping from the states to actions, learned from data

A heatmap illustrating the changes of the enzyme levels throughout the optimization process for a particular product (succinic acid), using the proposed MARL algorithm is presented in Fig 2.

Fig 2

Illustration of strain design optimization using the MARL approach for succinic acid.

The first 6 rows correspond to the enzyme levels and in the last row the response (product exchange flux * growth) is presented. Each column represents the enzyme levels design and the corresponding response which has been found in the iteration mentioned at the bottom.

Learning algorithm

In constructing the learning algorithm, we aim to find an optimal policy π* which maps the observed states to the promising actions. The basic task for an agent at round t is to find an action a, to likely improve the reward r, based on the observed state s of the system.

To implement this, we maintain a history , of the state-action-reward triples (s, a, r) obtained in the previous rounds, stored in matrices S, A and R, which will be used to learn the policy π for the next round.

To realize the prediction of next action, the policy function, we construct an optimization problem which aims to learn the relationship between the states and the actions.

We apply the Maximum Margin Regression (MMR) framework [24], which builds on the Support Vector Machine (SVM) and is able to predict not only binary outputs but arbitrary responses given by vectors of a Hilbert space. As a heritage from the SVM, the MMR can also use kernel represented input vectors to express potentially nonlinear relations as well. Since the MMR allows to learn vector outputs with internal structure, therefore the potential interdependence between the output components can be expressed. The MMR is the simplest known approach to implement these kind of regression problems. The alternatives, e.g. Structured SVM, see several implementations in [25], and the Input Output Kernel Regression, [26], have significantly higher computational complexity, thus their use in an iterative, Reinforcement Learning framework could be too expensive. The computational complexity of solving the optimization problem realizing an MMR based learner is the same as the complexity of the SVM applied to binary classification.

MMR learns the policy through a linear operator W: H → H, where H and H are feature spaces (Hilbert spaces) corresponding to the states and actions, respectively. Feature maps : S → H, and : A → H are used to map the states and actions to their respective feature spaces.

The predicted action by the learned policy in state s is given by where the inner product 〈(a), W (s)〉 can be interpreted as a predicted reward of action a in state s. When the feature map ψ is a surjection, the above maximization has the alternative solution π(s) = ψ−1(W (s)), where ψ−1: H → A is the pre-image of ψ.

The learning is guided by a constraint calling for all triplets (s, a, r) in training data the predicted reward to be lower bounded by a margin, a strictly monotonically increasing function of the observed reward:

The type of constraint encourages the model to use the triplets with a high reward as support vectors and thus improves the model’s performance in the high-reward regime.

The full optimization problem is formulated in the following way:

The objective minimizes a sum of a regularizer (squared Frobenius norm of W) and an error term consisting of sum of slack variables for the training data. The slacks allow outliers, triplets that fail to achieve a required margin. For the margin, we applied log-transformed rewards: M(r) = log2(r + 1). The hyper-parameter C controls the trade-off between regularization and error.

Like with SVM, efficient non-linear modelling can be achieved by the use of non-linear kernel functions, which circumvent the need for constructing the potentially high-dimensional feature maps explicitly. We apply Gaussian kernels and for the state and action spaces, as well as a joint kernel K(s, a; s′, a′) = K(s, s′)K(a, a′).

Using the above kernels, a Lagrangian dual problem of (1) can be derived, where the predicted action in state s is expressed as as an approximate pre-image solution of Gaussian prediction over actions.

The solution to the above problem is computed by solving the the dual optimization problem [24]: where K = K ∘ K is the element-wise product of the kernels matrices of the states and the actions.

Multi-agent algorithm

The learning algorithm described above learns sequentially, from one experiment per DBTL round. However, modern wet-lab arrangements typically allow parallel execution of experiments either by using multi-well plates or micro-bioreactors. In general, the parallel execution of experiments have a lower overhead, due to the lower human labour cost. However, from machine learning point of view, fully parallel execution does not yield optimal results as most learning happens between the DBTL rounds. MARL employs several co-operating agents that learn from each other and together explore the solution space.

MARL setups can be generally defined using centralized or decentralized schemes in either/both of training and execution phases [16]. A centralized training scheme enables the agents to learn from each other, which can improve sample efficiency [27, 28], critical in the applications where environment interaction (e.g. Build and Test phases in DBTL) is expensive. However, decentralized learning allows the agents to adapt to specific regimes of the state space, and thus, obtains better overall performance through diversity among the agents [29]. Various structures for mixing centralized and decentralized schemes have been studied, due to their diverse advantages which can be helpful in complex real-world tasks [30-32].

In the following, we show how the strain optimization task can be tackled with MARL approach, where several agents provide strain recommendations to be executed in parallel within a single DBTL round.

Our MARL algorithm (Algorithm 1) uses a mixture of centralized and decentralized training, realized by grouping the agents so that within group g the agents are learning a joint policy π based on the history H of the group (centralized part). Across groups the agents learn independently, which allows each agent group to specialize to a specific region of the state space (decentralized part). All agents share the same actions space A and state space S. However, each agents occupies its own state at any given iteration.

The algorithm begins with an initialization phase where an initial data set is generated as a warm-up to the MARL algorithm and initial enzyme levels are chosen (see section below for details).

In each iteration of the MARL algorithm t, an agent x in group g observes the state , and executes an action based on the group policy π, which transitions the system to the next state , and the reward is obtained.

The parameters of the joint policy π for each group g, are updated in each iteration using the Eq (3), using the groups history as the training data, where the matrices of states, actions and rewards are built by concatenating the individual matrices of the agents in the group.

To avoid the the agents of each group to become too homogeneous, the actions are perturbed to increase their diversity before their execution (see below for details). In addition, every few iterations, the worst performing agent is replaced with a copy of randomly selected other agent.

Initialization

Since the algorithm uses the experiments that have been collected previously as training examples, an empty memory in the beginning produces arbitrary actions that can be avoided by including the warming up phase to initiate the memory. The warm-up period and controllable factors initialization is implemented as described below:

Warm-up: In this study, as in [7], we have used Latin hyper-cube sampling (LHS) [33] to ensure having diverse initial experiments representing the input space variability. In LHS, each dimension’s range is divided into non-overlapping, equally probable intervals. Hence the input space is divided into equiprobable subspaces. Each sample should be taken from a hyper-row and hyper-column, which does not contain any of the previous samples. We build up a history of the data obtained from τ performed interactions with the environment for each agent x, which will be later used to predict the future actions.

Choosing initial enzyme levels: For the first round (t = τ + 1), the enzyme levels of the first agent are initiated using the median of the the enzyme levels corresponding to the upper quantile of the observed responses in the warm-up. Enzyme levels of the other agents are initialized randomly to ensure exploration of different sub-spaces.

Algorithm 1 Multi agent Max-Margin Strain Design Optimization Algorithm

Input: τ: warm-up iterations, X: agents, : agent groups, S: state space, A: action space

Output: e: enzyme levels

t = τ + 1

Repeat:

For j = 1…n: # loop over the groups of agents

∀x ∈ g # each agent makes an experiment

∀x ∈ g # update group history

π = MMR() # update group policy using MMR (3)

∀ x ∈ g # predict actions:

= Perturb_agents # diversify the actions of the group

, # compute new enzyme levels for each agent

If t mod k = 0 # every k iterations

# worst agent

# substitute with a random agent

t = t+1

Until convergence

Perturbation scheme

One of the main purpose of applying multi-agents in the learning procedure is to increase the explorative capacity of the learner. This requires an algorithm in which the divergence of the actions among the agents is sufficiently large. The basic idea is to increase the volume of the parallelotope spanned by the action vectors. That volume can be expressed by the determinant of the Gram matrix, consisting of the inner product between the action vectors, and the greater determinant implies greater volume. The maximum volume corresponds to case where those actions are orthogonal, and the minimum is achieved when they point into the same direction. In the presented algorithm the perturbation vectors for each group (g) of the agents are created by the Gram-Schmidt orthogonalization of the actions, then those perturbations are scaled by a fixed perturbation coefficient (c = 0.8) and added to the original action vectors. Finally the length of the actions are restored to the original one. It can be proved that the volume corresponding to the Gram matrix of the actions increases by this algorithm. A summary of this procedure in presented in algorithm 2.

Algorithm 2 Perturb_agents

Input: actions: a1, …, a, perturbation coefficient c

Output: Perturbed actions:

For i = 1…sz:

#action perturbation

#Normalizing to the original length

Enzyme level tuning in k-ecoli457

In this work, we have used simulations from k-ecoli457 [34] genome-scale kinetic model of E. coli metabolism as a surrogate for wet-lab experimentation. The k-ecoli457 model has been validated on fluxomic data corresponding to wild-type and 25 mutant strains (21 genetically perturbed strains with glucose as the carbon substrate and four with different carbon substrates under aerobic conditions) growing under different conditions on various substrates. As the authors of this model reported in their study, the Pearson correlation coefficient between the experimental data and the product yield predicted using this model for 320 engineered strains is notably higher than using flux balance analysis [2] or minimization of metabolic adjustment [35] for predicting the product yield. This comprehensive kinetic model with 457 reactions, 337 metabolites, and 295 substrate-level regulatory interactions seems to be an appropriate candidate for strain optimization simulation.

In k-ecoli457 model, change of the enzyme levels is implemented by varying the total pool of the normalized enzymes from wild-type () down to 10-fold down-regulation () and up to 10-fold up-regulation ().

For the purpose of the MARL algorithm, enzyme levels are transformed via , which induces nearly symmetric ranges for of (−1, 0), (1, 0.82], −1 and 0 for down-regulation, up-regulation, knock-out and no regulation, respectively. Inverse of the transformation function f (), is used to transform back the enzyme level for interacting with the simulator (environment).

Choice of controllable variables and products

In Table 1, the list of the perturbed enzymes for three investigated products acetate, ethanol and succinic acid is presented, chosen according to [34] and [36]. The strain designs for improved production are sought for using a target quantity [production x growth] similar to [37], and referred to as response from here on, to avoid solutions severely compromising the viability of cells. The response for each product on the carbon source glucose is computed as rr * rr, where the first term denotes the specific production rate of the target compound and the second term denotes the growth. All the simulations were performed under aerobic conditions with glucose as a non-depleting carbon source (i.e. external concentration independent of utilization) [34]. Enzyme level perturbation has been performed by varying (increasing or decreasing) the enzyme levels within the range [0.1, 10], in which [0.1, 1] corresponds to 10-fold down-regulation up to wild-type level range and [1, 10] corresponds to wild-type level up to 10-fold up-regulation (more information in [34]). We have discarded gene deletion due to causing abrupt (non-continuous) changes in the yield space. The modified enzyme levels set in each iteration of the algorithm is queried in the simulator, to check for the obtained response improvement over the best found response yet. In addition of the response improvement, we are interested to evaluate the stability of the obtained response in the response region. This is mainly crucial, as in the real-world scenario the error-prone measurements might land the experiment in a close neighborhood of the computed enzyme levels rather than in a exact point of in vivo optimality in the enzyme level space.

Table 1

List of the investigated products of interest and the corresponding enzymes whose levels were optimized, selected based on [34] and [36].

product	target reaction	perturbed enzymes
succinic acid [36]	SUCOS	ICL, PPC, ENO, FUM, ICDHyr, SUCCt2_2pp_ex_H
acetate [34]	EX_ac(e)	RPI, PFL, GLCptspp_ex_exi, HSK, THRS, IPPS, IPMD, OMCDC, IPPMIb, IPPMIa, HSDy, ASPK, LDH_D
ethanol [34]	PDH	PFL, FRD2, RPI, GLCptspp_ex_exi, HSK, LDH_D, ASPK, THRS, IPPS, IPMD, OMCDC, IPPMIb, IPPMIa

In this study, we assume running 20 parallel experiments (4 technical replicates times 5 biological replicates) and 40 sequential iterations search performed in the enzyme level space.

Compared strain optimization methods

To asses the proposed MARL method for strain design optimization, we carried out a comparison against Bayesian Optimization with Gaussian process as surrogate function (BO-GP) [38, 39] and also random search (RAND), which has been shown theoretically and empirically to be more efficient than grid search for hyper-parameter optimization [40].

For the all three algorithms, we have used the same examples for the warm up period. In the experiments, four parallel agents are considered by the MARL algorithm, and to match the setup, in the BO-GP and RAND experiments, four parallel experiments are executed. We have used an existing Bayesian Optimization implementation in SHERPA [41], for BO-GP.

All the reported experiments relate to the average of the best found response obtained by iteration 40—a point where the improvement of the response variable has slowed down for all methods—across four agents and five runs of varied seeds (resembling biological replicates).

Results

We report results on evaluating of the proposed MARL method both in simulation data arising from the k-ecoli457 model, as well as on experimental data arising from [8].

In the simulation experiments, the MARL approach is used to optimize the k-ecoli457 model’s production of acetate, ethanol or succinic acid multiplied by growth (i.e. response). All the reported responses are relative to their corresponding wild-type cases, that is, improvement obtained for each product with respect to its corresponding wild-type case. Besides the response optimization, we also study the stability of the solutions and the noise tolerance of algorithms.

For the evaluation of the proposed method on experimental data, we used a combinatorial yeast S. cerevisiae strain library provided in [8]. In this problem, the aim was improving L-tryptophan production by optimizing the expression levels of five genes PCK1, TAL1, TKL1, CDC19 and PFK1.

Response improvement

The response improvements obtained for the aforementioned products in k-ecoli457, using the studied algorithms, MARL, BO-GP and RAND, are presented in this section.

In Fig 3, we have illustrated the median (solid lines) and 25th to 75th percentiles (shaded areas) of the obtained best found response improvement of the MARL approach with respect to BO-GP and RAND throughout the iterations, for the three products, acetate, ethanol and succinic acid. MARL approach reached in earlier iterations than BO-GP and RAND improvement in both acetate and succinic acid, and it outperformed the two other approaches in ethanol as early as iteration 13.

Fig 3

The development of the median (solid lines) and 25th to the 75th percentile (shaded areas) of the response (growth*production, C-mmol biomass/(g CDW h)*mmol/(g CDW h)) in k-ecoli457 model using MARL, BO-GP and RAND, in acetate (top), ethanol (middle) and succinic acid (bottom) production.

Statistical stability of the optimal strain designs

Here we set out to to study how sensitive the final strain design obtained by the algorithm is to the situation that the actual enzyme levels realized are different from the designed ones. For quantifying the stability of a particular strain design, we computed the variation in the response (product exchange flux * growth) of the final strain design within a neighborhood generated by small perturbations of strains. Specifically, we sampled randomly within 0.5 fold change of enzyme levels of the final strain design and checked the target response in each sampled point. We used as a measure of strain stability the relative standard deviation (or coefficient of variation) of the yields, where σ and μ denote standard deviation and mean of the neighbors’ responses, respectively.

We ran MARL for 40 iterations, and considered the strain design with the best response within those iterations as the final strain design. We computed in total 20 trajectories per product, consisting of four agents, and five repeated trials per agent. We used neighborhood size of ten for computing the stability measure (RSD). RL appears to converge to more stable strain designs than RAND meaning that the close neighbourhood of the recommended strains by MARL are less likely to result in a drastic different response (Table 2). Compared to BO-GP, more stable strain designs are found by MARL for acetate and ethanol, and for succinic acid the stability of the strain designs obtained with BO-GP and MARL are comparable. However, the final succinic acid response with BO-GP is significantly lower than the response obtained with either MARL or BO-GP, which lessens the worth of the good stability.

Table 2

Mean and standard deviations of the strain stability measure (RSD) obtained with algorithms.

For each product and each method mean and standard deviation among 20 computed RSDs are shown.

product	method	mean±std
acetate	RL	30.49±16.21
	BO-GP	45.86±24.44
	RAND	49.67±14.37
ethanol	RL	23.59±21.18
	BO-GP	71.33±55.58
	RAND	51.36±28.38
succinic acid	RL	34.11±25.31
	BO-GP	33.47±22.53
	RAND	57.62±31.87

Production improvement in the presence of noise

Using the strain designs in practice is subject to noise arising from different sources (measurement errors, biological variability, inaccuracy of implementing enzyme level changes). Since the k-ecoli457 model is deterministic, it lacks all of these sources of variability and, consequently, might make the algorithm’s performance look unrealistically good compared to the real-life situation. To evaluate the robustness of the proposed MARL strategy under noise, we performed simulations related to MARL approach in a setup where noise was added to the input or output variables of the k-ecoli457 model. We experimented with noise levels of m% (m ∈ [10, 20, 30]) noise in either states, actions or response (product exchange flux * growth).

To simulate noisy states (concentrations of metabolites and enzymes in k-ecoli457 output), the original state vector was added to a Gaussian noise signal denoted as X ∼ N(μ, σ2) which is normally distributed with mean μ = 0 and standard deviation equal to m% of the median of the values in the state vector. Similarly, the action vector (changes on the enzyme levels input to k-ecoli457) is summed up with a Gaussian noise signal with mean μ = 0 and standard deviation equal to m% of the minimum absolute value in the action vector. The obtained growth-coupled production (response) is also perturbed by a Gaussian noise signal with mean μ = 0 and standard deviation equal to m% of the original growth-coupled production. The average percentage of median max-response decrements at iteration 40 over three studied products was computed for the three investigated noise levels for the case when noise is added to the action (blue bars), response (orange bars), states (green bars) and also to the all three elements (pink bars), with respect to no noise scenario. In none of the simulations a significant adverse effect was observed, which reveals the robustness of the proposed MARL method in the presence of noisy measurements (Fig 4).

Fig 4

Average percentage of median response decrements at iteration 40 over 3 studied products, with respect to the baseline with no noise.

Evaluation of strain design optimization using experimental data

In this section we report the evaluation of our proposed MARL approach on experimental data from a combinatorial yeast S. cerevisiae strain library provided by Zhang et al. (2020) [8] (Supplementary Data 3). The combinatorial strain library was aimed at improving L-tryptophan production by optimizing the expression levels of five genes PCK1, TAL1, TKL1, CDC19 and PFK1. As a proxy of L-tryptophan production a green-fluorescent protein (GFP) synthesis rate readout as a reporter of an engineered tryptophan biosensor had been determined. This comprehensive data set allowed us to evaluate designs proposed by the MARL algorithm, by selecting each time the closest matching design from the library.

To this end, the action vector (a) consists of 6 levels of genes strength coded as {0, 0.2, 0.4, 0.6, 0.8, 1}, respectively correspond to the lowest strength up to the highest of each gene. The state vector (s) consists of all the information available for each strain (the entire row for each strain of Supplementary Data 3 in [8]). The reward (r) is defined as the change of the GFP synthesis rate, which is the optimization objective in this problem.

The setting of the algorithm is kept as previously mentioned. At each iteration, for each recommendation given by the algorithm, we search for its nearest neighbor strain (element-wise subtraction for each coded genes strengths is less than 0.2) in the library to use it as the strain that has to be queried. In the case that no close strain is found, zeros is assigned to the GFP rate of that experiment. This assumption is not destructive in our experiments as we are looking at the best found GFP synthesis rate up to the current iteration, hence this choice will only miss the strains that can potentially produce higher GFP synthesis rates but are not provided in the available data set. We use this approach as a proxy for the lab experimentation.

We compare the proposed MARL algorithm with RAND and BO-GP as explained in the previous experiments. In Fig 5 we illustrate the average and 25th to 75th percentiles (across four agents and 20 times running the algorithms) of the best found GFP synthesis rates using MARL, BO-GP an RAND. Note that, for this experiment, we have not used the recommendation strains for training. However, to asses whether the algorithm finds the strains proposed by [8], we compare each recommendation with all the strains to monitor their corresponding GFP synthesis rates. As we can see that the proposed MARL approach outperforms BO-GP and RAND in terms of the average obtained GFP synthesis rate (MFI/h). Also we observe in this plot, MARL tends to find the highest GFP synthesis rates (blue shaded area).

Fig 5

Mean and 25th to 75th percentile (shaded areas) of the best GFP synthesis rate (MFI/h), using MARL, BO-GP and RAND.

Discussion

In this paper we have put forward a MARL framework for strain design optimization. The proposed MARL algorithm is able to find metabolic enzyme levels that would optimize a given target such as desired compound production. MARL does not assume any stoichiometry or metabolic kinetics to be known but learns those dependencies and beyond between the given target and the enzyme levels from strain characterization data. Optimizing the enzyme levels has been found useful in improving production [42] and in wet-lab the enzyme levels can be adjusted by using e.g. different promoters, numbers of gene copies, or synthetic degrons. Though wide promoter libraries and other options for tuning the enzyme levels are available for many hosts, desired levels may not be achieved exactly. For this reason it is important that the designs are stable around the optimal enzyme levels. In our simulations using the k-ecoli457 model as a proxy for the wet-lab, our MARL approach in general showed better stability of the enzyme level designs upon small perturbations than the alternative BO-GP algorithm. The actual enzyme levels in engineered strains can be quantified by sampling cultivations and using targeted mass spectrometry based proteomics approaches. Neither noise in the enzyme levels (i.e. input data) nor in the target (i.e. output data) substantially compromised the performance of our MARL algorithm, instead a gradual (linear) deterioration of obtained response levels was observed upon increasing noise levels.

We demonstrated the versatility of the algorithm by optimizing L-tryptophan production in S. cerevisiae using the combinatorial library data by Zhang et al. (2020) [8]. Our results show that MARL could optimize the L-tryptophan production (response variable being the GFP synthesis rate in response to L-tryptophan biosensor) effectively, achieving 95% of the optimum mean response within 12 iterations, while BO-GP algorithm not achieving this level within 19 iterations. Thus, MARL would be suitable for optimizing enzyme levels when integrated to an experimental DBTL cycle.

For the practical use of MARL integrated to an experimental DBTL cycle, we note that when the experiments were replaced with the k-ecoli457 model simulations, around 10-15 iterations of the MARL algorithm (i.e. corresponding to wetlab strain improvement and characterization iterations) were required to reach a median response that is within 10% of the optimum median response in case of improving ethanol or succinic acid production, while for improving acetate production more than thirty iteration were needed. This number depends, in part, on the number of enzymes that are being simultaneously optimized which here was between 6-13 enzymes. Focusing on a smaller set of enzymes is likely to lead to lower number of required samples to reach convergence, but at the same time, usually a lower optimum response will be exhibited. Another direction is to increase the parallelism by using more agents per iteration (i.e. corresponding to parallel modified strains in wetlab). This is also likely to decrease the number of iterations required, however, the cost per iteration will be higher, and the overall number of strain designs required will be higher.

MARL could also be used for optimizing metabolic gene expression levels. In that case, we expect that the amount of samples required by the algorithm is higher than in the enzyme level optimization, due to the more indirect link to metabolic fluxes. Similarly, response variables alternative to the product exchange flux by growth we restricted our attention to in the simulation studies or the biosensor’s GFP reporter readout in the L-tryptophan case (i.e. Zhang et al. (2020) [8] are directly compatible with MARL.

Current limitation of our MARL method is that the strain design optimization is restricted to the given target enzymes though the number of target enzymes can be varied. However, target enzymes whose levels should be tuned for optimizing production have previously been successfully identified using a genome-scale metabolic model simulations [8]. If suitable prior knowledge is not available, we believe using genome-scale metabolic models for the target enzyme identification for MARL is a strong approach as it relies on the ultimate limits of metabolic states, mass conservation and thermodynamics. Genome-scale metabolic model simulations could also be used for proposing enzyme deletions leading to overproduction previously successfully predicted [43-45] before strain design optimization using MARL. However, the gene deletions are limited to the very few pathways that compete with production but are not essential for sufficient cell growth. Further development of MARL could also be using fluxes of modified strains predicted with genome-scale metabolic model simulations (with established methods such as e.g. MOMA [35], ROOM [46]) as an input for MARL together with experimental strain characterization data. Finally, optimizing the actual enzyme levels for maximizing production calls for a data driven approach like our MARL method due to the complexity of cellular regulation beyond known stoichiometric dependencies and kinetic constraints.

Overall, our results show that MARL is a promising approach for guiding the strain optimization beyond mechanistic knowledge, and has potential to contribute to advancing novel microbial strains reaching industrial production processes.

25 Aug 2021

Dear Dr. sabzevari,

Thank you very much for submitting your manuscript "Strain design optimization using reinforcement Learning" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

The originality of your work was well received all around, but a major issue has been recognized relating to the practical application of your approach: you don't provide biological validation of the method and it seems that it could be very hard to implement in practice. In addition several reviewers commented on the lack of discussion of your results in the biological context. In order for the paper to be accepted in a future version it is critical that you address these points clearly. Of course we would like to see your reply (and corresponding changes) relating to the other points raised by the reviewers too. However the experimental feasibility of the method must be addressed, otherwise the manuscript fails to fit in the publication scope of this journal.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Pedro Mendes, PhD

Associate Editor

PLOS Computational Biology

Daniel Beard

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors recognise that strain optimization is a combinatorial search problem, and that RL is a useful approach in such cases (provided the data – real or simulated - are available). Of course in the hands of DeepMnd RL has had considerable success e.g. in learning to play Go, and so on. They (wisely) choose a model-free approach, not least since there are no directly predictive relations between the observables (i) concentration and (ii) flux (both are variables).

I set out thinking that this would be very exciting, but TBH it promised rather more than it delivered. I think the gushing tone of the abstract needs to be toned down to recognise the limitations more clearly: (i) only three somewhat unexciting fluxes were measured, (ii) it required FORTY iterations, which is a huge number for an experimental program, and (iii) many of the explanations were close to non-existent. It probably reflects the mathematical strengths of the first author, but these cannot come – in PlosBiol at least – at the expense of explaining what the findings actually meant. In consequence, it failed to convince me that I might seek to adopt their strategy, which is what I would wish it to have done. (Or alternatively that it was a nice idea but not in fact worth adopting in this way, a point nicely made in Harford T: How to

Make the World Add Up: ten rules for thinking differently about numbers. London: Bridge Stree Press, 2020.)

Abstract: “We demonstrate the method’s capabilities in comprehensive experiments using the genome-scale kinetic model of E. coli, k-ecoli457, as a surrogate for in vivo cell behaviour”. These were simulations, not experiments. Call them so. Comprehensive is a claim, not a fact.

I am not qualified other than to take the maths on pp 3-4 on trust. The background of the authors suggests this is reasonable.

The multi-agent aspect is nice as it recognises that this is effectively how experiments are performed in a DBTL cycle.

The authors might care to contrast their approach with BOED as in

Foster A, Jankowiak M, Bingham E, Horsfall P, Teh YW, Rainforth T, Goodman N: Variational Bayesian Optimal Experimental Design. arXiv 2019:1903.05480v05483.

Vanlier J, Tiemann CA, Hilbers PAJ, van Riel NAW: A Bayesian approach to targeted experiment design. Bioinformatics 2012; 28:1136-1142.

Fig 2. Numbers in heat map unreadable. Also try to choose non red-green for those who are colour blind. https://colorbrewer2.org/#type=diverging&scheme=BrBG&n=3 gives suggestions.

“As we observe in this figure, the algoritmh [sic] tunes the enzyme levels in generally smooth fashion while the larger changes correspond to replacing the worst agent with another”. Maybe so, but I cannot see this and so it needs a much better explanation. Tell me what I am supposed to look at and what I am supposed to draw from it.

I don’t understand the units in Figs 2 and 3 – you can’t have more moles out than went in…

Yield improvement p8. 40 iteration is not that realistic for real DBTL programs. It would be good to show the time course of the median per iteration.

Fig 3. Why is succinate different from ethanol and acetate?

Table 2. Strain stability measure. Needs far more explanation – what am I supposed to infer? Also I doubt that the second DP has much meaning, or even the first…

Fig 4. Again you tell us what you did but not what we are supposed to make of the results.

“The enzyme levels can be quantified using targeted mass spectrometry based proteomics approaches found useful in optimizing production”. I doubt it. Transcriptomics typically provide a much better surrogate (Machado D, Herrgård M: Systematic evaluation of methods for integration of transcriptomic data into constraint-based models of metabolism. PLoS Comput Biol 2014; 10:e1003580).

Reviewer #2: This is a good paper, the proposed reinforcement learning (RL) method is novel to the field of strain engineering, well explained and sound. The benchmark with other methods is also fair and clear. The introduction of noise is fully relevant especially where dealing with strain engineering. Yet, I find the method a bit limited in scope. First, one needs a strain ODE kinetics models and very few are available. Secondly, the method is applied only to a predefined set of enzymes and that information is not always available. Finally, to experimentally implement the method one needs to regulate the expression levels of enzymes (via promoters, RBSs..) and this is not a straightforward task where model simulations always agree with experimental results.

I would advise the authors to add in the discussion the current limitations of their method and discus if and how a similar method could be developed and used with genome scale models (GEMs) and steady state dynamics as there are plenty of such models, methods and results in the literature. The authors could also comment and if and how the actions could be simplified to propose set of genes to be knocked out. Note that in that latter case and when using GEMs, results could be compared with MILP solvers like OptKnock.

Minor comments:

Page 2: In addition to references 6-8 which make use of machine learning but not RL there are few papers explicitly making use of RL in the context of bioproduction and it is worth mentioning these //doi.org/10.1016/j.jprocont.2018.07.013 , //doi.org/10.1021/acssynbio.9b00447, //doi.org/10.1016/j.compchemeng.2019.106649 , //doi.org/10.1371/journal.pcbi.1007783 .

Figure 1. It would be wise to define the symbols used in the caption. Some can be found in the text but not all. Would also be wise to indicate where in the flow chart the policy is learned and used.

Page 4. The ML engine to build the policy is MMR (SVM based), could the author motivate this choice as other methods could have been used?

Page 5. The authors make use of a mixed centralized and decentralized training and defined groups where RL is carried out. Although within a group, exploration seems to be favored, I do not see the RL exploitation vs. exploration strategy being used nor discussed. Could the authors comment on that?

Table 1. The target reaction needs to be defined.

Figure. 2. The process is iterated 40 times but the Figure shows only 13 rows, what are they? Would also be good to state that the numbers and colors correspond to enzyme level.

Figure 4. The Figure shows yield decrement when noise is introduced separately to action, yield and state, has any test been carried out where noise was introduced simultaneously on the 3 elements?

In addition there are few typos which could be cleared with a spell checker.

Reviewer #3: This is an excellent paper outlining an important new approach to strain design optimization. I have two main areas of feedback:

1) The specific settings/decisions for tuning the MARL parameters are not justified. How was it decided for example the number of iterations and trajectories? Would slightly different options yield improved performance? Would one set of parameters be universally best for any problem, or would users of this method need to re-evaluate this for each new context? This should be explained.

2) Substantial editing is needed to improve the clarity of the manuscript. For example: on line 239 "algorithm" is spelled wrong and the word "a" is missing. The sentence starting on line 254 is grammatically incorrect.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

22 Mar 2022

Submitted filename: response.pdf

Click here for additional data file.

26 Apr 2022

Dear Dr. sabzevari,

Thank you very much for submitting your manuscript "Strain design optimization using reinforcement Learning" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations.

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Pedro Mendes, PhD

Associate Editor

PLOS Computational Biology

Daniel Beard

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #2: I have read the revised version and agree with the modifications the authors have introduced. The authors have answered well all the points I raised with the exception of one.

The comment was “I would advise the authors to add in the discussion the current limitations of their method and discus if and how a similar method could be developed and used with genome scale models (GEMs) and steady state dynamics as there are plenty of such models, methods and results in the literature. The authors could also comment and if and how the actions could be simplified to propose set of genes to be knocked out. Note that in that latter case and when using GEMs, results could be compared with MILP solvers like OptKnock.“

As a reply to this comment the authors are presenting the difficulties and solutions to modulate enzyme expression level (a point that was raised earlier in the initial review). Could the above comment be addressed?

Reviewer #3: Thank you for your updates, this addresses my concerns and comments.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #2: No

Reviewer #3: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

References:

Review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript.

If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

5 May 2022

Submitted filename: response_r2.pdf

Click here for additional data file.

6 May 2022

Dear Dr. sabzevari,

We are pleased to inform you that your manuscript 'Strain design optimization using reinforcement Learning' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Pedro Mendes, PhD

Associate Editor

PLOS Computational Biology

Daniel Beard

Deputy Editor

PLOS Computational Biology

***********************************************************

26 May 2022

PCOMPBIOL-D-21-01222R2

Strain design optimization using reinforcement learning

Dear Dr Sabzevari,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Anita Estes

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol