Modular, robust, and extendible multicellular circuit design in yeast
Abstract
Division of labor between cells is ubiquitous in biology but the use of multicellular consortia for engineering applications is only beginning to be explored. A significant advantage of multicellular circuits is their potential to be modular with respect to composition but this claim has not yet been extensively tested using experiments and quantitative modeling. Here, we construct a library of 24 yeast strains capable of sending, receiving or responding to three molecular signals, characterize them experimentally and build quantitative models of their input-output relationships. We then compose these strains into two- and three-strain cascades as well as a four-strain bistable switch and show that experimentally measured consortia dynamics can be predicted from the models of the constituent parts. To further explore the achievable range of behaviors, we perform a fully automated computational search over all two-, three-, and four-strain consortia to identify combinations that realize target behaviors including logic gates, band-pass filters, and time pulses. Strain combinations that are predicted to map onto a target behavior are further computationally optimized and then experimentally tested. Experiments closely track computational predictions. The high reliability of these model descriptions further strengthens the feasibility and highlights the potential for distributed computing in synthetic biology.
Editor's evaluation
This paper used multiple strains to build gene circuits and demonstrate the modular composition of strain circuits with an automated design strategy to achieve a target behavior from a large space of possible functional circuit architectures. This paper provides synthetic biologists with an alternative solution for the problems of scalability, robustness, and modularity.
https://doi.org/10.7554/eLife.74540.sa0Introduction
The leading paradigm for genetic circuit design is to combine biological parts in a delicate balance within the same cell (Ellis et al., 2009; Kosuri et al., 2013; Ottoz et al., 2014). This approach has resulted in increasingly large genetic circuits that realize functions such as logic gates and circuits (Moon et al., 2012; Bonnet et al., 2013; Nielsen et al., 2016; Gander et al., 2017), time pulses (Gao et al., 2018; Guo and Murray, 2019), incoherent feed-forward loops (Entus et al., 2007; Ellis et al., 2009), bistable switches (Chen et al., 2012; Huang et al., 2012; Yang et al., 2019; Barbier et al., 2020; Grant et al., 2020), or oscillators (Elowitz and Leibler, 2000; Tigges et al., 2009; Tigges et al., 2010). Albeit very successful, single cell circuit engineering has limited scalability and robustness because parts cannot be reused, the genetic burden on the cell grows with circuit size (Nikolados et al., 2019), and retroactivity (Del Vecchio et al., 2008) and component crosstalk (Del Vecchio, 2015) interferes with expected behavior. Furthermore, competition over shared gene expression resources makes unrelated circuits unintentionally coupled (Zhang et al., 2021).
An alternative approach is to generate complex behaviors using consortia of cells wherein different cell types perform distinct functions and communicate with each other through chemical signals to realize more complex behaviors. Such division of labor is common in nature and interest has recently emerged in engineering multi-cellular synthetic systems (Brenner et al., 2008). Although developed later, synthetic consortia have caught up with single-cell circuits in terms of the complexity of functions that have been realized, generating behaviors such as multicellular time pulses (Basu et al., 2004), oscillations (Chen et al., 2015; Danino et al., 2010), and logic gates (Regot et al., 2011; Tamsir et al., 2011), bioproduction (Egbert et al., 2016), or circuits that use quorum sensing to define social interactions (Kong et al., 2018; Balagaddé et al., 2008; Weber et al., 2007). Furthermore, the emergence of novel orthogonal cross-species signaling molecules (Billerbeck et al., 2018; Du et al., 2020) and genetic channel-selector devices to multiplex individual signals (Sexton and Tabor, 2020) has opened the path for engineering more complex and precise multicellular behaviors.
To date, synthetic multicellular systems were largely designed with specific target behaviors in mind, rather than optimizing modularity of components to make them usable in a large number of contexts. Still, mathematical modeling has demonstrated that multi-cellular computation should easily access a larger space of behaviors: for instance, just three different cell populations can, in theory, generate up to 100 different logic functions (Regot et al., 2011) and even bimodality and cell-synchronization (Thurley et al., 2018). Furthermore, the intuitive modular nature of cell-to-cell communication should provide a useful tool to rationally design synthetic circuits with predetermined performance. Rational design significantly speeds up circuit assembly (Chen et al., 2012; Chen et al., 2020) and allows design of global functions from local behaviors (Salis et al., 2009; Carothers et al., 2011). As a matter of fact, mathematical models have been successful to bridge the gap between individual processes and collective behaviors when applied to ecological interactions between distinct populations of communicating cells (Shou et al., 2007; Momeni et al., 2013, Egbert et al., 2016). However, an example of a model-driven strain selection for multicellular circuit design is currently missing.
Here, we propose a large vocabulary of yeast strains that use chemical signals for cell-to-cell communication and that can be modularly combined to realize a large number of functions (Figure 1A). Each strain senses one or two inputs and produces a single output. The transfer function relating inputs and outputs can be either activating or repressing. The output is a fluorescent protein which can be used to read out circuit behavior, another signaling molecule which can be used to connect strains or an enzyme that sequesters or degrades a signal thus disrupting communication. We experimentally characterize each strain, model their dynamics using differential equations and then use these models to predict the behavior of strain combinations (Figure 1B). Using a library of 24 strains, we rationally design complex multicellular behaviors, such as bandpass filters, negative and positive feedbacks, time pulses, logic gates, and bistable switches. These behaviors are common in nature and they underlie cell decisions on metabolism (for instance, responding to environmental signals, [Lee et al., 2002], utilization of carbon sources, [Brandman et al., 2005]) or cell differentiation (Kueh et al., 2016; Duddu et al., 2020), and have been repeatedly connected to multicellular organization (signal gradients in development) and signaling (auxin pulses in plants). The modularity of composition achieved by our multicellular circuits allow us to design these non-trivial behaviors based on model simulations alone. This qualitative and quantitative precision demonstrated in both time and steady-state experiments further extend circuit design automation (Nielsen et al., 2016; Chen et al., 2020).
Results
Engineering yeast strains for signal sensing, synthesis and depletion
As a first step toward the construction of a vocabulary of modular yeast strains, we selected a set of signals to enable cell-to-cell communication and then optimized chassis strains for signal sensing, synthesis, and depletion. We selected the plant hormone auxin (we will also refer to it as IAA, indole-3-acetic acid) and the yeast hormone α-factor as signaling molecules and additionally used the mammalian hormone β-estradiol as an inducer.
The α-factor pathway has been subject of extensive studies and mechanisms for sensing through the surface receptor STE2, synthesis by the MFα1 gene and degradation through the BAR1 protease are well understood and have been engineered to generate a wide range of behaviors (Youk and Lim, 2014; Groves et al., 2016; Shaw et al., 2019). To have a baseline strain that does not interfere with signaling, we knocked out the native BAR1 gene as in Youk and Lim, 2014 to prevent α-factor degradation. Moreover, to account for growth-arrest induced by α-factor, which affects gene expression on a large scale, we knocked out FAR1 a protein that contributes to arresting the cell cycle at G1Chang and Chang and Herskowitz, 1990, and constitutively expressed POG1, a protein that promotes growth-arrest recovery (Leza and Elion, 1999). Unexpectedly, we detected no α-factor output from strains that both sense and secrete α-factor. We suspected that this was caused by the surface receptor STE2 internally binding to α-factor. It has been shown that, upon binding to α-factor, STE2 undergoes endocytosis and then shares the secretory pathway of α-factor itself (Schandel and Jenness, 1994). We solved this problem by overexpressing STE2 on a pGPD promoter Sun et al., 2012 in these strains, aiming to have some protein copies escaping this interaction.
Next, we turned to the optimization of components for Auxin sensing, synthesis, and elimination. In prior work from our group, we (Khakhar et al., 2016; Pierre-Jerome et al., 2014) developed an auxin-responsive transcription factor using a chimeric dCas9-Aux/IAA protein regulation. In the presence of Auxin, the Aux/IAA degron part of the protein binds to the auxin signaling F-box protein (a modified TIR1 for this study) and acts as part of an E3 ubiquitin ligase to catalyze ubiquitination and degradation of the Aux/IAA-degron-containing protein. Similarly, we constructed a synthetic auxin synthesis pathway in yeast (Khakhar et al., 2016). We demonstrated conversion of the precursor molecule IAM (indole-3-acetamide) into IAA through expression of the IaaH gene in yeast. Here, we amplified Auxin secretion and thus effectively the strength of the signal produced, through integration of the auxin-efflux pump ABCB1/PGP1 from A. thaliana into our yeast strains, as previously reported in Geisler et al., 2005. We measured a significant increase in the auxin-synthesis yield as measured by a neighboring IAA-detecting cell when PGP1 was integrated (Figure 1—figure supplement 1).
Unlike for α-factor, a depletion mechanism for Auxin had not yet been reported in yeast. To create an IAA depletion mechanism, we thus selected the plant protein GH3.3 that has been shown to conjugate IAA to aspartic acid, forming the signaling-inactive IAA-Asp. GH3.3 is part of a family of proteins that encode IAA-amino synthetases, which have been reported to control auxin homeostasis (Staswick et al., 2005). To test whether GH3.3 or related proteins could be used to inactivate Auxin in yeast, we first expressed codon-optimized versions of GH3.3 and GH3.6 from A. thaliana and C. papaya from a highly expressed pGPD promoter. We then tested the IAA to IAA-Asp conversion rate using mass spectrometry and found that GH3.3 from A. thaliana had the higher efficiency (Figure 1—figure supplement 1). Finally, we tested that IAA-Asp does not activate the IAA-mediated degradation pathway in yeast, by adding 10 µM of IAA-Asp in an auxin-sensor yeast culture and detecting no variation in fluorescence (Figure 1—figure supplement 1).
Establishing a vocabulary of parts for cell to cell communication
Having established conditions for efficient signal sensing, synthesis and inhibition, we combined signals with activating and repressive transfer functions to create a vocabulary of strains. Transfer functions that use α-factor as input are mediated by the transcription factor (STE2) that either directly induces expression of the gene of interest (activation) or induces expression of a repressor that inhibits the output (repression). Similarly, β-estradiol binds and activates a transcription factor (ZEV4, McIsaac et al., 2013) that either directly activates gene expression (activation) or induces expression of a repressor and inhibits output synthesis (repression). For both signaling molecules, we chose dCas9 fused with the repressor domain Mxi1 (Gander et al., 2017) as repressor (Figure 1—figure supplement 1 for the full pathways).
To induce activation or repression with auxin, we build on the same auxin-mediated degradation pathway used for auxin sensing above. Specifically, activation results from degradation of a repressor, while repression results from degradation of an activator. We chose a dCas9-Mxi1-auxin degron (Khakhar et al., 2016) fusion as a repressor and used a dCas9-VP64-auxin degron fusion as the activator (Figure 1—figure supplement 1).
We conducted extensive pathway optimization to increase the separation between high and low expression levels and IAA sensitivity, using a mechanistic model to explore the parameter space and guide the genetic engineering (Figure 1—figure supplement 1). We adopted the model proposed in Pierre-Jerome et al., 2014, where each parameter easily translates to a biological function, and performed parameter sensitivity analysis with respect to fold change between the baseline and the Aux/IAA-induced fully repressed state. We then selected the top six highest scoring parameter perturbations, designed circuit variants that reflected those changes, and tested them resulting in a good agreement with our predictions (Figure 1—figure supplement 1). Guided by the model, we combined five of the tested circuit variants to increase the fold change from 1.3-fold (sensor from Khakhar et al., 2016) to 3.1-fold: we used this final circuit for all the repressive strains, swapping the fluorescent reporter gene with the output gene for non-sensor strains. For the auxin activating strains, we used a similar model-driven approach to rationally design the activating pathway novel to this paper, obtaining a threefold-change activation.
Combining the three different input signals, the activation/repression circuits and output secretion, we built and tested all the possible combinatorial designs presented in Figure 1A. The strains that sense β-estradiol and repress expression of BAR1 and GH3.3, the strain that senses α-factor and represses expression of GH3.3, and the strain that senses IAA and represses BAR1 expression are shown in Figure 1—figure supplement 1 and not used below, since their response was too slow to produce meaningful results. Of the remaining 24 strains, 6 sensor strains express GFP in response to the three input signals (three activators and three repressors), 12 strains (six activators and six repressors) synthesize a signaling molecule, and 6 strains act as signal attenuators (expressing BAR1 or GH3.3). Four of these 24 strains sense and secrete the same signaling molecule (α-factor or IAA, ‘positive’ or ‘negative’ feedback strains). Finally, two strains express repressors of their own input (α-factor expressing BAR1 and IAA expressing GH3), also describing a negative feedback topology (Figure 1—figure supplement 1).
Sensor strain characterization
For sensor strain characterization, we collected time series data for eight different input concentrations. Input concentrations were selected to fully cover the sensor dynamic range for model fitting. We normalized fluorescence data by cell size and took the mean of the histogram as in Groves et al., 2016 and then subtracted background fluorescence. Each measurement in the figure is an average of three experimental repeats (error bars representing the standard deviation).
With a scalable and modular system in mind, we fitted a set of three ordinary differential equations (ODEs) with eight parameters for each strain to describe input sensing (x1), signal processing (x2) and fluorescence output synthesis (x3) (Figure 1B, Model 1). Signal processing (activation or repression) is modeled with a simple Hill function (Model 1b), which naturally incorporates signal saturation. Input sensing and output synthesis are modeled as linear ODEs (Model 1 a and c). A constant term in the last ODE accounts for background promoter activity. For instance, for the auxin sensor in Figure 1B, x represents the TIR1-auxin complex concentration, x2, the dCas9-VP64-auxin degron population, and x3, GFP concentration. These simple models capture the system dynamics, with the benefit of being easy to fit and analytically approachable. Parameters were fitted independently for each experimental repeat to obtain a mean and a standard deviation for the Hill coefficient. We also separately fitted a Hill curve to an average of the three experimental repeats and the resulting Hill coefficient () is reported in the figures and used for simulations (Figure 1B). The sensor strains range in sensitivity depending on the input. The β-estradiol sensors respond to inputs concentration ranging from 0.1 to 100 nM with an of 0.5 ± 0.0 nM for repression and 12.6 ± 0.9 nM for activation. α-factor sensors are sensitive to input concentrations ranging from 1 to 500 nM range with an nM for activation and 89.0 ± 6.6 nM for repression. Finally the IAA sensor is least sensitive and responds to inputs ranging from 5 nM-10 uM withe an nM for activation and 276.5 ± 8.8 nM for repression. The Hill coefficients for these sensors vary between 0.8 (3-repeat interval: 0.8 ± 0.0, α-factor repression) and 1.0 (1.2 ± 0.02, α-factor activation), 2.2 (2.1 ± 0.3, β-estradiol repression) and 1.3 (1.25 ± 0.0, β-estradiol activation), 1.0 (1.0 ± 0.4, auxin repression) and 0.8 (, auxin activation) consistent with previously reported values for β-estradiol repression (dCas9-Mxi1 confers stronger and more consistent repression than dCas9 alone, leading to ultrasensitivity, Gander et al., 2017) and α-factor activation (Groves et al., 2016). Finally, the sensors achieve an ON/OFF fold-change that ranges between a minimum of 3 (IAA activating GFP) and a maximum of 42 (β-estradiol activating GFP). We tested signal orthogonality by adding pairwise combinations of inducers at saturation concentration to single-input repressive sensors (Figure 1—figure supplement 1). Variations across the treatments are contained within 11% of the nominal value for signal concentrations within the range used for circuit experiments.
In addition to single-input sensors, we constructed two sensor strains that respond to both α-factor and IAA. In the first strain, α-factor induces fluorescence expression and IAA represses it, while in the second strain, IAA induces fluorescence expression and α-factor represses it (Figure 1E).
We collected data points for eight different input combinations taken at five time points up to saturation. The output signal is monotonic with respect to each input and the input functional range is similar to the one measured for the correspondent strains that respond to only one of the two inputs.
To model these strains, we used a simple extension of our previously introduced models with five ODEs (Figure 1E, Model 2): two ODEs model input sensing (a), two model input processing (b), and one ODE (c) combines the signals according to their activating or repressing nature and defines output synthesis. These models fit the experimental data even when the output is non-monotonic over time.
Assembling modular, tunable and easily-extendable circuits using cell-to-cell communication
To determine the potential to build biological circuits using cell-to-cell communication, we experimentally tested if communication occurs between strains that secrete an output and their corresponding sensor strains. For example, a strain that produces auxin in response to α-factor sensing was grown in co-culture with a sensor that switches off GFP expression in response to auxin (Figure 1B). The two strain populations were added at the same concentration and the fluorescent output was measured at steady state (10 hours after induction). The experimental data is consistent with a negative α-factor sensor as expected (the more α-factor, the lower the fluorescent signal, as seen in Figure 2A).
Since the core genetic circuit of this sender strain is identical to the α-factor repressing sensor, we tested if the mathematical model previously fit to the sensor data preserves its predictive power. We re-fit only the 3 output parameters to account for the fact that the output is now auxin rather than GFP. The output of the sender cell model was used directly as input of the sensor cell model. To test the hypothesis that strains that have common input/processing parts share parameters and model structure, we collected two datasets on two different experiments composed of four data points each: we used one for fitting (yellow dots, first experiment), and the other for validation (orange dots, second experiment, Figure 2A). A time-series at EC50 input concentration was also used for fitting (not present in the figure).
Intuitively, we expected that higher initial sender cell concentration would result in an overall higher concentration of their output signal over the same time scale. Most importantly, we wanted to know if this output ‘gain’ can be predicted by our models so that we could use it to tune circuit behavior. We modeled this effect with a factor multiplied by the output signal, where is the fold-change with respect to the standard initial concentration: we represented this gain using the same iconography used in electronic circuits (Figure 2B). Our predictions matched closely with data collected using 5 X and 10 X the original concentration of sender cells (green and blue lines) using the same strain co-culture as in the previous panel (orange line). Notice that the fold-change, as we just defined it above, is not equivalent to the ratio between the strains. Hence, difference in growth rates are not relevant for model predictions provided that there is no dilution or competition for resources, which is minimal in these experiments (the strains are kept in log-phase throughout the experiments and their concentrations are kept below saturation).
We further tested if we could modulate the concentration of signaling molecules by removing it from the system through BAR1 and GH3.3-expressing strains (Figure 2C). We modeled signal degradation as a first-order Hill repressing function, where the output of the sender cell acts as a negative regulator. As before, only the three output parameters of the sender strain were fitted using receiver cell fluorescence data. Finally, we tested all the sender-receiver pairs in our vocabulary with the exception of those generating positive or negative feedback. All activating strains function within the sensor range of their receiver strains (Figure 2D), and the models correctly fit or predict the data. On the other hand, the output of repressing strains (Figure 2E) did not fully cover the input range of their receiver strains, as seen in the limited fold-change of the fluorescent signal. Models suggested that increasing the sender strain population to 10-fold its original value would produce a more noticeable response, which we successfully verified experimentally (purple dots and lines, Figure 2D and E).
After the two-strain combinations, we also verified the predictive power of our model on two 3-strain chains with different topologies and different strain stoichiometries (Figure 2F). Here too, the simple model strategy we outlined earlier captured the overall dynamics even when different strain concentrations were used. These results support the hypothesis that multicellular circuits behave like a sum of their individual parts (modularity), are easy to modulate (whether through altering initial strain concentrations or signal degradation) and can be extended to longer chains.
Increasing nonlinear response using external positive feedback
Thus far, we explored ways to simulate and design multicellular circuits with tunable gains to obtain monotonic, quasi-linear dynamic systems with a single equilibrium point. To extend the range of observable behaviors and generate non-linear responses to the inputs, we used positive feedback strains that sense and secrete the same signaling molecule (Figure 3A).
To highlight the increase of nonlinear response, on top of fitting models to the positive feedback strains (as in Figure 2A), we also re-fitted the sensor strains to estimate a new Hill coefficient (much as in Shaw et al., 2019). We then plotted these fitted curves, the data points and the two Hill coefficients for both the feedback system and the nominal response (in this case, the sensor alone). In both cases, the positive feedback increased the nonlinearity of the response (from 0.8 ± 0.0 to 1.5 ± 0.1 and from 1.0 ± 0.4 to 1.2 ± 0.3 for the α-factor and the auxin case respectively). We also considered positive feedback circuits that operate through double repression (Figure 3B). In this case, we tested topologies where either α-factor represses BAR1 synthesis, or auxin represses GH3.3 synthesis. Ideally, at low signaling molecule concentration, signal degradation through BAR1 or GH3.3 is predominant so there is no fluorescence response in the receiver cells. But at high input concentration, the degradation pathway is switched off and the input is free to reach the receiver cells. As in Figure 3A, we measured the Hill coefficients of these circuits and reported an increase in nonlinear response (from 0.8 ± 0.0 to 1.0 ± 0.0 and from 1.0 ± 0.4 to 1.4 ± 0.6 for α-factor and auxin circuits respectively).
Constructing a multicellular bistable switch
We next turned to the construction of a bistable switch circuit. We opted to use a mutual-repression topology to generate bistability (Gardner et al., 2000; Oyarzún and Chaves, 2015). For a first design, we combined the strain that senses α-factor and represses auxin synthesis with the strain that senses auxin and represses α-factor synthesis. Here, the main state variables are signaling molecule concentration in the media rather than the internal state of the cells (specifically, high auxin/low α-factor and low auxin/high α-factor). However, both model simulation and experimental data showed that this circuit can generate only a single equilibrium, independently of the strain stoichiometry (Figure 3—figure supplement 3). This result is not unexpected given the low Hill coefficients of the two repressive circuits (0.8 and 1, respectively). Rather than re-designing the gene regulatory circuits internal to these strains to be more suitable for a bistable switch, we decided to take advantage of the composable nature of the system and incorporate more strains to the architecture.
As shown in Figure 3B, combining strains that work cooperatively increases the Hill coefficients. Hence, to boost non-linearity, we added two strains to the mix that induce signal degradation: a strain that senses α-factor and synthesizes GH3 and one that senses auxin and synthesizes BAR1 (Figure 3C). The resulting system can still be seen as two modules that repress each other’s activity: α-factor lowers auxin concentration (through pathway repression and GH3 expression) and auxin reduces α-factor concentration (through pathway repression and BAR1 expression). We studied this new system with a steady-state model for auxin and α-factor concentrations (Figure 3C, Model 3 and Appendix 1 for model derivation). The model suggests that one of the two Hill coefficients (n1) is higher than 2, a necessary condition for the existence of more than one equilibrium.
We investigated the existence and properties of the equilibria while varying the individual concentrations of the four strains in the circuit: these variations are captured by the four parameters in the model, representing the gains of the four strains in the circuit as explained earlier. Using the model, we identified a range of concentrations predicted to result in multiple equilibria (for the analysis, see Appendix 1). For further investigation, we picked a set of concentrations that maximize the distance between equilibria such that the equilibria are robust to small variations in strain concentrations (red diamond, Figure 3D). This solution generates two nullclines that intersect three times, resulting in 2 stable (red crosses) and 1 unstable (black empty circle) equilibrium as expected (Figure 3E).
We tested this model-guided design experimentally (Figure 3F). Strains were mixed according to the chosen concentrations and an auxin negative sensor strain was added. Upon reaching log phase (Time 0 in Figure 3F), samples were diluted at regular time intervals to prevent the cells from saturating (see the Materials and methods section). To alternate between the states, we exogenously added first auxin (at 3 hr) and then α-factor (at 15 hr) and let dilution reduce their concentration to below detectable levels for our sensor (at 12 and 24 hr, respectively).
Model 3 simulation and data largely agree, although there is a lag between the two. We hypothesize that the lag is due to a signaling delay between the strains that was not fully captured by the models in this five-strain circuit (four strains for bistability plus one auxin sensor). More importantly, the ‘high’ equilibrium is stable in the first 3 hr of the experiment, but seems to become unstable after 25 hr, as shown by a negative trend in the data toward the ‘low’ equilibrium. Modeling suggests that a small increase in the K4 gain (strain: IAA expressing BAR1) would lead the system to a single ‘low’ equilibrium (Figure 3—figure supplement 3). This could occur if the corresponding strain grows slightly faster than the other strains (although, in practice, we could not detect any growth difference between the single-input/single-output strains using an ANOVA test, Figure 3—figure supplement 3). An alternative explanation is that metabolic load affects the growth rates depending on strain circuit activity (Sadeghpour et al., 2017; Zhang et al., 2020). Accordingly, active strains grow more slowly than inactive ones, affecting their intended concentrations over time, which could result in leaving the bistability region. To test this hypothesis, we introduced metabolic load to our models (Appendix 1). Following Sadeghpour et al., 2017, we assumed metabolic load to be linearly proportional to (normalized) gene expression This dependency is captured by a parameter , where represents no impact and represents a high impact on growth rate (no growth at maximum gene expression). Simulations show that a behavior qualitatively similar to the loss of stability after 25 hours is possible for high metabolic load (), although at a later time than the experimental data (∼36 hrs, Figure 3—figure supplement 3) according to this model representation. Moreover, since growth rates are statistically identical for active and inactive strains (Figure 3—figure supplement 3), we concluded that metabolic load is unlikely to significantly affect our bistable circuit. This result is consistent with previous work supporting robustness of the repressive circuit topology to growth feedback (Zhang et al., 2020).
Automated design of strain circuits to generate logic gates
In order to expand our target behaviors, we developed an automated approach to select circuits using the twenty-four strains introduced above. Specifically, we simulated all possible strain combinations for networks of size 2, 3, and 4 strains using Models 1 and 2. For each network, we simulated the system response over 12 hours to all possible combinations of α-factor (in the discretized range of nM), auxin ( nM), and β-estr ( nM). The ranges were chosen to properly sample the operational range of the strains (see Figure 1C, D and E).
Next, we screened our simulation space for steady-state behaviors whose profiles resemble AND, OR, NAND, or NOR logic gates. For each strain combination, we restricted the steady-state behavior space to simulations obtained using combinations of the highest input concentrations (α-factor = 1000 uM, auxin = 5000 µM, β-estr = 100 nM) or no input to match the logic table input entries. This selection resulted in 120 combinations from the 2-node networks, 560 from the 3-node networks, and 1,820 from the 4-node networks. After normalization to allow for comparison between different sensors (Appendix 2), each of these combinations was scored according to how well they match one of the logic gate truth tables. If the predicted output for one of the expected OFF states was higher than the output for one of the expected ON states then the metric would return a value below 1 labeling that strain combination to be a poor logic gate realization. On the other hand, values above 1.0 imply a match between the output vector and the target truth table: the higher the value, the better the separation between the ON and OFF state.
We first applied this method to identify circuit topologies that generate AND gate profiles (Figure 4A). Each network topology is represented as a diamond, color-coded according to the network size (blue for 2-node, red for 3-node and yellow for 4-node networks), and divided in 4 groups according to the sensor strain reporter. Each network was scored, and its value reported on the x-axis.
Next, we selected all the network topologies that scored higher than 1.0 and performed an optimization step. Optimization aimed to maximize the target metric using strain concentration (gains) as optimization parameters (Figure 4A, left panel). Finally, we picked the highest scoring topology for experimental testing. The best AND gate was a 4-node network topology that used α-factor and β-estradiol as inputs and a negative auxin sensor to determine the output.
This strain combination included two strains repressing auxin output in the presence of α-factor and β-estradiol respectively. These strains by themselves should ideally generate an AND gate in concert with the negative auxin sensor. The other two strains improve performance of the AND gate, likely by reducing the effect of leaky auxin synthesis from the α-factor-sensing/auxin-repressing strain. The experimental realization of this circuit shows separation between the ON and OFF states. In fact, the data (blue bars) seem to slightly outperform the predictions (Figure 4A, right panel).
We repeated this same procedure for NOR, NAND, and OR gates (Figure 4B). The optimal NOR gate is also a four-node network, while the optimal NAND gate is a three-node network, and both use the same auxin sensor as the optimal AND gate. The bar separating ON and OFF states as defined for the AND gate holds for all the gates sharing the same sensor. The optimal OR gate is the naive realization with two strains synthesizing auxin from different inputs and an activating an auxin sensor. All the gate profiles are close to their predicted values from simulation, showing the high degree of modularity of our vocabulary of strains, even for complex systems with internal feedback as the NOR gate architecture. For each gate, we also tested the optimal realization of the strain combinations for the two remaining network sizes with similarly positive results (Figure 4—figure supplement 4).
Identification of circuit designs for time pulses and band pass filters
We next extended our automated design strategy to circuits that either generate time pulses or that acts as band pass filters on the input signal concentrations. Starting from the simulation dataset of all possible strain combinations, we selected all those that displayed non-monotonic behaviors (having at least one local maximum/minimum) as a function of time or as a function of the steady state input concentration. We defined a non-monotonicity metric as the distance between the local maximum/minimum and the maximum/minimum between the initial value and the final value of the series. Higher values of the metric hint to more pronounced non-monotonic behaviors, while 0 implies that no local minimum/maximum is present (Figure 5A). We then selected the top six candidates (one for each sensor type in Figure 1B, D and E, excluding the β-estradiol ones), and performed optimization to maximize the metric using cell concentrations and input concentration as parameters. Finally, we experimentally tested the top two topologies overall for both time pulses and steady-state band-pass filters. It is worth noticing that only 0.52% of all possible topologies across network sizes generated non-monotonic behaviors in time (of about 106 in total, Figure 5A) and only 0.32% at steady state (of about 105 in total, Figure 5B). Hence, a ‘brute force’ experimental approach to test all possible strain combinations would be evidently out of reach.
The highest performing time pulse topology (Figure 5C) is induced by α-factor that activates both fluorescence expression and auxin synthesis. In turn, auxin induces BAR1 production, which degrades the exogenous α-factor signal: unsurprisingly, this is an incoherent feed-forward loop, type 1 (as in Mangan and Alon, 2003). The optimal ‘reversed’ time pulse, that is, a dip in the output at intermediate times (Figure 5D), responds to α-factor induction and implements a modified incoherent feed-forward loop type 3, where α-factor both represses and activates fluorescences. Model predictions suggested that three different α-factor concentrations would generate this behavior at different capacities (0.1 nM, 1 nM and 5 nM α-factor).
According to the models, both these nonmonotonic behaviors are a consequence of delay between an activator and an inhibitor pathway.
The low-pass concentration filter (Figure 5E) uses very similar components but the sensor with opposite topology. In this case, α-factor activates fluorescence expression but also auxin synthesis, which then activates its own expression in a positive feedback. This topology is also an incoherent feed-forward loop (type 1) and results in a band-pass filter output. Finally, the band-stop filter (Figure 5F) uses both β-estr and α-factor as inputs, and the circuit topology presents a combination of negative feedback and an incoherent feed-forward loop (type 3). The circuit operates as a ‘reversed’ band-pass filter over α-factor concentration, while β-estr is kept constant and only adds a baseline of α-factor and auxin through the two top strains in the circuit.
The experimental data qualitatively agreed with the predictions, although we notice a time delay for the first two circuits and a shift in baseline for the latter two. The time delay is caused by underestimation of sensor strain activation due to difficulty in isolating the sensor strain from the other strains through gating of the fluorescence data (Figure 5C) and possibly by an under-modeled delay of auxin synthesis (Figure 5D). The baseline shift is consistent with a higher concentration of auxin in the system in both cases. As a matter of fact, these two experiments (Figure 5E F) ended at a higher cell concentration than usual (running time of 12 and 14 hr, respectively), and it is known that yeast natively synthesizes auxin at saturation (Rao et al., 2010).
Discussion
Here, we demonstrated the potential of synthetic multicellular circuits to generate a wide range of behaviors starting from simple activating or repressing individual strains. Instead of implementing complex circuits in isogenic populations, we designed simple monotonic circuits in different strains and allowed them to communicate using just two signaling molecules (α-factor and auxin) or to affect their environment by attenuating those signals (through BAR1 and GH3.3). These simple constructs alone were capable of recapitulating many behaviors previously realized with synthetic gene circuits such as bistability, band-pass filters, pulses, and logic gates. While these behaviors were previously realized with fewer strains, it is important to point out that none of the strains in this study were designed with specific behaviors in mind other than monotonic activation/repression. Strain composition here plays the role of genetic fine-tuning, with a significantly lower experimental cost while maintaining high predictability. Strains with higher nonlinear response would simplify the topologies presented in this study but would not demonstrate the property of composition as much. Moreover, our last results on automated identification of circuit topologies that realize a target behavior (Figure 5) hint that the space of possible functional circuit architectures is larger than we explored.
We demonstrated through an extensive use of mathematical models, that these synthetic multicellular circuits are modular, easy-to-tune and extendable. Modularity is achieved through cell-cell communication that avoids cross-talk, and is demonstrated by combining input and output of our simple models. Multi-cellular circuits are tuned using different cell concentrations or positive-feedback architectures. Finally, we realized that we could tune the circuit behavior by extending or shortening the length of the signaling chain by adding or removing intermediate strains. We exploited this property to build the two time-pulses (Figure 5, A and B): a slower repression or activation dynamic allows for the opposite signal to operate first, generating the non-monotonic behaviors we observed. We imagine this property will gain more practical applications when the number of signaling molecules increases, which is the current major limitation in our system.
We reported that our most complex circuits display some discrepancies between simulations and experimental data, whose explanation is not always clear. It is possible that the delay between sending and receiving is not adequately modeled, or that factors that were not modeled, such as the effect of metabolic load on cellular population distribution (Blanchard et al., 2018), are affecting those circuits. Future work should account for these factors at the modeling steps, for instance using distributions to describe input-output relationships Thurley et al., 2018 or more complex models.
Finally, we leveraged the mathematical description to define an automated method to design behaviors according to performance specifications. Computationally, the method simulates all the possible strain combinations given a fixed number of nodes, and then scores them according to how well they reproduce the behavior of interest. As it is defined, the method is not easily scalable because the number of total simulations increases exponentially with the number of strains and does not account for differences in the initial strain concentrations (also defined as ’gains’ in this study). However, we found no catastrophic failures in the experimental validations, such as a qualitatively different behavior, owing to the modularity of our system. Future efforts should be directed towards more efficient ways to simulate the networks, for instance by training a neural network on the current simulation sets to predict the output of interactions.
Materials and methods
Construction of yeast strains
Request a detailed protocolYeast transformations were carried out using a standard lithium acetate protocol used by Gander et al., 2017. Yeast cells were made competent by growing 50 ml cultures in rich media to log growth phase, then spinning down the cells and washing with sterile deionized water. Next, linearized DNA, salmon sperm donor DNA, 50% polyethylene glycol and 1 M LiOAc were combined with cell pellet and the mixture was heat shocked at 42° for 15 min. The cells were then spun down, supernatant was removed and they were resuspended in sterile deionized water and then plated on selective agar media. Transformations were done into Saccharomyces cerevisiae strain MATa W303 - 1A.
Cytometry
Request a detailed protocolFluorescence intensity was measured with a BD Accuri C6 flow cytometer equipped with a CSampler plate adapter at room temperature using excitation wavelengths of 488 and 640 nm and an emission detection filter at 533 nm (FL1 channel). A total of 10,000 (20,000 for the bistable switch samples) events above a 400,000 FSC-H threshold (to exclude debris) were recorded for each sample using the Accuri C6 CFlow Sampler software at fast flow rate (66 µL/min, 22 μ core). Cytometry data were exported as FCS 3.0 files and processed using a custom Python script to obtain the mean FL1-A value for each data point.
Data collection for sensor strains
Request a detailed protocolSynthetic complete growth medium was used to grow the cells overnight from glycerol stock, while 300 µM IAM (Indole-3-acetamide) was added in all the experimental medium (also synthetic complete). All yeast cultures in all experiments of this study were grown in a 30° shaker incubator at 250 rpm in 14 ml Corning Falcon polypropylene round-bottom tubes (cat #: 352059) unless otherwise stated. Experiments involving time course data were taken during log phase via the following preparation: 16 hr of overnight growth in the synthetic complete medium in a shaker incubator followed by dilution to 30 events/μL into fresh, room-temperature medium. After 10 hr of growth at 30°, we performed a new dilution to 30 events/μL in 3 ml of medium, added the inducers from dilution aliquots of 100 µM stocks, and started collecting 100 µL samples without replacement for measurements periodically until the completion of the experiment.Ten thousand events were collected for each condition. Experimental replicates are intended as biological replicates of the same overnight sample (replicates conducted on the same day) or the same glycerol stock (replicates conducted on different days). Fluorescence data was labeled as outlier and discarded if the negative control differed with previous measurements to an evident amount.
Data collection for co-culture experiments
Request a detailed protocolSample preparation was conducted as described above with each strain in a separate test tube. At the start of the experiment, each strain concentration was first measured and then strains were combined at the concentration specified by the experiment. We considered the concentration of 30 events/μL as the baseline concentration for all the experiments and all concentrations are expressed as multiples of this reference concentration. Once the samples were added together at the desired concentrations, inducers were added and measurements were taken periodically until the completion of the experiment as described for experiments with sensor strains. For the duration of the experiment, the samples were kept in a 30° shaker at 250 rpm. Experimental replicates are intended as biological replicates of the same overnight sample (replicates conducted on the same day) or the same glycerol stock (replicates conducted on different days). Fluorescence data was labeled as outlier and discarded if the negative control differed with previous measurements to an evident amount.
Data collection for bistable switch data
Request a detailed protocolThe sample preparation was conducted as described above, each strain in individual test tubes. Then, each of the five strains in the bistable switch were combined together at the concentration as outlined in the main text in a 3 ml test tube kept in a 30°C shaker at 250 rpm for the duration of the experiment. Every 40 min, we performed a dilution with fresh media with 300 µM concentration of IAM, through manual pipetting. Samples of 120 µL were collected from the sample without replacement. Alternatively, samples were collected from the dilution discard when available for the duration of the experiment approximately every 3 hr. Experimental replicates are intended as biological replicates of the same overnight sample.
Model fitting procedure and simulations
Request a detailed protocolModel parameters were estimated using Matlab fminsearch function to minimize the L2-norm of the difference between observations and simulations. For sensor strains in Figure 1, each parameter was estimated three times on three different experimental repeats to identify mean and standard deviation. Then, the parameters used for all the simulations in this study were estimated by fitting the average of the three measurements. Model parameters for the non-sensor strains were estimated by minimizing the L2-norm of the difference between the simulations and the average of three experimental repeats as data points. Models were simulated using the Matlab code 15s function. All models parameter values are contained in Appendix 3.
Codes and data availability
Request a detailed protocolCodes and data are available at https://github.com/Alby86/MulticellularYeast, (Carignano, 2021 copy archived at swh:1:rev:7dc5d2f016054123df8a2bdbdd9543a06e9be63d).
Appendix 1
Model derivation for the bistable switch
Starting from the steady state expressions with general parameters:
One obtains that the steady-state quantities of signaling molecules are:
where IAA0 and α0 represents exogenous concentrations added to the mix, and K1, K2, K3, and K4 are the initial concentrations of each cell type (where K = 1 is the standard value).
Substituting the above expressions (1)-(4) in (5) and (6) and re-labelling the parameters, one obtains:
If one simplifies the expression above considering only the highest and lowest terms, one obtains:
Finally, if one relabels the parameters using and , obtains the same expression as in Figure 3C, Model 3.
Expressions (7) and (8) have then been fitted to simulation data of the correspondent strain combinations for different α-factor and IAA concentration to estimate the values of the parameters. Accordingly, the new exponents are: and .
Equation | |||||||
---|---|---|---|---|---|---|---|
IAA | 101.90 | 2.003 | 4.79e + 03 | 6.16e + 06 | 494.41 | 3.21e-12 | 17.59 |
Equation | |||||||
---|---|---|---|---|---|---|---|
α | 428.28 | 1.78 | 1.65e + 04 | 3.28 | 1.406e + 03 | 3.59e + 04 | 233.68 |
To study the existence of two stable solutions, we solved system (9) + (10) numerically with different values of the parameters ( and identified the area when multiple solutions arise. In general, bistability arise for:
1. low concentration of the IAA-repressing strain () and of the BAR1-synthesis strain ()
2. high concentration of the α-factor-repressing strain () and the GH3-synthesis strain ()
In the bistability region, higher values are paired with lower values in their admissible range (or vice-versa) and higher values are paired with lower values in their admissible range (or vice-versa). The simulation in Figure 3D in the main manuscript is generated by selecting =0.2 and , while varying the values of between 0.1 and 30, and in the range between 0.1 and 1.
To simulate the evolution of α-factor and IAA concentrations over time, we extended the steady state equations (9) and (10) above to an ODE description by fitting the following to a time series simulation of the four strains on a 1:1:1:1 ratio:
where and . These estimated models were used to simulate the α-factor and IAA concentrations that were then fed to the IAA repressing GFP sensor model in Figure 1B: the resulting simulation is shown in Figure 3F overlapping the experimental data.
As suggested by one of the reviewers, we considered the possibility that metabolic weight affects cell growth, which in turn could affect the stability of the equilibria. As in [3], the metabolic weight is proportional to the functions and and decreases the degradation rates and (higher metabolic rate implies slower dilution rate) and the cell population ratios (in this context, the parameters and ). Assuming equal exponential growth for all the strains (parameters estimated in Figure 3—figure supplement 3), we obtained:
where and are the same functions introduced before but normalized to 1, and determines the impact of the metabolic load on the growth rates. Notice that for , we obtain back (11) and (12). The parameters and are now a function of time as follows:
Simulations with varying values of from 0 to 1 (Figure 3—figure supplement 2) shows that high metabolic load affects the stability of the equilibria, similarly to what is shown in [3].
Appendix 2
Normalization procedure for logic gate vectors
Each simulation ran 3 differential equations for each node, plus pooling parameters that represent the overall concentration of auxin and alpha-factor accounting for exogenous inputs, strain secretion, and BAR1 or GH3 synthesis, the two latter quantities being dependent on the strain selection. To automate the simulations, we represented each strain as a vector of parameters: the first 9 entries being the differential equation parameters, followed by one-hot encoding vector of size 3 to represent the sensed input ([beta-estr, alpha-factor, auxin]), and a one-hot encoding vector of size 4 to represent the secreted output ([alpha-factor, auxin, BAR1, GH3]). The last entry of the vector is a 0 for repressing strains and 1 for activating strains. Along the selected strains, we simulated the 8 sensor strains (Figure 1C, D and E) as system outputs.
We then grouped the simulations according to the input entries of the logic tables obtained from the 3 possible permutations of the inputs ([beta,alpha],[beta,auxin],[auxin,alpha]). These groups are shaped as a vector with 4 entries that correspond to the input states ([0 0], [1 0], [0 1], [1 1]). Within each group, we normalized each vector by their minimum value and subtracted 1 (the minimum becoming 0). This normalization scheme allowed comparison between different topologies accounting for differences in the sensor strains. At the end of this process, we had 120 vectors resulting from the 2-node networks, 560 from the 3-node networks, and 1,820 from the 4-node networks.
Appendix 3
All model parameters
Parameter fit values Table describing parameter estimated for Models 1 in Figure 1. Parameters were estimated using Matlab fminsearch function to minimize the L2 norm of the difference between observation and simulation. Each parameter was estimated three times on three different experimental repeats. Parameter mean and standard deviation are reported in the table against strain number and function using the ->/activation and -|/repression formalism. The parameters used for all the simulation in this study were estimated by fitting the average measurement directly, and it is presented in brackets below the mean ± std values.
Parameter | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|
K1 | δ1 | K2 | ψ | N | δ2 | K3 | δ3 | B | ||
Strain # | 1: IAA-|GFP | 1.32e + 06 ± 5.79e + 05(1.76e + 06) | 3.41e + 06 ± 4.25e+05(3.13e + 02) | 7.08e + 05 ± 2.75e + 05(8.29e + 05) | 60.52 ± 84.57(256.87) | 0.59 ± 0.24(0.89) | 2.46e + 05 ± 2.51e + 04(4.41e + 5) | 2.50e05 ±3.52e + 04(2.12e + 05) | 5.39 ± 1.81e + 94(7.46e + 04) | 864.63 ± 166.77(867.98) |
2: IAA-|2xGFP | 0.65 ± 0.27(0.69) | 326.27 ± 101.52(416.02) | 201.66 ± 59.58(277.36) | 2.07 ± 0.51(1.79) | 1.01 ± 0.30(1.011) | 10.99 ± 1.63(11.40) | 0.0019 ± 7.45e-04(0.0011) | 0.498 ± 0.019(0.49) | 0.0044 ± 0.0023(0.0049) | |
3: IAA->GFP | 9.64e + 04 ± 2.39e + 04(8.95e + 04) | 0.092 ± 0.086(0.082) | 2.53e + 07 ± 8.37e + 06(1.73e + 07) | 1.097e + 07 ± 1.835e + 06(1.24e + 07) | 0.83 ± 0.0059(0.836) | 1.66e + 07 ± 2.50e + 06(1.88e + 07) | 2.66e + 04 ± 4.95e + 03(3.15e + 04) | 2.25e + 06 ± 4.47e + 05(1.66e + 06) | 1.98e + 04 ± 5.49e + 03(1.46e + 04) | |
4: alpha-|GFP | 5.73 ± 2.35(6.05) | 51.20 ± 18.34(59.54) | 4.46 ± 4.15(10.61) | 0.72 ± 0.17(0.67) | 0.80 ± 0.0015(0.80) | 0.79 ± 0.23(0.82) | 0.0026 ± 0.0028(3.66e-04) | 1.89 ± 1.59(1.10) | 0.0054 ± 0.0045(0.0032) | |
5: alpha->Venus | 7.49e + 05 ± 1.23e + 06(1.06e + 03) | 3.19e + 05 ± 6.36e + 05(0.245) | 3.06e + 06 ± 3.6e + 06(3.47e + 05) | 8.54e + 05 ± 4.96e + 05(7.08e + 05) | 1.60 ± 0.88(1.04) | 1.77e + 06 ± 2.58e + 06(1.56e + 05) | 8.036e + 04 ± 1.23e + 05(1.62e + 04) | 8.20e + 06 ± 1.19e + 07(2.15e + 06) | 2.36e + 04 ± 3.46e + 04(5.96e + 03) | |
6: beta->GFP | 48.74 ± 10.42(50.66) | 35.04 ± 10.09(25.86) | 8.77 ± 2.41(11.00) | 38.14 ± 12.70(57.12) | 1.26 ± 0.031(1.26) | 87.98 ± 16.44(110.43) | 0.093 ± 0.02(0.089) | 0.16 ± 0.0035(0.16) | 2.13e-04 ± 3.68e-05(2.155e-04) | |
7: beta-|GFP | 1.18 ± 0.6(0.89) | 0.18 ± 0.13(0.17) | 243.78 ± 62.08(174.32) | 9.16 ± 1.93(6.65) | 2.12 ± 0.26(2.18) | 0.58 ± 0.097(0.56) | 1.82e-04 ± 7.53e-05(1.5e-04) | 0.63 ± 0.073(0.55) | 6.32e-04 ± 1.97e-04(3.77e-04) | |
Units | Molecule NucVol–1nM–1hour–1 | hour–1 | Molecule NucVol–1 hour–1 | Molecule NucVol–1 | dimensionless | NucVol Molecule–1 hour–1 | Fluorescence Arbitrary Units (FAU) hour–1 | NucVol Molecule–1 hour–1 | FAU Molecule NucVol–1 hour–1 | |
Description | Ligand binding affinity of membrane protein(α)/co-factor (IAA)/ TF (β) | Dissociation/degradation/dilution rate of ligand-binder complex | Maximum/minimum transcription rate of TF activator/repressor | Ligand-binder complex concentration producing half occupation | Hill-coefficient | Dissociation rate of DNA-TF | Fluorescence protein accumulation over time normalized to cell size | Fluorescence protein degradation/dilution rate over time | Baseline fluorescence protein accumulation over time normalized to cell size |
Table presenting the parameters estimated for Model 1 for Single-Input/Single-Output strains in Figures 2 and 3. In the first column, the strain number and function is reported; the second column reports the sensor strain architecture that was used as baseline but with different output pathway; the parameters here reported refers uniquely to the output pathway: the remaining parameters are identical to the ones estimated for the baseline strain. Each parameter was fitted to data representing the average of three experimental repeats.
Parameter | |||||
---|---|---|---|---|---|
Baseline strain # | K3 | δ3 | b | ||
Strain # | 10: alpha->IAA | 5 | 566.24 | 0.575 | 55.83 |
11: beta->BAR1 | 6 | 26.038 | 1.92e-06 | 0.35 | |
12: beta->IAA | 6 | 3.48e + 03 | 0.16 | 0.21 | |
13: beta->alpha | 6 | 121.60 | 0.062 | 0.14 | |
14: IAA->alpha | 3 | 2.285 | 0.28 | 0.74 | |
15: beta->GH3 | 6 | 109.96 | 36.71 | 1.60e-04 | |
16: IAA->BAR1 | 3 | 1.89 | 1.83e-13 | 0.365 | |
17: alpha->GH3 | 5 | 582.32 | 368.67 | 1.17e-14 | |
18: alpha-|IAA | 4 | 4.09 | 8.74e-11 | 47.78 | |
19: beta-|IAA | 7 | 2.024e + 11 | 4.29e + 10 | 1.85e + 09 | |
20: beta-|alpha | 7 | 0.077 | 1.77 | 0.18 | |
21: IAA-|alpha | 1 | 471.73 | 0.41 | 1.34 | |
22: alpha->alpha | 5 | 419.52 | 2.10e + 04 | 2.32e + 04 | |
23: IAA->IAA | 3 | 775.05 | 0.84 | 780.90 | |
24: alpha-|BAR1 | 4 | 0.0014 | 3.74e + 07 | 1.096e + 08 | |
25: IAA-|GH3 | 1 | 225.51 | 42.46 | 3.14e-06 | |
Units | nM hour–1 | NucVol Molecule–1 hour–1 | nM Molecule NucVol–1 hour–1 | ||
Description | Output protein accumulation over time | Output protein degradation/dilution rate | Baseline output protein accumulation over time |
Table reporting parameter values for the two Multi-Input/Single-Output (MISO) strains in Figure 1E. The parameters were estimated by minimizing the L2 norm of the difference between the simulation of Model 2 (Figure 1) and the average of three experimental repeats.
Strain # | |||||
---|---|---|---|---|---|
8: IAA->GFP|-alpha | 9: alpha->GFP|-IAA | Units | Description | ||
Parameter | K1 | 0.0036 | 3.14E-06 | Molecule NucVol–1nM–1 hour–1 | Ligand binding affinity of membrane protein(α)/co-factor (IAA)/ TF (β) |
δ1 | 0.0021 | 6.74E-04 | hour–1 | Dissociation/degradation/dilution rate of ligand-binder complex | |
K2 | 6.92E-05 | 1.34E-19 | Molecule NucVol–1nM–1hour–1 | Ligand binding affinity of membrane protein(α)/co-factor (IAA)/ TF (β) | |
δ2 | 0.0047 | 7.46E-05 | hour–1 | Dissociation/degradation/dilution rate of ligand-binder complex | |
K3 | 0.059 | 298.86 | Molecule NucVol–1nM–1hour–1 | Maximum/minimum transcription rate of TF activator/repressor | |
ψ1 | 0.3 | 1.3 | Molecule NucVol–1 | Ligand-binder complex concentration producing half occupation | |
n1 | 0.31 | 0.003 | dimensionless | Hill coefficient | |
δ3 | 3.235 | 574.73 | NucVol Molecule–1 hour–1 | Dissociation rate of DNA-TF | |
K4 | 74 | 0.012 | Molecule NucVol–1nM–1hour–1 | Maximum/minimum transcription rate of TF activator/repressor | |
ψ2 | 1.93 | 0.13 | Molecule NucVol–1 | Ligand-binder complex concentration producing half occupation | |
n2 | 0.85 | 0.6 | dimensionless | Hill coefficient | |
δ4 | 58.34 | 0.13 | NucVol Molecule–1 hour–1 | Dissociation rate of DNA-TF | |
K5 | 0.015 | 0.015 | FAU hour–1 NucVol Molecule–1 | Background fluorescence protein accumulation over time normalized to cell size | |
φ1 | 0.32 | 0.043 | FAU hour–1 | Fluorescence protein synthesis by the activator over time normalized to cell size | |
φ2 | 0.87 | 2.09E + 03 | dimensionless | Repressor molecule number needed to degrade one activating molecule | |
δ5 | 0.9 | 0.52 | hour–1 | Fluorescence protein degradation/dilution rate over time |
Data availability
Figure 1 - Source Data 1, Figure 2 - Source Data 1, Figure 3 - Source Data 1, Figure 4 - Source Data 1, Figure 5 - Source Data 1 contain the numerical data used to generate the figures.
References
-
A synthetic Escherichia coli predator-prey ecosystemMolecular Systems Biology 4:187.https://doi.org/10.1038/msb.2008.24
-
Circuit-host coupling induces multifaceted behavioral modulations of a gene switchBiophysical Journal 114:737–746.https://doi.org/10.1016/j.bpj.2017.12.010
-
Amplifying genetic logic gatesScience (New York, N.Y.) 340:599–603.https://doi.org/10.1126/science.1232758
-
Interlinked fast and slow positive feedback loops drive reliable cell decisionsScience (New York, N.Y.) 310:496–498.https://doi.org/10.1126/science.1113834
-
Engineering microbial consortia: a new frontier in synthetic biologyTrends in Biotechnology 26:483–489.https://doi.org/10.1016/j.tibtech.2008.05.004
-
SoftwareMulticellularYeast, version swh:1:rev:7dc5d2f016054123df8a2bdbdd9543a06e9be63dSoftware Heritage.
-
Model-driven engineering of rna devices to quantitatively program gene expressionScience (New York, N.Y.) 334:1716–1719.https://doi.org/10.1126/science.1212209
-
Automated design of genetic toggle switches with predetermined bistabilityACS Synthetic Biology 1:284–290.https://doi.org/10.1021/sb300027y
-
Emergent genetic oscillations in a synthetic microbial consortiumScience (New York, N.Y.) 349:986–989.https://doi.org/10.1126/science.aaa3794
-
Genetic circuit design automation for yeastNature Microbiology 5:1349–1360.https://doi.org/10.1038/s41564-020-0757-2
-
Modular cell biology: retroactivity and insulationMolecular Systems Biology 4:161.https://doi.org/10.1038/msb4100204
-
Modularity, context-dependence, and insulation in engineered biological circuitsTrends in Biotechnology 33:111–119.https://doi.org/10.1016/j.tibtech.2014.11.009
-
Multi-stability in cellular differentiation enabled by a network of three mutually repressing master regulatorsJournal of the Royal Society, Interface 17:20200631.https://doi.org/10.1098/rsif.2020.0631
-
Design and implementation of three incoherent feed-forward motif based biological concentration sensorsSystems and Synthetic Biology 1:119–128.https://doi.org/10.1007/s11693-007-9008-6
-
Digital logic circuits in yeast with CRISPR-dCas9 NOR gatesNature Communications 8:15459.https://doi.org/10.1038/ncomms15459
-
Programmable protein circuits in living cellsScience (New York, N.Y.) 361:1252–1258.https://doi.org/10.1126/science.aat5062
-
Interpretation of morphogen gradients by a synthetic bistable circuitNature Communications 11:5545.https://doi.org/10.1038/s41467-020-19098-w
-
A genetic bistable switch utilizing nonlinear protein degradationJournal of Biological Engineering 6:1–13.https://doi.org/10.1186/1754-1611-6-9
-
Designing microbial consortia with defined social interactionsNature Chemical Biology 14:821–829.https://doi.org/10.1038/s41589-018-0091-7
-
Transcriptional regulatory networks in Saccharomyces cerevisiaeScience (New York, N.Y.) 298:799–804.https://doi.org/10.1126/science.1075090
-
Expression-level optimization of a multi-enzyme pathway in the absence of a high-throughput assayNucleic Acids Research 41:10668–10678.https://doi.org/10.1093/nar/gkt809
-
Genetic circuit design automationScience (New York, N.Y.) 352:aac7341.https://doi.org/10.1126/science.aac7341
-
Growth defects and loss-of-function in synthetic gene circuitsACS Synthetic Biology 8:1231–1240.https://doi.org/10.1021/acssynbio.8b00531
-
Design of a bistable switch to control cellular uptakeJournal of the Royal Society, Interface 12:20150618.https://doi.org/10.1098/rsif.2015.0618
-
Bistability and oscillations in co-repressive synthetic microbial consortiaQuantitative Biology (Beijing, China) 5:55–66.https://doi.org/10.1007/s40484-017-0100-y
-
Automated design of synthetic ribosome binding sites to control protein expressionNature Biotechnology 27:946–950.https://doi.org/10.1038/nbt.1568
-
Direct evidence for ligand-induced internalization of the yeast alpha-factor pheromone receptorMolecular and Cellular Biology 14:7245–7255.https://doi.org/10.1128/mcb.14.11.7245-7255.1994
-
Multiplexing cell-cell communicationMolecular Systems Biology 16:e9618.https://doi.org/10.15252/msb.20209618
-
Cloning and characterization of a panel of constitutive promoters for applications in pathway engineering in Saccharomyces cerevisiaeBiotechnology and Bioengineering 109:2082–2092.https://doi.org/10.1002/bit.24481
-
A synthetic low-frequency mammalian oscillatorNucleic Acids Research 38:2702–2711.https://doi.org/10.1093/nar/gkq121
-
Synthetic bistability and differentiation in yeastACS Synthetic Biology 8:929–936.https://doi.org/10.1021/acssynbio.8b00524
-
Secreting and sensing the same molecule allows cells to achieve versatile social behaviorsScience (New York, N.Y.) 343:1242782.https://doi.org/10.1126/science.1242782
-
Topology-dependent interference of synthetic gene circuit function by growth feedbackNature Chemical Biology 16:695–701.https://doi.org/10.1038/s41589-020-0509-x
Article and author information
Author details
Funding
Office of Naval Research (N00014-16-1-3189)
- Alberto Carignano
- Georg Seelig
- Eric Klavins
National Science Foundation (1807132)
- Alberto Carignano
- Eric Klavins
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We are grateful to Cameron Cordray, Samer Halabiya, and Klavins lab technicians for technical support and to Mitchell Szeto and Alex Carr for help with the Python scripts.
Copyright
© 2022, Carignano et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 1,755
- views
-
- 295
- downloads
-
- 7
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Computational and Systems Biology
- Physics of Living Systems
Planar cell polarity (PCP) – tissue-scale alignment of the direction of asymmetric localization of proteins at the cell-cell interface – is essential for embryonic development and physiological functions. Abnormalities in PCP can result in developmental imperfections, including neural tube closure defects and misaligned hair follicles. Decoding the mechanisms responsible for PCP establishment and maintenance remains a fundamental open question. While the roles of various molecules – broadly classified into ‘global’ and ‘local’ modules – have been well-studied, their necessity and sufficiency in explaining PCP and connecting their perturbations to experimentally observed patterns have not been examined. Here, we develop a minimal model that captures the proposed features of PCP establishment – a global tissue-level gradient and local asymmetric distribution of protein complexes. The proposed model suggests that while polarity can emerge without a gradient, the gradient not only acts as a global cue but also increases the robustness of PCP against stochastic perturbations. We also recapitulated and quantified the experimentally observed features of swirling patterns and domineering non-autonomy, using only three free model parameters - rate of protein binding to membrane, the concentration of PCP proteins, and the gradient steepness. We explain how self-stabilizing asymmetric protein localizations in the presence of tissue-level gradient can lead to robust PCP patterns and reveal minimal design principles for a polarized system.
-
- Computational and Systems Biology
- Neuroscience
The basolateral amygdala (BLA) is a key site where fear learning takes place through synaptic plasticity. Rodent research shows prominent low theta (~3–6 Hz), high theta (~6–12 Hz), and gamma (>30 Hz) rhythms in the BLA local field potential recordings. However, it is not understood what role these rhythms play in supporting the plasticity. Here, we create a biophysically detailed model of the BLA circuit to show that several classes of interneurons (PV, SOM, and VIP) in the BLA can be critically involved in producing the rhythms; these rhythms promote the formation of a dedicated fear circuit shaped through spike-timing-dependent plasticity. Each class of interneurons is necessary for the plasticity. We find that the low theta rhythm is a biomarker of successful fear conditioning. The model makes use of interneurons commonly found in the cortex and, hence, may apply to a wide variety of associative learning situations.