The cell free protein synthesis (CFPS) reactions were carried out using the myTXTL Sigma 70 Master Mix (Arbor Biosciences) in 1.5 mL Eppendorf tubes. The working volume of all the reactions was 12 μL, composed of the Sigma 70 Master Mix (9 μL) and the plasmids (3 μL total): P70a-deGFP (5 nM) for the single-gene system; P70a-deGFP-ssrA (8 nM), P70a-S28 (1.5 nM), and P28a-cI-ssrA (1 nM) for the negative feedback circuit. The plasmids were bought in lyophilized form (Arbor Biosciences) and purified using QIAprep Spin Miniprep Kit (Qiagen) using cell lines DH5-Alpha (for P28a-cI-ssrA) or KL740 (for P70a-deGFP, P70a-deGFP-ssrA, and P70a-S28). The CFPS reactions were incubated at 29°C.

Following each CFPS run, the total RNA was extracted from 1 μL of the reaction mixture using PureLink RNA Mini Kit (Thermo Fisher Scientific) and stored at −80°C. The quantitative RT-PCR reactions were done using Applied Biosystems™ TaqMan™ RNA-to-CT™ 1-Step Kit and Custom TaqMan Gene Expression Assays (Thermo Fisher Scientific). An mRNA standard curve was used to determine absolute mRNA concentrations for each of the samples. The mRNA standards were prepared as follows: separate CFPS reactions for 5 nM of plasmids (P70a-S28, P70a-deGFP, and P70a-deGFP-ssrA) were carried out for 2 h. Total RNA was extracted using the full reaction volume. This was followed by the removal of 16S and 23S rRNA using the MICROBExpress™Bacterial mRNA Enrichment Kit (Life Technologies Corporation). Lastly, the MEGAclear™Kit (Life Technologies Corporation) was used to further purify the mRNA. The mRNA concentrations were determined using the Qubit™RNA assay kit (ThermoFisher Scientific). At least three technical replicates were performed for each standard. The concentration of cI-ssrA mRNA was quantified using the deGFP-ssrA mRNA standard. Green fluorescent protein (deGFP) fluorescence was measured using the Varioskan Lux plate reader at 488 nm (excitation) and 535 nm (emission). At the end of the CFPS run, 3 μL of the reaction mixture was diluted in 33 μL phosphate buffered saline (PBS) and stored at −80°C. The fluorescence was measured in triplicate with 10 μL each of this mixture. For all measurements, at least three biological replicates were performed.

The two genetic circuits (C1 and C2) used in this study were based upon the bacterial sigma factor regulatory system (Figure 1). Sigma factor 70 (σ₇₀), endogenously present in the extract, was the primary driver of each circuit. In C1, σ₇₀ induced green fluorescent protein (deGFP) expression was explored in the absence of additional regulators or protein degradation (Figure 1A). In C2, σ₇₀ induced the expression of sigma factor 28 (σ₂₈) and deGFP-ssrA (Figure 1B). Sigma 28 induced the expression of the lambda phage repressor protein cI-ssrA, which was under the σ₂₈ responsive P28 promoter. The cI-ssrA protein repressed the P70a promoter, thereby down-regulating σ₂₈ and deGFP-ssrA transcription (Marshall and Noireaux, 2018). Simultaneously, the C-terminal ssrA degradation tags present on the deGFP and cI proteins were recognized by the endogenous ClpXP protease system in the cell free extract, thereby promoting the degradation of these proteins into peptide fragments (Flynn et al., 2003; Garamella et al., 2016). In addition, messenger RNAs (mRNAs) were always subject to degradation due to the presence of degradation enzymes in the extract (Karzbrun et al., 2011; Garamella et al., 2016). Taken together, the interactions of the components manifested in an accumulation of deGFP protein for C1, and a pulse signal of deGFP-ssrA in C2. Studying C1 allowed us to estimate parameters governing the interaction of σ₇₀ with the P70a promoter. Whereas, the C2 allowed us to characterize the interaction of σ₂₈ with the P28 promoter, the strength of the transcriptional repression by cI-ssrA, and the kinetics of protein degradation by the endogenous ClpXP protease system. Finally, both circuits tested the effective model formulation for the transcription and translation rates.

Consider a cell free synthetic circuit composed of the genes G=1,2,…,N. Each gene in the circuit is described by two differential equations, one for mRNA (m_j) and a second for the corresponding protein (p_j): (1)ṁj=rX,juj(…)-θm,jmj j=1,2,…,N (2)ṗj=rL,jwj(…)-θp,jpj The term r_X,ju_j(…) in the mRNA balance, which denotes the regulated rate of transcription for gene j, is the product of a kinetic limit r_X,j (nM h⁻¹) and a transcription control function 0 ≤ u_j(…) ≤ 1 (dimensionless). Similarly, the rate of translation of mRNA j, denoted by r_L,jw_j, is also the product of the kinetic limit of translation (nM h⁻¹) and a translational control term 0 ≤ w_j(…) ≤ 1 (dimensionless). Lastly, θ_{⋆, j} denotes the first-order rate constant (h⁻¹) governing degradation of protein and mRNA. The model equations, encoded in the Julia programming language (Bezanson et al., 2017), were automatically generated using the JuGRN tool². The model equations were solved numerically using the Rosenbrock23 routine of the DifferentialEquations.jl package (Rackauckas and Nie, 2017).

Model parameters were estimated from published studies, were specified by experimental conditions (Table 1) or were estimated by minimizing the squared difference between model simulations and messenger RNA (mRNA), or protein measurements generated in this study. For the P70-deGFP model (C1), 11 parameters were estimated, while 33 parameters were estimated for the negative feedback model (C2).

The minimization problem to estimate the unknown model parameters was structured as a multiobjective optimization problem in which each measured mRNA or protein trajectory was treated as a separate objective. The minimization problem was solved using the Pareto Optimal Ensemble Technique in the Julia programming language (JuPOETs) (Bassen et al., 2017). JuPOETs is a multiobjective optimization approach which integrates simulated annealing with Pareto optimality to estimate parameter values on or near the optimal tradeoff surface between N potentially competing objectives (squared difference between model simulations and experimental measurements). JuPOETs minimized a problem of the form: (11)minkEj=∑i=1Tj(M^ij-xij(k))2  j=1,2,…,N subject to (12)ẋ=f(x,k) (13)L≤k≤U (14)x(to)=xo where Equation (12) denotes the model equations, Equation (13) denotes the parameter bounds, and Equation (14) denotes the initial conditions. The objective function(s) Ej measured the squared difference between model simulations and experiment j (either a protein or mRNA trajectory). The symbol M^ij denotes an experimental observation at time index i from experiment j, while the symbol x_ij denotes the model simulation output at time index i from experiment j. The quantity i denotes the sampled time-index and Tj denotes the number of time points for experiment j. For the P70-deGFP model (C1), E1 corresponded to mRNA deGFP, while E2 corresponded to the deGFP protein concentration. On the other hand, for the negative feedback model (C2), E1 corresponded to mRNA deGFP-ssrA, E2 to mRNA σ₂₈, E3 to mRNA cI-ssrA and E4 to the deGFP-ssrA protein concentration. Lastly, we penalized accumulation of the cI-ssrA protein (unmeasured) reaching unrealistically high levels with a term of the form: E5 = C × max(0, x_cI − U_cI) where C denotes a penalty parameter (C = 1 × 10⁵), x_cI denotes the maximum simulated cI-ssrA protein concentration, and U_cI denotes a concentration upper bound (U_cI = 100μM). This bound was chosen to be approximately five-fold higher than the protein levels observed in an uninhibited circuit (e.g., C1).

The lower and upper bounds for unknown model parameters were established from previously published studies, or from previous model analysis; parameter values estimated for the P70-deGFP model were also used to establish ranges for the negative feedback model. JuPOETs searched over ΔG_i, K_L, and τ_L,1/2 values directly, while other unknown parameter values took the form of corrections to order of magnitude characteristic literature estimates. For example, we set the mRNA degradation rate constant (θ_m) to a characteristic value taken from literature. Then, the degradation constant for any particular mRNA was represented as: θ_m,i = α_iθ_m, where α_i was an unknown (but bounded) modifier. In this way, we guaranteed the parameter search (and the resulting estimated parameters) were within a specified range of literature values. We used this procedure for all degradation constants (both mRNA and protein) and all time constants (for both transcription and translation). The baseline parameter values are given in Table 1. JuPOETs was run for 20 generations for both models, and all parameters sets with Pareto rank less than or equal to two were collected for each generation. The JuPOETs parameter estimation routine is encoded in the sa_poets_estimate.jl script in the model repositories.

JuPOETs uses a simulated annealing approach to generate candidate parameter solutions whose Pareto rank is then evaluated; ranks below a threshold are kept while higher rank solutions are discarded. Thus, all the advantages (and disadvantages) associated with simulated annealing have been inherited by JuPOETs; for example, the time required to generate a family of low rank solutions will be significantly longer than a derivative based approach. Beyond these specific performance issues, a unique pathology of JuPOETs is the use of Pareto rank as a surrogate for training error. JuPOETs attempts to find low rank solutions, but rank is a relative measure of the quality of a solution. Thus, during the early iterations, low rank solutions often have large errors. As the iteration count increases the approach tends to find low error solutions with low rank, however, for complex models the rate of convergence to these low rank low error solutions is slow. To address this concern, we periodically switch to single objective mode where we minimize the total training error (summation of all objective functions) instead of finding low rank solutions. The best solutions from single objective mode can then be used to restart the multi-objective calculation. This hybrid approach, which was used in this study, has previously been shown to increase the rate of finding low rank and low error solutions (see Bassen et al., 2017).

Morris sensitivity analysis was used to understand which model parameters were sensitive (Morris, 1991). The Morris method is a global method that computes an elementary effect value for each parameter by sampling a model performance function, in this case the area under the curve for each model species in their respective timeplots, over a range of values for each parameter; the mean of elementary effects measures the direct effect of a particular parameter on the performance function, while the variance of each elementary effect indicates whether the effects are non-linear or the result of interactions with other parameters (large variance suggests connectedness or non-linearity). The Morris sensitivity coefficients were computed using the DiffEqSensitivity.jl package (Rackauckas and Nie, 2017). The parameter ranges were established by calculating the minimum and the maximum value for each parameter in the parameter ensemble generated by JuPOETs. Each range was then subdivided into 10,000 samples for the sensitivity calculation. Elementary effect values were then calculated one at a time by measuring the perturbation in the AUC on changing one parameter, where the AUC was calculated by solving the set of ODEs for each change. In order to calculate the mean and variance, the top 1000 perturbations with the highest spread in parameter values were used. In total, the model was evaluated 10000n times, where n is the number of parameters varied. The Morris sensitivity coefficients are calculated using the compute_sensitivity_coefficients.jl script in the model repositories.

The effective biophysical transcription and translation model captured σ₇₀ induced deGFP expression at the mRNA and protein level within the experimental error for C1 (Figure 2). JuPOETs produced an ensemble (N = 140) of the 11 unknown model parameters which captured the transcription of mRNA (Figure 2A) and the translation of deGFP protein (Figure 2B). The mean and standard deviation of key parameters is given in Table 2. The deGFP mRNA reached its steady state concentration of approximately 580 nM within 2 h, and stayed at this level for the remainder of the reaction. Thus, the cell free reaction maintained continuous transcriptional activity with an average mRNA lifetime of 27 min; Garamella et al. (2016) reported a similar lifetime of 17–18 min. On the other hand, deGFP protein concentration increased more slowly, and began to saturate between 8 and 10 h at approximately 15 μM. Given there was negligible protein degradation (the mean deGFP half-life was estimated to be ~11 days, which was similar to the value of 6 days estimated by Horvath et al., albeit in a different cell free system, Horvath et al., 2020). The saturating protein concentration suggested that the translational capacity of the cell free system decreased over the course of the reaction. The decrease in translational capacity, which could stem from several sources, was captured in the simulations using a monotonically decreasing translation capacity state variable ϵ, and the translational control variable w(…). In particular, the mean half-life of translational capacity was estimated to be τ_L,1/2 ~ 4 h in the C1 experiments. Taken together, JuPOETs produced an ensemble of model parameters that captured the experimental training data. Next, we considered which C1 model parameters were important to the model performance using Morris sensitivity analysis, a global sensitivity analysis method.

The importance of C1 model parameters was quantified using Morris sensitivity analysis (Figure 2B). The Morris method computes the influence of each parameter, known as the elementary effect, on a model performance function. The mean of elementary effects measures the direct effect of a particular parameter, while the variance indicates whether the effects are non-linear or the result of interactions with other parameters (large variance suggests non-linearity). The performance function was defined as the integrated area under the curve (AUC) for each mRNA and protein species in their respective timeplots, calculated for each parameter value range. The Morris sensitivity measures (mean and variance) were binned into categories based upon their relative magnitudes, from no influence (white) to high influence (black). Only four parameters (translation saturation coefficient K_L, translational capacity half-life τ_L,1/2, translation time constant, and protein degradation constant) influenced the protein level. On the other hand, six parameters influenced both mRNA and protein abundance; all six of these parameters were either directly or indirectly associated with transcription. Thus, these parameters influenced the production or stability of mRNA which in turn influenced the protein level. In particular, the mRNA degradation constant and the cooperativity of σ₇₀ in the P70a promoter function had the largest direct effect and variance. Surprisingly, the ΔG of σ₇₀/RNAP/promoter configuration was the least influential of the six parameters and had a small elementary effect variance. Taken together, Morris sensitivity analysis of the C1 model parameters highlighted the hierarchical structure of the transcriptional and translational model, suggesting experimentally tunable parameters such as mRNA stability were globally important. Next, we used the ensemble of P70a, time constant and degradation parameters estimated for C1 to constrain the analysis of C2.

The effective biophysical transcription and translation model captured the deGFP-ssrA expression dynamics in the negative feedback circuit C2 (Figure 3A). JuPOETs produced an ensemble (N = 498) of the 33 unknown model parameters which captured transcription and translation dynamics for σ₂₈, cI-ssrA and deGFP-ssrA. The mean and standard deviation of key parameters is given in Table 3. Compared with the estimated parameters for C1, the C2 model had almost a two fold change in the half life of translation and the translation saturation coefficient. Similarly, there were variations in the values of the transcription and translation time constants for the two systems. However, for both circuits, the small values of the transcription and translation time constants qualitatively suggested elongation limited reactions; the exception was σ₂₈ translation which was closer to initiation limited. Unlike C1, the mRNA expression pattern for σ₂₈ and deGFP-ssrA both showed an initial spike, to a concentration similar with the previous pseudo steady state, before the cI-ssrA regulator protein could be expressed. However, once cI-ssrA began to accumulate, the concentrations of the regulated mRNAs dropped by approximately an order of magnitude compared with the unregulated case. Again, as shown with C1, the regulated mRNA concentrations reached an approximate steady-state. This further confirmed continuous transcription and mRNA degradation throughout the cell free reaction. The mean estimated mRNA lifetime for cI-ssrA and deGFP were similar (approximately 16 min), while the degradation of σ₂₈ mRNA was predicted to be slower (mean mRNA lifetime was estimated to be approximately 30 min). Lastly, the mean peak degradation rate for GFP was approximately 47 nM/min, while the mean peak cI-ssrA degradation rate was predicted to be approximately 63 nM/min; both of these degradation rate estimates were consistent with the previously reported range of 15–150 nM/min measured by Garamella et al. (2016).

The secondary effect of cI-ssrA repression was visible in the cI-ssrA mRNA expression pattern. The expression of cI-ssrA was induced by σ₂₈, however, σ₂₈ expression was repressed by cI-ssrA, thereby completing a negative feedback loop. Initially, before appreciable levels of cI-ssrA had been translated, the cI-ssrA transcription rate was maximum (approximately 200 nM/h). However, the transcription rate decreased to approximately 12 nM/h after 2 h and remained constant for the remainder of the cell free reaction. Similarly, transcription rates for σ₂₈ (approximately 1,200 nM/h) and deGFP-ssrA (approximately 750 nM/h) were initially at a maximum due to the presence of endogenous σ₇₀, but then quickly dropped as cI protein accumulated. Protein synthesis followed a similar trend, with the translation rates for σ₂₈ and deGFP-ssrA initially present at their maximum values before quickly dropping. After 1 h, deGFP levels reached a peak and decayed due to the ClpXP-mediated degradation, whereas σ₂₈ protein levels continued to slowly rise at a steady rate (approximately 15 nM/h). The C2 model also predicted the expected lag present during the initial phase of cI-ssrA protein synthesis due to the need for σ₂₈ protein to reach appreciable levels. Moreover, the combination of high cI-ssrA mRNA abundance (expressed because σ₂₈ does not have a degradation tag) and ClpXP-mediated degradation led to the saturation of the cI-ssrA protein concentration. However, the cI-ssrA protein concentration could not be verified because we did not have an experimental measurement for this species. Taken together, the effective model simulated cell free expression dynamics for C2. Next, we considered which C2 model parameters were important using Morris sensitivity analysis.

Morris sensitivity analysis of the negative feedback circuit C2 stratified the parameters into locally and globally important groups (Figure 3B). The influence of 33 parameters was computed using the AUC of each mRNA and protein species as the performance function. The Morris sensitivity metrics (mean and variance) were binned into categories based upon their relative magnitudes, from no influence (white) to high influence (black). Some parameters affected only their respective mRNA or protein target, whereas others had widespread effects. For example, the time constant (tc) modifiers, stability of deGFP-ssrA protein and mRNA, and the binding dissociation constant (K) and cooperativity parameter (n) of cI-ssrA and σ₇₀ for the deGFP-ssrA promoter affected only the values of deGFP-ssrA protein and mRNA. On the other hand, the tc, stability, K and n parameters for σ₇₀, σ₂₈, or cI-ssrA influenced mRNA and protein expression globally. The σ₇₀ and σ₂₈ proteins acted as inducers or repressors for more than one gene product: σ₇₀ induced both deGFP-ssrA and σ₂₈, and cI-ssrA protein repressed both of these genes. Degradation constants (denoted as stability) affected the half-lives of the transcribed messages or the translated proteins in the mixture, while the time constant modifiers changed the time required to form the open gene complex (or translationally active complex). Dissociation and cooperativity constants affected the binding interactions of the inducer (or repressor in the case of cI-ssrA) in the promoter control function. Varying these parameters, therefore, had a strong effect on their respective targets. Similarly, the translation saturation and its half-life, which captured the depletion in the translation activity over the course of the reaction, not only affected protein levels but also mRNA levels. This is because these parameters tuned the rate of formation of cI-ssrA, which in turn affected the mRNA levels of its gene targets. Given that cI-ssrA was the main regulator (repressor) of the circuit, the parameters that dictated the levels of cI-ssrA mRNA and protein had a global effect. We also observed high sensitivity variance for several parameters, in particular involving cI-ssrA. For example, the time constant modifiers for cI-ssrA mRNA and protein had a two-pronged effect. On the one hand, they positively influenced the transcription/translation rates of the gene and mRNA products, directly increasing the cI-ssrA protein. On the other hand, increased cI-ssrA expression reduced the level of σ₂₈, in turn reducing the cI-ssrA levels over time. Taken together, Morris sensitivity analysis of the C2 model stratified that parameters into local and globally important groups, with the parameters governing the synthesis rates of the cI-ssrA mRNA and protein being globally important. The sensitivity analysis also gave insight into the organization of the circuit, suggesting cI to be highly connected within the circuit.