Engineering a Model Cell for Rational Tuning of GPCR Signaling

1G protein-coupled receptors (GPCRs) are widely represented in most lifeforms and comprise the largest family of signaling proteins in humans, with over 800 members detecting structurally diverse agonists. Their abundance and ubiquity to all cell types makes them one of the most important signaling pathway classes in healthcare but also one of the most complex. Multiple types of G protein-based signaling are seen and the downstream signal transduction to activate gene expression is typically complex and intertwined with other pathways. The nature of signal transduction through the pathway also depends on many different factors, including the stoichiometry of the signaling proteins, the presence of inherent feedback mechanisms, and even cellular history. Altogether, this makes it difficult to delineate receptor and signaling properties simply from measuring the activation of downstream targets. It also makes it a major challenge to predict how changes in the levels of pathway components, for example due to different environments or mutations, can affect the performance of a given signaling pathway.

2One of the most studied examples of eukaryotic GPCR signaling is the pheromone response pathway of Saccharomyces cerevisiae, having been the focus of significant efforts from systems biology to model its actions via quantification of its behavior. To understand this pathway, researchers have parsed the contributions of numerous studies that have perturbed the dose-response and dynamics of the native system by changing growth conditions, by protein mutagenesis, via traditional gene overexpression and knockout methods, and more recently using optogenetics. While these efforts have helped to build our best picture of the events required for the transduction of signal from agonist to gene activation, inability to control the whole pathway in these experiments has meant that a complete system for exploring the dose-response relationship has not yet been achieved.

3In silico approaches typically model a system by concentrating only on the key components and varying important parameters of these such as their expression levels, while removing other non-key interactions from consideration. With advanced genome engineering and synthetic biology tools available, it now becomes possible to take an equivalent modeling approach in vivo, removing any non-essential interactions via gene knockout and finely tuning the expression of the key components using promoter libraries. This engineering approach—known as refactoring—makes a system easier to study by removing all non-essential natural regulation and feedback, thus enabling the system to be more efficiently tuned and directly measured. Effectively this generates cells streamlined for improved understanding of pathways and systems, while also making these cells more straightforward to utilize in downstream applications.

4Here, we used genome engineering to construct a heavily modified yeast suitable as an in vivo model for tuning GPCR signaling. By removing non-essential components, native transcriptional feedback regulation, and all connections to the mating response, we built a model strain retaining only the core signaling elements. In conjunction with a mathematical model, we used promoter libraries to vary the key components in this simplified, refactored pathway and uncovered principles for tuning the sensitivity, basal activity, and signal amplitude of the dose-response curve via expression level. This new knowledge provides us with a rational approach for tuning signaling characteristics and, as we demonstrate, enables us to quickly reprogram yeast to sense and measure a variety of different inputs, either in single-cell systems or community-based systems.

5Glucose sensing and the pheromone response pathway are the two native GPCR signaling pathways in S. cerevisiae, and the latter has long been the go-to choice for coupling heterologous GPCRs to yeast gene expression or for building systems for evolving GPCRs to desired targets. Core to this pathway is an extensively studied mitogen-activated protein kinase (MAPK) signaling cascade that functions with its own intrinsic feedback to maintain a robust input-output relationship in varying conditions. Because this MAPK cascade displays a graded, linear response with respect to dose and can be considered as a black-box processing unit in transduction through the pathway, we chose to make this natural system the core from which we build and tune GPCR signaling pathways.

6Keeping the five genes of the MAPK cascade fixed, we set out to generate a model strain for our work by first removing all other GPCR pathway-related genes from S. cerevisiae (Figure 1). This required making precise changes at 18 genomic loci in BY4741 yeast, generating our model strain, yWS1922, via nine rounds of CRISPR/Cas9-mediated editing (Figures S1A–S1C). Genomic changes were validated at each round by PCR and locus sequencing followed by long-read nanopore sequencing of the final strain (Figures S1D and S1E).

7During strain construction, we added addressable 24 bp targets in place of the open reading frames (ORFs) of all deleted genes. These allow for rapid and markerless (re)insertion of native or heterologous ORFs into these locations at high efficiency by CRISPR-aided multiplex integration (Figures S1F–S1H). For stable single-copy addition of further genes, three highly characterized landing pads were also introduced that interface with the MoClo Yeast Toolkit (YTK) modular cloning system, which enables rapid multigene construction from high-characterized parts. These changes were designed to facilitate rapid exploration of the effects of altering individual components of the pheromone response pathway.

8To determine how the model strain performed with all non-essential components removed, the native GPCR (STE2), Gα (GPA1), Gβ (STE4), Gγ (STE18), and pheromone responsive transcription factor (STE12) genes were restored at their natural loci to generate the “Quasi-WT” strain (Figure S1H). The α-factor dose-response of this strain was then compared to BY4741 yeast using the pheromone response FUS1 promoter driving sfGFP expression. As expected, a substantial increase in sensitivity and signal output was observed upon minimizing the signaling pathway, largely due to the removal of the native regulator of G protein signaling, Sst2 (Figure S2A).

9Previous work has shown that the sensitivity of the yeast pheromone response pathway can be changed by altering the receptor number in the absence of Sst2. Basal activity from constitutive receptor activity can also be reduced by overexpressing Gα, as this acts as a negative regulator of signaling since Gβγ propagates the response to the MAPK cascade. Therefore, to explore how sensitivity and basal activity could be varied by altering expression of the GPCR, Gα, and Gβγ, we built a mathematical model derived from a previously described cubic ternary complex model, designing this to capture G protein signaling in our minimized pheromone response pathway (Figure 2A). We then systematically probed the response of the pathway in this model by individually altering the initial GPCR, Gα, and Gβγ concentrations, while keeping all other components fixed (Figures 2D and 2E, model).

10While this model demonstrated a clear monotonic relationship between receptor number, sensitivity, and maximum signal, as previously shown by, the relationship between Gα and Gβγ levels and pathway response was non-monotonic. At lower concentrations of Gα, constitutive activation of the pathway was observed due to increased free Gβγ. This, in combination with the incorrect trafficking and instability of Gβγ, due to low levels of Gα, leads to a lower maximum-fold change in activation. At higher Gα concentrations, free Gβγ is rapidly sequestered, also leading to a decrease in the maximum fold change by acting as a “sponge” to signaling. A reverse behavior was seen when varying the initial concentration of Gβγ, with low concentrations decreasing maximum signaling and high concentrations leading to increased basal activity, albeit with a lower maximum signal, due to incorrect trafficking of free Gβγ, and therefore the Ste5 scaffold. Taken together, the model predicts an optimal level of Gα and Gβγ expression where all three members of the heterotrimeric G protein are required to be in balance to give a high-fold change in signal upon activation.

11Next, using the model strain and modular cloning system, we experimentally validated the findings of the mathematical model using a minimized pheromone response pathway. This was constructed by refactoring the GPCR, Gα, Gβγ, and transcription factor so that these are expressed from constitutive promoters that we determined to give approximately the same expression of the native promoters in their uninduced state (Figures S2B–S2G). For the Gβγ component, we expressed both subunits together as a bicistronic protein (Ste4-2A-Ste18) where the self-cleaving 2A peptide releases equimolar amounts of the two proteins. The FUS1 promoter was used to express sfGFP as the response reporter (Figure S2H), enabling pathway activation to be measured by flow cytometry. With this strain, we then individually varied the expression of the GPCR, Gα, and Gβγ by changing their promoter strengths while keeping all other components fixed. This gave in vivo results that qualitatively matched the model, demonstrating tunability of sensitivity by changing GPCR expression and revealing the predicted optimal levels of Gα and Gβγ that give the peak response (Figures 2D and 2E, experimental; Figure S3). As basal activity could be reduced by altering levels of either Gα or Gβγ, we decided to concentrate on Gα tuning, using a strain with fixed expression of Gβγ (yWS677) for all further experiments.

12Next, we sought to modulate maximum pathway output by varying the expression of the pheromone-responsive transcription factor, Ste12, as levels of transcription factor in a system often dictate gene expression strength. However, increasing the expression of Ste12 led to poor cell growth, presumably by triggering high basal activity of mating response genes (Figures S4A and S4B). We therefore needed a new means to tune pathway-mediated expression without changing Ste12 concentrations.

13By fusing the pheromone-responsive domain of Ste12 (PRD; 216-688) to the full-length LexA bacterial repressor protein, we generated a synthetic transcription factor (sTF; LexA-PRD) able to target the pathway output to a library of modular synthetic promoters containing LexA operator sequences (LexO) (Figure 3C). We could then vary the recruitment of the LexA-PRD transcription factor by altering the number of operator sequences in the promoter upstream activating sequence (UAS), enabling us to modulate the maximum output of the response over a 3-fold range without compromising the tightness of the OFF state (Figure 3D). Changing the identity of the core region of the promoter to alter the transcription initiation rate offers a further approach to tuning the output (Figure 3E). Externally tuning the activation of the synthetic promoters with chemical inducers was also made possible by fusing the PRD to ligand-inducible DNA binding domains (DBDs) (Figures 3F–3K). Finally, we validated that the pheromone response pathway and the downstream mating response were now decoupled by this design, by demonstrating that native Ste12-regulated genes were no longer transcriptionally activated by pathway activity (Figure S4C).

14Together, our mathematical and in vivo models reveal that although GPCR dose-response is the output of a complex system, variability can be achieved by altering the promoter identify for just three components (receptor, Gα, and reporter), offering a simple approach to rationally tune the sensitivity, basal activity, and signal output of a GPCR signaling pathway. To demonstrate this in practice, we next set out to optimize the α-factor response of our minimized response pathway through iterative refactoring of these components. Our starting strain, with constitutive expression of the three components set at native levels (Figures 3B and 4A, design 1), performed poorly compared to the Quasi-WT strain, which we hypothesized was due to removal of regulated expression of pathway components. Indeed, Ste2, Gpa1, and Ste12 all have Ste12-activated promoters and so benefit from positive feedback in their native setting. Guided by our model, we were able to restore performance by engineered reintroduction of Ste12-mediated feedback (Figures S4D–S4G), however, we did not move forward with this design because it would couple sensitivity to pathway output.

15Instead, to improve on the performance of the design 1 strain, we increased Gpa1 levels using the PGK1 promoter to reduce the basal activity of the response to effectively zero (Figures 3B and 4A, design 2). The dose-response sensitivity could then also be boosted by increasing Ste2 expression via the strongest promoter available (pCCW12). This version (Figures 4A and 4B, design 3) now approached the sensitivity of the Quasi-WT strain for the α-factor inducer and the physical limits of the receptor (Figure S4H). Finally, the output strength of the response could be maximized by linking the pathway to activate the best performing synthetic promoter, LexO(6x)-pLEU2m, via the LexA-PRD sTF (Figures 4A and 4B, design 4). The resulting pathway was highly optimized in comparison to the Quasi-WT strain with improved operational range, tightness, dynamic range, and maximum output. Notably, this is achieved without the need for feedback regulation of signaling components and also offers a pathway decoupled from the >100 genes usually upregulated in the mating response. Furthermore, the engineering we had performed to the design 4 strain had come with no penalty to the growth rate when compared to wild-type or Quasi-WT yeast (Figures S4I and S4J).

16Understanding how the pathway dose-response can be shifted in this manner advances our basic knowledge of how component level changes effect signal transduction. Alongside this, it also offers direct applications for synthetic biology, where reprogramming cells to receive specific signals and respond in a desired manner is a core goal. Indeed, GPCRs represent the ideal sensory module for eukaryotic synthetic biology because they are responsive to a plethora of ligands and stimuli, often operate with high specificity, and naturally have modularity written into their signaling architecture.

17With this in mind, we set out to establish our model strain as a host cell for rationally engineering yeast as sensors that detect diverse inputs via heterologous GPCRs (see Figure S5 for a description of the final toolkit). As an initial demonstration we took the Mam2 receptor from Schizosaccharomyces pombe, which detects a 23-amino acid peptide pheromone called P-factor. Using the optimized tuning levels determined for design 4, we generated a P-factor-sensing strain in less than a week. We then compared the response of this sensor strain to the response observed in its native context, as reported by (Figure 4D). The Mam2 sensor strain behaved almost exactly as in S. pombe, achieving an identical potency (EC50) to P-factor. Furthermore, the Mam2 sensor strain displayed no detectable basal activity and exhibited a 180-fold change in signal after activation, suggesting the optimization we had performed would also be suitable for other GPCRs.

18Efforts to create sensor strains over the last two decades have coupled heterologous GPCRs to the yeast pheromone response pathway with varying success. Now with the approach described here, the sensitivity, basal activity, and response output of cells sensing via heterologous GPCRs can be rationally tuned. However, one further important characteristic of a sensor—its operational range (the Hill slope of its dose-response curve)—is more difficult to adjust, as it is determined largely by the ligand-binding properties of the receptor. Some receptors will confer a narrow range switch-like behavior, only requiring a small increase in signal to trigger maximum output (i.e., a digital response), whereas others will give a wide operational range where there is a proportional relationship between signal and output (i.e., a linear response). For sensor applications, a linear response is typically required, whereas the digital response is more desirable for point-of-care and gene circuit applications.

19Given these considerations, we next set-out to solve how the operational range can also be tuned using an engineering approach, so that the Hill slope can be reduced to expand an operational range, or the Hill slope can be increased to narrow an operational range. For this we first built two new sensor strains, both sensing medically relevant metabolites by having human GPCRs coupled to our refactored yeast pathway. The chosen receptors were the adenosine-responsive A2BR receptor, previously shown to give a digital-like response in yeast, and the melatonin-responsive MTNR1A receptor, previously shown to give a linear-like response in yeast (Figures S6A–S6D).

20While previous efforts have tuned the yeast pheromone response Hill slope by overlaying synthetic feedback loops into the MAPK cascade, feedback loops on the components tuned in our approach proved ineffective in changing the dose-response curves (Figures S6E–S6L), likely because our orthogonal sTFs do not have the autoregulatory feedback that Ste12 has via its native promoter. Without this avenue, we instead choose a different tactic, tuning the Hill slope by creating engineered communities of cells that change the average response at the population level.

21First, to create a population that linearizes the steep response curve of our adenosine-sensing cells, we took inspiration from a strategy employed by previous artificial sensor systems, where receptors with different sensitivities are combined and their average response determines the output. We used rational tuning of GPCR levels to create two new strains with increased and decreased sensitivity to adenosine, and then tuned the output promoters so that their maximum outputs match (Figure 5A), using total fluorescence of each culture to measure GFP output (Figures S7A and S7B). We then co-cultured the three-sensing strains in a 1:1:1 ratio to create a consortium whose average response integrates the signal from all cells to give an extended operational range. This almost halved the Hill slope of the response while maintaining a similar potency, yielding an operational range 50-fold greater than the initial response (Figures 5B and 5C).

22Narrowing the operational range of the melatonin-responsive MTNR1A sensor strains required more complex engineering as the Hill slope of a response can only be increased via mechanisms such as cooperativity, sequestration, or positive feedback. As before, we utilized a community-based approach, but here using cell-to-cell communication to enable feedback at the population level. A two-cell system was designed where the first cell acts as an amplifier, sensing via MTNR1A and responding by secreting α-factor from the reintroduced gene. The second cell senses α-factor and responds with reporter gene (sfGFP) expression, and also secretes constitutive levels of the α-factor degrading protease, Bar1, to create a threshold for activation (Figures 5D and 5E).

23Computational modeling revealed that fine-tuning Bar1 levels was important for a generating a steep Hill slope while maintaining tightness and high dynamic range, albeit with some loss of sensitivity (Figure S7C). To determine the optimal expression of Bar1, we tested a range of promoters and measured the dose-response of the two-cell system. The RPL18B promoter emerged as ideal choice for driving Bar1 expression as it eliminated basal activity in the two-cell amplifier-reporter system while also maintaining high levels of maximum signal (Figure S7D).

24Following the creation of the α-factor-detecting reporter cells and the tuning of Bar1 levels secreted by these, we engineered the digital response in two steps using our toolkit (Figures S7F–S7H). Amplifier cells were first created by linking the output of MTNR1A sensor strain to the production of peptide pheromone, α-factor. The sensitivity of these cells was then tuned by adjusting MTNR1A receptor expression in the amplifier strain, so that the potency of the response in the two-cell system matched that of the single-cell version. When co-cultured in a 1:1 ratio, the final tuned two-cell system maintained the same potency but now provided a dose-response curve with a 2.3-fold increase in Hill slope and more than 200-fold decrease in the operational range (Figures 5F and 5G).

25The ability to rationally tune all dose-response properties for GPCR-signaling enables yeast strains to be optimized as sensors appropriate for different types of applications. To demonstrate this, we tackled two different example cases: metabolite quantification, where specificity and a wide dynamic range are desired, and point-of-care pathogen detection, where a sensitive digital-like response is needed.

26recently reported production of melatonin from S. cerevisiae by constructing a biosynthetic pathway that converts L-tryptophan into melatonin via three non-native intermediates. We thus sought to create sensor strains optimized for measuring melatonin production from these engineered yeast, exploiting the fact that our MTNR1A sensor strains demonstrate exquisite specificity for melatonin over its precursors (Figure 6A).

27To enable quantification of melatonin from our strains we engineered a two-strain consortia to widen the operational range to match the range of melatonin concentrations reported by from their production strains (Figure 6B). The resulting curve was linear over four orders of magnitude. We harvested supernatant from 88 different melatonin producer strains from the previous study, and in parallel, measured melatonin levels with our consortia of sensor strains and by liquid chromatography-mass spectrometry (LC-MS) (Figures 6C and 6D). The calculated concentrations of melatonin determined by the two methods agreed strongly, with the yeast sensor arguably more accurate over the large span of concentrations, due to the known limitations of the linear range in LC-MS measurement.

28For pathogen detection, we next set out to tune a point-of-care sensor for Paracoccidioides brasiliensis, recently developed by. This yeast GPCR-based sensor, utilizing the Ste2 homolog from P. brasiliensis (PbSte2), detects environmental levels of pheromone peptide released by haploid P. brasiliensis cells (PbPeptide) as a means of specifically detecting this human pathogen in complex samples (Figure 6E). This receptor exhibits a wide operational range with nanomolar sensitivity in S. cerevisiae, but a digital “yes/no” response would be more appropriate for point-of-care use, as would greater sensitivity. To optimize these two properties, we applied our amplifier-receiver two strain approach to the PbSte2 sensor (Figures 6F–6H). Increasing receptor levels with the CCW12 promoter and incorporating signal amplification created a consortia that now detected in the picomolar range, reporting with a digital-like response due to an operational range now narrowed by 1,000-fold.

29In this study, we used genome engineering and synthetic biology tools to refactor a minimized GPCR signaling pathway so that it could be rationally tuned both in silico and in vivo. This revealed three key design principles for tuning the signaling pathway dose-response curve: (1) sensitivity is increased by increasing GPCR numbers, (2) basal activity is reduced by finding the G protein expression optima, and (3) amplitude (pathway output) is tunable via synthetic transcription factors and engineered promoters. With these principles we were able to engineer yeast strains for desired performance as sensors for peptide inducers, and for primary and secondary metabolites.

30Uncovering these principles was achieved here by the close connection between in silico modeling and experimental validation. This was made possible by refactoring yeast into a strain that effectively acts as an in vivo model, with component expression levels that can be individually varied and performance simple to measure via fluorescent output. The in silico model allowed us to first understand where component level tuning has the greatest impact on the signaling response, with experimental work then confirming this and identifying the promoters that will achieve these levels in subsequent engineering.

31As well as confirming past observations on signaling stoichiometry, the model also revealed a complex, non-monotonic relationship between G protein subunit levels and pathway activity. An unexpected decrease in the response maximum output was predicted by the model for low Gα or high Gβγ levels, with this confirmed by experiments. By modeling different mechanisms, we were able to attribute this reduced activity to incorrect trafficking of Gβγ in the absence of free Gα. This insight told us that fine-tuned expression of G protein subunits would be essential for signal transduction that has minimal basal activity but maintains a high output when activated.

32Varying signaling component stoichiometries in our in silico and in vivo work also demonstrated that the GPCR-mediated response is not merely defined by the receptor’s intrinsic properties (i.e., ligand affinity); instead, it is a function of the properties of all components in the signaling pathway and particularly their relative levels. This fact has important consequences. It explains how signal transduction behavior could be significantly altered by a change in component levels, whether due to a change in environmental conditions or due to altered expression and protein turnover, for example, in different tissues. Indeed, via this mechanism cells can have different sensitivities and activation thresholds for the same agonist while expressing identical receptors. Importantly, this fact also underlines why non-coding genetic variation, such as mutation in promoter regions, has to also be considered alongside protein polymorphisms when assessing how genetic variation links to health and to the efficacy of treatments. Already, receptor variation in humans is recognized as a major cause of GPCR-targeting drugs being ineffective in many individuals, and it is possible that non-coding mutations that alter pathway stoichiometries could further explain such cases.

33We anticipate that the tuning principles uncovered here in yeast will also be relevant for GPCR signaling in all eukaryotes, however, it is worth re-stressing the large diversity in the type and structure of downstream signaling pathways paired with GPCRs in different organisms and cell types. The next steps for our approach will therefore be to use equivalent tools to refactor a canonical mammalian GPCR pathway so that its components can be tuned and assessed in isolation to the point where the dose-response to the agonist can be set as desired. This would also accelerate applications in pharmacology and healthcare that utilize GPCRs, such as in cell-based theranostics where cells are engineered to detect and act upon receiving defined cues within the human body.

34While mammalian systems remain the go-to choice for studying GPCR activity, the genetic amenability, ease of use, and low cost make yeast an ideal organism for many sensor applications. As demonstrated here, our model strain offers immediate applications for engineering yeast to sense its environment, whether as single strains or in engineered consortia. Already efforts in synthetic biology have used engineered yeast sensors as medical diagnostics, for pathogen detection, and as a tool for accelerating metabolic engineering. In all these applications, it is desirable for the user to have control over the response to input and the magnitude of gene expression it triggers. Full ability to tune sensors, as shown here, allows engineering for desired detection windows, and could be used in further work to define thresholds for activation (e.g., for directed evolution) or for matching the input-output levels when cells are engineered to detect and act or to communicate in connected systems.

35A current limitation of using yeast as sensors is that most medically relevant GPCRs do not port directly into S. cerevisiae without requiring optimization of expression, membrane translocation, and pathway coupling. Co-expression of mammalian accessory proteins and the humanization of the yeast membrane have also been shown to improve porting of receptors from mammalian species, and at the last count at least 50 different foreign GPCRs have been demonstrated to function in yeast (Table S1). These past successes provide a guide for those looking to generate yeast sensors for ligands and metabolites. However, without further experimental testing there are no guarantees that they will behave as required. In particular, the specificity of the receptor for the intended ligand is a major consideration as many receptors are promiscuous in what they bind. Another consideration is whether the ligand of interest can easily access a GPCR on the yeast cell membrane, especially given that yeast has a primitive cell wall. Fortunately for extracellular metabolite sensing, the porous structure of the cell wall is less of a concern, because it allows the free passage of molecules as large as 800 kDa.

36Our model strain now offers a new starting point for the many applications of GPCR-sensing enabling further work in a systematic, plug-and-play manner. Our overall strategy of simplifying and refactoring a natural pathway to first understand it and then rationally expand it should also be applicable to other systems in and beyond signal transduction. With the accelerating capabilities of genome engineering and synthetic biology in all organisms, it is likely that we will see the creation of equivalent in vivo model strains to rationally explore and exploit the key features and parameters of other important biological systems.

37Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Tom Ellis (t.ellis@imperial.ac.uk).

38NEB® Turbo Competent E. coli was used for all cloning experiments. Selection and growth of E. coli was performed in Lysogeny Broth (LB) medium at 37°C with aeration. With the exception of generating competent cells, the LB medium was supplemented with appropriate antibiotics (ampicillin 100 μg/mL, chloramphenicol 34 μg/mL, or kanamycin 50 μg/mL).

39For a list of all yeast strains used in this study see Key Resources Table. Apart from BY4741, all strains used in this study are a derivative of the yWS677 model strain (sst2Δ0 far1Δ0 bar1Δ0 ste2Δ0 ste12Δ0 gpa1Δ0 ste3Δ0 mf(alpha)1Δ0 mf(alpha)2Δ0 mfa1Δ0 mfa2Δ0 gpr1Δ0 gpa2Δ0), which is itself a derivative of BY4741 (MATɑ his3Δ1 leu2Δ0 met15Δ0 ura3Δ0). The yWS677 model strain was generated using iterative markerless CRISPR/Cas9 genome engineering (see METHOD DETAILS). Further editing of yWS677 to create the additional strains was performed in a single step using markerless CRISPR/Cas9 genome engineering.

40Yeast extract peptone dextrose (YPD) was used for culturing cells in preparation for transformation: 1% (w/v) Bacto Yeast Extract (Merck), 2% (w/v) Bacto Peptone (Merck), 2% glucose (VWR). Cells were cultured at 30°C shaking at 250 rpm.

41Selection of yeast transformants was performed on synthetic complete (SC) dropout agar medium: 2% (w/v) glucose (VWR), 0.67% (w/v) Yeast Nitrogen Base without amino acids (Sigma), 0.14% (w/v) Yeast Synthetic Drop-out Medium Supplements without histidine, leucine, tryptophan, and uracil (Sigma) supplemented with 20 mg/L tryptophan (Sigma), and 20 g/L bacteriological agar (VWR). Depending on the required selection, SC dropout media was supplemented with 20 mg/L uracil (Sigma), 100 mg/L leucine (Sigma), and 20 mg/L histidine (Sigma). Cells were grown at 30°C static.

42All liquid experiments were performed in synthetic complete (SC) medium with 2% (w/v) glucose (VWR), 0.67% (w/v) Yeast Nitrogen Base without amino acids (Sigma), 0.14% (w/v) Yeast Synthetic Drop-out Medium Supplements without histidine, leucine, tryptophan, and uracil (Sigma), 20 mg/L uracil (Sigma), 100 mg/L leucine (Sigma), 20 mg/L histidine (Sigma), and 20 mg/mL tryptophan (Sigma). Unless otherwise stated, all yeast strains were cultured in 500 μL of SC medium and grown in 2.2 mL 96 deep-well plates at 30°C in an Infors HT Multitron, shaking at 700 rpm.

43Chemically competent cells were created following the TSS protocol for KCM transformations. A colony of E. coli was grown to saturation overnight in 10 mL of LB and then split into two 2 L baffled flasks with 500 mL of LB. The culture was grown for 2-3 h to OD600 ∼1.0, chilled on ice to stop growth, split between 50 mL conical tubes, and centrifuged at 4000 rpm at 4°C for 10 minutes. The supernatant was then discarded, and the cell pellets resuspended by aspiration in ice-cold TSS (85 mL LB, 10 g PEG-3350, 5 mL DMSO, and 2 mL 1 M MgCl2). 200 μL of the cell suspension was then aliquoted into PCR reaction tubes, flash frozen on dry ice, and put into a −80°C freezer for long term storage. To transform the DNA, 50 μL of 5 x KCM (500 mM KCl, 150 mM CaCl2, 250 mM MgCl2) was added to 200 μL of the competent cell prep after 10 minutes of thawing on ice. 50 μL of the competent cell-KCM cocktail was then added to DNA and transferred to a thermocycler with the following protocol: 4°C for 10 minutes, 42°C for 1 minute, 4°C for 1 minute, and then 30-60 minutes recovery at 37°C. Cells were then plated on solid LB medium supplemented with the appropriate antibiotics.

44Chemically competent yeast cells were created following the lithium acetate protocol. Yeast colonies were grown to saturation overnight in YPD. The following morning the cells were diluted 1:100 in 15 mL of fresh YPD in a 50 mL conical tube and grown for 4-6 h to OD600 0.8-1.0. Cells were pelleted and washed once with 10 mL 0.1 M lithium acetate (LiOAc) (Sigma). Cells were then resuspended in 0.1 M LiOAc to a total volume of 100 μL/transformation. 100 μL of cell suspension was then distributed into 1.5 mL reaction tubes and pelleted. Cells were resuspended in 64 μL of DNA/salmon sperm DNA mixture (10 μL of boiled salmon sperm DNA (Invitrogen) + DNA + ddH2O), and then mixed with 294 μL of PEG/LiOAc mixture (260 μL 50% (w/v) PEG-3350 (Sigma) + 36 μL 1 M LiOAc). The yeast transformation mixture was then heat-shocked at 42°C for 40 mins, pelleted, resuspended in 200 μL 5 mM CaCl2 and plated onto the appropriate synthetic dropout medium.

45All genomic edits were performed via CRISPR/Cas9, using a two-plasmid system, consisting of a gRNA expression plasmid (pWS082) and a CRISPR/Cas9 expression plasmid with the choice of three different selection markers (pWS158, URA3; pWS171, LEU2; pWS172, HIS3). This system supplies Cas9 and gRNAs on two different plasmids, which are first linearized and then gap repair with each other in yeast to create single plasmid which contains both components. For multiplexing edits, multiple gRNA fragments can be introduced into yeast simultaneously. Individual gRNAs were expressed using the format used in the YTK system, whereby a tRNA promoter drives the expression of a tRNA and HDV ribozyme fused to the 5′ of the gRNA followed by the SNR52 terminator. Cas9 is then expressed using the PGK1 promoter and terminator. To prepare the two plasmids for transformation, the CRISPR/Cas9 and gRNA expression plasmids were first digested with BsmBI or EcoRV, respectively, and the following size fragments were gel purified: pWS082, 1022 bp; pWS158, 10051 bp; pWS171, 10909 bp; pWS172, 10102 bp.

46New gRNA targets were designed in Benchling (http://www.benchling.com) using the CRISPR tool. gRNA sequences were then created using two annealed 26 bp oligonucleotides and assembled into the gRNA expression vector using a BsmbI Golden Gate assembly using the (GACT) overhang at the 5′ and (GTTT) at the 3′. An additional (TT) was included between the 5′ overhang and gRNA sequence to complete the HDV ribozyme sequence:5′ GACTTTNNNNNNNNNNNNNNNNNNNN 3′3′ AANNNNNNNNNNNNNNNNNNNNCAAA 5′

475′ GACTTTNNNNNNNNNNNNNNNNNNNN 3′

483′ AANNNNNNNNNNNNNNNNNNNNCAAA 5′

49Oligonucleotides were assembled into the pWS082 gRNA entry vector using the following protocol: Oligonucleotides were first resuspended at 100 μM in H2O. Each oligonucleotide was then treated separately with T4 polynucleotide kinase (PNK) in the following reaction: 1 μL oligonucleotide (100 μM), 1 μL 10 x T4 DNA ligase buffer (Promega), 0.5 μL T4 PNK (NEB), and 7.5 μL H2O. The mixture was then incubated at 37°C for 1 hour. The 10 μL reactions for both oligonucleotides in the fragment pair were then added together and brought to a total volume of 200 μL in H2O (10 μL oligo (sense) + 10 μL oligo (antisense) + 180 μL H2O). The oligonucleotides were then annealed under slowly decreasing temperatures using the following program: 96°C for 6 minutes followed by 0.1°C/s ramp down to 20°C, and then hold at 20°C. The resulting fragment was then ligated into the gRNA expression vector using Golden Gate assembly, which was prepared as follows: 0.1 μL of pWS082 (50 fmol/ μL), 1 μL of the small fragment, 1 μL T4 DNA ligase buffer (Promega), 0.5 μL T7 DNA Ligase (NEB), 0.5 μL BsmBI (NEB), and water to bring the final volume to 10 μL. Reaction mixtures were then incubated in a thermocycler using the following program: (42°C for 2 min, 16°C for 5 min) x 10 cycles, followed by a final digestion step of 60°C for 10 min, and then heat inactivation at 80°C for 10 min. The entire reaction mixture was then transformed directly into E. coli and plated on LB medium with ampicillin (100 μg/mL). For a list of oligonucleotides used to create the gRNA sequences in this study, see Table S2.

50To edit the targeted regions, donor DNA was introduced into yeast alongside the CRISPR DNA to facilitate homology-directed repair at the double-strand break. Donor DNA was created by first cloning the sequence into the pYTK001 part entry vector. Donor DNA consisted of 500 bp arms of homology flanking a unique 24 bp sequence containing a new CRISPR/Cas9 target (landing pad), which were designed using the Benchling CRISPR tool to have an off-target score of 100 and a high on-target score. The plasmid was then sequence verified and primers were designed to amplify the donor DNA from the plasmid. The 1024 bp PCR amplicon was then gel purified to generate the donor DNA for transformation. For more information on donor DNA design see Figure S1A. For a list of landing pad sequences at the edited loci in the yWS677 and yWS1922 model strains, see below.

51Reagents for CRIPSR markerless genome editing were prepared as follows: 50 ng of linear CRISPR/Cas9 plasmid, 500 ng of each linearized gRNA expression plasmid, and 1000 ng of each donor DNA. The DNA was then transformed directly into yeast and plated on the appropriate selective media. Edits were validated by colony PCR followed by sanger sequencing of the amplified genomic region (for a list of primers used in this study, see Table S2). To rapidly iterate between successive edits when generating the model strain, yWS677, a marker cycling protocol was used, where two deletions or changes were performed per round of editing. Plasmid curing was skipped, and the marker for selecting the CRISPR/Cas9 plasmid was cycled between 3 markers (URA3, LEU2, and HIS3) for each iteration. The final edit was performed with CRISPR/Cas9 plasmid containing the URA3 selection marker, which was counter selected using 5-FoA to cure the yeast of CRISPR machinery (Figures S1B and S1C). The absence of all CRISPR/Cas9 plasmids was validated by colony PCR and replica plating on selective media.

52DNA was isolated from yWS677 for Nanopore sequencing using the 100/G Genomic-Tip kit (QIAGEN), sheared to 20 kb using a g-TUBE (Covaris) and prepared for sequencing using a Ligation Sequencing Kit 1D2 R9.5 (Oxford Nanopore Technologies). The genomic DNA was then run on an R9.5 flow cell using a MinION Mk 1B (Oxford Nanopore Technologies). A standard 48h sequencing run was performed using the MinKnow 1.5.5 software using local basecalling. Reads were exported directly to fastq using MinKNOW. Canu (v1.5) was used to correct raw reads (http://www.canu.readthedocs.io) and smartdenovo (http://www.github.com/ruanjue/smartdenovo) was used to de novo assemble the reads into contiguous sequences (contigs) using default flags. Resulting contigs were compared to a WT reference genome (S288C, SGD) using lastdb/lastal (http://www.last.cbrc.jp) and viewed on integrative genome viewer (IGV) (http://www.software.broadinstitute.org) to inspect genomic changes. We then probed all discrepancies between the yWS677 genome and the S288C reference using a minimum alignment length of 100 bp.

53Single colonies were grown to saturation overnight in 3 mL YPD. The next day, the yeast cultures were back diluted to an OD600 of 0.175, and 200 μL was transferred to a 96-well clear, flat-bottom microplate (Corning). OD600 was the measured over 24 h by a Synergy HT Microplate Reader (BioTek) taking measurements every 15 minutes with shaking at 30°C in between readings. Growth rate per hour was calculated according to Equation 1, where t is time in hours.

54All plasmids within this study were created using the MoClo Yeast Toolkit (YTK) system. Additional sequences not included within the YTK system that were used within this study can be found in the Table S3. For a list of all plasmid constructs used in this study, see Table S4. Unless indicated, part sequences were either mutated or synthesized to remove or avoid all instances of the BsmBI, BsaI, BpiI, and NotI recognition sequences.

55Construction of all plasmid constructs in Table S4 was achieved using Golden Gate assembly. All parts were set to equimolar concentrations of 50 fmol/μL (50 nM) prior to experiments. Golden Gate reactions were prepared as follows: 0.1 μL of backbone vector, 0.5 μL of each plasmid, 1 μL T4 DNA ligase buffer (Promega), 0.5 μL T7 DNA Ligase (NEB), 0.5 μL restriction enzyme (BsaI or BsmBI) (NEB), and water to bring the final volume to 10 μL. Reaction mixtures were then incubated in a thermocycler using the following program: (42°C for 2 min, 16°C for 5 min) x 25 cycles, followed by a final digestion step of 60°C for 10 min, and then heat inactivation at 80°C for 10 min. The entire reaction mixture was then transformed directly into E. coli and plated on LB medium with the appropriate antibiotics.

56G protein C-terminal variants were created by substituting the GFP dropout cassette in the pWS936 Gpa1 C-terminal truncation vector for a small DNA fragment consisting of two annealed oligonucleotides. To create the small DNA fragment, oligonucleotides were first resuspended at 100 μM concentration in H2O. Each oligonucleotide was then treated separately with T4 polynucleotide kinase (PNK) in the following reaction: 1 μL oligonucleotide (100 μM), 1 μL 10 x T4 DNA ligase buffer (Promega), 0.5 μL T4 PNK (NEB), and 7.5 μL H2O. The mixture was then incubated at 37°C for 1 hour. The 10 μL reactions for both oligonucleotides in the fragment pair were then added together and brought to a total volume of 200 μL in H2O (10 μL oligo (sense) + 10 μL oligo (antisense) + 180 μL H2O). The oligonucleotides were then annealed under slowly decreasing temperatures using the following program: 96°C for 6 minutes followed by 0.1°C/s ramp down to 20°C, and then hold at 20°C. The resulting fragment was then ligated into the Gpa1 C-terminal truncation vector using Golden Gate assembly, which was prepared as follows: 0.1 μL of pWS936 (50 fmol/ μL), 1 μL of the small fragment, 1 μL T4 DNA ligase buffer (Promega), 0.5 μL T7 DNA Ligase (NEB), 0.5 μL BsmBI (NEB), and water to bring the final volume to 10 μL. Reaction mixtures were then incubated in a thermocycler using the following program: (42°C for 2 min, 16°C for 5 min) x 10 cycles, followed by a final digestion step of 60°C for 10 min, and then heat inactivation at 80°C for 10 min. The entire reaction mixture was then transformed directly into E. coli and plated on LB medium with chloramphenicol (34 μg/mL) For a list of oligonucleotides used to create the G protein library, see Table S2.

57The yWS677 model strain was prepared for multiplex integration of selectable plasmids by integrating landing pads (LPs) at the URA3, LEU2, and HO loci, conforming to the YTK integration plasmid format. Transient expression of Cas9 and the gRNAs targeting the landing pads was achieved by individually assembling Cas9 (pTDH3-Cas9-tTDH1) and the gRNAs (created using the YTK050 sgRNA dropout) into cassette plasmids without a yeast marker or yeast replicon. Double and triple integration of marker plasmids were performed with 100 and 200 ng of plasmid, respectively, with 100 ng of the Cas9 expression cassette and 200 ng of each gRNA expression cassette. All plasmids were first linearized by digestion before transformation using NotI-HF (NEB). Successful plasmid integration was selected for using synthetic drop-out media missing the appropriate supplements. Cas9 and gRNA expression was transient and quickly lost due to lack of selection or flanking homology to genome. Initially validated by colony PCR and then assumed thereafter.

58All quantitative PCR (qPCR) reactions were performed in an MasterCycler ep realplex 4 (Eppendorf) using SYBR FAST Universal qPCR Master Mix (Kapa Biosystems). For RNA purification, RNA was isolated from yeast culture grown to an OD600 of 1 ± 0.1 using a YeaStar RNA Kit (Zymo Research). RNA was quantified by nanodrop spectrophotometer (Thermo Fisher) and cDNA was generated from each RNA prep using a High Capacity cDNA Reverse Transcription Kit (Applied Biosystems). Each qPCR reaction contained 20 ng of cDNA. qPCR results were normalized to the housekeeping gene HTB2. All qPCR primers were designed manually using Benchling.

59All sensor strains were picked into 500 μL of synthetic complete media and grown in 2.2 mL 96 deep-well plates at 30°C in an Infors HT Multitron, shaking at 700 rpm overnight. The next day, saturated strains were then diluted 1:100 into fresh media. After 2 h of incubation the strains were induced with their respective ligands and incubated for a further 4 h. All ligands were dissolved in DMSO, and the final concentration of DMSO in all cultures was 1 %. For strains using the TetR-PRD or Z3E-PRD transcription factor, aTc and β-estradiol was added during the back dilution at time 0 h. To perform flow cytometry and plate reader measurements, 200 μL from each well was directly transferred to a 96-well clear, flat-bottom microplate (Corning).

60For monoclonal cell experiments, cell fluorescence was measured by an Attune NxT Flow Cytometer (Thermo Scientific) with the following settings for measuring sfGFP: FSC 300 V, SSC 350 V, BL1 500 V. Fluorescence data was collected from 10,000 cells for each experiment and analyzed using FlowJo software. For polyclonal cell experiments, cell fluorescence was measured by a Synergy HT Microplate Reader (BioTek) with the following settings for measuring sfGFP: excitation 485/20, emission 528/20, gain 80. Unnormalized, raw fluorescence readings from each well were used for data analysis.

61Samples for mass spectrometry and the MTNR1A sensor were prepared by centrifuging yeast cultures at 4000 rpm in a large desktop centrifuge for 10 minutes at 4°C and extracting the supernatant. Supernatant samples were kept on ice before running on the LC-MS or transferred to the −20°C for later use. No further sample preparation was performed on the supernatant sample before running on the LC-MS. Melatonin standards were kept in 100% DMSO before being diluted in spent media. Spent media was prepared from BY4741 in the same manner as the measured yeast.

62An LC-MS method was developed for the measurement of melatonin in media, using an Agilent 1290 UPLC and 6550 quadrupole – time-of-flight (Q-ToF) mass spectrometer with electrospray ionization (Santa Clara, CA). The UPLC column was an Agilent Zorbax Eclipse Plus C-18, 2.1 × 50mm and 1.8um particle size. The UPLC buffers were 0.1% formic acid in water and 0.1% formic acid in acetonitrile (v/v). The gradient elution method is detailed below.

63Quantitation was based on the LC retention time from melatonin standard solutions and the area of accurately measured diagnostic ions from the molecule, namely the protonated molecule, [M+H]+, along with an in-source fragment (see below). The solutions of a melatonin standard in media were used to generate calibration curves.

64To measure the melatonin from producer strains using the MTNR1A sensor strains, 50 μL of supernatant was added to an adjusted 450 μL volume of the sensor cells and run on the flow cytometer according to the ligand sensing protocol. Melatonin concentrations were then calculated from a standard curve.

65All presented dose-response fittings were generated in Prism 7 (GraphPad), which was used to determine the logEC50 and Hill slope. To determine all remaining properties of the dose-response curve, curve fitting was performed using Python (SciPy and Matplotlib) using the 4PL model (Equation 2), where x is the concentration, A is the minimum asymptote, B is the steepness, C is the inflection point and D is the maximum asymptote.

66To model biochemical reactions systems, it is customary to use a set of ordinary differential equations (ODEs) to describe changes in concentrations of biochemical species. The systems of equations will be defined as shown in Equation 3.Where in this case, are the concentrations of molecular species within the system, and are functions describing the molecular interactions.

67Numerical integration of the system of ODEs may allow us to derive the concentration of all molecular species at a certain time point. Simulations from numerical integration were performed using MATLAB R2017a offered by MathWorks.

68The aim to formulate a system of ODEs of the Ste2, Gpa1, and Ste4:Ste18 promoter library required a top-down approach via reproduction of similar trends against general experimental results. In terms of the Ste2 library, an increase in receptor concentration directly increases maximum response while retaining low constitutive signaling. An increase in [Gpa1] or [Gɑ] in the pathway should demonstrate an optimal response, where a specific concentration of [Gɑ] produces the maximum response while having constitutive signaling. Increases in [Ste4:Ste18] or [Gβγ] should also demonstrate an optimal response with some basal activity. Thus, it is important to design a detailed system of ODEs which differentiates the ligand-induced responses against the constitutive response.

69Moreover, since the refactored strain observe changes in response by modifying the receptor (Ste2), Gɑ subunit (Gpa1), and Gβγ subunit (Ste4:Ste18) concentration, a detailed model focusing on receptor/G protein interactions was used. Although there are models presenting full yeast pheromone ODE models in detail, these models are not be ideal for demonstrating the intricate changes in multiple constitutive activities and normal responses within the receptor/G protein interactions. Thus, we have chosen as our basis of the mathematical model our previously described cubic ternary complex model, which focuses upon the receptor and heterotrimeric G protein complex. In contrast to our previous work, our model output is the changes in free Gβγ subunits since this is how the S. cerevisiae pheromone-response pathway mediates downstream signaling. Upon release of the Gβγ dimer from Gɑ downstream signaling continues through interaction with Ste5 and recruitment of a MAPK cascade (Ste11, Ste7 and Fus3 respectively). It has been reported that the MAPK-cascade does not modify the Ste2/G protein dynamics and simply acts to amplify and transduce the response from the plasma membrane to the nucleus. As such, given our specific interest in the dynamics of the R/G protein and the fact that our refactored yeast cells are designed to be linearly correlated we have not included the MAPK in our model simulations with the exception of Ste5, which is considered as the downstream effector.

70Here, we formulate the ODE model dynamics, as we have described previously for refactored S. cerevisiae strain yWS1922 to recapitulate results obtained through varying expressions of Ste2 or Gpa1 through a promoter library (Figures 2D and 2E, Experimental). The basis of the mathematical model used is from a cubic ternary complex model, in which a receptor may be under its inactive conformation R or its active conformation R∗. Furthermore, a receptor may bind to G protein G, and activation of receptor may cause G to dissociate and undergo the G protein cycle. The model describes ligand binding, receptor activation, G protein binding and the G protein cycle. The scheme used for the reaction model is depicted in Methods S1A with detailed ODEs in Equation 4. The synthesis and degradation of the models were assumed to be slower compared to the signaling events and thus omitted for simplicity. Taking ligand concentration as a constant, our system consists of 8 receptor states and 8 non-receptor bound G protein states, with a total for the system of 16 ODEs. For model outputs, concentrations [Gβγ∗STE5] will be taken as subsequent downstream pathways are designed to be linearly correlated in the system. For the individual kinetic rate constants, the lowercase k was used, along with subscripts + and – to denote forward and backward reactions.

71Here we present the numerical results used to illustrate the Ste2 and Gpa1 promoter library model. Both concentration-response curves are created through MATLAB R2017a using the pre-installed ODE solver ode15s. ODEs were solved to obtain endpoint concentration of [Gβγ∗] after 1000 time points (Figures 2D and 2E, Model). Prior to the addition of ligand, the system was run for 1e8 time points to enable complex intermediates to obtain equilibrium. The simulated results demonstrate similar trends to the ‘wet’ experimental data (Figures 2D and 2E, Experimental – Ste2), where increases in Ste2 (R in the model) concentration shows an increase in maximum response. Conversely, for a fixed concentration of R increasing concentrations of Gpa1 (G in the model) demonstrates a non-monotonic relationship between G protein concentration and maximal signaling response (Figures 2D and 2E, Experimental – Gpa1). For fixed concentrations of R and Gɑ while increasing Ste4/Ste18 (Gβγ in the model), a non-monotonic relationship shows an optimal signaling response (Figures 2D and 2E, Experimental – Ste4:Ste18). The base parameter set, including initial species concentrations used for both simulations can be found below. Parameter values were qualitatively fit through experimental dose-response results of the Gpa1 and Ste2 promoter library, retaining source values when possible.

72Since the model allows G to dissociate in the absence of Ligand L or R∗ interaction, multiple levels of constitutive activity (elevated signaling in the absence of ligand) can be made possible. Along with the main pathway of G activation through LR∗ Gɑβγ, R∗Gɑβγ may also activate the G protein cycle in the absence of ligand. Gɑβγ can also dissociate and activate independently of R through rate constant kGRA-. To further reproduce constitutive features of the yeast pheromone cycle, inactive Gβγ may also be activated to be Gβγ∗ through kconact, which is a simplification of Gβγ interacting with downstream effectors. This detailed description of the constitutively active pathways allowed us to recapitulate the low constitutive signaling in the experimental results. To demonstrate how signaling decreases from an optimal point when increasing Gβγ, it has been anticipated that at very high concentrations, Gβγ is unable to be trafficked to the plasma membrane due to the lack of Gɑ, but still affecting the pathway by binding to downstream effectors prematurely. This has been described simply in the model by showing how Gβγ and Gβγ∗ compete for the same effector (Ste5). At low initial concentrations of Gβγ, little is available to bind to Ste5 so only the Gβγ that has been activated by interaction with the Ste2/Gpa1 complex (Gβγ∗) will actively bind to Ste5 (due to its higher affinity). However, when the initial concentrations of Gβγ are increased free Gβγ is able to compete with the Gβγ∗ and so sequesters the Ste5 species into an inactive complex so reducing the overall signaling maxima.

73It has been well documented that the intracellular concentrations both Ste2 and Gpa1 increase upon pheromone stimulation. This arises due to positive feedback (at the transcriptional level – Ste12) thereby increasing the amount of available Ste2/Gpa1 present at the shmooing tip. Since expression of all components in our refactored yeast cells are driven through constitutive promoters, not pheromone-inducible ones, we wondered what impact positive feedback would have upon the response output. We therefore constructed a model derived from the Ste2/Gpa1 model described previously (see Methods S1A) but incorporating reaction rates to enable the concentration of either/both Ste2/Gpa1 to increase upon production of free Gβγ dimer (our model output). We hypothesized that increasing both Ste2 and Gpa1 concentrations (R and G respectively in the mode) through a positive feedback loop would enable enhanced maximal signaling (increased EMax) production while retaining a low basal activity through suppression by Gpa1 (Note Gpa1 acts to sequester free Gβγ so blocking signaling).

74The formulation of the system of ODEs for the feedback model requires additional components to the previous system. Since the feedback loop is likely to involve multiple transcription/translation processes, our model required extension from the single ternary complex model to involve the whole refactored pheromone signaling pathway. In addition, long-term changes in the signaling components must also be considered. To retain a simplified model but that reflected dynamics at the Ste2/Gpa1 level but accounted for the MAPK cascade and gene expression we introduced two arbitrary molecular species (downstream of our readout free Gβγ): X1 and X2. Both can either be in an inactive or active state (denoted X or X∗ respectively). Once X2 is in its active state it acts as a promoter, directly increasing concentrations of Ste2, Gpa1, or both. The basic scheme of the system of ODEs is denoted in Methods S1B.

75The detailed system of ODEs can be found in Equation 5. In addition to the feedback loop, synthesis and degradation rates have also been considered in the parameter rates due to the longer timescale of the system. This relatively detailed system of ODEs consists of 8 receptor states, 5 non-receptor bound G protein states, and 4 downstream signaling states totaling to 17 ODEs. Model outputs would be concentrations [Gβγ].

76Here we present the numerical results to illustrate the Ste2 and Gpa1 feedback model. Final parameter values can be found below. ODE results shown are endpoint concentration of [Gβγ] after 1000 time points, with a 1e8 time point equilibration. The simulated results seem to show similar trends against the experimental data, with an increase in maximum while retaining a low basal activity to the system (Figures S4F and S4G).

77An application demonstrating the flexibility in the operational range of the yeast GPCR-based sensors was to use mixed populations of communicating strains to narrow the range, creating concentration-response curves close to a digital, ‘all or none’ response. The aim of this digital sensor model was to describe, in detail, the interactions between yeast cells in which their properties are modified so the product of one cell (ɑ-factor; Amplifier cell) becomes the agonist of another (Reporter cell). The Reporter cell expresses, under the constitutive RPL18B promoter, Bar1, a protease which degrades ɑ-factor. It has been suggested that the presence of Bar1 reduces ‘non-productive mating’ events in S. cerevisiae by rapidly degrading low levels of ɑ-factor. In essence, Bar1 works as a barrier between the ɑ-factor producing cell and the reporter cell, degrading the low levels of constitutively released ɑ-factor from the Amplifier cell. As concentrations of stimulating agonist (melatonin) increase, the proteolytic activity of Bar1 is saturated with increased ɑ-factor, allowing the activation of the reporter strain.

78Our previous modeling efforts have concentrated upon the dynamics of the R/Gɑβγ complex. However, to enable qualitative modeling of the digital sensor (effectively two entire pheromone signaling cascades) we decided to use a heavily reduced model for the R/Gɑβγ complex. We felt this important since our focus is on the interactions/behaviors of the two cells rather than solely the G protein cycle within a cell. Thus, we have chosen a reduced model which accurately describes overall in vivo responses for time- and dose-dependent effects in the system.

79To deviate from processes within the cell and focus on the interactions between cells, a reduced yeast pheromone pathway based on a model developed by was used. The structure of the reduced model consisting of 9 ODEs for each cell, modified to fit the experimental system through the addition of Bar1 interaction with ɑ-factor, along with a modification to set a maximum amount of possible product each cell can produce.

80The overall scheme used for the model is shown in Methods S1C. In this system, receptor states are simplified to two (R and R∗) and G protein states are simplified to three species (Goff, Gon, and GonEffector). The activated G protein will interact with an Effector protein (Ste5) which consequently activates a series of delay species to simulate the MAPK-like signaling cascade. prez3 is a term responsible for creating the final product of each cell by being directly converted to its corresponding product. It will act as a ‘cap’ by limiting maximum production within the concentration of prez3. Constitutive activity of each of the cell is achieved by directly adding GFP (the product generated from the reporter cell) to the system.

81The detailed system of ODEs can be found in Equation 6. Synthesis and degradation of the system are omitted for simplicity, apart from R∗ internalization. Our system includes the full pheromone response of two yeast cells, interactions between the two cells, and the constitutive protease activity of Bar1 consisting of 23 ODEs. Model outputs are either the product of the ɑ-factor producing cell (ɑfactor) or the product of the reporter cell (GFP), depending on the setup of the assay.

82260-minute endpoint readings of the refactored melatonin-sensitive yeast cell (MTNR1A sensor), Digital Feedback including the reporter cell without constitutive Bar1 activity, and Digital Feedback with Bar1 activity data were fitted simultaneously. Since most of the yeast pheromone system except the core MAPK cascade were refactored within the system, parameter values related to the receptor and its subsequent G protein activation rates were fitted using COPASI 4.16. Moreover, parameters and species introduced in the system for the first time were fitted against experimental data mentioned above. The set of finalized parameter values were obtained through parameter estimation algorithms provided in COPASI. To simulate and accurately fit incubation times in COPASI, constitutive fluorescence has been set as the initial concentrations of the corresponding products (3.1 nM for MTNR1A sensor, 10.4 nM for digital feedback without Bar1, and 4.4 nM for digital feedback with Bar1). Two optimization algorithms were tested: the Hooke and Jeeves algorithm and the evolutionary programming algorithm. Both methods reliably converged to the same minimum within the parameter space, which were thus considered as a global minimum.

83Here we present the parameter values and fitted graphs of the model, which can be found below. Ligand and Receptor specific parameters were refitted to the experimental time course data, as well as Bar1 interaction and initial concentration. Initial concentrations of Receptor, Goff, and Effector were set consistently at 16 μM for both cells.

84Final figures and fitting results created through the model are shown in Methods S1D. Fitting results demonstrate a lower EMax, higher Hill slope and increased potency upon introduction of Bar1, which was also observed in the ‘wet’ experimental data (Methods S1D–S1H). To provide an estimate of quality of the fitting of the model to the biological data we used residual analysis (Methods S1E), by which each residual is taken from the deviation of an observed value against the predicted value. Residual values are defined in Equation 7.Where a dataset has n values marked y1,…,yn (collectively defined as yi), the predicted values marked f1,…fn (collectively fi). All residual plots suggest an unbiased distribution with an equal deviation from 0. The residual plot digital feedback without Bar1 suggest a good fitting result with an unbiased distribution and a consistent variance throughout. Although heteroscedasticity may be observed in MTNR1A sensor and the digital feedback with Bar1, differences in variation within each ligand concentration suggests the natural tendency for variation to occur under the system with higher ligand concentration. Nevertheless, the residual plots suggest an unbiased distribution of experimental results. R2 values also suggest and give some information about the goodness of fit, defined in Equation 8.With as the mean of the dataset. Along with a R2 value of 97.93%, 99.90%, and 99.68% for the MTNR1A sensor, the digital feedback without Bar1, and digital feedback with Bar1 respectively, the non-linear regression suggests a good fit with the experimental data.

85Using the model, we simulate the effects of increasing concentrations of Bar1 (Figure S7C). At low concentrations of Bar1 constitutive signaling is decreased while obtaining similar or slightly lower concentrations close to the EMax. These biological responses suggest that the reporter cell has a ‘capped’ response where from a defined concentration, the reporter cell will not exceed its maximum capacity to produce GFP. This interpretation enables the reporter cell to experience ɑfactor saturation and even with a relatively small concentration of Bar1 it will not be able to degrade ɑfactor sufficiently at higher [L] (Figure S7D). With the introduction of Bar1, basal activity can be reduced while retaining the maximum possible response of the system.

86Ultrasensitivity, or the ability to create a digitalized response from a graduated input within a steady state response, is a well-studied field in systems biology and is said to be possible through various motifs such as multistep, feed-forward cooperativity, zero-order, and inhibitor ultrasensitivity. Within this ODE model, it has been designed to have two different mechanisms of ultrasensivity: through the amplification of the reporter cell and inhibition through Bar1. The MAPK cascade is a well known pathway where ultrasensitivity can arise from multistep-phosphorylation, positive feedback, and the cooperativity obtained through scaffolding protein, which is referred to as Ste5 in yeast. Bar1 interaction is also considered as a source of ultrasensitivity, where a competitive inhibitor will essentially compete against the Ste2 receptor.

87Statistical tests of all experiments was performed using GraphPad Prism version 7 or MathWorks MATLAB version R2017a and are detailed within the legend of each figure. For flow cytometry experiments, fluorescence data was collected from 10,000 cells for each datapoint and analyzed using FlowJo software, representing the average fluorescence of the population as the geometric mean. In all figures, the data points represent the mean ± SD. Curves fitted to all dose-response data was fitted in Prism 7 using the Nonlinear Regression: Variable slope (four parameter) curve fitting. One-way analysis of variance (ANOVA) was used to determine signifiance between the mean growth rates of the wild-type, model, and sensor yeast strains. Residual analysis was used to provide an estimate of quality of the fitting of the Cell-Cell Digital model to the biological data, by which each residual is taken from the deviation of an observed value against the predicted value.

88All part plasmids described in this paper are available at Addgene (Addgene ID: 123024-123065). Corrected nanopore sequencing reads for the yWS677 model strain are available at Sequence Read Archive (SRA) (SRA accession: PRJNA516326). MATLAB code for the mathematical models in this paper are available at BioModels Database (BioModels ID: Cubic Ternary Complex Model, MODEL1901300001; Cell-To-Cell Digital Sensor Model, MODEL1901300002).