Mathematical model of the IRMA synthetic network IRMA was developed as a testbed synthetic network in yeast for the design and validation of reverse engineering and modelling approaches [1] . IRMA consists of ve genes, CBF1, ASH1, SWI5, GAL4 and GAL80 and its topology comprises both transcriptional and protein-protein regulation mechanisms ( [1]). Figure 2 in the main text shows a schematic diagram of the regulatory interactions among the ve genes. The network topology comprises a positive feedback loop from CBF1 back to itself via GAL4 and SWI5, and a negative feedback loop from CBF1 back to itself via ASH1. A further regulation is present between SWI5 and GAL80 via the GAL10 promoter bound by GAL4. Transcription of network genes can be controlled by the presence of Galactose (GAL) in the growth medium, whose presence inhibits transcription of SWI5 from the GAL10 promoter. To capture the dynamics of the network a hybrid model (Supplementary Figure S1 ) approximating the dynamics in Glucose (F 1 ) and Galactose (F 2 ), has been readapted from [1] . Both the vector elds F 1 and F 2 share the same model structure as well as most of the parameters (v 3 ,k 4 andγ need a specic argumentation) as reported below: where are the system states. We used Hill functions to model transcription rates from promoters; the multiple regulation on CBF1 is modelled by the product of two Hill functions (AND regulation). A time delay τ is present in the equation for x 1 modelling the transcription of CBF1, which is aected by a 100 minute-long time delay due to the sequential recruitment of chromatin-modifying complexes to the HO promoter (which follows binding of SWI5 and other transcription factors) [2] . A list of all model parameters can be found in Supplementary Table S1 in [1] . Note that the model is hybrid as parametersv 3 ,k 4 andγ switch between two dierent sets of values depending on the carbon source (Galactose or Glucose). Control algorithm: design and implementation In this system, the GAL1 promoter drives the expression of the Gal1-Gfp fusion protein. Thus it can be described as a single input-single output (SISO) dynamical system. The input u(t) describes the presence of galactose or glucose in the growth medium. The output y(t) is the measured average level of uorescence of the Gal1-Gfp protein in the cell population. The control objective is a set-point regulation, where the cell population is required to express, over several generations, a xed amount of uorescence (control reference r(t)). To this end, we designed a simple control algorithm based on a Proportional-Integral regulator (PI) whose outputû(t) is a function of the control error e(t) (the mismatch between the desired and current output of the system e(t) = r(t) − y(t)) dened as:û Since cells can consume either galactose or glucose in a mutually exclusive manner, the continuous signalû(t) has to be decoded into a discrete signal denoting either galactose (u(t) = 2) or glucose (u(t) = 0). The above constraint allows an analogy with the problems faced in the design of feedback control strategies for power electronic circuits [3] . Here, switches and SCRs (silicon controlled rectiers) can 2 only be turned on or o, some output is typically measured or estimated and, particularly in industrial applications, compensating noise and external disturbances is of utmost importance. The simplest and most widely used control technique in this context is to use the PI regulator coupled to a PWM (Pulse Width Modulation) control strategy with a classic limited integral anti-windup [4] . This is also the strategy we adopted to control the cell population. In the simplest feedback implementation of the PI-PWM, a sawtooth signal is compared withû(t) in order to modulate the width of a rectangular pulse train, which is then used as control input (see [5] for further details). Namely, let be the sawtooth signal; then To control the GAL1 promoter we chose the sawtooth wave parameters as follows: α = 0, β = 2 and T = 5min. The gains of the PI controller, namely K p = 6 and K i = 0.3, were tuned following the Cohen-Coon strategy ( [6]) using a linear transfer function describing the GAL1 promoter previously derived in [7] : with parameters µ = 0.28, d = 68 and Θ = 134. IRMA can be modeled as an input-output system where the input u models the presence/absence of Galactose and the output y is the concentration of one of its genes, namely Cbf1 (x 1 ). Note that the input acts nonlinearly on the dynamics of the network as the presence of Galactose changes the values of all the Galactose-dependent parameters (namelyv 3 ,k 4 andγ). Moreover, as soon as Glucose is administered to the cells, these stop responding to Galactose, even if it is still present in the medium. Therefore, the control input (interpreted as Galactose concentration of 2 w/v% in the total volume of uid reaching yeasts) is restricted to be either ON (u = 1) or OFF (u = 0). The system output y = x 1 cannot be measured directly as a concentration. Instead, the cells were engineered so that CBF1p is fused with a GFP, the green uorescent protein [1] . In this way, higher concentrations of Cbf1p are associated to higher levels of uorescence. Moreover, as outlined in [8] , the amount of estimated uorescence can be directly related to the actual number of molecules. From a control perspective, the gene network model is, therefore, a highly nonlinear, hybrid, time-delayed dynamical system of the form: T , µ is the vector of parameter values in Glucose (u = 0) andμ is the vector of parameter values in Galactose (u = 1). Hybrid systems are often used to model gene networks (e.g. see [911]), where it is quite common to observe threshold dependent and switch-like activation or inhibition functions governing the dynamics of protein-protein or protein-gene interactions. From a control viewpoint, the task is to regulate the expression of CBF1 to some desired value, say r(t) by modulating the control input u(t) (the carbon source being administered to the cell). In what follows, we will ignore phenomena like Zeno behaviors (which would imply the gene network to switch on and o an unbounded number times in a nite and bounded length of time). To solve this problem, the challenge is to design an eective, yet simple, control strategy able to cope with the many unavoidable 3 constraints, which characterize biological systems (i.e. high levels of noise, incomplete knowledge and unmodelled dynamics). Cells are living organisms and they resist external actions very eectively. Moreover, it is important to minimize the need of knowing in great detail the biological process to be controlled, in order to make the control strategy as general as possible and applicable to any process of interest. The main design constraint that needs to be taken into account when synthesizing a controller for a biological system such as IRMA is that the cells tend to consume Glucose as their primary carbon source, even when galactose is present, since they obtain energy from it at a lesser energetical cost [12] . Moreover GLU-genes repress GAL genes as outlined in [13] ; therefore our control input is binary: either galactose or glucose, but not both. As in the case of the GAL1 promoter, we used the PI-PWM control strategy for the set-point control task. The sawtooth wave parameters for the PWM were set to α = 0, β = 10E − 5 and T = 5 min. In this way, the amplitude of the signal is 10% of the Cbf1 level at steady state in Glucose and the period is twice the settling time of Cbf1. With this choice, the constraint is satised of the control input being binary. In the case of the signal tracking control task, in order to compensate for the estimated delay of 100min in the CBF1 gene transcription, a strategy was implemented, which is inspired by the classical Smith Predictor scheme in [14] , as shown in Supplementary Figure S2 . The delay-free version of IRMA's model (block M in Supplementary Figure S2 ) is used to produce a numerical estimate of the anticipated system outputŷ. This is then delayed and subtracted from the quantied uorescence output, y, of the network so as to produce a signal which is ltered via a linear lter F and added back to y. The lter was designed empirically to suppress high frequency noise entering the control loop and is characterized y: The output y s of the delay compensation scheme is then compared to the desired output (y ref ) to obtain the error signal e. A Proportional-Integral controller takes e in input and computes the control signalû(t) with the gains K p = 175.6 and K I = 2.11. These gains were found with the Cohen-Coon method, as previously described for the Gal1 promoter; in the case of IRMA the parameters of the transfer function in eq. (8) were found to be µ = 0.0467, d = 146.85 and Θ = 667.62. A comparison between the step response of the linear approximation with the non linear model is shown in Supplementary Figure  S3 . The computed analog control signalû is then converted into binary values by means of the PWM modulator described above that feeds both the physical plant and the model used for delay compensation with the same binary input u, thus closing the feedback loop. The control strategy described in the previous section was implemented as a Finite State Automaton in MATLAB (Mathworks Inc.) for the set-point regulation of both GAL1 promoter and the IRMA network. The Finite State Automaton implementing the PI-PWM control algorithm is very simple: at each step (k) an image is acquired by the microscope, and the normalized uorescence signal is computed thanks to an image segmentation algorithm described in the next section. The uorescence signal y s (k) is compared against the reference signal y ref (k), to obtain the error e(k). The control input u(k) is then computed using the discrete-time implementation of the PI controller discussed in [4] .The control input u(k) is used to determine the duration of the pulse of Glucose or Galactose by means of the PWM strategy. The duration of each pulse corresponds to the time interval during which the syringe loaded with Galactose remains higher than the one containing Glucose (or vice-versa). At the next instant (k + 1) a new image is acquired and the feedback computation takes place. The error e(k + 1) is available for a new control iteration and each step is repeated again. In the case of the signal tracking control task for the IRMA network, the Finite State Automaton implementing the control algorithm is more complex, since in this case there is the need of simulating the IRMA model M in the Predictor block, necessary for the delay compensation. As represented in Supplementary Figure S4 , at the k-th step, e(k) is known, thus it is possible to calculate the control input u(k) to be applied to both the plant and the IRMA model M via the PI controller [4]. The control input is then used to compute the duration of the pulse of Glucose or Galactose by means of the PWM strategy. In-silico analysis of the control algorithm for the IRMA network Given the complexity of the IRMA model, prior to the experimental implementation of the control strategy in -vivo, we decided to test and validate the performance of the control scheme in -silico. In the numerical implementation, several alternative design options were evaluated. For example, it was found that using a nonlinear model for delay compensation (in the block labelled as 'M' in Supplementary Figure S2 ) gives a comparable performance to that obtained when using a linearized one (as suggested in [14] ). In addition, we performed numerical simulations of the control scheme against uncertainties in the delay estimation by considering the scenario where the delay τ is set to zero. The results are presented in Supplementary Figure S5 and Supplementary Figure S6 showing the ability of the control scheme to guarantee a reasonable performance also in this case. Finally, we explored the eects of simplifying the control strategy further by removing the prediction scheme used to compensate the delay. Although a simpler control strategy may decrease the overall performance, the reason for taking this step is that we prefer to sacrice precision in favour of a control strategy which can easily be used to control other gene networks and proteins, without requiring a dynamical model of the system to be controlled. In-silico results reported in Supplementary Figure S7 and Supplementary Figure S8 show the performance of such a simple PI/PWM control strategy when implemented to control gene expression levels in IRMA. The in silico experiments conrm that this strategy force eectively the Cbf1 concentration to track the desired reference signal when the delay τ = 0. The closed-loop control platform described above, employs an inverted epiuorescence Nikon-TI Eclipse microscope. The microscope is programmed to acquire two types of images: (a) a phase contrast image (phase-contrast) and (b) two uorescence images (one for the green spectrum for GFP and one in the red spectrum for Sulforhodamine B). The red dye Sulforhodamine B is added to the galactose medium and it is used to check for the proper administration of the control input. Once cells have been imaged, image analysis methods can be applied to estimate their uorescence [15] . To this end, we developed a custom image processing algorithm to fully exploit the peculiarities of our platform. [16] Our aim was to carry out both cell segmentation and cell-tracking, maximizing sensitivity as primary objective, and then rening the results by improving specicity. We based our segmentation method on global edge linking system. Yeast cells in phase contrast images occur in clustered, low intensity, convex and often quasi-circular shapes surrounded by a white halo (Supplementary Figure S9 C) . The contrast between the pixels belonging to the cells and the pixels belonging to the halos is usually so high that edge points can be detected by the evaluation of the magnitude of the gradient calculated in each point of the image (Supplementary Figure S9) . Due to the shape of yeast cell, edge points can be connected with the Circular Hough Transform (CHT) [17] . CHT can detect almost all cells within the image, even when cells edges overlap; unfortunately 5 the CHT algorithm is computationally expensive (both in memory and time) and shows limited specicity capabilities. Therefore, in order to limit the computational cost associated to the algorithm, we introduce a preprocessing stage of the image to select the regions in the image where cells are located ( Supplementary  Figure S9 B ). Such regions are identied via a mask obtained through morphological operations (aperture, closure and lling) and thresholding (see Supplementary Figure S9 A). False detections are reduced as a result of the segmentation of these regions. Moreover, since the area of the regions in the image containing cells is smaller than the area of the background, computational time is considerably reduced. Once cells have been located in the phase contrast image, a binary lter is used to detect only pixels of the GFP uorescence eld within cells. Let us dene the uorescence eld image I as:
In the view of the above, we will consider the setup documented in REF as the reference platform for the in vivo implementation of our model selection approach.
In-Vivo Real-Time Control of Protein Expression from Endogenous and Synthetic Gene Networks
PLoS Computational Biology
PLoS Computational Biology
Biology Computer Science Medicine