A tunable population timer in multicellular consortia
AA tunable multicellular timer in bacterial consortia
Carlos Toscano-Ochoa ∗ and Jordi Garcia-Ojalvo † Department of Experimental and Health Sciences, Universitat Pompeu Fabra,Barcelona Biomedical Research Park, Dr. Aiguader 88, 08003 Barcelona, Spain
Processing time-dependent information requires cells to quantify the durations of past regulatoryevents and program the time span of future signals. Such timer mechanisms are difficult to implementat the level of single cells, however, due to saturation in molecular components and stochasticityin the limited intracellular space. Multicellular implementations, on the other hand, outsourcesome of the components of information-processing circuits to the extracellular space, and therebymight escape those constraints. Here we develop a theoretical framework, based on a trilinearcoordinate representation, to study the collective behavior of a three-strain bacterial populationunder stationary conditions. This framework reveals that distributing different processes (in ourcase the production, detection and degradation of a time-encoding signal) across distinct bacterialstrains enables the robust implementation of a multicellular timer. Our analysis also shows thecircuit to be easily tunable by varying the relative frequencies of the bacterial strains composing theconsortium.
Keywords: biological time encoding, polybacterial communities, chemical wiring problem, distributed com-puting, synthetic consortia, quorum sensing
I. INTRODUCTION
One of the defining features of cells is their ability toprocess information, by way of sensing internal and ex-ternal signals and generating a response. This paradigm,while serving as an inspiration for the development ofsynthetic biology [2], has also guided our search for geneand protein circuits responsible for all kinds of cellularfunctions [1]. Frequently, studies of cellular circuits havefocused on situations in which the input signal is con-stant throughout time. However, cells usually live intime-dependent environments [29], and thus they face theneed to process information that varies in time in a com-plex manner [11]. An important instance of temporalinformation processing by cells is the quantification ofthe time interval that a given input signal has been act-ing upon a cell. Complementarily, cells sending signalsmight need to program the amount of time that the signalis active. In technological settings, both these functionsare performed by timers. It is thus necessary to establishhow timers are implemented in cells.Cell-intrinsic molecular timers have been identifiedin recent years, having been reported to regulate awide variety of cellular processes, including apoptosis[10, 13, 14, 19], cellular proliferation [32], cell-fate speci-fication [6, 24], and infection response [7]. These timersusually depend on intricate molecular processes [18] thatare difficult to tune and are sensitive to noise. From asystems-level perspective, feedback-based timer circuitshave been proposed in synthetic biology applications[8, 22], and have also been found to operate naturallythrough pulses [17, 26] and oscillations [4], and even to beinduced by noise [31]. These intracellular timer circuits ∗ Email address (corresponding author): [email protected] † Email address: [email protected] usually depend on an accumulating signal surpassing athreshold [16], and are thus limited by the amount of sig-nal molecules that cells can produce and store in their in-terior. This limitation disappears if the factor controllingthe timer is exported outside the cell and stored in theextracellular medium. Cells in such multicellular timerswould operate collectively, and consequently their func-tion would be less affected by noise [9]. In this paperwe introduce a minimal implementation of a multicellu-lar timer, using a novel theoretical framework to describethe collective dynamics of microbial consortia.Microbial consortia are mixed populations of microor-ganisms, in particular bacteria, that exhibit division oflabor among different species or strains. Extensive effortshave been devoted in recent years to engineer syntheticconsortia [3, 20], pushing the development of distributedcomputation systems [25, 28, 30]. Under this approach,individual circuit components such as logic gates or smallnetworks motifs are implemented in distinct cell subpop-ulations, and connected between them by using extracel-lular signals to transfer information. Cell-cell commu-nication in this context is usually based on the quorumsensing systems of gram-negative bacteria [21]. Thesesystems commonly consist in an enzyme that producesAHL (acil-homoserin-lactone), a family of chemical com-pounds with the ability to diffuse through cellular mem-branes, bind to a transcription factor and thereby regu-late gene expression. In the distributed computing sce-nario, the sender and receiver circuit modules that pro-duce and sense AHL, respectively, are placed in differ-ent bacterial strains [25, 30]. Even though this sender-receiver framework allows to implement systems with theability to self-regulate and process information collec-tively, they lack a component that actively degrades thechemical signal in the medium. This implies that whenthe sender strain is turned off, the chemical substancethat transfers the information may disappear very slowlyif it is relatively stable. This limits the ability of the cir- a r X i v : . [ q - b i o . M N ] N ov cuit to process time-dependent information. Therefore,the consortium should include a strain that degrades thechemical signal [5, 23, 27].Our microbial consortium constitutes a minimal in-formation transmission system that we call a chemicalwire. The system is composed of three different bacte-rial strains (Fig. 1A): the first one produces a chemicalcompound that mediates the communication (emitter, E,green cells in the figure), the second detects the signal (re-ceptor, R, blue cells), and the third actively degrades it(sink, S, red cells). In what follows, we establish a theo-retical framework to describe this wire architecture, andprovide mathematical definitions of what we call the wirespace. This framework is built with the aim to provideboth a practical tool for bioengineers and a systems-leveldescription of multicellular information transmission andprocessing. Using this approach, we show that bacterialconsortia can work as a timers, by remembering how longthe sender strain has been in the ON state. In turn, thetime during which the receiver strain is active can be pro-grammed by tuning the duration of the input pulse. Weconsider that the population is not structured in spaceand is well mixed (corresponding to a liquid culture),and that it operates in the stationary phase (i.e. thecells do not proliferate). The latter is not a limiting as-sumption, since it has been shown that bacteria such as Escherichia coli are able to respond to external stimulieven in stationary phase [12]. This establishes a less re-strictive environment for the implementation of consortiacircuits connected through chemical wires, since the rel-ative grow rates of the different strains involved will notplay a relevant role in the operation of the system. Anearlier proposal of a multicellular timer involved a grow-ing population of cells, but as a result it required theuse of a single cell type, and the output was dominatedby noise [15]. In what follows, we present a theoreticalframework that allows us to analyze the behavior of ourthree-strain multicellular system quantitatively.
II. WIRE DEFINITION AND MODEL FORCULTURES OF NON-GROWING CELLS
As discussed above, we define a wire as an entity withthree elements, corresponding to three bacterial strains:an emitter half-wire (EHW), a receptor half-wire (RHW)and a sink (Fig. 1B). We call the communicating signalbbit, which stands for biological bit, or biobit. In ourcase the bbit will be the quorum sensing molecule AHL.A bbit can take any value within the [0,1] interval de-pending on the input signal, and as we will see it canalso store input duration. The EHW component is givenby the minimum set of genes required to synthesize bbitat a constant rate. The RHW component comprises theminimum set of genes required to induce gene expressionfrom a promoter when bbit is present in the medium at aconcentration above/below a specific threshold (depend-ing on whether AHL activates or represses its target gene, respectively). Finally, the sink component consists of theminimum set of genes, constitutively expressed, requiredto degrade bbit. The bacterial strains containing eachwire component will be called EHW, RHW and SINKstrains, respectively. Our first goal is to analyze howthese strains communicate with each other, and how therelative frequencies between them affect this communi-cation.We assume that the production of bbit, with concen-tration B , is linearly proportional to the density of theEHW strain. In contrast, the sink circuit will be sup-posed to degrade bbit following a Michaelis-Menten satu-ration dynamics, proportional to the SINK strain density(Fig. 1A): dBdt = K E E − D m BK S + B S (1)In this equation E and S represent the concentrationof EHW and SINK cells, respectively, while below we willalso introduce the concentration R of RHW cells. Theseconcentrations can be given in cells/mL or directly inoptical density units (which could be more convenientwhen designing an experimental implementation of thesystem). For simplicity, our wire will be conceived instationary phase, a condition that still allows cells to re-spond to external stimuli [12]. This means that cells arealive but are not dividing nor growing, and there is nodilution of bbit. We also assume that intrinsic bbit degra-dation is negligible in the time scales considered here.It is convenient to define the relative frequency of eachstrain in the population: ε = E/P , σ = S/P , and ρ = R/P , where P is the total population density (whichis constant provided the population is in the stationaryphase). By definition, ε + σ + ρ = 1 . We now establish abijection between any possible bacterial consortium anda point inside an equilateral triangle (Fig. 1C). In thisrepresentation, each vertex in the triangle correspondsto one of the three strains, and each of the relative fre-quencies ε , σ and ρ is given by the distance between thepoint representing the consortium and the side of the tri-angle opposite the corresponding vertex . This is calleda trilinear coordinate representation, and it will be ourgraphical framework to explore the behavior of our bac-terial consortium. The intuition behind this geometri-cal representation is simple: the closer the system is toa vertex, the more represented the corresponding strainwill be in the consortium. Since the population is in thestationary phase, none of the bacterial strains change indensity, and the point describing the system will be fixedinside the triangle. We define the set of all possible con-sortia as the wire space ∆ (the interior of the triangle). The triangle must have equal sides to obey the constraint thatthe sum of the distances from any point lying within it to itsthree sides is constant (equal to 1 in our case, as we have seenabove).
Figure 1.
A minimal biochemical wire. (A) Wire structure. The emitter strain E emits bbits to the medium at a constantrate, proportional to the strain density; the receptor strain R senses the bbit signal intensity; and the sink strain S degradesbbits with a Michaelis-Menten dynamics proportional to the strain density. B stands for the bbit signal intensity, which willbe bbit concentration in most cases. (B) Wire elements. The emitter half-wire module (top) is the minimum set of genes thatemit bbits when the corresponding promoter is induced. The receptor half-wire module (middle) is the minimum set of genes(constitutively expressed from the P C promoter) that sense the bbit and induce gene expression from the promoter P RHW promoter when bbit concentration in the medium is above some threshold. The sink module (bottom) is the minimum set ofgenes (constitutively expressed from the P C promoter) that degrade bbit in the medium. In the figure, the red stop signs standfor transcription terminators. (C) Trilinear coordinate system. Each point within the triangle represents a unique combinationof emitter, receptor and sink. Any specific point within the wire space will correspondto a unique population, whose behavior we will explorein the following sections.
III. CHARACTERIZING THE BBIT BUILD-UP
We can calculate the fixed points of Eq. (1) to obtainthe explicit expression of bbit concentration at equilib-rium: B eq = K E K S εD m σ − K E ε (2) The pairs ( ε, σ ) that make the denominator in B eq equal to zero divide the wire space in two regions(see dashed line inside the triangle in the left panel ofFig. 2A). When D m σ > K E ε (region at the left ofthe line, covered by green dots), the bbit concentrationreaches the equilibrium value B eq defined by Eq. (2).We call this the convergence region (represented by C below). Solving the differential equation (1) in timewithin the convergence region shows that the bbit timetraces converge to their predicted equilibrium concentra-tion (green lines in the middle panel of Fig. 2A). Thiscorresponds to a situation in which the production ofthe bbit by the EHW strain is balanced by its degrada-tion by the SINK strain (bottom inset in the right sideof Fig. 2A). In contrast, when D m σ < K E ε (region atthe right of the dashed line in Fig. 2A, covered by ma-genta dots), the production of bbit by the EHW strainis always larger than (and thus cannot be balanced by)the saturated degradation by the SINK strain (top insetin the right side of Fig. 2A). This leads asymptoticallyto a linear growth of B in time (magenta lines in themiddle panel of Fig. 2A). We call the region in which D m σ < K E ε the divergence region D . These tworegions are separated by what is called the horizon ofdivergence (represented by ˆ H ). In set notation, the wirespace is the union of convergence and divergence regionsand the horizon of divergence: ∆ = C ∪D∪ ˆ H . The math-ematical characterization of these regions is summarizedin Table I.Within the convergence region, the equilibrium con-centration of bbit obeys the following lemma: Lemma 1 (Isolines Lemma) . All points along anystraight line that radiates from the top ρ vertex of thewire space triangle and hits its σε bottom side reach thesame B eq value.Proof. We begin by noting that the expression (2) of thebbit concentration at equilibrium can be rewritten as B eq = K E K S D m σ/ε − K E . Thus, B eq depends on the consor-tium composition only via the ratio σ/ε . In turn, geo-metric similarity shows in a straightforward manner thatthe ratio σ/ε is constant for all points belonging to anystraight line radiating from the ρ vertex of the triangle.To see this, consider two points along one of these radiat-ing lines, characterized respectively by strain abundancefractions (relative frequencies) ( σ , ε ) and ( σ , ε ) , asshown in Fig. 2B. Let h and h be distances betweenthese two points and the ρ vertex. The similarity be-tween the two triangles delimited by the ε - h sides in onecase, and by the σ - h sides in the other, leads to: ε ε = h h = σ σ = ⇒ σ ε = σ ε (3)Thus B eq, = B eq, , which proves the lemma. The greenlines and dots in the left and center panels of Fig. 2Cdepict some isolines and the behavior of B ( t ) along them.Within the divergence region, each culture-point canbe characterized by the asymptotic rate ˆ v at which B ( t ) grows in time. Considering that B ( t ) has grown well We note that in the divergence region, bbit production wouldeventually stop due to limited capacity of the medium, toxic-ity, lack of precursors, etc. We ignore these limitations in whatfollows. past K S , the degradation term in Eq. (1) becomes linear,leading to ˆ v = P ( K E ε − D m σ ) (4)We can ask which are the points in the divergence re-gion that have the same ˆ v values. Geometric considera-tions show that all points lying in a segment parallel tothe horizon of divergence will have the same asymptoticgrowth speed given by Eq. (4). These growth-rate iso-lines are shown in the right and center panels of Fig. 2C(purple lines and dots).There are simple rules to gain intuition about both theconvergence and divergence regions. Within the conver-gence region, for instance, the closer an isoline is to thehorizon of divergence, the greater its associated B eq valuewill be. And within the divergence region, the closer agrowth-rate isoline is to the ε vertex, the greater its ˆ v value will be. At the horizon of divergence, ˆ v = 0 . Thetrilinear coordinate representation introduced here pro-vides us with an intuitive understanding of this multi-cellular information-processing system. We now examinethe response of this system to a modulation of the activ-ity of the emitter strain, and argue that in a significantportion of the wire space this response can take the formof a timer. IV. WIRE ON/OFF STATES DEFINE DIGITALAND BUFFER REGIONS
So far we have not considered the response of the re-ceptor strain RHW to the bbit present in the medium.As discussed above (Fig. 1B), in that strain bbit inducesthe expression of an output protein, such as a fluorescentprotein. We assume that the transfer function relatingthe bbit concentration with its output follows an acti-vating sigmoidal dynamics, described for instance by aHill function (Fig. 3A, left). In what follows, we normal-ize the response so that it falls within the interval (0 , .We also define two arbitrary thresholds for the receptorresponse, T and T , corresponding to the bbit concen-trations at which the receptor cells reach 5% and 95% oftheir maximum response, respectively. We note at thispoint that in an actual implementation of this analysis,the response of the system would be measured at the pop-ulation level, not at the single-cell level. This allows us toignore the heterogeneity existing not only in the responseof individual receptor cells, but also in the bbit produc-tion rate of individual emitter cells. Since all emittercells contribute with the same communication molecule(the chosen bbit) to the medium, where the bbit signalis averaged, this model is not affected by single-cell het-erogeneity. In this way, our arbitrary thresholds T and T could be measured precisely in the experimental im-plementation of this system.The receptor thresholds defined above, T and T , al-low us to split the convergence domain into three re-gions according to the induction level of the receptor Figure 2.
Dynamical dissection of the wire space . (A) All the points lying in the convergence region (green dots inwire space, left) are associated with convergent bbit time-traces (center, green lines). This is due to the intersection betweenthe bbit production term and the bbit degradation term of Eq. (1) (right, bottom inset). Meanwhile, all points lying in thedivergence region (magenta dots in wire space, left), are associated with quasi-linear bbit growth (center, magenta lines). Thisis due to the lack of intersection between the bbit production term and the bbit degradation term of Eq. (1) (right, top inset).(B) Graphical support for the B eq Isolines Lemma: ( σ, ε ) values for two points along an isoline define two similar triangles. (C)Isolines: when cultures are placed on the same line radiating from the top vertex to the bottom side within the convergenceregion, they reach the same B eq value (green dots and lines). In the divergence region, lines parallel to the divergence horizonhave the same asymptotic growth rate of B ( t ) . Parameter values for all simulations shown here: K E = 0 . , K S = 5 , D m = 1 ,and P = 5 . Convergence Region Horizon of Divergence Divergence RegionNotation C ˆ H D
Definition εσ < D m K E εσ = D m K E εσ > D m K E Property B eq = K E K S εD m σ − K E ε ( σ, ε, ρ ) = ( k ˆ σ, k ˆ ε, − k ) ˆ v = P ( D m σ − K E ε ) Table I. Definition and properties of the convergence and divergence regions. (Fig. 3A, right): (i) the uninduced region, named U , inwhich B eq ≤ T so that the receptor cells do not produceany measurable output; (ii) the induced region, named I ,where B eq ≥ T corresponding to a saturated maximaloutput readout; and (iii) the variable region in between,named V , where receptor cells reach a variable inductionlevel with B eq ∈ ( T , T ) . Within the divergence do-main, receptor cells obviously reach eventually the max-imum level of induction upon constant activation of theemitter strain, but we will keep distinguishing between this region and the induced region I for reasons thatwill become clear in what follows. We also note thatin the right panel of Fig. 3A, region D is color-codedaccording to the asymptotic bbit accumulation rate ˆ v .So far we have characterized the wire when the emit-ter is constantly active. But wires have two states bydefinition: ON and OFF, corresponding to the emit-ter being activated (EHW ON ⇒ K ON E > ) or inacti-vated (EHW OFF ⇒ K OFF E ≈ ), respectively. Thesetwo wire states have different response regions ( U , V , I , Figure 3.
Response characteristics and emitter modulation . (A) Transfer function associated to the receptor’s response(left panel). Two bbit concentration thresholds are fixed to define (right panel) when the receptor is uninduced ( B eq ≤ T ,region U ), when it is fully induced ( B eq ≥ T , region I ), or when it has a variable induction level B eq ∈ ( T , T ) , region V .These three regions span the convergence domain, and are color coded according to the value of B eq . The rest of the wirespace is covered by the divergence domain D . Within that region, color coding corresponds to the asymptotic accumulationrate of B ( t ) . These regions structure the wire space in a way that informs us about the response of the wire given a set ofparameters. (B) When two states of the wire alternate (induced and uninduced emitter), two new regions can be defined. Thedigital region Φ is the intersection of the induced region I (where the receptor is fully induced) when the emitter is ON, withthe uninduced region U (where the receptor is uninduced) when the emitter is OFF. The buffer region Ψ is the intersection ofthe divergence region D (where bbit accumulates) when the emitter is induced, with the uninduced region U when the emitteris OFF. Parameter values for all simulations shown here: K ON E = 0 . , K OFF E = 0 . , K S = 5 , D m = 1 , and P = 5 . Hill functionparameters: K = 1 , with Hill coefficient n = 3 . and D ) because the parameter K E changes between thetwo states (Fig. 3B). If the difference between K ON E and K OFF E is large enough, interesting intersections betweenthe two states can emerge. The digital region Φ is theintersection between the U OFF region (uninduced region U when the emitter is OFF), and the I ON region (in-duced region I when the emitter is ON). In set notation, Φ = U OFF ∩ I ON . Within this region the consortium be-haves in a fully digital manner, transducing directly intothe response module/strain the input driving the emit-ter module/strain. A more interesting region is the bufferregion Ψ , which corresponds to the intersection between the U OFF region (uninduced region U when the emitteris OFF), and the D ON region (divergence region D whenthe emitter is ON). In set notation, Ψ = U OFF ∩ D ON .As we will see in what follows, within that region the re-sponse strain encodes the time during which the emitterstrain has been ON. Complementarily, the time duringwhich the response is active can be programmed by fix-ing the duration of the ON state of the input. In otherwords, within the Ψ region the consortium can behave asa timer. V. THE BUFFER REGION OPERATES AS ATIMER
The two regions described in the previous section havedifferent and useful behaviors. First we focus on the dig-ital region Φ . All culture-points within this region willbehave in a digital manner, meaning that when the emit-ter is ON, the receptor will be fully induced. In turn,when the emitter changes its state to OFF, then the re-ceptor changes its state to fully uninduced. A key pointof this digital region is that the time that it takes to tran-sition from the uninduced state to the induced state and,most importantly, the time needed for the reporter tobecome uninduced again after the emitter is turned OFFis constant anywhere within this region for an arbitraryculture-point within it, and independent of how long theemitter is ON (Fig. 4A). This is the crucial feature thatallows the system to behave in a digital manner (whichwill not be fulfilled in the buffer region, and we will seebelow). A potential limitation of the digital behavior ofthe system in this region is the time needed for the systemto react when the emitter is turned on. For the parame-ters chosen in Fig. 4A that time is negligible in compari-son with the duration of the pulse, but it is worth askinghow this turn-on time depends on the location of the con-sortium within the wire space. To address this question,we define the time t ON that it takes for the receptor toreach 95% of its corresponding B eq . This quantity can-not be calculated analytically, but it can be computedsystematically by simulating the behavior of the systemfor culture-points equally distributed over the wire space.The results are plotted in Fig. 4C, which shows that thisturn-on time is low for a large region of the wire space,diverging to infinity only when the culture is close to the ρ vertex (where the receptor strain dominates over theother two). In the bottom region of the wire space, theresponse is the fastest where the emitter strain domi-nates (i.e. most of the consortium is formed by emittercells). Supplementary Video 1 shows how the responsegets delayed eventually preventing receptor cells to getfully induced. This supplementary video also emphasizesthe behaviour of the wire for cultures lying on top of thesame isoline.Within the buffer region Ψ , all culture-points alsoreach maximum receptor induction when the emitter isturned ON, in a time that is fixed (i.e. is independentof the duration of the preceding OFF state) and given aswell by Fig. 4C as discussed above. However, in strikecontrast with the digital region, the time that the systemtakes to transition back from the induced to the unin-duced state when the emitter strain is turned OFF isnot fixed, but crucially depends on the duration of theprevious ON state. In other words, the longer the ONstate of the emitter, the longer it will take for the re-ceptor to transition from the induced to the uninducedstate (Fig. 4B). This prevents the use of the system as adigital encoder within this region, but enables a second,potentially more interesting, ability, namely measuring the time that the ON signal was present (or conversely,programming the time that the receptor is going to beinduced, by tuning the time during with the emitter isON). See Supplementary Video 2 to better appreciate thedifference between buffer and digital regions.To quantify this behavior, we define t IN as the durationof the input pulse, and introduce the magnitude t OFF tomeasure how long it takes for the wire receptor to tran-sition from the induced to the uninduced state, startingwhen the input signal disappears, and ending when thereceptor reaches again a response less or equal to 5%of its maximum response value at the uninduced state(see double-headed arrows in Fig. 4B). Figure 4D com-pares the behavior of this magnitude in the two regimes:whereas in the digital region t OFF is not affected when t IN changes (as we anticipated above), within the bufferregion the response pulse is extended a time t OFF thatincreases linearly in response to increasing t IN , therebyproviding the population with a clean control of the re-sponse time duration. The fact that t OFF remains fixedbut greater than zero when the culture is within the digi-tal region associates a lag time to that region, whereas inthe buffer region we can talk about buffer time. There-fore, biological wires can store information about thetime they have been in the ON state when the relativeratios of the different strains forming the consortium liewithin the buffer region. This establishes the underlyingprinciple of the collective timer behaviour.
VI. FILTER RESPONSE TO PERIODIC INPUTSWITHIN THE BUFFER REGION
So far we have characterized the response of our timerto a single pulse of emitter activity. However, in real-lifesituations such activity pulses will occur repeatedly alongtime. In order to analyze the timer functionality in theface of repeated stimulus presentations, we consider inwhat follows that the emitter activation takes the formof a periodic train of pulses. Such digital periodic signalscan be characterized by their period τ and their dutycycle dc , which corresponds to the fraction of the cycleduring which the signal is ON (Fig. 5A).We are interested in establishing the operation limits ofthe timer due to the interaction between consecutive ONpulses. To that end, we explore in Fig. 5B the response ofthe system in time across all the different regions of thewire space. Moving from left to right across the triangle,we first note that even before reaching the digital region(white region at the left) the system can already respondto the input pulses, although with a non-maximal out-put. As soon as we enter into the digital region Φ (greenregion in Fig. 5B) the response of the system surpassesthe T threshold, so that the response becomes fully in-duced, at least during part of the input pulse. Withinthat region, as we move further to the right, the turn-offlag time t OFF starts already to increase. This increasebecomes more evident when the system enters the buffer
Figure 4.
Digital response versus timer function. (A) Response of the consortium to an input signal with two differentdurations when the culture is in the digital region. (B) Corresponding behavior in the buffer region. (C) Response time t ON across the wire space. Darker colors correspond to a slower responses. (D) Dependence of the turn-off time t OFF on the inputduration ( t IN ) for both the digital (green) and buffer (magenta) regions. Parameters used for these simulation are the same asin figure 3. region Ψ (magenta region in Fig. 5B), until finally theresponse to two consecutive pulses overlap and the signalis completely filtered out, rendering the timer unusable.For the period and duty cycle values used in these sim-ulations, the system starts to filter before reaching theend of the buffer region, but this will depend on the sys-tem parameters. Figure 5C shows the location of this filter horizon for varying values of the duty cycle and afixed period. The horizon gets closer to the triangle side ρε , beyond the buffer region, for small duty cycles, andcan extend into the digital region (and thus away fromthe region where the timer operates) for large duty cyclevalues. This behavior imposes a limit in the wire spacedomain where the timer can function. This domain canbe visualized in Fig. 5D, which represents the regions inwhich the signal is either always filtered or never filteredirrespective of the specific values of the relative frequen-cies ε , σ and ρ , depending on the value of the ε/σ ratioand of the duty cycle d c . The figure shows that as theduty cycle increases, the ratio of emitters with respect to that of sink cells should decrease in order to avoid fil-tering the signal. For large enough duty cycles, no partof the buffer region is free from the filtering effect, andthus the timer cannot operate. This plot thus serves as adesign guide for multicellular timers in the case in whichthe system must respond to repetitive pulses. In real-lifesituations, both the period and the duty cycle will varyin time, in which case the consortium will move acrossthe phase diagram. See supplementary video 3 to visual-ize how a culture-point lying within the digital or bufferregions processes a periodic digital signal. VII. DISCUSSION
Here we have proposed a multicellular architecturebased on three bacterial strains, an emitter, a receptorand a sink, that communicate with each other via a bio-logical bit (bbit). Together, these three components con-stitute what can be considered a chemical wire. The sink
Figure 5.
Filtering periodic pulses. (A) Graphical definition of duty cycle ( d c ): given a digital periodic signal acting asinput signal, the duty cycle is the percentage of the period ( τ ) that is in the ON state. (B) Simulations of the receptor statefor different culture locations within the wire space (color plot) for a fixed period of 10 time units and a 30% duty cycle of thepulse train. Solid black lines show where the T and T are crossed in the simulations. The figure on the left represents theinput signal. The bottom figure shows the wire space with the digital (green) and buffer (magenta) regions. The vertical blackdashed line on the color plot shows the ratio ε/σ at which the input signal is filtered (filter horizon). (C) Filter horizons fora fixed period of 20 time units and different duty cycles ranging from 10% up to 90%. Pink shading emphasizes the locationof the buffer region, and shows how the different filter horizons cover this region. (D) Phase diagram showing the regions inwhich the periodic pulse train is filtered in the plane formed by ε/σ and the duty cycle, for a fixed period of 20 time units.The digital and buffer regions are delimited by light grey dotted lines. In the region between the two orange lines the filteringresponse depends on the specific values of the relative frequencies ε , σ and ρ . Parameter values for all plots remain the sameas in figure 3. strain degrades the bbit in the culture, allowing the sys-tem to adapt its response to the input signal in the fastestway possible. In the absence of that strain, the commu-nication molecule can only disappear from the culture ifthe media is washed or if the molecule is degraded on itsown. Controlling the bbit degradation allows fine tun-ing of the system dynamics, unlocking its responsivenessto a dynamical input. To gain intuition about how thecollective behavior of these three strains in culture, wehave proposed a conceptual representation using trilin-ear coordinates to visualize the entire wire space. This theoretical and representation framework is intended toprovide a simple way to interpret the capabilities of achemical wire implemented with our proposed architec-ture. Ideally, in an experimental setting, just measuringthe constants K OFF E , K ON E , K S , D m and the parame-ters of the receptor response, K m and its Hill coefficient,would allow us to predict the range of relative frequenciesfor each strain in the media at which our chemical wirewill work, and its performance.Using the approach outlined above, we have identifieda domain in the wire space, which we call the buffer re-0gion, that allows the wire to store information about thetime during which the emitter has been in the ON state.This represents an extra level of information, beyond thecurrent state of the emitter, that is stored directly onthe wire, operating as a timer. Our theoretical and com-putational results show that this biological function isrobust not only to the specific values of the relative fre-quencies of the three strains for which the system canbehave as a timer, but also in its response to repetitivestimulation pulses. In particular, our characterization ofthe consortium response to a periodic train of pulses re-vealed that the signal is filtered out only for large enoughduty cycles, and when the relative frequency of emitters islarge enough in comparison with that of sink cells. In therest of the parameter space, the consortium performs ro-bustly as a multicellular timer, which can be easily tunedby varying the relative frequencies of the three strains,thereby moving freely across wire space.Our model is based on several simplifying assumptions,including that the input signal is directly transduced intobbit production by the emitter strain, and that bbit isdirectly converted into output by the receptor strain.These assumptions are not limiting, since our model can be considered a starting point for models of increasingcomplexity and richness. In that way, our model setslower bounds for the wire to work properly. Specifically,if our model predicts that a wire will not have a digitalor buffer region, increasing the descriptive capability ofthe model by adding more complexity will not fix thatissue. On the other hand, if the model claims that eitherdigital or buffer regions exist, a more descriptive modelcould be used to improve that prediction quantitatively,but the results would not change qualitatively.The analysis presented above has been phrased for clar-ity from a design perspective, but it is also applicable tonatural multicellular circuits. The type of cellular pro-cesses described here (production, detection and degra-dation of a diffusible molecule) are very generic and existnaturally in cells. Here we show that combining theseprocesses in different cell types forming a multicellularpopulation could give rise to a robust and cellular timer.Thus, beyond their obvious implications for syntheticbiology applications, our results also constitute a fur-ther instance of how cell-cell communication can lead tocollective functionality in cellular populations that goesfar beyond the capabilities of gene and protein circuitswithin individual cells. [1] Alon, U., 2019. An introduction to systems biology: de-sign principles of biological circuits. CRC press.[2] Andrianantoandro, E., Basu, S., Karig, D.K., Weiss,R., 2006. Synthetic biology: new engineering rules foran emerging discipline. Molecular Systems Biology 2,2006.0028. doi:doi:10.1038/msb4100073.[3] Brenner, K., You, L., Arnold, F.H., 2008. Engineer-ing microbial consortia: a new frontier in synthetic bi-ology. Trends in Biotechnology 26, 483 – 489. doi:doi:10.1016/j.tibtech.2008.05.004.[4] Cai, H., Katoh-Kurasawa, M., Muramoto, T., San-thanam, B., Long, Y., Li, L., Ueda, M., Iglesias, P.A.,Shaulsky, G., Devreotes, P.N., 2014. Nucleocytoplas-mic shuttling of a gata transcription factor functions asa development timer. Science 343, 1249531. doi:doi:10.1126/science.1249531.[5] Danino, T., Mondragón-Palomino, O., Tsimring, L.,Hasty, J., 2010. A synchronized quorum of genetic clocks.Nature 463, 326–330. doi:doi:10.1038/nature08753.[6] Durand, B., Raff, M., 2000. A cell-intrinsictimer that operates during oligodendrocyte develop-ment. BioEssays 22, 64–71. doi:doi:10.1002/(SICI)1521-1878(200001)22:1<64::AID-BIES11>3.0.CO;2-Q.[7] Eckert, B., Martin, A., Balbach, J., Schmid, F.X., 2005.Prolyl isomerization as a molecular timer in phage infec-tion. Nature Structural & Molecular Biology 12, 619–623.doi:doi:10.1038/nsmb946.[8] Ellis, T., Wang, X., Collins, J.J., 2009. Diversity-based,model-guided construction of synthetic gene networkswith predicted functions. Nature Biotechnology 27, 465–471. doi:doi:10.1038/nbt.1536.[9] Enright, J., 1980. Temporal precision in circa-dian systems: a reliable neuronal clock from unre- liable components? Science 209, 1542. doi:doi:10.1126/science.7433976.[10] Fullstone, G., Bauer, T.L., Guttà, C., Salvucci, M.,Prehn, J.H.M., Rehm, M., 2020. The apoptosome molec-ular timer synergises with xiap to suppress apoptosis ex-ecution and contributes to prognosticating survival incolorectal cancer. Cell Death & Differentiation doi:doi:10.1038/s41418-020-0545-9.[11] Gabalda-Sagarra, M., Carey, L.B., Garcia-Ojalvo, J.,2018. Recurrence-based information processing in generegulatory networks. Chaos 28, 106313. doi:doi:10.1063/1.5039861.[12] Gefen, O., Fridman, O., Ronin, I., Balaban, N.Q., 2014.Direct observation of single stationary-phase bacteria re-veals a surprisingly long period of constant protein pro-duction activity. Proceedings of the National Academy ofSciences 111, 556–561. doi:doi:10.1073/pnas.1314114111.[13] Gerecht, K., Margiola, S., Müller, M.M., 2020.p53 deamidation as a molecular timer for celldeath. Biophysical Journal 118, 485a. doi:doi:10.1016/j.bpj.2019.11.2687.[14] Hou, W.S., Van Parijs, L., 2004. A bcl-2-dependentmolecular timer regulates the lifespan and immunogenic-ity of dendritic cells. Nature Immunology 5, 583–589.doi:doi:10.1038/ni1071.[15] Koseska, A., Zaikin, A., Kurths, J., Garcia-Ojalvo, J.,2009. Timing cellular decision making under noise viacell–cell communication. PLOS ONE 4, e4872. doi:doi:10.1371/journal.pone.0004872.[16] Levine, J.H., Elowitz, M.B., 2014. Polyphasic feedbackenables tunable cellular timers. Current Biology 24,R994–R995. doi:doi:10.1016/j.cub.2014.08.030.[17] Levine, J.H., Fontes, M.E., Dworkin, J., Elowitz,1