Modular, robust, and extendible multicellular circuit design in yeast

  1. Alberto Carignano
  2. Dai Hua Chen
  3. Cannon Mallory
  4. R Clay Wright
  5. Georg Seelig  Is a corresponding author
  6. Eric Klavins  Is a corresponding author
  1. Department of Electrical and Computer Engineering, University of Washington, United States
  2. Department of Biological Systems Engineering, Virginia Tech, United States
  3. Paul G Allen School of Computer Science and Engineering, University of Washington, United States

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.sa0

Introduction

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).

Figure 1 with 7 supplements see all
Modular components for engineering multi-cellular signaling circuits.

(A) Three input signals, two transfer functions and five outputs are used to assemble 24 distinct strains. (B) Differential equations are used to model and predict behavior of all strains and strain combinations in the paper. In this example, a model (Model 1) and symbolic representation for a fluorescent reporter repressed by high auxin concentration are shown. (C) Left: Time course fluorescence data for different auxin concentrations for the sensor strain shown in (B). Full lines are model simulations. Right: End point fluorescence data is shown as a function of auxin concentration color-matching the time series on the left. The steady-state simulation is shown in orange. (D) Symbolic representation, steady state data and model fit for all other single input reporter strains. (E) Modeling framework (Model 2) and fluorescence kinetics data for two different two-input reporter strains. Error bars represents the s.d. of three biological replicates.

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 (n) 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 EC50 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 EC506.0±0.4 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 EC50964.8±93.8 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 (0.8±0.0, 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).

Multi-strain signaling cascades.

(A) Symbolic representation of a two-strain cascade that uses α-factor as the input and auxin as an intermediate signal. To model multi-strain cascades, models for individual strains are concatenated and only the last differential equation of the first model is fit: the GFP output (dashed yellow semi-circle) is substituted with IAA (solid blue semi-circle). End-point fluorescence data and fit are shown for the example cascade. (B) By varying the concentration of the upstream strain, the strength of the signal seen by the downstream strain can be predictably controlled. This change in signal is modeled with a single parameter K. Right: data and model predictions for three experiments with the same strains but varying concentrations of the upstream strain. (C) Symbolic representation, model and data for a two strain cascade wherein the upstream strain removes a signal rather than secreting it. (D) Symbolic representation, model and data for all two strain cascades where the upstream input is activating the production of the intermediate signal. (E) Symbolic representation, model and data for all two strain cascades where the upstream input is repressing the production of the intermediate signal. (F) Symbolic representation, model and data for three-layer signaling cascades. Error bars represents the s.d. of three biological replicates.

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 K multiplied by the output signal, where K 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).

Figure 3 with 3 supplements see all
Feedback strains and bistability.

(A) Symbolic representation, model and data for two self-activating positive feedback circuits using alpha factor (left) and auxin (right) as the activating signal. (B) Symbolic representation, model and data for two positive feedback circuits that act through double repression. (C) Architecture and model for a bistable switch assembled from four separate strains. (D) Area of bistability as a function of strain concentrations, where the dark blue color represents one stable equilibrium and color shades report predicted fold-change between ON and OFF states. The coordinates of the red diamond are the concentrations chosen for experimental testing. (E) For the chosen solution, the nullclines show the IAA and alpha-factor concentration for the three equilibria (× = stable, ο = unstable) (F) Experimental data for the bistable switch. Full line is Model 3 simulation, and exogenous auxin and alpha-factor concentrations are reported as a function of dilution over time,depicted in the figure as jagged blue and red lines (the peaks are normalized to simplify their graphical representation). Error bars represents the s.d. of three biological replicates.

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 K 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 0ρ1, where ρ=0 represents no impact and ρ=1 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 (ρ0.5), 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 [0 >1 >5 >10 >50 >200 >1000] nM), auxin ([0 >100 >500 >1000 >5000] nM), and β-estr ([0 >1 >5 >10 >100] 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 [01] 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.

Figure 4 with 4 supplements see all
Model-generated implementation of Boolean logic gates.

(A) We used an automated search algorithm to screen all possible strain combinations up to size four (and including exactly one sensor strain) for their ability to realize a set of logic functions. In the example, each strain combination is scored according to how close the combination is to realizing an ideal Boolean AND gate. Each colored diamond corresponds to one strain combination. Circuits consisting of two strains (+ a sensor) are shown in blue, three-strain gates are shown in red, and four-strain gates in yellow. All strains in these simulations are at equal stoichiometry. Higher scores indicate more AND-like behavior. The top-performing strains (score >1, teal box) are selected for further computational optimization of strain stoichiometry. Optimized strain combinations have higher AND scores. Optimized circuits consisting of two strains are shown in light blue, three-strain gates are shown in green, and four-strain gates in purple. Simulations are separated according to which sensor strain is used. A specific high-scoring four-strain combination was chosen for experimental testing (red box). The experimental data (blue bars) show good quantitative agreement with the predictions (orange bars). (B) Similar optimization procedures to the ones shown in (A) were used to identify strain combinations that realize NOR, NAND, and OR logic functions. Example implementations, model predictions and experimental data are shown for all three logic functions. Error bars represents the s.d. of three biological replicates.

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.

Model generated implementations of analog functions.

(A) An automated search algorithm was used to screen all possible strain combinations up to size four (plus one sensor strain) for non-monotonic behaviors. Specifically, we set up the search to find combinations that result in a pulse as a function of time or as a function of concentration (B). Top row: We define a non-monotonicity metric and rank all combinations according to that score. Bar graphs show the total number of possible strain combinations (blue) and the percentage that show non-linear behaviors. Results are organized by size of the strain combination and the type of target behavior. (C) Strain stoichiometries are optimized for the most promising strains to obtain more extreme maxima or minima. Left: a diagram of the highest scoring circuit for time pulses selected for testing. Center: behavior prediction after circuit stoichiometry optimization. Right: Experimental data for the time pulse circuits is shown. (D) Strain combination, model prediction, and experimental data for a system that results in a dip in fluorescence as a function of time rather than a peak. (E) Strain combination, model prediction and experimental data for a system that generates a peak at intermediate auxin concentration thus realizing a band-pass filter. (F) Strain combination, model prediction and experimental data for a system that generates a dip at intermediate alpha-factor concentrations. This system realizes a ‘band-stop’ filter, which is a combination of a low and a high-pass filter. Error bars represents the s.d. of three biological replicates.

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 protocol

Yeast 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 protocol

Fluorescence 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 protocol

Synthetic 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 protocol

Sample 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 protocol

The 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 13 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 protocol

Model 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 protocol

Codes 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:

(1) αIAA:xss=x0+ a0φ0+αn1
(2) αGH3:yss=y0+ a1αn2φ1+αn2
(3) IAAα:zss=z0+ a2φ2+IAAn3
(4) IAABAR1:wss=w0+ a4IAAn4φ4+IAAn4

One obtains that the steady-state quantities of signaling molecules are:

(5) IAA=IAA0+K1xss1+K3yss
(6) α=α0+K2zss1+K4wss

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:

(7) IAA=a1IAA0+a2IAA0αn1+a3IAA0αn2+IAA0αn1+n2+K1c1+c2αn1+c3αn2+c4αn1+n2b0+b1K3+αn1c1+c2K3+αn2d1+d2K3+αn1+n21+eK3
(8) α=f1α0+f2α0IAAn3+f3α0IAAn4+α0IAAn3+n4+K2h1+h2IAAn3+h3IAAn4+h4IAAn3+n4g0+g1K4+IAAn3j1+j2K4+IAAn4l1+l2K4+IAAn3+n41+mK4

If one simplifies the expression above considering only the highest and lowest terms, one obtains:

(9) IAA=fIAAα=a1IAA0+IAA0αn1+n2+K1c1+c4αn1+n2b0+b1K3+αn1+n21+eK3
(10) α=fαIAA=f1α0+α0IAAn3+n4+K2h1+h4IAAn3+n4g0+g1K4+IAAn3+n41+mK4

Finally, if one relabels the parameters using n~1=n1+n2 and n~2=n3+n4 , 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: n~1=2.003 and n~2=1.78.

Equationa1n~1c1c4b0b1e
IAA101.902.0034.79e + 036.16e + 06494.413.21e-1217.59
Equationf1n~2h1h4g0g1m
α428.281.781.65e + 043.281.406e + 033.59e + 04233.68

To study the existence of two stable solutions, we solved system (9) + (10) numerically with different values of the parameters Ki (0.1Ki50) and identified the area when multiple solutions arise. In general, bistability arise for:

1. low concentration of the IAA-repressing strain (0.1K11) and of the BAR1-synthesis strain (0.1K41)

2. high concentration of the α-factor-repressing strain (10K250) and the GH3-synthesis strain (1K310)

In the bistability region, higher K1 values are paired with lower K4 values in their admissible range (or vice-versa) and higher K2 values are paired with lower K3 values in their admissible range (or vice-versa). The simulation in Figure 3D in the main manuscript is generated by selecting K1 =0.2 and K2=25, while varying the values of K3 between 0.1 and 30, and K4 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:

(11) IAA˙=δIAAfIAA(α)-βIAAIAA
(12) α˙=δαfα(IAA)-βαα

where δIAA=βIAA=0.01 and δα=βα=0.49. 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 fIAA(α) and fα(IAA) and decreases the degradation rates βIAA and βα (higher metabolic rate implies slower dilution rate) and the cell population ratios (in this context, the parameters K1, K2, K3, and K4). Assuming equal exponential growth for all the strains (parameters estimated in Figure 3—figure supplement 3), we obtained:

(13) βIAA= βIAA(1ρfIAA(α))
(14) βα~= βα(1-ρfα~(IAA))

where fIAA~(α) and fα~(IAA) are the same functions introduced before but normalized to 1, and 0ρ1 determines the impact of the metabolic load on the growth rates. Notice that for ρ=0, we obtain back (11) and (12). The parameters K1, K2, K3, and K4 are now a function of time as follows:

(15) Ki~t=Ki e-ρfIAA~(α)t i=1,3
(16) Ki~t=Ki e-ρfα~(IAA)t i=2,4

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δ1K2ψNδ2K3δ3B
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-|2xGFP0.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->GFP9.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-|GFP5.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->Venus7.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->GFP48.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-|GFP1.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)
UnitsMolecule NucVol–1nM–1hour–1hour–1Molecule NucVol–1 hour–1Molecule NucVol–1dimensionlessNucVol Molecule–1 hour–1Fluorescence Arbitrary Units (FAU) hour–1NucVol Molecule–1 hour–1FAU Molecule NucVol–1 hour–1
DescriptionLigand binding affinity of membrane protein(α)/co-factor (IAA)/ TF (β)Dissociation/degradation/dilution rate of ligand-binder complexMaximum/minimum transcription rate of TF activator/repressorLigand-binder complex concentration producing half occupationHill-coefficientDissociation rate of DNA-TFFluorescence protein accumulation over time normalized to cell sizeFluorescence protein degradation/dilution rate over timeBaseline 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δ3b
Strain #10: alpha->IAA5566.240.57555.83
11: beta->BAR1626.0381.92e-060.35
12: beta->IAA63.48e + 030.160.21
13: beta->alpha6121.600.0620.14
14: IAA->alpha32.2850.280.74
15: beta->GH36109.9636.711.60e-04
16: IAA->BAR131.891.83e-130.365
17: alpha->GH35582.32368.671.17e-14
18: alpha-|IAA44.098.74e-1147.78
19: beta-|IAA72.024e + 114.29e + 101.85e + 09
20: beta-|alpha70.0771.770.18
21: IAA-|alpha
1471.730.411.34
22: alpha->alpha5419.522.10e + 042.32e + 04
23: IAA->IAA3775.050.84780.90
24: alpha-|BAR140.00143.74e + 071.096e + 08
25: IAA-|GH31225.5142.463.14e-06
UnitsnM hour–1NucVol Molecule–1 hour–1nM Molecule NucVol–1 hour–1
DescriptionOutput protein accumulation over timeOutput protein degradation/dilution rateBaseline 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|-alpha9: alpha->GFP|-IAAUnitsDescription
ParameterK10.00363.14E-06Molecule NucVol–1nM–1 hour–1Ligand binding affinity of membrane protein(α)/co-factor (IAA)/ TF (β)
δ10.00216.74E-04hour–1Dissociation/degradation/dilution rate of ligand-binder complex
K26.92E-051.34E-19Molecule NucVol–1nM–1hour–1Ligand binding affinity of membrane protein(α)/co-factor (IAA)/ TF (β)
δ20.00477.46E-05hour–1Dissociation/degradation/dilution rate of ligand-binder complex
K30.059298.86Molecule NucVol–1nM–1hour–1Maximum/minimum transcription rate of TF activator/repressor
ψ10.31.3Molecule NucVol–1Ligand-binder complex concentration producing half occupation
n10.310.003dimensionlessHill coefficient
δ33.235574.73NucVol Molecule–1 hour–1Dissociation rate of DNA-TF
K4740.012Molecule NucVol–1nM–1hour–1Maximum/minimum transcription rate of TF activator/repressor
ψ21.930.13Molecule NucVol–1Ligand-binder complex concentration producing half occupation
n20.850.6dimensionlessHill coefficient
δ458.340.13NucVol Molecule–1 hour–1Dissociation rate of DNA-TF
K50.0150.015FAU hour–1 NucVol Molecule–1Background fluorescence protein accumulation over time normalized to cell size
φ10.320.043FAU hour–1Fluorescence protein synthesis by the activator over time normalized to cell size
φ20.872.09E + 03dimensionlessRepressor molecule number needed to degrade one activating molecule
δ50.90.52hour–1Fluorescence 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

Decision letter

  1. Xiaojun Tian
    Reviewing Editor; Arizona State University, School of Biological and Health Systems Engineering, United States
  2. Naama Barkai
    Senior Editor; Weizmann Institute of Science, Israel

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Modular, robust and extendible multicellular circuit design in yeast" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen Naama Barkai as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1. One of the major concerns is that I am not sure if and how the authors consider the cell growth and growth feedback into their design and modeling? The growth rates could heavily depend on the burden caused by the circuits, and host cell growth could dilute the gene expression. Will a small difference in growth rates lead to the imbalanced ratio of the strain or even loss of one of the strains in the long term? Especially in the experiment with the concentration of the upstream strain increased to10-fold. If appliable, the authors could discuss the relevance of growth rate and growth feedback in the design of multistrain circuits when they talk about the discrepancies between the simulations and experimental data. Here are some relevant publications (PMID: 29414718, 32251409)

2. The multicellular circuits in the currents look much more complicated than the original design. For example, the original toggle switch only needs two genes. In contrast, now four strains are needed in the authors' multistrain toggle switch. Similar to the logic gates. So it will add complexity and limits the scalability and robustness. What we really need is the simplicity and modularity of the circuit design is still maintained for the multistrain circuit. An in-depth discussion is expected in the revised manuscript.

3. Looks like the increased nonlinearity with external positive feedback was not used at all for the multistrain circuit designs. If so, what is the purpose of this section? Would that be used for the bistable switch? Maybe two strains instead of four is enough to make a multistrain bistable switch.

4. In the introduction, I suggest the author discuss the potential loss of modularity in the single-strain circuit but which can be solved with multiple strains. Some relevant publications (PMID: 33558556).

5. While the modeling framework can generally be applicable to all the strains, it would be better for the audience to understand if the authors specify what x1 and χ2 represent in the circuits.

6. What is the potential mechanism of ultrasensitivity in the strain in Figure 1D 2nd panel where n=2.2 was found? Would using the 1st-order Hill repressing function add or decrease some of the nonlinearity if you compare activation cases where a linear function is used?

All the reference of figures needs 'Figure ' or 'Figure' before figure number/panel

Figure 2A, what is the yellow semicircle? Should it be removed?

What is the reason why there are only 3 data points in Figure 2D-panel 4?

How did the author decide which data points are used for fitting while the rest for validation?

7. The area of bistability in Figure 3D needs to be defined.

Page 9, typo, the toggle switch circuit should be four-strain, not five-strain.

Add unit for the growth rate in Figure S10.

8. Orthogonality or no cross-talk between signaling molecules should be supported by either references or the first-hand data in this work. Please add this information when combing signals was mentioned at the first time in line 130. This data could further support the claim in line 406 that modularity is achieved through cell-cell communication that avoids cross-talk.

9. Line 309, based on the estimated growth rates in SI Figure 10, growth rates are actually statistically different by using t-test. Why did the authors say "we could not detect any growth difference…" here?

10. Line 389, what do the authors mean by "cytometrically isolating the sensor strain from other strains"? Do you mean data gating or cell sorting? The cytometer Accuri C6 cannot do cell sorting.

11. Line 418, it is unclear what individual cell difference that the authors refer to? Please be explicit about the difference.

12. Line 426, similarly, it is unclear what differences are in the initial strain populations? Please be explicit about the difference.

13. Consider changing the chemical formula H2O (lines 435 and 438) to sterile deionized water if applicable.

14. Line 439, please specify as yeast Saccharomyces cerevisiae strain MATa W303-1A.

15. Line 441, when using CSampler plate adapter, what is the temperature of the plate exposed to during the measurement?

16. For others to reproduce the results, please provide the details.

17. Line 444, what is the flow rate of cytometer?

18. Line 451, in 16 hrs of overnight growth, what is the rpm setting of the shaker? Are the shaker settings the same across all the experiments? What is the vessel, a tube or a plate, and what catalog number from which company?

19. Line 452, in 10 hrs of growth, what is the vessel? Is it the same rpm setting? After dilution to 3 mL medium, what is the vessel, a tube or a plate?

20. Line 453, what is the stock concentration of the inducers? 1000x stock?

21. Line454, when taking 100 uL samples periodically, do the author replenish 100 uL fresh medium after each sampling?

22. Line478, did the authors replenish fresh medium after sampling 120 uL every 3 hours?

23. Line 33. This sentence "Albeit very successful, this approach shows its limitations when it comes to scalability and robustness." is redundant.

All the numbering of the figures in the main text did not start with "Figure." Please add it.

24. Line 123, A. thaliana and C. papaya should be italic.

25. Line 155, please indicate where readers could find the data of the auxin activating strains which have a 3 fold-change activation.

26. The labeling of figure panels which use uppercase A,B,C… is not consistent with the labeling in the caption which use lowercase a, b, c….

27. SI Figure 8 and SI Figure 11 are better re-named as SI Note 1 and SI Note 2 because they are not figures.

28. Line 292, should "a set concentrations" be written as "a set of concentrations"?

29. Line 296, should Figure 3E be Figure 3F?

30. It is much clear to draw the bistable switch and the fifth reporter strain in Figure 3F.

https://doi.org/10.7554/eLife.74540.sa1

Author response

Essential revisions:

1. One of the major concerns is that I am not sure if and how the authors consider the cell growth and growth feedback into their design and modeling? The growth rates could heavily depend on the burden caused by the circuits, and host cell growth could dilute the gene expression. Will a small difference in growth rates lead to the imbalanced ratio of the strain or even loss of one of the strains in the long term? Especially in the experiment with the concentration of the upstream strain increased to10-fold. If appliable, the authors could discuss the relevance of growth rate and growth feedback in the design of multistrain circuits when they talk about the discrepancies between the simulations and experimental data. Here are some relevant publications (PMID: 29414718, 32251409)

We thank the reviewer for their support and detailed feedback.

We added references to the publications suggested by the author. Furthermore, we expanded our analysis to address the effect of growth feedback into our design.

We believe growth feedback is most relevant in the context of the bistable switch, since it is our most complex circuit and involves feedback regulation. We incorporated growth feedback in our original model as proposed in [Sadeghpour et al., 2017; DOI: 10.1007/s40484-017-0100-y]. Simulations with varying metabolic weight show that the system can indeed become unstable and could provide a useful explanation for some discrepancies between the experimental data and ideal model that we previously reported. We added these new results to the main text and as Figure 3—figure supplement 2.

We edited the paper as follows (line 329-338):

“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 behavior0 ≤ ρ ≤ 1 ρ = 0 ρ = 1 qualitatively similar to the loss of stability after 25 hours is possible for high metabolic load (supplement 2) according to this model representation. ρ ≥ 0. 5 ), although at a later time than the experimental data (~36hrs, Figure 3-figure. 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)).”

Also, line 454-456 in the conclusions now reads:

“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.”

2. The multicellular circuits in the currents look much more complicated than the original design. For example, the original toggle switch only needs two genes. In contrast, now four strains are needed in the authors' multistrain toggle switch. Similar to the logic gates. So it will add complexity and limits the scalability and robustness. What we really need is the simplicity and modularity of the circuit design is still maintained for the multistrain circuit. An in-depth discussion is expected in the revised manuscript.

We thank the reviewers for bringing up the important point of scalability and robustness. The goal of our multicellular approach was to prove that a set of strains that individually can only monotonically activate or repress gene expression are sufficient to reproduce complex behaviors through composition. In the example of the bistable switch, the selected strains were not individually designed with that behavior in mind: higher Hill coefficient would be an obvious improvement for instance and could be fairly easily achieved through strain engineering. However, thanks to the composable nature of multicellular circuits, we still are able to demonstrate bistability without the need to engineer any tailor-made gene circuit.

Our goal is to show that distributed computation can be an alternative to genetic fine tuning. ‘Multicellular tuning’ has the advantage of being experimentally simpler and cheaper, while maintaining high predictive power. We incorporated these important points to the result section and the discussion.

We edited the paper as follows (line 292-295):

“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.”

And in the conclusion (line 433-439):

“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.”

3. Looks like the increased nonlinearity with external positive feedback was not used at all for the multistrain circuit designs. If so, what is the purpose of this section? Would that be used for the bistable switch? Maybe two strains instead of four is enough to make a multistrain bistable switch.

Although those strains were not directly used for the multistrain circuits, they are a proof of concept that strain composition can tune the level of nonlinearity of the response. This concept is the cornerstone of the bistable switch topology. As the reviewer points out, it is likely that, by using these positive feedback strains, we could generate a bistable system. However, while positive feedback can switch between and OFF state to an ON state, it cannot, by itself, switch back. Alternatively, one can genetically incorporate a repressive element to the external positive feedback strains, and it is possible that this strategy would reduce the number of strains in the bistable switch to two. However, we wanted to demonstrate that strain composition achieves the same results as genetically refining the strains.

4. In the introduction, I suggest the author discuss the potential loss of modularity in the single-strain circuit but which can be solved with multiple strains. Some relevant publications (PMID: 33558556).

We thank the reviewers for pointing out this reference. We have added it to the introduction (line 37-38):

“Furthermore, competition over shared gene expression resources makes unrelated circuits unintentionally coupled Zhang et al. (2021).”

5. While the modeling framework can generally be applicable to all the strains, it would be better for the audience to understand if the authors specify what x1 and χ2 represent in the circuits.

We thank the reviewers for pointing out this lack of clarity in our presentation. We have added this clarification to the text in line 179-181:

“For instance, for the auxin sensor in Figure 1B, x1 represents the TIR1-auxin complex concentration, χ2, the dCas9-VP64-auxin degron population and x3, GFP concentration.”

6. What is the potential mechanism of ultrasensitivity in the strain in Figure 1D 2nd panel where n=2.2 was found? Would using the 1st-order Hill repressing function add or decrease some of the nonlinearity if you compare activation cases where a linear function is used?

We believe that this ultrasensitivity is caused by the chromatin remodeling effect of Mxi1, which is genetically fused to dCas9 in our circuit, as suggested by [Gander et al., 2017]. The value of n=2.2 is consistent with what was estimated in that paper using a 1st-order Hill repressing function, suggesting that ultrasensitivity is not an artifact of our modeling strategy. In the manuscript, we edited lines 195-196:

“(dCas9-Mxi1 confers stronger and more consistent repression than dCas9 alone, leading to ultrasensitivity, Gander et al. (2017))”.

All the reference of figures needs 'Figure ' or 'Figure' before figure number/panel.

We thank the reviewers for pointing out this inconsistency in the text. We have edited the text accordingly.

Figure 2A, what is the yellow semicircle? Should it be removed?

We thank the reviewers for pointing out this lack of clarity in our presentation. We edited the caption of the figure to clarify the purpose of the yellow semicircle. The yellow semicircle represents the previous output (fluorescence), now swapped for auxin synthesis (blue semicircle). The strains are genetically identical with the exception of the output gene swap. We edited the caption of Figure 2A. It now says:

“To model multi-strain cascades, models for individual strains are concatenated and only the last differential equation of the first model is fit: the GFP output (dashed yellow semi-circle) is substituted with IAA (solid blue semi-circle).”

What is the reason why there are only 3 data points in Figure 2D-panel 4?

We thank the reviewers for pointing out this lack of clarity in our presentation. All strain datasets collected for Figure 2 comprise of at least 8 data points per strain. In Figure 2D-panel4, there are actually 3 data points for each of the 4 IAA concentrations, which makes for 12 data points. Hence, that dataset is in line with the other strains. The concentrations were selected to cover the operational range of the β-estr activating strains and IAA repressing strains.

How did the author decide which data points are used for fitting while the rest for validation?

We thank the reviewers for pointing this out. The two datasets come from two different experiments. We added a reference to the main text in lines 231-233:

“…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)…”

7. The area of bistability in Figure 3D needs to be defined.

We thank the reviewers for pointing this out. We have added a description of the bistability region in Appendix 1, second page, starting from the first paragraph.

Page 9, typo, the toggle switch circuit should be four-strain, not five-strain.

We thank the reviewers for pointing this out. We have added a description of the bistability region in Appendix 1, second page, starting from the first paragraph.

Add unit for the growth rate in Figure S10.

We thank the reviewers for pointing this out. We have edited the text (added min^{-1}) of Figure S10 accordingly.

8. Orthogonality or no cross-talk between signaling molecules should be supported by either references or the first-hand data in this work. Please add this information when combing signals was mentioned at the first time in line 130. This data could further support the claim in line 406 that modularity is achieved through cell-cell communication that avoids cross-talk.

We thank the reviewers for bringing up the important point of signal orthogonality. We added Figure 1—figure supplement 7 that shows a high degree of orthogonality between the 3 different signals across the repressive sensors. The signals were combined and compared to the individual signals and little difference was detected. The figure is also referenced in the main text in lines 198-202 as:

“We tested signal orthogonality by adding pairwise combinations of inducers at saturation concentration to single-input repressive sensors (Figure 1—figure supplement 7). Variations across the treatments are contained within 11% of the nominal value for signal concentrations within the range used for circuit experiments.

A maximum increase of 24% above the baseline is observed for the IAA-sensor in a regime where both α factor and auxin are at or above saturating concentrations (1uM α-factor and 10uM IAA). However, this combination is not used in any circuit experiment.”

9. Line 309, based on the estimated growth rates in SI Figure 10, growth rates are actually statistically different by using t-test. Why did the authors say "we could not detect any growth difference…" here?

We thank the reviewers for pointing out the lack of statistical validation.

We performed and added to Figure 3—figure supplement 3 an ANOVA test that does not find the growth rates to be statistically different across strains with a 90% confidence. Pairwise t-tests also show no growth difference between the strains that sense-and-secrete signaling molecules, but did find differences if comparing with the sensor strain with 95% confidence. However, changes in the sensor strain concentration do not affect the existence of two stable equilibria in the system. According to our simulations, the dilution scheme of the bistable switch does not result in the sensor strain taking over the co-culture during the duration of the experiment since the difference is minimal. Furthermore, it is not obvious that statistical tests would be meaningful with a low sample size (n=3) as in this case. We added this analysis to Figure 3—figure supplement 3 and referenced it in the main text, line 328: “…we could not detect any growth difference between the single-input/single-output strains using an ANOVA test, Figure 3—figure supplement 3.”

10. Line 389, what do the authors mean by "cytometrically isolating the sensor strain from other strains"? Do you mean data gating or cell sorting? The cytometer Accuri C6 cannot do cell sorting.

We thank the reviewers for pointing out this lack of clarity in the text. We have edited the text, and now lines 419-420 read: “…difficulty in isolating the sensor strain from the other strains through gating of the fluorescence data.”

11. Line 418, it is unclear what individual cell difference that the authors refer to? Please be explicit about the difference.

We thank the reviewers for pointing out this lack of clarity in the text. We have edited the text, and now lines 454-455 read: “…factors that were not modeled, such as the effect of metabolic load on cellular population distribution…”.

12. Line 426, similarly, it is unclear what differences are in the initial strain populations? Please be explicit about the difference.

We thank the reviewers for pointing out this lack of clarity in the text. We have edited the text, and now lines 464-465 read: “…does not account for differences in the initial strain concentrations (also defined as ’gains’ in this study)…”.

13. Consider changing the chemical formula H2O (lines 435 and 438) to sterile deionized water if applicable.

We thank the reviewers for pointing out this lack of clarity in the text. Lines 473 and 476 (former lines 435 and 438) now read: “sterile deionized water”.

14. Line 439, please specify as yeast Saccharomyces cerevisiae strain MATa W303-1A.

We thank the reviewers for pointing out this lack of clarity in the text. Line 477 (former line 439) now reads: “Saccharomyces cerevisiae strain MATa W303-1A”.

15. Line 441, when using CSampler plate adapter, what is the temperature of the plate exposed to during the measurement?

We thank the reviewers for pointing out this lack of clarity in the text. We have added to the text at line 480: “at room temperature”.

16. For others to reproduce the results, please provide the details.

Details are included below for specific questions.

17. Line 444, what is the flow rate of cytometer?

We thank the reviewers for pointing out this lack of clarity in the text. We have added line 483: “at fast flow rate (66L/min, 22 core)”.

18. Line 451, in 16 hrs of overnight growth, what is the rpm setting of the shaker? Are the shaker settings the same across all the experiments? What is the vessel, a tube or a plate, and what catalog number from which company?

We thank the reviewers for pointing out this lack of detail in the text. We have added line 489-491: “All yeast cultures in all experiments of this study were grown in a 30C shaker incubator at 250rpm in 14 ml Corning Falcon polypropylene round-bottom tubes (cat #: 352059) unless otherwise stated.”

19. Line 452, in 10 hrs of growth, what is the vessel? Is it the same rpm setting? After dilution to 3 mL medium, what is the vessel, a tube or a plate?

We thank the reviewers for pointing out this lack of clarity in the text. We have edited the text as described in the previous point.

20. Line 453, what is the stock concentration of the inducers? 1000x stock?

We thank the reviewers for pointing out this lack of clarity in the text. We added to lines 494-495: “from dilution aliquots of 100M stocks”.

21. Line454, when taking 100 uL samples periodically, do the author replenish 100 uL fresh medium after each sampling?

We thank the reviewers for pointing out this lack of clarity in the text. We added line 495 that reads: “without replacement”.

22. Line478, did the authors replenish fresh medium after sampling 120 uL every 3 hours?

We thank the reviewers for pointing out this lack of clarity in the text. We have edited lines 519-520 as follow: “…from the sample without replacement. Alternatively, samples were collected from the dilution discard when available…”.

23. Line 33. This sentence "Albeit very successful, this approach shows its limitations when it comes to scalability and robustness." is redundant.

We thank the reviewers for pointing out this redundancy in the text. We have removed the redundant line from the text.

All the numbering of the figures in the main text did not start with "Figure." Please add it.

24. Line 123, A. thaliana and C. papaya should be italic.

We thank the reviewers for pointing out this lack of clarity in the text. We have edited the text accordingly and now the numbering starts with ‘Figure’. Also, the name of the two species are now in italic.

25. Line 155, please indicate where readers could find the data of the auxin activating strains which have a 3 fold-change activation.

In the main manuscript, at line 533, we reported that all data can be retrieved on the github repository (https://github.com/Alby86/MulticellularYeast.git). The data for that strain are also available in the AllData.xlsx spreadsheet.

26. The labeling of figure panels which use uppercase A,B,C… is not consistent with the labeling in the caption which use lowercase a, b, c….

We thank the reviewers for pointing out this inconsistency in the text. We have edited the text accordingly and it is now consistently using uppercase letters.

27. SI Figure 8 and SI Figure 11 are better re-named as SI Note 1 and SI Note 2 because they are not figures.

We thank the reviewers for pointing out this inconsistency in the text. We have relabelled them as Appendix 1 and Appendix 2.

28. Line 292, should "a set concentrations" be written as "a set of concentrations"?

We thank the reviewers for pointing out this typo. Line 309 now reads: “a set of concentrations”.

29. Line 296, should Figure 3E be Figure 3F?

We thank the reviewers for pointing out this mistake. We now reference Figure 3F.

30. It is much clear to draw the bistable switch and the fifth reporter strain in Figure 3F.

We thank the reviewers for pointing out this lack of clarity in the text. We have modified Figure 3F to incorporate the reporter strain.

https://doi.org/10.7554/eLife.74540.sa2

Article and author information

Author details

  1. Alberto Carignano

    Department of Electrical and Computer Engineering, University of Washington, Seattle, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing - original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3306-9365
  2. Dai Hua Chen

    Department of Electrical and Computer Engineering, University of Washington, Seattle, United States
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  3. Cannon Mallory

    Department of Electrical and Computer Engineering, University of Washington, Seattle, United States
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  4. R Clay Wright

    Department of Biological Systems Engineering, Virginia Tech, Blacksburg, United States
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
  5. Georg Seelig

    1. Department of Electrical and Computer Engineering, University of Washington, Seattle, United States
    2. Paul G Allen School of Computer Science and Engineering, University of Washington, Seattle, United States
    Contribution
    Funding acquisition, Supervision, Writing – review and editing
    For correspondence
    gseelig@uw.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3163-8782
  6. Eric Klavins

    Department of Electrical and Computer Engineering, University of Washington, Seattle, United States
    Contribution
    Conceptualization, Funding acquisition, Supervision, Writing – review and editing
    For correspondence
    klavins@uw.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3805-5117

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.

Senior Editor

  1. Naama Barkai, Weizmann Institute of Science, Israel

Reviewing Editor

  1. Xiaojun Tian, Arizona State University, School of Biological and Health Systems Engineering, United States

Publication history

  1. Received: October 8, 2021
  2. Preprint posted: October 14, 2021 (view preprint)
  3. Accepted: March 20, 2022
  4. Accepted Manuscript published: March 21, 2022 (version 1)
  5. Accepted Manuscript updated: March 22, 2022 (version 2)
  6. Version of Record published: April 11, 2022 (version 3)

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,034
    Page views
  • 229
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

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)

  1. Alberto Carignano
  2. Dai Hua Chen
  3. Cannon Mallory
  4. R Clay Wright
  5. Georg Seelig
  6. Eric Klavins
(2022)
Modular, robust, and extendible multicellular circuit design in yeast
eLife 11:e74540.
https://doi.org/10.7554/eLife.74540

Further reading

    1. Computational and Systems Biology
    2. Genetics and Genomics
    Stella Tamana, Maria Xenophontos ... Petros Kountouris
    Research Article

    Haemoglobinopathies are the commonest monogenic diseases worldwide and are caused by variants in the globin gene clusters. With over 2400 variants detected to date, their interpretation using the ACMG/AMP guidelines is challenging and computational evidence can provide valuable input about their functional annotation. While many in silico predictors have already been developed, their performance varies for different genes and diseases. In this study, we evaluate 31 in silico predictors using a dataset of 1627 variants in HBA1, HBA2, and HBB. By varying the decision threshold for each tool, we analyse their performance (a) as binary classifiers of pathogenicity, and (b) by using different non-overlapping pathogenic and benign thresholds for their optimal use in the ACMG/AMP framework. Our results show that CADD, Eigen-PC, and REVEL are the overall top performers, with the former reaching moderate strength level for pathogenic prediction. Eigen-PC and REVEL achieve the highest accuracies for missense variants, while CADD is also a reliable predictor of non-missense variants. Moreover, SpliceAI is the top performing splicing predictor, reaching strong level of evidence, while GERP++ and phyloP are the most accurate conservation tools. This study provides evidence about the optimal use of computational tools in globin gene clusters under the ACMG/AMP framework.

    1. Computational and Systems Biology
    2. Physics of Living Systems
    Mario García-Navarrete, Merisa Avdovic ... Krzysztof Wabnik
    Research Article Updated

    Cells convert electrical signals into chemical outputs to facilitate the active transport of information across larger distances. This electrical-to-chemical conversion requires a tightly regulated expression of ion channels. Alterations of ion channel expression provide landmarks of numerous pathological diseases, such as cardiac arrhythmia, epilepsy, or cancer. Although the activity of ion channels can be locally regulated by external light or chemical stimulus, it remains challenging to coordinate the expression of ion channels on extended spatial–temporal scales. Here, we engineered yeast Saccharomyces cerevisiae to read and convert chemical concentrations into a dynamic potassium channel expression. A synthetic dual-feedback circuit controls the expression of engineered potassium channels through phytohormones auxin and salicylate to produce a macroscopically coordinated pulses of the plasma membrane potential. Our study provides a compact experimental model to control electrical activity through gene expression in eukaryotic cell populations setting grounds for various cellular engineering, synthetic biology, and potential therapeutic applications.