Supporting information for “theoretical investigation of a genetic switch for metabolic adaptation”

1 Department of Chemistry and Applied Biosciences, ETH Zurich, 8093 Zurich, Switzerland 2 Department of Physics, California Institute of Technology, Pasadena, CA, USA 3 Division of Biology and Biological Engineering, California Institute of Technology, Pasadena, CA, USA 4 Howard Hughes Medical Institute, University of Washington, Seattle, WA, USA 5 Department of Systems Biology, Harvard Medical School, Boston, MA, USA 6 Department of Physics and Astronomy, University of Southern California, Los Angeles, CA, USA

Results for different models of transcription. The parameters are the same as the ones that were used for all other phase portraits in this report and that are listed in a table in the main part. In (A), the phase portrait after addition of non-zero transcription rates to all other polymerase-bound states is shown (third assumption lifted). The new parameter values are r I = 10 −2 r p for the transcription rate when only one XapR is bound and r 0 = 10 −3 r p for the basal rate. In (B), the phase portrait after introducing different XapR dissociation constants can be seen (second assumption lifted). The new parameter values are [XapR] R,1 = 10 and [XapR] R,2 = 0.1. Table 1. States and weights for a more precise model of transcription.

XapR
Pol. Transcription Weight none no The first column shows which XapR sites are occupied in this state, the second one whether polymerase is bound to the promoter, and the third the transcription rate. In the weights, the notation is as follows: w X1 is the weight for XapR binding to the first binding site on the promoter, w X2 is that to the second site, w p is the weight of polymerase binding to the promoter, e ∆ X1 is the interaction energy between one XapR on the first binding site and the polymerase that is bound, e ∆ X2 is the same for the second binding site, and finally, e ∆ coop is the interaction energy of the two XapR at both binding sites as in the main part.
Furthermore, we do not expect our approximation to have a large influence: the interaction energy will make the doubly bound state dominate, and the polymerase-XapR interaction should be weak in the singly bound state, because very weak activation is observed in the experiment when one of the XapR binding sites is removed (see section E in this supporting information text). If e ∆ X1 and e ∆ X2 are small, we have w p w X1 e ∆ X1 ≈ w p w X2 e ∆ X2 ≈ w p w X . If, in addition, the transcription rates for the singly bound states are much smaller than those for the doubly bound state, we can neglect the singly bound states in the nominator of p active (which, as argued in the previous paragraph "Simplistic modeling of transcription", is a good assumption). Thereby, we can obtain the same form of p active that we are working with in the main part by including the interaction between polymerase and XapR into coop .
Working with only one protein species. The proteins are transcribed together, and thus, their rate of transcription and mRNA decay should be the same. The protein decay is assumed to be dominated by cell division, making the decay rates roughly equal as well. What could differ is the translation rate. Changes in this rate rescale the protein synthesis rate by the same factor and thus, in steady-state, change the total number of proteins. The only other place where the proteins appear in the ODE's is in their respective Michaelis-Menten term, and the latter is proportional to k cat [E] 0 . Hence, different numbers of XapA and XapB can be accounted for by changing their respective values of k cat .
Exponential protein decay. The protein decay should be dominated by dilution through cell division, which is a discrete process. However, the cell grows before division, resulting in a continuous dilution in concentration. Furthermore, replication leads to a varying number of gene copies and thereby also a varying transcription rate. Trying to capture all these mechanistic details would lead to an overcomplicated model. Treating the protein decay as a Poisson process has produced correct results in past research, which is why we expect this assumption to be reasonable.
Kinetics of XapB, NupC, and NupG. All statements in this paragraph are based on data from [3]. Xanthosine entering the cell consists of two steps. In a first transport step, xanthosine passes the outer membrane through porins like Tsx, OmpF and OmpC. This does not seem to be a rate-limiting step. Then, XapB actively transports it across the inner membrane, powered by the proton motive force. There are two other active nucleoside transporters that seem to be transporting xanthosine with a very low affinity: NupC and NupG. It was found that ΔnupC ΔnupG strains cannot grow on xanthosine. Hence, these seem to be necessary to "start" the system by importing small amounts of xanthosine that then activate xapB, and there appears to be no other significant way in which xanthosine can enter the cell.
As mentioned in the main text, we assume the kinetic scheme of all three transporters to be similar to that of lac permease (as it is described in [4]), which is discussed in more detail in the next subsection. One transporter of this type can do two things to change the intracellular xanthosine concentration [x]: it can actively transport one xanthosine into the cell, powered by the proton gradient, or it can transport one out of the cell, against the proton gradient. For net transport, a proton and a substrate need to bind on one side of the membrane and detach from the transporter on the other side of the membrane before the transporter goes back to the first side again (see next subsection). Because of the proton gradient, influx is overall thermodynamically more likely than export at low intracellular xanthosine concentrations, which leads to a net influx. For much higher intra-than extracellular xanthosine concentrations, the difference in the chemical potential of xanthosine across the membrane can dominate that of the protons and there is net efflux.
We model influx and efflux separately and use Michaelis-Menten kinetics for XapB (this is discussed in more detail in the next subsection). For the turnover rate and Michaelis constant for influx we write k b,i and K b,i , respectively, and k b,e , K b,e for efflux. These parameters implicitly include the proton gradient.
For NupC and NupG, we approximate the transport reactions as linearly dependent on the substrate. This should be appropriate, since their affinity for xanthosine is so low that it cannot even be measured properly [3], which makes K M much larger than c or [x] will ever be. [E] 0 is denoted by k nup for influx, and the efflux rate is written as ξk nup . For the XapB kinetics, no such approximation can be made, because the xanthosine concentration in the dynamic system can range from far below the respective K M value to far above.
The kinetics of lac permease. All transporters in our system are driven by the proton gradient [3], and we assume their kinetic scheme to be like that of lac permease. This scheme is shown schematically in Fig 2. Even though it includes several steps, the essential part is substrate binding, conformational change of the transporter, and substrate release on the other side of the membrane. We assume the proton gradient across the membrane to be constant and proton binding and unbinding to the transporter to be fast. Hence, the latter does not need to be included explicitly in the kinetic model (it is implicitly part of the turnover rate).
Substrate release should be much faster than conformational change, which is analogous to ignoring the "EP" state in Michaelis-Menten (enzyme with product bound, here the conformationally changed transporter with substrate still bound). Also, it should be reasonable to neglect the immediate backwards reaction (before proton release), and thus, Michaelis-Menten kinetics is obtained. Data from [5] and [6] actually show that, at a physiologically reasonable proton gradient, both, influx and efflux, can be described by the Michaelis-Menten equation very well.
What is neglected in this modeling is the fact that the transporter needs to release the proton and undergo a conformational change again before being able to conduct the same reaction again. However, by using experimentally determined apparent values to estimate the kinetic constants, this is already included in these parameter values. Another simplification we make is to treat influx and efflux as two independent reactions instead of thinking of it as one reversible kinetic process. This seems reasonable, though, because the transporter needs to lose its proton between reactions in order to obtain a net transport. If the transporter has for example done steps 1-6 in Fig 2, it can from there on really enter two entirely different reactions: it can do its conformational change in its empty state (steps 7 and 8) to then transport a substrate again (steps 1-6) or it can bind a new proton and a new substrate directly and transport them back (steps 6-1). The former would give influx and the latter efflux, so the two of them can be modeled as separate reactions.
Passive transport across the membrane. Most small molecules diffuse passively through the membrane, but at physiological pH, xanthosine is mostly negatively charged and charged substances do not significantly diffuse across the membrane on their own. Nevertheless, from the pK a of xanthosine we find that at physiological pH, about 1.5% Kinetic scheme for transport by proton-driven pumps. The mechanism of lac permease as for example described in [4] is used as a model for xanthosine transport. The enzyme binds a proton and a substrate (steps 1 and 2), undergoes conformational change (steps 3 and 4), and then releases substrate and proton on the other side of the membrane (steps 5 and 6). When the enzyme is empty, it can change its conformation again to repeat the transport (steps 7 and 8). Instead, it could also bind a new proton and substrate and transport them in the opposite direction (steps 6, 5, 4, 3, 2, 1). Alternatively, it could also keep the proton and just cycle back and forth, exchanging substrates from outside and inside the cell (steps 5, 4, 3, 2). of it is uncharged and could pass the membrane. This could be modeled as follows: Here, J is the flux through the membrane, which follows directly from Fick's first law (if steady-state diffusion is assumed). The other new parameters are the membrane area A and the membrane permeability coefficient p of the diffusing substrate. The diffusion rate k d = pA can be estimated with the help of Overton's rule. In our model, we find that this passive diffusion alone (i.e. without NupC and NupG) is enough to obtain bistability in the deterministic analysis as well as bimodality in the stochastic simulations. However, experiments show that ΔnupC ΔnupG strains cannot grow on xanthosine [3]. This suggests that passive diffusion alone is not enough to make the cells switch. Experiments with cells in minimal media containing some glucose and a xanthosine concentration in the bimodal regime show that the switching time can be several generations. The same is observed in our simulations. The xanthosine concentration that was used for the double knockout experiments, where no glucose was present, is just above the switching threshold, which is when switching is slow. For this reason, we suspect that without NupC and NupG and with only xanthosine as an energy resource, the cells starve before they have switched. This probably is the reason why passive diffusion is, in reality, not enough to observe switching. NupC and NupG can probably keep the cells going at very low intracellular xanthosine concentrations. Because of the active transport, they can accumulate a higher intracellular than extracellular concentration of xanthosine.
Since NupC and NupG transport seems to be much more significant than diffusion, we do not include the latter explicitly in our model. However, because it takes the form given in Eq 1, it can easily be accounted for implicitly by slightly changing the values of k nup and ξ.
Other mechanisms for xanthosine degradation. There probably are some ways in which xanthosine is degraded, other than by XapA. Otherwise, initially uninduced cells that are transferred to only xanthosine as an energy source might not survive until they have switched. As mentioned previously, it was found that they do [3], so they need to be able to metabolize xanthosine somehow. Nevertheless, these degradation pathways cannot have a big effect: as argued later on when estimating the Nup parameters, the intracellular xanthosine concentration in the uninduced state must be a bit larger, but not too much larger, than the extracellular one. Thus, degradation of xanthosine has to be small because import through Nup is weak. Due to this, it should be appropriate to linearize the equations describing these effects. Then, they can simply be accounted for in the term for Nup transport by making slight changes in the values of ξ and k nup . It could, however, also be that the basal level of XapA is enough to obtain the small degradation effect without any other pathways.

B Parameter estimation
In this section we present a detailed estimation of the values of all nondimensional parameters. Two dimensionful parameters appear in most nondimensional ones, namely γ p and K a . For this reason, we start by estimating these.
The protein decay rate γ p : Protein decay is generally dominated by dilution through cell growth and division. The average doubling time of E. coli is 20-30 minutes. We take this as an estimate for the protein half life, which leads us to a decay rate of Note that these upper and lower bounds correspond to half-lifes between 15 and 60 minutes, which should be reasonable. For the following nondimensionalization, we work with a fixed value of γ p = 5 · 10 −4 s −1 .
The Michaelis-Menten parameters for XapA (K a , k α ): The Michaelis-Menten parameters of XapA kinetics have been measured in two independent publications [7,8]. Both have found similar values and we conclude K a ≈ 6 ± 3 · 10 −5 M and k a ≈ 10 −1±0.8 s −1 . This gives k α ≈ 10 2±0.8 . We choose to work with K a = 5 · 10 −5 M in the following.
The extracellular xanthosine concentration [c] a : In the experiment, switching occurred if cells were grown in a solution with a concentration of xanthosine of roughly a few mM [3]. According to these experimental observations, the interesting range of xanthosine should be [c] a ∈ [0, 10 3 ] (nondimensionalized), which corresponds to c ∈ [0, 5 · 10 −2 ] M. However, we are careful because there are issues with the solubility of xanthosine.
The MWC parameters K χA , K IA , and ∆ x : The three MWC parameters are the energy barrier between the active and inactive state of XapR and the equilibrium constants of xanthosine binding to XapR in the two states. In Fig 3, the induction curve as a function of the nondimensional parameters is shown. We assume that evolution has tuned the induction curve such that less than 10% of XapR is active in the absence of xanthosine, and more than 90% is active for [x] approaching infinity. This is a relatively weak assumption, otherwise XapR would not really be fulfilling its job. It helps us by putting a constraint on K IA and ∆ x : A further constraint is found from K IA = KxI KxA = e β∆E , where ∆E is the energy difference between xanthosine binding to the two states. It is unlikely that this energy difference is more than 7 k B T , which corresponds to roughly three strong hydrogen bonds. This gives us an upper bound K IA < 10 3 and thereby ∆ x < 10 1 . For K IA = 10 2 , the upper bound is ∆ x < 7. This gives us a range of K IA ≈ 10 1 − 10 3 and ∆ x ≈ 2 − 10 or, more precisely, ∆ x ≈ 2 to ln 1 9 (K IA ) 2 ≈ 2 (ln (K IA ) − 1) < 12. Compared with other systems, these values do seem reasonable. To find K χA we furthermore assume that the steep part of the induction curve lies in the biologically relevant xanthosine regime, which estimated experimentally, see the paragraph about [c] a above. Again, this is a rather weak assumption. As can be seen from Fig 3, this assumption sets the value range of K χA . The line for ∆ x = 5 in Fig 3 is at the experimentally expected position in the xanthosine regime. It becomes clear from the figure that the position of the steep range is strongly affected by ∆ x , and a change of 1 in the latter leads to roughly a change of 10 1 in the former. Thus, we estimate K χA ≈ 10 2±1 · 10 ∆ x −5 and thereby finish our estimation of the MWC parameters.
The XapR concentration [XapR] R : It is experimentally known that the XapR copy number is very low, on the order of 10 −8 M (10 1 molecules per cell) [9]. We estimate the XapR -DNA dissociation constant from the values for lac that were found in [2] and obtain 10 −8±2 M, where we have added one additional order of magnitude as insecurity to the upper and lower bound. This gives The XapR cooperative energy ∆ coop : Any cooperative energy should come mainly from an interaction between the two XapR molecules. It would be surprising if it was larger than 10 k B T , which corresponds to roughly 4 strong hydrogen bonds. Hence, we estimate ∆ coop ≈ 0 − 10.
The transcription rate parameter ρ m : The transcription rate r m This estimate comes from the assumption that this rate should be roughly 10 1 − 10 2 times larger than a normal basal rate (of other genes), for which common values are around 10 −3±1 nM s −1 [10]. Because experiments show that compared to the activated promoter, the basal transcription rate is almost zero for xapAB (probably due to the weak polymerase binding site), this actually means a 10 2 -to 10 3 -fold increase in the transcription rate due to activation by XapR.
Keep in mind that the rate includes not only the effect of XapR on the pure speed of transcription, bt also any effect of XapR on the polymerase dissociation constant. Still, we do not expect higher values such as 10 0 nM s −1 since this is roughly the rate for transcription of rRNA and the latter should be significantly faster. To further validate our estimate, we can compare it to measured values for lac, which should be roughly the same because the gene function is similar. The transcription rate for lac in the absence of any repression that was found in [11] corresponds well with our estimate.
The Michaelis-Menten parameters for XapB (k β,i , K β,i , k β,e , K β,e ): To estimate the parameters for both influx and efflux, we need data from experiments where these two processes have been measured independently. Lactose permease, which we expect to function in a similar way as XapB, is comparatively well studied such that these rates have been measured in multiple ways and for different conditions. Some data from several publications is collected in [6] and we will base our estimates on this and the corresponding references. Comparing the effective rates and Michaelis constants of lactose permease and the nucleoside transporters NupC and NupG shows that their orders of magnitude are similar: from [12,13] we find k cat ≈ 10 1 s −1 , K M ≈ 10 1 − 10 2 µM for LacY and from [3,14,15] k cat ≈ 10 2 s −1 , K M ≈ 10 0 − 10 1 µM for pyrimidines through NupC and NupG. Note that especially for Nup, these are only rough measurements, and it is known that NupC behaves different with purines (thus also with xanthosine). Given these values, we expect lactose permease to be similar enough to nucleoside transporters so that we can use it as a starting point for our estimation of the XapB kinetic values.
We use the measured influx values that are summarized in [6] for a membrane potential of 100 mV and obtain k b,i ≈ 10 1±1 s −1 and K b,i ≈ 10 −4±2 M. Here, the uncertainty is estimated from the distribution of values of many kinds of enzymes that is presented in [10]. The bounds also reflect the difference between LacY and Nup that was found above. For the dimensionless quantities, this gives k β,i ≈ 10 4±1 and K β,i ≈ 10 1±2 .
For the efflux parameters, we proceed similarly and obtain k b,e ≈ 10 0±2 s −1 and K b,e ≈ 10 −3±3 M, therefore k β,e ≈ 10 3±2 and K β,e ≈ 10 2±3 . Because of disagreement in literature, difficulty in measuring these parameters, and a lack of knowledge about what value range should be expected, we put larger error bounds than for the influx.
The non-specific transport parameters k η , ξ: Estimating the rates of the Nup transporters for xanthosine is difficult since there are no experimental values that could be used as a reference. Because the rate cannot even be measured [3], it must be much lower than that of XapB. The kinetic parameters should, however, be such that the steady-state intracellular xanthosine concentration is significantly larger than the extracellular one. Otherwise, there would be no relevant difference to passive diffusion, which was shown to be insufficient for switching [3]. This essentially means ξ < 1, which on the other hand implies that the Nup rate is larger than or at least equal to the diffusion rate: If Nup transport were slower than diffusion, the latter would rule and ξ would be close to 1 (because this is the ξ value of diffusion).
In addition, we know from lac permease that influx is not much faster than efflux, which implies that ξ also cannot be too much smaller than 1. We guess ξ ≈ 0.8 ± 0.1. The number of Nup transporters per cell is roughly 10 2 − 10 3 [9]. We expect kcat KM to be at least two or three orders of magnitude lower for Nup than for saturated XapB, where we have k b K b ≈ 10 5±3 M −1 s −1 . These arguments set an estimated upper bound on the non-specific transport rate: To find a lower bound, we estimate the diffusion rate. As described in Eq 1 in the previous section, the diffusion rate can be written as 10 −2 · p · A, where p is the permeability coefficient of the membrane for xanthosine, A is the total cell surface area, and 10 −2 is the fraction of uncharged xanthosine at physiological pH. From the permeabilities of various substances that are collected in [10] and from the measured permeability values for purine nucleosides in [16], we estimate a permeability coefficient ofp = 10 −8±3 m s −1 . We write the tilde because this value is given in a convention wherẽ p · A = dx dt with x being the number of xanthosine molecules and not the concentration. Since we work in concentrations, we divide by an estimated cell volume of 10 −18 m 3 and, with a cell surface area of 10 −11 m 2 , find a diffusion rate of roughly 10 −3±3 s −1 . This sets the lower bound of k nup .

C Additional plots and explanations of the results
The 2D set of equations. If the mRNA concentration is set to steady-state as discussed in the main text, the dynamical system is reduced to the two equations shown below: Here, several of the parameters from the 3D system could now be combined into one by defining ρ . . = Increasing [c] a . Fig 5 shows a family of nullclines for increasing extracellular xanthosine concentration. It shows the transition from monostability, to bistability, and back to monostability again.
Including transcriptional bursting. To test the robustness of our results to larger stochasticity (which is expected in the real biological system), we explored how the addition of burstiness in transcription altered our stochastic simulations. We assumed a geometric distribution for the burst size s, i.e. p(s|θ) = θ · (1 − θ) s−1 with mean θ −1 = 4, which is comparable to or slightly larger than reported values [17]. We then truncate this distribution at a burst size of 10, meaning we did not consider larger bursts, leading to an actual mean of around 3.2. The transcription rate was rescaled such that the mean number of mRNA's produced remained the same.
We observe the same qualitative behavior but with stronger fluctuations around the mean and larger variation in the time to reach the steady-state mean. As an example, Fig 6A shows the time evolution of one run of the simulation with bursting for the same conditions (i.e. parameters and starting values) as the time evolution that was shown in the main text. In addition, the distribution of adaptation times can be seen in Fig 6B. Fig 7 shows the distribution obtained for the same and an additional condition (i.e. parameters and starting values) as for the distributions shown in the main text. The variance in the timescales and distributions is larger than without bursting. Furthermore, switching already happens for lower [c] a : in 7C, most cells are in the high expression state, whereas in the corresponding distribution in the main text, there were many uninduced cells, i.e. the first peak was higher. To obtain a bimodal distribution like in the main text, [c] a now needs to be set to 17.5 instead of 18.5 as shown in 7B. This was to be expected, since the larger fluctuations make it easier for the system to jump across the middle fixed point.
Still, there are no big qualitative differences to the results without bursting. Thus, the assumptions we made that lead to less stochasticity than in the real biological

D Chemical master equation
Here we provide some further details on our stochastic simulations by writing the chemical master equation for our model. The state of our system is described by the number of mRNA m, protein p, and xanthosine molecules x, and we write P (m, p, x) for the probability of that state. The chemical master equation for the system is: Here, all the parameters of course have to be measured in counts, not in concentrations.
As can be seen from the equation, none of the propensities depends explicitly on time and hence, we can use the hybrid algorithm between classical Gillespie and τ-leaping that we relied on. Note that in the case with bursts that is discussed above, the transcription term r m p active (x)P (m − 1, p, x) is changed to r m p active (x)

E Experimental materials and methods
Construction of ∆xapABR deletion strain. Cells containing plasmid pSIM6 were diluted 1:100 in 50 mL LB media and grown to an OD600 of ≈ 0.4 at 30 • C. The culture was immediately placed in a water bath shaker at 43 • C for 15 minutes and then cooled in an ice bath for 10 minutes. Cells were then spun down for 10 minutes (4, 000 g, 4 • C) and resuspended on ice in 50 mL of chilled water. This was repeated three times before resuspending in 200 µL of chilled water to generate competent cells. Homologous primer extension sequences were obtained from Baba et al. [18] and used to generate linear DNA containing a kanamycin resistance gene insert by PCR (amplified from plasmid pKD4), using primers that contained homology for the region on the chromosome to be deleted [19]. Electroporation was performed using 1 µL purified PCR product (about 100 ng DNA), mixed with 50 µL cells. Cells were resuspended in 750 µL SOC media and outgrowth was performed on a shaker at 30 • C for 90-120 minutes. Cells were then selected on an LB-agar plate containing kanamycin (30 µg/mL) and grown overnight at 30 • C. The deletions were confirmed by both colony PCR and sequencing. After confirmation, the deletion was transferred to a clean strain through P1 transduction and selection on kanamycin. Removal of kanamycin resistance was achieved by transforming cells with plasmid pCP20 (expressing FLP recombinase, carbenicillin resistance, on a temperature-sensitive replicon), with cells grown at 30 • C. Lastly, several colonies were then growth colony-purified once non-selectively at 43 • C and then tested for loss of antibiotic resistance for both kanamycin and carbenicillin.
Construction of ∆xapAB strain. Further experiments are underway in which precise control of the xapR copy number are needed and the preliminary data from these experiments are shown in the right panel of Fig. 2A of the main text. This strain was constructed by creating a XapR-mCherry fusion protein expressed from a tetR regulated promoter. This fusion construct was integrated into the ycbN locus of the ∆xapABR strain making an effective ∆xapAB strain. A pZS3-PN25-TetR expression plasmid was also inserted into this genetic background which permits titratable control over the xapR-mCherry copy number. For the data presented in this work, expression was maximized using 10 ng/mL of the TetR inducer anhydro-tetracycline. Further information of a similar construct (LacI-mCherry) can be found in [20] and [21].