Modular, robust, and extendible multicellular circuit design in yeast

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.


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        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 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 . 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 1figure 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    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 (x 1 ), signal processing (x 2 ) and fluorescence output synthesis (x 3 ) ( 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, x 2 , the dCas9-VP64-auxin degron population, and x 3 , 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 EC50 ∼ 6.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 EC50 ∼ 964.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 foldchange 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 cellto-cell communication To determine the potential to build biological circuits using cell-to-cell communication, we experimentally tested if communication occurs between strains that secrete an output and their corresponding sensor strains. For example, a strain that produces auxin in response to α-factor sensing was grown in co-culture with a sensor that switches off GFP expression in response to auxin ( Figure 1B). The two strain populations were added at the same concentration and the fluorescent output was measured at steady state (10 hours after induction). The experimental data is consistent with a negative α-factor sensor as expected (the more α-factor, the lower the fluorescent signal, as seen in Figure 2A).
Since the core genetic circuit of this sender strain is identical to the α-factor repressing sensor, we tested if the mathematical model previously fit to the sensor data preserves its predictive power. We re-fit only the 3 output parameters to account for the fact that the output is now auxin rather than GFP. The output of the sender cell model was used directly as input of the sensor cell model. To test the hypothesis that strains that have common input/processing parts share parameters and model structure, we collected two datasets on two different experiments composed of four data points each: we used one for fitting (yellow dots, first experiment), and the other for validation (orange dots, second experiment, Figure 2A). A time-series at EC50 input concentration was also used for fitting (not present in the figure).
Intuitively, we expected that higher initial sender cell concentration would result in an overall higher concentration of their output signal over the same time scale. Most importantly, we wanted to know if this output 'gain' can be predicted by our models so that we could use it to tune circuit behavior. We modeled this effect with a factor 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).
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).
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.
The online version of this article includes the following source data for figure 2: Source data 1. The data used for plotting Figure 2.

Figure 2 continued
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 (n 1 ) 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 K 4 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;. 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 .

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. 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.
The online version of this article includes the following source data and figure supplement(s) for figure 4: Source data 1. The data used for plotting Figure 4 and supplements.
Figure supplement 1. Optimal circuits, simulation, and experimental realization of the AND gate for 2, 3, and 4-node topologies.

Figure 4 continued on next page
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 10 6 in total, Figure 5A) and only 0.32% at steady state (of about 10 5 in total, Figure 5B). Hence, a 'brute force' experimental approach to test all possible strain combinations would be evidently out of reach.
The highest performing time pulse topology ( Figure 5C) is induced by α-factor that activates both fluorescence expression and auxin synthesis. In turn, auxin induces BAR1 production, which degrades the exogenous α-factor signal: unsurprisingly, this is an incoherent feed-forward loop, type 1 (as in Mangan and Alon, 2003). The optimal 'reversed' time pulse, that is, a dip in the output at intermediate times ( Figure 5D), responds to α-factor induction and implements a modified incoherent feed-forward loop type 3, where α-factor both represses and activates fluorescences. Model predictions suggested that three different α-factor concentrations would generate this behavior at different capacities (0.1 nM, 1 nM and 5 nM α-factor).
According to the models, both these nonmonotonic behaviors are a consequence of delay between an activator and an inhibitor pathway.
The low-pass concentration filter ( Figure 5E) uses very similar components but the sensor with opposite topology. In this case, α-factor activates fluorescence expression but also auxin synthesis,    which then activates its own expression in a positive feedback. This topology is also an incoherent feedforward 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 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 online version of this article includes the following source data for figure 5: Source data 1. The data used for plotting Figure 5.  (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.

Construction of yeast strains
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
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
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
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
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 1 3 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
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.

Appendix 1 Model derivation for the bistable switch
Starting from the steady state expressions with general parameters: One obtains that the steady-state quantities of signaling molecules are: where IAA 0 and α 0 represents exogenous concentrations added to the mix, and K 1 , K 2 , K 3 , and K 4 are the initial concentrations of each cell type (where K = 1 is the standard value).
To simulate the evolution of α-factor and IAA concentrations over time, we extended the steady state equations (9) and (10) above to an ODE description by fitting the following to a time series simulation of the four strains on a 1:1:1:1 ratio: where δ 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 f IAA ( α ) 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 K 1 , K 2 , K 3 , and K 4 ). Assuming equal exponential growth for all the strains (parameters estimated in Figure 3-figure supplement 3), we obtained: where 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 K 1 , K 2 , K 3 , and K 4 are now a function of time as follows: Simulations with varying values of ρ from 0 to 1 (Figure 3-figure supplement 2) shows that high metabolic load affects the stability of the equilibria, similarly to what is shown in [3].

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