Supplemental Information For Chapter IV: The Physiological Adaptability of a Simple Genetic Circuit

Published as …

---

A version of this chapter is currently under review. A preprint is released as Chure, G, Kaczmarek, Z. A., Phillips, R. Physiological Adaptability and Parametric Versatility in a Simple Genetic Circuit. bioRxiv 2019. DOI: 10.1101/2019.12.19.878462. G.C. and R.P. designed experiments and developed theoretical models. G.C. and Z.A.K. collected and analyzed data. G.C. and R.P. wrote the paper.

Non-parametric Inference of Growth Rates

In this section, we discuss the measurement of the bacterial growth curves as well as our strategy for estimating the growth rates for all experimental conditions in this work.

Experimental Growth Curves

As is described in the Materials and Methods section of Chapter 4, we measured the growth of E. coli strains constitutively expressing YFP using a BioTek Cytation 5 96-well plate reader generously provided by the lab of Prof. David Van Valen at Caltech. Briefly, cells were grown overnight in a nutrient rich LB medium to saturation and were subsequently diluted 1000 fold into the appropriate growth medium. Once these diluted cultures reached an OD_600nm of  ≈ 0.3, the cells were again diluted 100 fold into fresh medium preheated to the appropriate temperature. Aliquots of 300 μL of this dilution were then transferred to a 96-well plate leaving two-rows on all sides of the plate filled with sterile media to serve as blank measurements. Once prepared, the plate was transferred to the plate reader and OD_600nm measurements were measured every ≈ 5 - 10 min for 12 to 18 hours. Between measurements, the plates were agitated with a linear shaking mode to avoid sedimentation of the culture. A series of technical replicates of the growth curve in glycerol supplemented medium at 37^∘ C is shown in Fig. 1 (A).

Inference of Maximum Growth Rate

The phenomenon of collective bacterial growth has been the subject of intense research for the better part of a century (Schaechter, Maaløe, and Kjeldgaard 1958; Jun et al. 2018) yielding many parametric descriptions of the bulk growth rates, each with varying degrees of detail (Jun et al. 2018; Allen and Waclaw 2018). For the scope of this work, we are not particularly interested in estimating numerous parameters that describe the phenomenology of the growth curves. Rather, we are interested in knowing a single quantity – the maximum growth rate – and the degree to which it is tuned across the different experimental conditions. To avoid forcing the bacterial growth curves into a parametric form, we treated the observed growth curves using Gaussian process modeling using as implemented in the FitDeriv software described in Ref. (Swain et al. 2016). This method permits an estimation of the most-likely OD_600nm value at each point in time given knowledge of the adjacent measurements. The weighting given to all points in the series of measurements is defined by the covariance kernel function and we direct the reader to Ref. (Swain et al. 2016) for a more detailed discussion on this kernel choice and overall implementation of Gaussian process modeling of time-series data.

     As this approach estimates the probability of a given OD_600nm at each time point, we can compute the mean and standard deviation. Fig. 1 (B) show the raw measurements and the mean estimated value in dark green and light green, respectively. With the high temporal resolution of the OD_600nm measurements, modeling the entire growth curve becomes a computationally intensive task. Furthermore, as we are interested in only the maximum growth rate, there is no benefit in analyzing the latter portion of the experiment where growth slows and the population reaches saturation. We therefore manually restricted each analyzed growth curve to a region capturing early and mid exponential phase growth, illustrated by the shaded region in Fig. 1 (A).

     With a smooth description of the OD_600nm measurements as a function of time with an appropriate measure of uncertainty, we can easily compute the time derivative which is equivalent to the growth rate. A representative time derivative is shown in Fig. 1 (C). Here, the dark green curve is the mean value of the time derivative and the shaded region is ± one standard deviation. The maximum value of this inferred derivative is the reported maximum growth rate of that experimental condition.

Approximating Cell Volume

In Fig. 4.2 of Chapter 4, we make reference to the volume of the cells grown in various conditions. Here, we illustrate how we approximated this estimate using measurements of the individual cell segmentation masks.

      Estimation of bacterial cell volume and its dependence on the total growth rate has been the target of numerous quantitative studies using a variety of methods including microscopy (Pilizota and Shaevitz 2012, 2014; Taheri-Araghi et al. 2015; Schmidt et al. 2016; Schaechter, Maaløe, and Kjeldgaard 1958) and microfluidics (Kubitschek and Friske 1986), revealing a fascinating phenomenology of growth at the single cell level.

     Despite the high precision and extensive calibration of these methods, it is not uncommon to have different methods yield different estimates, indicating that it is not a trivial measurement to make. In the present work, we sought to estimate the cell volume and compare it to the well-established empirical results of the field of bacterial physiology to ensure that our experimental protocol does not alter the physiology beyond expectations. As the bulk of this work is performed using single-cell microscopy, we chose to infer the approximate cell volume from the segmentation masks produced by the SuperSegger MATLAB software (Stylianidou et al. 2016) which reported the cell length and width in units of pixels which can be converted to meaningful units given knowledge of the camera interpixel distance.

     We approximated each segmented cell as a cylinder of length a and radius r capped on each end by hemispheres with a radius r. With these measurements in hand, the total cell volume was computed as
$$ V_\text{cell} =\pi r^2 \left(a + \frac{4}{3}r\right). \qquad(1)$$
The output of the SuperSegger segmentation process is an individual matrix for each cell with a variety of fluorescence statistics and information regarding the cell shape. Of the latter category, the software reports in pixels the total length ℓ and width w of the cell segmentation mask. Given these measurements, we computed the radius r of the spherocylinder as
$$ r = \frac{w}{2} \qquad(2)$$
and the cylinder length a as

a = ℓ − w.   (3)

     Fig. 2 shows the validity of modeling the segmentation masks as a spherocylinder in two dimensions. Here, the thin colored lines are the contours of a collection of segmentation masks and the black solid lines are spherocylinders using the average length and width of the segmentation masks, as calculated by Eq. 2 and Eq. 3. It appears that this simple approximation is reasonable for the purposes of this work.

Counting Repressors

In this section, we expand upon the theoretical and experimental implementation of the fluorescence calibration method derived in Rosenfeld et al. (2005). We cover several experimental data validation steps as well as details regarding the parameter inference. Finally, we comment on the presence of a systematic error in the repressor counts due to continued asynchronous division between sample preparation and imaging.

Theoretical Background of the Binomial Partitioning Method

A key component of this work is the direct measurement of the repressor copy number in each growth condition using fluorescence microscopy. To translate between absolute fluorescence and protein copy number, we must be able to estimate the average brightness of a single fluorophore or, in other words, determine a calibration factor α that permits translation from copy number to intensity or vice versa. Several methods have been used over the past decade to estimate this factor, such as measuring single-molecule photobleaching steps (Garcia et al. 2011a; Bialecka-Fornal et al. 2012), measurement of in vivo photobleaching rates (Nayak and Rutenberg 2011; Kim et al. 2016), and through measuring the partitioning of fluorescent molecules between sibling cells after cell division (Rosenfeld et al. 2005, 2006; Brewster et al. 2014). In this work, we used the latter method to estimate the brightness of a single LacI-mCherry dimer. Here, we derive a simple expression which allows the determination of α from measurements of the fluorescence intensities of a collection of sibling cells.

     In the absence of measurement error, the fluorescence intensity of a given cell is proportional to the total number of fluorescent proteins N_prot by some factor α,
I_cell = αN_prot.   (4)
Assuming that no fluorophores are produced or degraded over the course of the cell cycle, the fluorescent proteins will be partitioned into the two siblings cells such that the intensity of each sibling can be computed as
I₁ = αN₁ ; I₂ = α(N _tot − N₁),   (5)
where N _tot is the total number of proteins in the parent cell and N₁ and N₂ correspond to the number of proteins in sibling cells 1 and 2, respectively. This explicitly states that fluorescence is conserved upon a division,
I_tot = I₁ + I₂.   (6)
As the observed intensity is directly proportional to the number of proteins per cell, measuring the variance in intensity between sibling cells provides some information as to how many proteins were there to begin with. We can compute these fluctuations as the squared intensity difference between the two siblings as
⟨(I₁−I₂)²⟩ = ⟨(2I₁−I _tot)²⟩.   (7)
We can relate Eq. 7 in terms of the number of proteins using Eq. 4 as
⟨(I₁−I₂)²⟩ = 4α²⟨N₁²⟩ − 4α²⟨N₁⟩N _tot + α²N _tot²,   (8)
where the squared fluctuations are now cast in terms of the first and second moment of the probability distribution for N₁.

     Without any active partitioning of the proteins into the sibling cells, one can model the probability distribution g(N₁) of finding N₁ proteins in sibling cell 1 as a binomial distribution,
$$ g(N_1\vert\, N_\text{ tot}, p) = \frac{N_\text{ tot}!}{N_1!(N_\text{ tot} - N_1)!}p^{N_1}(1 - p)^{N_\text{ tot} - N_1}, \qquad(9)$$
where p is the probability of a protein being partitioned into one sibling over the other. With a probability distribution for N₁ in hand, we can begin to simplify Eq. 8. Recall that the mean and variance of a binomial distribution are
⟨N₁⟩ = N _totp,   (10)
and
⟨N₁²⟩ − ⟨N₁⟩² = N _totp(1 − p),   (11)
respectively. With knowledge of the mean and variance, we can solve for the second moment as
⟨N₁²⟩ = N _totp(1 − p) + N _tot²p² = N _totp(1−p+N _totp).   (12)
By plugging Eq. 10 and Eq. 12 into our expression for the fluctuations (Eq. 8), we arrive at
$$ \langle \left( I_1 - I_2\right)^2 \rangle = 4\alpha^2\left[(N_\text{ tot}p[1 - p - N_\text{ tot}p]) - N_\text{ tot}^2\left(p + \frac{1}{4}\right)\right], \qquad(13)$$
which is now defined in terms of the total number of proteins present in the parent cell. Assuming that the proteins are equally partitioned p = 1/2, Eq. 13 reduces to
⟨(I₁−I₂)²⟩ = α²N _tot = αI _tot.   (14)
Invoking our assertion that fluorescence is conserved (Eq. 6), I_tot is equivalent to the sum total fluorescence of the siblings,
⟨(I₁−I₂)²⟩ = α(I₁ + I₂).   (15)
Thus, given snapshots of cell intensities and information of their lineage, one can compute how many arbitrary fluorescence units correspond to a single fluorescent protein.

Cell Husbandry and Time-Lapse Microscopy

The fluorescence calibration method was first described and implemented by Rosenfeld et al. (2005) followed by a more in-depth approach on the statistical inference of the calibration factor in Rosenfeld et al. (2006). In both of these works, the partitioning of a fluorescent protein was tracked across many generations from a single parent cell, permitting inference of a calibration factor from a single lineage. Brewster et al. (2014) applied this method in a slightly different manner by quantifying the fluorescence across a large number of single division events. Thus, rather than examining the partitioning of fluorescence down many branches of a single family tree, it was estimated from an array of single division events where the fluorescence intensity of the parent cell was variable. In the present work, we take a similar approach to that of Brewster et al. (2014), and examine the partitioning of fluorescence among a large number of independent cell divisions.

     A typical experimental work-flow is shown in Fig. 3. For each experiment, the strains were grown in varying concentrations of ATC to tune the expression of the repressor. Once the cells had reached exponential phase growth (OD_600nm ≈ 0.3), the cells were harvested and prepared for imaging. This involved two separate sample handling procedures, one for preparing samples for lineage tracking and estimation of the calibration factor and another for taking snapshots of cells from each ATC induction condition for the calculation of fold-change.

     To prepare cells for the calibration factor measurement, a 100 μL aliquot of each ATC induction condition was combined and mixed in a 1.5 mL centrifuge tube. This cell mixture was then pelleted at 13000 × g for 1 – 2 minutes. The supernatant containing ATC was then aspirated and the pellet was resuspended in an equal volume of sterile growth media without ATC. This washing procedure was repeated three times to ensure that any residual ATC had been removed from the culture and that expression of LacI-mCherry had ceased. Once washed and resuspended, the cells were diluted ten fold into sterile M9 medium and then imaged on a rigid agarose substrate. Depending on the precise growth condition, a variety of positions were imaged for 1.5 to 4 hours with a phase contrast image acquired every 5 to 15 minutes to facilitate lineage tracking. On the final image of the experiment, an mCherry fluorescence image was acquired of every position. The experiments were then transferred to a computing cluster and the images were computationally analyzed, as described in the next section.

     During the washing steps, the remaining ATC induced samples were prepared for snapshot imaging to determine the repressor copy number and fold-change in gene expression for each ATC induction condition. Without mixing the induction conditions together, each ATC induced sample was diluted 1:10 into sterile M9 minimal medium and vigorously mixed. Once mixed, a small aliquot of the samples were deposited onto rigid agarose substrates for later imaging. While this step of the experiment was relatively simple, the total preparation procedure typically lasted between 30 and 60 minutes. As is discussed later, the continued growth of the asynchronously growing culture upon dilution into the sterile medium results in a systematic error in the calculation of the repressor copy number.

Lineage Tracking and Fluorescence Quantification

Segmentation and lineage tracking of both the fluorescence snap shots and time-lapse growth images were performed using the SuperSegger v1.03 (Stylianidou et al. 2016) software using MATLAB R2017B (MathWorks, Inc). The result of this segmentation is a list of matrices for each unique imaged position with identifying data for each segmented cell such as an assigned ID number, the ID of the sibling cell, the ID of the parent cell, and various statistics. These files were then analyzed using Python 3.7. All scripts and software used to perform this analysis can be found on the associated paper website and GitHub repository.

     Using the ID numbers assigned to each cell in a given position, we matched all sibling pairs present in the last frame of the growth movie when the final mCherry fluorescence image was acquired. These cells were then filtered to exclude segmentation artifacts (such as exceptionally large or small cells) as well as any cells which the SuperSegger software identified as having an error in segmentation. Given the large number of cells tracked in a given experiment, we could not manually correct these segmentation artifacts, even though it is possible using the software. To err on the side of caution, we did not consider these edge cases in our analysis.

     With sibling cells identified, we performed a series of validation checks on the data to ensure that both the experiment and analysis behaved as expected. Three validation checks are illustrated in Fig. 4. To make sure that the computational pairing of sibling cells was correct, we examined the intensity distributions of each sibling pair. If siblings were being paired solely on their lineage history and not by other features (such as size, fluorescence, etc), one would expect the fluorescence distributions between the two sibling cells to be identical. Fig. 4(A) shows the nearly identical intensity distributions of all siblings from a single experiment, indicating that the pairing of siblings is independent of their fluorescence.

     Furthermore, we examined the partitioning of the fluorescence between the siblings. In Eq. 10, we defined the mean number of proteins inherited by one sibling. This can be easily translated into the language of intensity as
⟨I₁⟩ = α(I₁+I₂)p,   (16)
where we assume that fluorescence is conserved during a division event such that I _tot = I₁ + I₂. As described previously, we make the simplifying assumption that partitioning of the fluorophores between the two sibling cells is fair, meaning that p = 0.5. We can see if this approximation is valid by computing the fractional intensity of each sibling cell as
$$ p = \frac{\langle I_1 \rangle}{\langle I_1 + I_2 \rangle}. \qquad(17)$$
Fig. 4(B) shows the distribution of the fractional intensity for each sibling pair. The distribution is approximately symmetric about 0.5, indicating that siblings are correctly paired and that the partitioning of fluorescence is approximately equal between siblings. Furthermore, we see no correlation between the cell volume immediately after division and the observed fractional intensity. This suggests that the probability of partitioning to one sibling or the other is not dependent on the cell size. An assumption (backed by experimental measurements (Garcia et al. 2011a; Phillips et al. 2019)) in our thermodynamic model is that all repressors in a given cell are bound to the chromosome, either specifically or nonspecifically. As the chromosome is duplicated and partitioned into the two siblings without fail, our assumption of repressor adsorption implies that partitioning should be independent of the size of the respective sibling. The collection of these validation statistics give us confidence that both the experimentation and the analysis are properly implemented and are not introducing bias into our estimation of the calibration factor.

Statistical Inference of the Fluorescence Calibration Factor

As is outlined in the Materials & Methods section of the Chapter 4, we took a Bayesian approach towards our inference of the calibration factor given fluorescence measurements of sibling cells. Here, we expand in detail on this statistical model and its implementation.

     To estimate the calibration factor α from a set of lineage measurements, we assume that fluctuations in intensity resulting from measurement noise is negligible compared to that resulting from binomial partitioning of the repressors upon cell division. In the absence of measurement noise, the intensity of a given cell I_cell can be directly related to the total number of fluorophores N through a scaling factor α, such that
I_cell = αN.   (18)
Assuming no fluorophores are produced or degraded over the course of a division cycle, the fluorescence of the parent cell before division is equal to the sum of the intensities of the sibling cells,
I_parent = I₁ + I₂ = α(N₁+N₂),   (19)
where subscripts 1 and 2 correspond to arbitrary labels of the two sibling cells. We are ultimately interested in knowing the probability of a given value of α which, using Bayes’ theorem, can be written as
$$ g(\alpha\, \vert\, I_1, I_2) = \frac{f(I_1, I_2\, \vert\, \alpha) g(\alpha)} {f(I_1, I_2)}, \qquad(20)$$
where we have used g and f to denote probability densities over parameters and data, respectively. The first quantity in the numerator f(I₁, I₂ | α) describes the likelihood of observing the data I₁, I₂ given a value for the calibration factor α. The term g(α) captures all prior knowledge we have about what the calibration factor could be, remaining ignorant of the collected data. The denominator f(I₁, I₂) is the likelihood of observing our data I₁, I₂ irrespective of the calibration factor and is a loose measure of how well our statistical model describes the data. As it is difficult to assign a functional form to this term and serves as a multiplicative constant, it can be neglected for the purposes of this work.

     Knowing that the two observed sibling cell intensities are related, conditional probability allows us to rewrite the likelihood as
f(I₁ I₂ | α) = f(I₁ | I₂, α)f(I₂ | α),   (21)
where f(I₂ | α) describes the likelihood of observing I₂ given a value of α. As I₂ can take any value with equal probability, this term can be treated as a constant. Through change of variables and noting that I₂ = αN₂, we can cast the likelihood f(I₁ | I₂, α) in terms of the number of proteins N₁ and N₂ as
$$ f(I_1\,\vert\, I_2, \alpha) = f(N_1\,\vert\,N_2, \alpha)\bigg\vert {d N_1 \over dI_1}\bigg\vert= \frac{1}{\alpha}f(N_1\,\vert\,N_2,\alpha). \qquad(22)$$
Given that the proteins are binomially distributed between the sibling cells with a probability p and that the intensity is proportional to the number of fluorophores, this likelihood becomes
$$ f(I_1\,\vert\,I_2, \alpha, p) = \frac{1}{\alpha}{\left(\frac{I_1 + I_2}{ \alpha}\right)! \over \left({I_1 \over \alpha}\right)!\left({I_2 \over \alpha}\right)!}p^(1 - p)^{I_2 \over \alpha}. \qquad(23)$$
However, the quantity I/α is not exact, making calculation of its factorial undefined. These factorials can therefore be approximated by a gamma function as n! = Γ(n + 1).

     Assuming that partitioning of a protein between the two sibling cells is a fair process (p = 1/2), Eq. 23 can be generalized to a set of lineage measurements [I₁, I₂] as
$$ f([I_1]\,\vert\,[I_2], \alpha) = {1 \over \alpha^k} \prod\limits_{i}^k {\Gamma\left({I_1 + I_2 \over \alpha} + 1\right) \over \Gamma\left({I_1 \over \alpha} + 1\right)\Gamma\left({I_2 \over \alpha} + 1 \right)}2^{-{I_1 + I_2 \over \alpha}}, \qquad(24)$$
where k is the number of division events observed.

     With a likelihood in place, we can now assign a functional form to the prior distribution for the calibration factor g(α). Though ignorant of data, the experimental design is such that imaging of a typical highly-expressing cell will occupy 2/3 of the dynamic range of the camera. We can assume that it is more likely that the calibration factor will be closer to 0 a.u. than the bit depth of the camera (4095 a.u.) or larger. We also know that it is physically impossible for the fluorophore to be less than 0 a.u., providing a hard lower-bound on its value. We can therefore impose a weakly informative prior distribution as a half normal distribution,
$$ g(\alpha) = \sqrt{2 \over \pi\sigma^2}\exp\left[{-\alpha^2 \over 2\sigma^2}\right]\,;\, \forall \alpha > 0, \qquad(25)$$
where the standard deviation is large, for example, σ = 500 a.u. / fluorophore. We evaluated the posterior distribution using Markov chain Monte Carlo (MCMC) as is implemented in the Stan probabilistic programming language (Carpenter et al. 2017). The .stan file associated with this model along with the Python code used to execute it can be accessed on the paper website and GitHub repository. Fig. 5 shows the posterior probability distributions of the calibration factor estimated for several biological replicates of the glucose growth condition at 37^∘ C. For each posterior distribution, the mean and standard deviation was used as the calibration factor and uncertainty for the corresponding data set.

Correcting for Systematic Experimental Error

While determination of the calibration factor relies on time-resolved measurement of fluorescence partitioning, we computed the repressor copy number and fold-change in gene expression from still snapshots of each ATC induction condition and two control samples, as is illustrated in Fig. 3. Given these snapshots, individual cells were segmented again using the SuperSegger software in MATLAB R2017b. The result of this analysis is an array of single-cell measurements of the YFP and mCherry fluorescence intensities. With these values and a calibration factor estimated by Eq. 24 and Eq. 25, we can compute the estimated number of repressors per cell in every condition.

     However, as outlined previously and in Fig. 3, the sample preparation steps for these experiments involve several steps which require careful manual labor. This results in an approximately 30 to 60 minute delay from when production of the LacI-mCherry construct is halted by removal of ATC to actual imaging on the microscope. During this time, the diluted cultures are asynchronously growing, meaning that the cells of the culture are at different steps in the cell cycle. Thus, at any point in time, a subset of cells will be on the precipice of undergoing division, partitioning the cytosolic milieu between the two progeny. As LacI-mCherry is no longer being produced, the cells that divided during the dwell time from cell harvesting to imaging will have reduced the number of repressors by a factor of 2 on average. This principle of continued cell division is shown Fig. 6.

     How does this partitioning affect our calculation and interpretation of the fold-change in gene expression? Like the LacI-mCherry fusion, the YFP reporter proteins are also partitioned between the progeny after a division event such that, on average, the total YFP signal of the newborn cells is one-half that of the parent cell. As the maturation time of the mCherry and YFP variants used in this work are relatively long in E. coli (Balleza, Kim, and Cluzel 2018; Nagai et al. 2002), we can make the assumption that any newly-expressed YFP molecules after cells have divided are not yet visible in our experiments. Thus, the fold-change in gene expression of the average parent cell can be calculated given knowledge of the average expression of the progeny.

     The fold-change in gene expression is a relative measurement to a control which is constitutively expressing YFP. As the latter control sample is also asynchronously dividing, the measured YFP intensity of the newborn constitutively expressing cells is on average 1/2 that of the parent cell. Therefore, the fold-change in gene expression can be calculated as
$$ \langle \text{ fold-change}_{parent}\rangle = \frac{2 \times \langle I^\text{(YFP)}_\text{ newborn}(R > 0) \rangle}{2 \times \langle I^\text{(YFP)}_\text{ newborn}(R = 0)\rangle} = \frac{\langle I^\text{(YFP)}_\text{ newborn}(R > 0) \rangle}{\langle I^\text{(YFP)}_\text{ newborn}(R = 0)\rangle}. \qquad(26)$$
Thus, when calculating the fold-change in gene expression one does not need to correct for any cell division that occurs between sample harvesting and imaging as it is a relative measurement. However, the determination of the repressor copy number is a direct measurement and requires a consideration of unknown division events.

     To address this source of systematic experimental error, we examined the cell length distributions of all segmented cells from the snapshots as well as the distribution of newborn cell lengths from the time-lapse measurements. Fig. 7 shows that a significant portion of the cells from the snapshots (colored distributions) overlap with the distribution of newborn cell lengths (grey distributions). For each condition, we partitioned the cells from the snapshots into three bins based on their lengths – “small” cells had a cell length less than 2.5 μm, “medium” cells had lengths between 2.5 μm and 3.5 μm, and “large” cells being longer than 3.5 μm. These thresholds were chosen manually by examining the newborn cell-size distributions and are shown as red vertical lines in Fig. 7. Under this partitioning, we consider all “small” cells to have divided between cessation of LacI-mCherry production and imaging, “medium”-length cells to have a mixture of long newborn cells (from the tail of the newborn cell length distribution) and cells that have not divided, and cells in the "long" group to be composed entirely of cells which did not undergo a division over the course of sample preparation.

     Given this coarse delineation of cell age by length, we examined how correction factors could be applied to correct for the undesired systematic error due to dilution of repressors. We took the data collected in this work and compared the results to the fold-change in gene expression reported in the literature for the same regulatory architecture. Without correcting for undesired cell division, the observed fold-change in gene expression falls below the prediction and does not overlap with data from the literature (Fig. 8(A), light purple). Using the uncorrected measurements, we estimated the DNA binding energy to be Δε_R ≈  − 15 k_BT which does not agree with the value for the O2 operator reported in Garcia et al. (2011a) or with the inferred DNA binding energies from the other data sources (Fig. 8 (B)).

      These results emphasize the need to correct for undesired dilution of repressors through cell division during the sample preparation period, and we now consider several different manners of applying this correction. We first consider the extreme case where all cells of the culture underwent an undesired division after LacI-mCherry production was halted. This means that the average repressor copy number measured from all cells is off by a factor of 2. The result of applying a factor of 2 correction to all measurements can be seen in Fig. 8 as dark red points. Upon applying this correction, we find that the observed fold-change in gene expression agrees with the prediction and data from other sources in the literature. The estimated DNA binding energy Δε_R from these data is also in agreement with other data sources (Fig. 8(B), dark red). This result suggests that over the course of sample preparation, a non-negligible fraction of the diluted culture undergoes a division event before being imaged.

     We now begin to relax assumptions as to what fraction of the measured cells underwent a division event before imaging. As described above, in drawing distinctions between “small,” “medium,” and “large” cells, we assume the latter represent cells which did not undergo a division between the harvesting and imaging of the samples. Thus, the repressor counts of these cells should require no correction. The white-faced points in Fig. 8(A) shows the fold-change in gene expression of only the large cell fraction, which falls within error of the theoretical prediction. Furthermore, the inferred DNA binding energy falls within error of that inferred from data of Garcia et al. (2011a) and that inferred from data assuming all cells underwent a division event Fig. 8(B), though it does not fall within error of the binding energy reported in Garcia et al. (2011a).

     The most realistic approach that can be taken to avoid using only the “large” bin of cells is to assume that all cells with a length ℓ < 2.5 μm have undergone a division, requiring a two-fold correction to their average repressor copy number. The result of this approach can be seen in Fig. 8 as purple points. The inferred DNA binding energy from this correction approach falls within error of that inferred from only the large cells white-faced points in 8(B) as well as overlapping with the estimate treating all cells as having undergone a division (red).

     There are several experimental techniques that could be implemented to avoid needing to apply a correction factor as described here. In Brewster et al. (2014), the fold-change in gene expression was measured by tracking the production rate of a YFP reporter before and after a single cell division event coupled with direct measurement of the repressor copy number using the same binomial partitioning method. This implementation required an extensive degree of manual curation of segmentation as well as correcting for photobleaching of the reporter, which in itself is a non-trivial correction (Garcia et al. 2011a). The experimental approach presented here sacrifices a direct measure of the repressor copy number for each cell via the binomial partitioning method, but permits the much higher throughput needed to assay the variety of environmental conditions. Ultimately, the inferred DNA binding energy for all of the scenarios described above agree within 1 k_BT, a value smaller than the natural variation in the DNA binding energies of the three native lacO, which is  ≈ 6 k_BT. For the purposes of this work, we erred on the side of caution and only used the cells deemed “large” for the measurements reported in Chapter 4 and the remaining sections of this chapter.

Parameter Estimation of DNA Binding Energies and Comparison Across Carbon Source

In Chapter 4, we conclude that the biophysical parameters defining the fold-change input-output function are unperturbed between different carbon sources. This conclusion is reached primarily by comparing how well the fold-change and the free energy shift ΔF is predicted using the parameter values determined in glucose supported medium at 37^∘ C. In this section, we redetermine the DNA binding energy for each carbon source condition and test its ability to predict the fold-change of the other conditions.

Reparameterizing the Fold-Change Input-Output Function

As described previously, the fold-change in gene expression is defined by the total repressor copy number R, the energetic difference between the active and inactive states of the repressor Δε_AI, and the binding energy of the repressor to the DNA Δε_RA. Using fluorescence microscopy, we can directly measure the average repressor copy number per cell, reducing the number of variable parameters to only the energetic terms.

Estimating both Δε_RA and Δε_AI simultaneously is fraught with difficulty as the parameters are highly degenerate (Razo-Mejia et al. 2018). We can avoid this degeneracy by reparameterizing the fold-change input-output function as
$$ \text{fold-change} =\left( 1 + \frac{R}{N_{NS}}e^{-\beta\epsilon}\right)^{-1}, \qquad(27)$$
where ϵ is the effective energetic parameter
ϵ = Δε_R − k_BTlog (1+e^{− βΔε_AI}).   (28)
Thus, to further elucidate any changes to the parameter values due to changing the carbon source, we can infer the best-fit value of ϵ for each condition and explore how well it predicts the fold-change in other conditions.

Statistical Inference of ϵ

We are interested in the probability distribution of the parameter ϵ given knowledge of the repressor copy number R and a collection of fold-change measurements fc. This quantity can be calculated via by Bayes’ theorem as
$$ g(\epsilon\,\vert\,R, \mathbf{fc}) = \frac{f(\mathbf{fc}\,\vert\, R, \epsilon)g(\epsilon)}{f(\mathbf{fc})}, \qquad(29)$$
where g and f represent probability densities of parameters and data, respectively. In this context, the denominator term f(fc) serves only as normalization constant and can be neglected. As is described in detail in Refs. (Razo-Mejia et al. 2018; Chure et al. 2019), we assume that a given set of fold-change measurements are normally distributed about the theoretical value μ defined by the fold-change function. The likelihood can be mathematically defined as
$$ f(\mathbf{fc}\,\vert\,\epsilon, R) = \frac{1}{\left(2\pi\sigma^2\right)^{N/2}}\prod\limits_i^N \exp\left[\frac{-\left(\text{fc}_i - \mu(\epsilon, R)\right)^2}{2\sigma^2}\right], \qquad(30)$$
where N is the total number of fold-change measurements, and σ is the standard deviation of the observations about the true mean and is another parameter that must be included in the estimation. As the fold-change in gene expression in this work covers several orders of magnitude (from  ≈ 10^− 3 − 10⁰), it is better to condition the parameters on the log fold-change rather than linear scaling, translating Eq. 30 to
$$ f(\mathbf{fc}^*\,\vert\, \epsilon, R) = \frac{1}{\left(2\pi\sigma\right)^{N/2}}\prod\limits_i^N\exp\left[-\frac{\left(\text{fc}^*_i - \mu^*(\epsilon, R)\right)^2}{2\sigma^2}\right], \qquad(31)$$
where fc^* and μ^*(ϵ, R) are the transformations
fc^* = log (fc)
and
$$\mu^*(\epsilon, R) = -\log\left(1 + \frac{R}{N_{NS}}e^{-\beta\epsilon}\right), \qquad(32)$$
respectively.

     With a likelihood in hand, we can now turn towards defining functional forms for the prior distributions g(σ) and g(ϵ). For these definitions, we can turn to those used in Chure et al. (2019),
g(ϵ) ∼ Normal(μ =  − 12, σ = 6)   (33)
and
g(σ) ∼ HalfNormal(ϕ = 0.1),   (34)
where we introduce the shorthand notation of “Normal” and “HalfNormal.” Combining Eq. 31 and Eq. 34 - Eq. 33 yields the complete posterior distribution for estimating the DNA binding energy for each carbon source medium. The complete posterior distribution was sampled using Markov chain Monte Carlo in the Stan probabilistic programming language (Carpenter et al. 2017).

     The sampled posterior distributions for ϵ and σ for each carbon source condition are shown in Fig. 9 and are summarized in Table 8.1. The posterior distributions of ϵ across the conditions are approximately equal with highly overlapping 95% credible regions. The predictive capacity of each estimate of ϵ is shown in Fig. 10 where all fold-change measurements fall upon the theoretical prediction regardless of which carbon source that the parameter value was conditioned upon. With this analysis, we can say with quantitative confidence that the biophysical parameters are indifferent to the physiological changes resulting from variation in carbon quality.

Summarized parameter estimates of ϵ and σ given a single growth condition. Reported as median and upper/lower bounds of 95% credible region.

Growth Condition	Parameter	Value
Glucose, 37∘ C	ϵ	− 14.5 − 0.3 + 0.2 kBT
	σ	0.3 − 0.1 + 0.1
Glycerol, 37∘ C	ϵ	− 14.6 − 0.1 + 0.1 kBT
	σ	0.39 − 0.08 + 0.10
Acetate, 37∘ C	ϵ	− 14.1 − 0.3 + 0.2 kBT
	σ	0.35 − 0.1 + 0.1

Statistical Inference of Entropic Costs

In the main text, we describe how a simple rescaling of the energetic parameters Δε_R and Δε_AI is not sufficient to describe the fold-change in gene expression when the growth temperature is changed from 37^∘ C. In this section, we describe the inference of hidden entropic parameters to phenomenologically describe the temperature dependence of the fold-change in gene expression.

Definition of Hidden Entropic Costs

The values of the energetic parameters Δε_R and Δε_AI were determined in a glucose supplemented medium held at 37^∘ C which we denote as T_ref. A null model to describe temperature dependence of these parameters is to rescale them to the changed temperature T_exp as
$$ \Delta\varepsilon^* = \frac{T_{ref}}{T_{exp}}\Delta\varepsilon, \qquad(35)$$
where Δε^* is either Δε_R or Δε_AI. However, we found that this null model was not sufficient to describe the fold-change in gene expression, prompting the formulation of a new phenomenological description.

      Our thermodynamic model for the fold-change in gene expression coarse-grains the regulatory architecture to a two-state model, meaning many of the rich features of regulation such as vibrational entropy, the material properties of DNA, and the occupancy of the repressor to the DNA are swept into the effective energetic parameters. As temperature was never perturbed when this model was developed, modeling these features was not necessary. However, we must now return to these features to consider what may be affected.

     Without assigning a specific mechanism, we can say that there is a temperature-dependent entropic parameter that was neglected in the estimation of the energetic parameters in Garcia et al. (2011a) and Razo-Mejia et al. (2018). In this case, the inferred energetic parameter Δε^* is composed of enthalpic (ΔH) and entropic (ΔS) parameters,
Δε^* = ΔH − TΔS.
For a set of fold-change measurements at a temperature T_exp, we are interested in estimating values for ΔH and ΔS for each energetic parameter. Given measurements from Refs. (Garcia et al. 2011a; Razo-Mejia et al. 2018), we know at 37^∘C what Δε_RA and Δε_AI are inferred to be, placing a constraint on the possible values of ΔH and ΔS,
ΔH_R = Δε_R + T_refΔS_R = T_refΔS_R − 13.9 k_BT,   (36)
and
ΔH_AI = Δε_AI + T_refΔS_AI = T_refΔS_AI + 4.5 k_BT   (37)
for the DNA binding energy and allosteric state energy difference, respectively.

Statistical Inference of ΔS_R and ΔS_AI

Given the constraints from Eq. 36 and Eq. 37, we are interested in inferring the entropic parameters ΔS_R and ΔS_AI given literature values for Δε_RA and Δε_AI and the set of fold-change measurements fc at a given temperature T_exp. The posterior probability distribution for the entropic parameters can be enumerated via Bayes’ theorem as
$$ g(\Delta S_R, \Delta S_{AI}\,\vert\, \mathbf{fc}) = \frac{f(\mathbf{fc}\,\vert\,\Delta S_R, \Delta S_{AI})g(\Delta S_R, \Delta S_{AI})}{f(\mathbf{fc})}, \qquad(38)$$
where g and f are used to denote probability densities over parameters and data, respectively. As we have done elsewhere in this SI text, we treat f(fc) as a normalization constant and neglect it in our estimation of g(ΔS_R, ΔS_AI | fc). Additionally, as is discussed in detail earlier in this chapter, we consider the log  fold-change measurements to be normally distributed about a mean fc^* defined by the fold-change input-output function and a standard deviation σ. Thus, the likelihood for the fold-change in gene expression is
$$ \begin{aligned} f(\mathbf{fc}\,\vert\,\Delta S_R, \Delta S_{AI}, \sigma, T_{exp}) &= \frac{1}{\left(2\pi\sigma^2\right)^{N/2}}\times \\ &\prod\limits_i^N \exp\left[-\frac{\left[\log \text{fc}_i - \log \text{fc}^*(\Delta S_R, \Delta S_{AI}, T_{exp})\right]^2}{2\sigma^2}\right]. \end{aligned} \qquad(39)$$

     In calculating the mean fc^*, the effective energetic parameters Δε_R^* and Δε_AI^* can be defined and constrained using Eq. 36 and Eq. 37 as
Δε_R^* = ΔH_R − T_expΔS_R = ΔS_R(T_ref − T_exp) − 13.9 k_BT_ref,   (40)
and
Δε_AI^* = ΔH_AI − T_expΔS_AI = ΔS_AI(T_ref − T_exp) + 4.5 k_BT_ref.   (41)

     As the enthalpic parameters are calculated directly from the constraints of Eq. 36 and Eq. 37, we must only estimate three parameters, ΔS_R, ΔS_AI, and σ, each of which need a functional form for the prior distribution.

     A priori, we know that both ΔS_R and ΔS_AI must be small because T_ref and T_exp are defined in K. As these entropic parameters can be either positive or negative, we can define the prior distributions g(ΔS_R) and g(ΔS_AI) as a normal distribution centered at zero with a small standard deviation,
$$ g(\Delta S_R) \sim \text{\itshape Normal}(\mu=0, \sigma=0.1), \qquad(42)$$
and
$$ g(\Delta S_{AI}) \sim \text{\itshape Normal}(\mu=0, \sigma=0.1). \qquad(43)$$
The standard deviation σ can be defined as a half-normal distribution centered at 0 with a small standard deviation ϕ,
$$ g(\sigma) \sim \text{\itshape HalfNormal}(\phi=0.1). \qquad(44)$$

     With the priors and likelihood functions in hand, we sampled the posterior distribution using Markov chain Monte Carlo as implemented in the Stan probabilistic programming language (Carpenter et al. 2017). We performed three different estimations – one inferring the parameters using only data at 32^∘ C, one using only data from 42^∘ C, and one using data sets from both temperatures pooled together.

     The sampling results can be seen in Fig. 11. The estimation of ΔS_R is distinct for each condition whereas the sampling for ΔS_AI is the same for all conditions and is centered at about 0. The latter suggests that the value of Δε_AI determined at 37^∘ C is not dependent on temperature within the resolution of our experiments. The difference between the estimated value of ΔS_R between temperatures suggests that there is another component of the temperature dependence that is not captured by the inclusion of a single entropic parameter. Fig. 12 shows that estimating ΔS from one temperature is not sufficient to predict the fold-change in gene expression at another temperature. The addition of the entropic parameter leads to better fit of the 32^∘ C condition than the simple rescaling of the energy as described by Eq. 35 (Fig. 12, dashed line), but poorly predicts the behavior at 42^∘C. Performing the inference on the combined 32^∘ C and 42^∘ C data strikes a middle ground between the predictions resulting from the two temperatures alone Fig. 12, grey shaded region.

Entropy as a Function of Temperature

In the previous section, we made an approximation of the energetic parameters Δε_R and Δε_AI to be defined by an enthalpic and entropic term, both of which being independent of temperature. However, entropy can be (and in many cases is) dependent on the system temperature, often in a non-trivial manner.

     To explore the effects of a temperature-dependent entropy in our prediction of the fold-change, we perform a Taylor expansion of Eq. ¿eq:utds? about the entropic parameter with respect to temperature keeping only the first order term such that
ΔS = ΔS₀ + ΔS₁T,   (45)
where ΔS₀ is a constant, temperature-independent entropic term and ΔS₁ is the entropic contribution per degree Kelvin. With this simple relationship enumerated, we can now define the temperature-dependent effective free energy parameters Δε_R^* and Δε_AI^* as
Δε_R^* = (S_0R + S_1RT_exp)(T_ref − T_exp) − 13.9 k_BT,   (46)
and
Δε_AI^* = (S_0AI + S_1AIT_exp)(T_ref − T_exp) + 4.5 k_BT,   (47)
respectively, again relying on the constraints defined by Eq. 40 and Eq. 41.

      Using a similar inferential approach as described in the previous section, we sample the posterior distribution of these parameters using Markov chain Monte Carlo and compute the fold-change and shift in free energy for each temperature. As seen in Fig. 13, there is a negligible improvement in the description of the data by including this temperature dependent entropic parameter.

     These results together suggest that our understanding of temperature dependence in this regulatory architecture is incomplete and requires further research from both theoretical and experimental standpoints.

Media Recipes and Bacterial Strains

The primary interest in varying the available carbon source in growth media was to modulate the quality of the carbon rather than the quantity. With this in mind, we developed the various growth media to contain the same net number of carbon atoms per cell. The standard reference was 0.5% (w/v) glucose (Garcia et al. 2011a), which results in 10⁸ carbon atoms per 10^− 15 L. The base recipe is given in Table 8.2. The bacterial strains used in this work are given in Table 8.3.

M9 minimal medium recipe for each carbon-supplemented medium.

Ingredient	[Concentration]	Volume	[Final Concentration]
ddH2O	-	773 mL	-
CaCl2	1M	100 μL	100 μM
MgSO4	1M	2mL	2mM
M9 Salts	5X	200 mL	-
(BD Medical, Cat. No. 248510)
Carbon Source			108 C / fL
Glucose	20% (w/v)	25 mL	0.5% (w/v)
Glycerol	20% (w/v)	25 mL	0.5% (w/v)
Acetate	20% (w/v)	25 mL	0.5% (w/v)

Bacterial strains used in various physiological states.

Genotype	Plasmid	Notes
MG1655:ΔlacZYA; intC<>4*CFP	–	Autofluorescence control
MG1655:ΔlacZYA; intC<>4*CFP
galK<>25-O2+11-YFP	–	Constitutive expression control
MG1655:ΔlacZYA;intC<>4*CFP
galK<>25-O2+11-YFP
ybcN<>1-lacI(Δ 353-363)-mCherry	pZS3PN25-tetR	Strain with ATC inducible lacI-mCherry

References