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


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.


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.


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


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.


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, (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.

IAA101.902.0034.79e + 036.16e + 06494.413.21e-1217.59
α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:

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

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.

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


Article and author information

Author details

  1. Alberto Carignano

    Department of Electrical and Computer Engineering, University of Washington, Seattle, United States
    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
    Competing interests
    No competing interests declared
  3. Cannon Mallory

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

    Department of Biological Systems Engineering, Virginia Tech, Blacksburg, United States
    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
    Funding acquisition, Supervision, Writing – review and editing
    For correspondence
    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
    Conceptualization, Funding acquisition, Supervision, Writing – review and editing
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3805-5117


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.


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.

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


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


  • 1,653
  • 287
  • 5

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

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
Modular, robust, and extendible multicellular circuit design in yeast
eLife 11:e74540.

Share this article

Further reading

    1. Cell Biology
    2. Computational and Systems Biology
    N Suhas Jagannathan, Javier Yu Peng Koh ... Lisa Tucker-Kellogg
    Research Article

    Bats have unique characteristics compared to other mammals, including increased longevity and higher resistance to cancer and infectious disease. While previous studies have analyzed the metabolic requirements for flight, it is still unclear how bat metabolism supports these unique features, and no study has integrated metabolomics, transcriptomics, and proteomics to characterize bat metabolism. In this work, we performed a multi-omics data analysis using a computational model of metabolic fluxes to identify fundamental differences in central metabolism between primary lung fibroblast cell lines from the black flying fox fruit bat (Pteropus alecto) and human. Bat cells showed higher expression levels of Complex I components of electron transport chain (ETC), but, remarkably, a lower rate of oxygen consumption. Computational modeling interpreted these results as indicating that Complex II activity may be low or reversed, similar to an ischemic state. An ischemic-like state of bats was also supported by decreased levels of central metabolites and increased ratios of succinate to fumarate in bat cells. Ischemic states tend to produce reactive oxygen species (ROS), which would be incompatible with the longevity of bats. However, bat cells had higher antioxidant reservoirs (higher total glutathione and higher ratio of NADPH to NADP) despite higher mitochondrial ROS levels. In addition, bat cells were more resistant to glucose deprivation and had increased resistance to ferroptosis, one of the characteristics of which is oxidative stress. Thus, our studies revealed distinct differences in the ETC regulation and metabolic stress responses between human and bat cells.

    1. Computational and Systems Biology
    2. Neuroscience
    Sara Ibañez, Nilapratim Sengupta ... Christina M Weaver
    Research Article

    Normal aging leads to myelin alterations in the rhesus monkey dorsolateral prefrontal cortex (dlPFC), which are positively correlated with degree of cognitive impairment. It is hypothesized that remyelination with shorter and thinner myelin sheaths partially compensates for myelin degradation, but computational modeling has not yet explored these two phenomena together systematically. Here, we used a two-pronged modeling approach to determine how age-related myelin changes affect a core cognitive function: spatial working memory. First, we built a multicompartment pyramidal neuron model fit to monkey dlPFC empirical data, with an axon including myelinated segments having paranodes, juxtaparanodes, internodes, and tight junctions. This model was used to quantify conduction velocity (CV) changes and action potential (AP) failures after demyelination and subsequent remyelination. Next, we incorporated the single neuron results into a spiking neural network model of working memory. While complete remyelination nearly recovered axonal transmission and network function to unperturbed levels, our models predict that biologically plausible levels of myelin dystrophy, if uncompensated by other factors, can account for substantial working memory impairment with aging. The present computational study unites empirical data from ultrastructure up to behavior during normal aging, and has broader implications for many demyelinating conditions, such as multiple sclerosis or schizophrenia.