Through The Intramolecular Grapevine: Cellular Decision Making Via Allosteric Transcription factors

Published as …

---

A version of this chapter originally appeared as Razo-Mejia, M., Barnes, S.L., Belliveau, N.M., Chure, G., Einav, T.*, Lewis, M., and Phillips, R. (2018). Tuning transcriptional regulation through signaling: A predictive theory of allosteric induction. Cell Systems 6, 456-469.e10. DOI:https://doi.org/10.1016/j.cels.2018.02.004. M.R.M, S.L.B, N.M.B, G.C., and T.E. contributed equally to this work from the theoretical underpinnings to the experimental design and execution. M.R.M, S.L.B, N.M.B, G.C, T.E., and R.P. wrote the paper. M.L. provided extensive guidance and advice.

Abstract

Allosteric regulation is found across all domains of life, yet we still lack simple, predictive theories that directly link the experimentally tunable parameters of a system to its input-output response. To that end, we present a general theory of allosteric transcriptional regulation using the Monod-Wyman- Changeux model. We rigorously test this model using the ubiquitous simple repression motif in bacteria by first predicting the behavior of strains that span a large range of repressor copy numbers and DNA binding strengths, and then constructing and measuring their response. Our model not only accurately captures the induction profiles of these strains, but also enables us to derive analytic expressions for key properties such as the dynamic range and [EC₅₀]. Finally, we derive an expression for the free energy of allosteric repressors that enables us to collapse our experimental data onto a single master curve that captures the diverse phenomenology of the induction profiles.

Introduction

     Understanding how organisms sense and respond to changes in their environment has long been a central theme of biological inquiry. At the cellular level, this interaction is mediated by a diverse collection of molecular signaling pathways. A pervasive mechanism of signaling in these pathways is allosteric regulation, in which the binding of a ligand induces a conformational change in some target molecule, triggering a signaling cascade (Lindsley and Rutter 2006). One of the most important examples of such signaling is offered by transcriptional regulation, where a transcription factors’ propensity to bind to DNA will be altered upon binding to an allosteric effector.

     Despite the ubiquity of allostery, we largely lack a formal, rigorous, and generalizable framework for studying its effects across the broad variety of contexts in which it appears. A key example of this is transcriptional regulation, in which allosteric transcription factors can be induced or corepressed by binding to a ligand. An allosteric transcription factor can adopt multiple conformational states, each of which has its own affinity for the ligand and for its DNA target site. In vitro studies have rigorously quantified the equilibria of different conformational states for allosteric transcription factors and measured the affinities of these states to the ligand (Harman 2001; Lanfranco et al. 2017). In spite of these experimental observations, the lack of a coherent quantitative model for allosteric transcriptional regulation has made it impossible to predict the behavior of even a simple genetic circuit across a range of regulatory parameters, physiological states of the organism, and evolutionary isoforms of the regulatory sequences.

     The ability to predict circuit behavior robustly— that is, across both broad ranges of parameters and regulatory architectures —is important for multiple reasons. First, in the context of a specific gene, accurate prediction demonstrates that all components relevant to the genes’ behavior have been identified and characterized to sufficient quantitative precision. Second, in the context of genetic circuits in general, robust prediction validates the model that generated the prediction. Possessing a validated model also has implications for future work. For example, when we have sufficient confidence in the model, a single data set can be used to accurately extrapolate a system’s behavior in other conditions. Moreover, there is an essential distinction between a predictive model, which is used to predict a system’s behavior given a set of input variables, and a retroactive model, which is used to describe the behavior of data that has already been obtained. We note that even some of the most careful and rigorous analysis of transcriptional regulation often entails only a retroactive reflection on a single experiment. This raises the fear that each regulatory architecture may require a unique analysis that cannot carry over to other systems, a worry that is exacerbated by the prevalent use of phenomenological functions (e.g. Hill functions or ratios of polynomials) that can analyze a single data set, but cannot be used to extrapolate a system’s behavior in other conditions (Setty et al. 2003; Poelwijk et al. 2011; Vilar and Saiz 2013; Rogers et al. 2015; Rohlhill, Sandoval, and Papoutsakis 2017).

     This work explores what happens when theory takes center stage, namely, we first write down the equations governing a system and describe its expected behavior across a wide array of experimental conditions, and only then do we set out to experimentally confirm these results. Building upon previous work (Garcia et al. 2011a; Brewster et al. 2014; Weinert et al. 2014) and the work of Monod, Wyman, and Changeux (Monod, Wyman, and Changeux 1965), we present a statistical mechanical rendering of allostery in the context of induction and corepression (shown schematically in Fig. 1 and henceforth referred to as the MWC model), and use it as the basis of parameter-free predictions which we then test experimentally. More specifically, we study the simple repression motif – a widespread bacterial genetic regulatory architecture in which binding of a transcription factor occludes binding of an RNA polymerase, thereby inhibiting transcription initiation. The MWC model stipulates that an allosteric protein fluctuates between two distinct conformations – an active and inactive state – in thermodynamic equilibrium (Monod, Wyman, and Changeux 1965). During induction, for example, effector binding increases the probability that a repressor will be in the inactive state, weakening its ability to bind to the promoter and resulting in increased expression. To test the predictions of our model across a wide range of operator binding strengths and repressor copy numbers, we design a genetic construct in Escherichia coli in which the binding probability of a repressor regulates gene expression of a fluorescent reporter.

     In total, the work presented here demonstrates that one extremely compact set of parameters can be applied self-consistently and predictively to different regulatory situations including simple repression on the chromosome, cases in which decoy binding sites for repressor are put on plasmids, cases in which multiple genes compete for the same regulatory machinery, cases involving multiple binding sites for repressor leading to DNA looping, and induction by signaling (Garcia et al. 2011a, 2011a; Brewster, Jones, and Phillips 2012; Brewster et al. 2014; Boedicker et al. 2013; Boedicker, Garcia, and Phillips 2013). Thus, rather than viewing the behavior of each circuit as giving rise to its own unique input-output response, the MWC model provides a means to characterize these seemingly diverse behaviors using a single unified framework governed by a small set of parameters.

Theoretical Model

Inducible Transcriptional Repression Via the MWC Model of Allostery

     We begin by considering a simple repression genetic architecture in which the binding of an allosteric repressor occludes the binding of RNA polymerase (RNAP) to the DNA (Ackers and Johnson 1982; Buchler, Gerland, and Hwa 2003). When an effector molecule (hereafter referred to as an “inducer” for the case of induction) binds to the repressor, it shifts the repressor’s allosteric equilibrium towards the inactive state as specified by the MWC model (Monod, Wyman, and Changeux 1965). This causes the repressor to bind more weakly to the operator, increasing the probability of RNAP binding the promoter which ultimately leads to gene expression. Simple repression motifs in the absence of inducer have been previously characterized by an equilibrium model where the probability of each state of repressor and RNAP promoter occupancy is dictated by the Boltzmann distribution (Ackers and Johnson 1982; Buchler, Gerland, and Hwa 2003; Vilar and Leibler 2003; Bintu, Buchler, Garcia, Gerland, Hwa, Kondev, Kuhlman, et al. 2005; Garcia et al. 2011a; Brewster et al. 2014) (we note that non-equilibrium models of simple repression have been shown to have the same functional form that we derive below (Phillips 2015)). We extend these models to consider allostery by accounting for the equilibrium state of the repressor through the MWC model.

     Thermodynamic models of gene expression begin by enumerating all possible states of the promoter and their corresponding statistical weights. As shown in Fig. 2 (A), the promoter can either be empty, occupied by RNAP, or occupied by either an active or inactive repressor. The probability of binding to the promoter will be affected by the protein copy number, which we denote as P for RNAP, R_A for active repressor, and R_I for inactive repressor. We note that repressors fluctuate between the active and inactive conformation in thermodynamic equilibrium, such that R_A and R_I will, on average, remain constant for a given inducer concentration (Monod, Wyman, and Changeux 1965). We assign the repressor a different DNA binding affinity in the active and inactive state. In addition to the specific binding sites at the promoter, we assume that there are N_NS non-specific binding sites elsewhere (i.e. on parts of the genome outside the simple repression architecture) where the RNAP or the repressor can bind. All specific binding energies are measured relative to the average non-specific binding energy. Thus, Δε_P represents the energy difference between the specific and non-specific binding for RNAP to the DNA. Likewise, Δε_RA and Δε_RI represent the difference in specific and non-specific binding energies for repressor in the active or inactive state, respectively.

     Thermodynamic models of transcription (Ackers and Johnson 1982; Buchler, Gerland, and Hwa 2003; Vilar and Leibler 2003; Bintu, Buchler, Garcia, Gerland, Hwa, Kondev, and Phillips 2005; Bintu, Buchler, Garcia, Gerland, Hwa, Kondev, Kuhlman, et al. 2005; Kuhlman et al. 2007; Daber, Sochor, and Lewis 2011; Garcia et al. 2011a; Brewster et al. 2014; Weinert et al. 2014) posit that gene expression is proportional to the probability that the RNAP is bound to the promoter p_bound, which is given by
$$ p_\text{bound} = \frac{\frac{P}{N_{NS}}e^{-\beta \Delta\varepsilon_{P}}}{1+\frac{R_A}{N_{NS}}e^{-\beta \Delta\varepsilon_{RA}}+\frac{R_I}{N_{NS}}e^{-\beta \Delta\varepsilon_{RI}}+\frac{P}{N_{NS}}e^{-\beta\Delta\varepsilon_{P}}}, \qquad(1)$$
with β = 1/k_BT, where k_B is the Boltzmann constant and T is the temperature of the system. As k_BT is the natural unit of energy at the molecular length scale, we treat the products βΔε_j as single parameters within our model. Measuring p_bound directly is fraught with experimental difficulties, as determining the exact proportionality between expression and p_bound is not straightforward. Instead, we measure the fold-change in gene expression due to the presence of the repressor. We define fold-change as the ratio of gene expression in the presence of repressor relative to expression in the absence of repressor (i.e. constitutive expression), namely,
$$ \text{fold-change} \equiv \frac{p_\text{bound}(R > 0)}{p_\text{bound}(R = 0)}. \qquad(2)$$

     We can simplify this expression using two well-justified approximations: (1) (P/N_NS)e^{− βΔε_P}≪ 1 implying that the RNAP binds weakly to the promoter (N_NS = 4.6 × 10⁶, P ≈ 10³ (Klumpp and Hwa 2008), Δε_P ≈  − 2  to  − 5 k_BT (Brewster, Jones, and Phillips 2012), so that (P/N_NS)e^{− βΔε_P} ≈ 0.01) and (2) (R_I/N_NS)e^{− βΔε_RI} ≪ 1 + (R_A/N_NS)e^{− βΔε_RA} which reflects our assumption that the inactive repressor binds weakly to the promoter of interest. Using these approximations, the fold-change reduces to the form
$$ \text{fold-change} \approx \left(1+\frac{R_A}{N_{NS}}e^{-\beta \Delta\varepsilon_{RA}}\right)^{-1} \equiv \left( 1+p_\text{act}(c) \frac{R}{N_{NS}}e^{-\beta\Delta\varepsilon_{RA}} \right)^{-1}, \qquad(3)$$
where in the last step we have introduced the fraction p_act(c) of repressors in the active state given a concentration c of inducer, such that R_A(c) = p_act(c)R. Since inducer binding shifts the repressors from the active to the inactive state, p_act(c) grows smaller as c increases.

     We use the MWC model to compute the probability p_act(c) that a repressor with n inducer binding sites will be active. The value of p_act(c) is given by the sum of the weights of the active repressor states divided by the sum of the weights of all possible repressor states (see Fig. 2 (B)), namely,
$$ p_\text{act}(c)=\frac{\left(1+\frac{c}{K_A}\right)^n}{\left(1+\frac{c}{K_A}\right)^n+e^{-\beta \Delta \varepsilon_{AI} }\left(1+\frac{c}{K_I}\right)^n}, \qquad(4)$$
where K_A and K_I represent the dissociation constant between the inducer and repressor in the active and inactive states, respectively, and Δε_AI = ε_I − ε_A is the free energy difference between a repressor in the inactive and active state (the quantity e^{− βΔε_AI} is sometimes denoted by L (Monod, Wyman, and Changeux 1965; Marzen, Garcia, and Phillips 2013) or K_RR* (Daber, Sochor, and Lewis 2011)). In this equation, c/K_A and c/K_I represent the change in free energy when an inducer binds to a repressor in the active or inactive state, respectively, while e^{− βΔε_AI} represents the change in free energy when the repressor changes from the active to inactive state in the absence of inducer. Thus, a repressor which favors the active state in the absence of inducer (Δε_AI > 0) will be driven towards the inactive state upon inducer binding when K_I < K_A. The specific case of a repressor dimer with n = 2 inducer binding sites is shown in Fig. 2 (B).

    Substituting p_act(c) from Eq. 4 into Eq. 3 yields the general formula for induction of a simple repression regulatory architecture (Phillips 2015), namely,
$$ \text{fold-change} = \left(1+\frac{\left(1+\frac{c}{K_A}\right)^n}{\left(1+\frac{c}{K_A}\right)^n+e^{-\beta \Delta \varepsilon_{AI} }\left(1+\frac{c}{K_I}\right)^n}\frac{R}{N_{NS}}e^{-\beta \Delta\varepsilon_{RA}} \right)^{-1}. \qquad(5)$$

While we have used the specific case of simple repression with induction to craft this model, the same mathematics describe the case of corepression in which binding of an allosteric effector stabilizes the active state of the repressor and decreases gene expression (see Fig. 1). Interestingly, we shift from induction (governed by K_I < K_A) to corepression (K_I > K_A) as the ligand transitions from preferentially binding to the inactive repressor state to stabilizing the active state. Furthermore, this general approach can be used to describe a variety of other motifs such as activation, multiple repressor binding sites, and combinations of activator and repressor binding sites (Bintu, Buchler, Garcia, Gerland, Hwa, Kondev, Kuhlman, et al. 2005; Brewster et al. 2014; Weinert et al. 2014).

    The formula presented in Eq. 5 enables us to make precise quantitative statements about induction profiles. Motivated by the broad range of predictions implied by Eq. 5, we designed a series of experiments using the lac system in E. coli to tune the control parameters for a simple repression genetic circuit. As discussed in Fig. 1 (C), previous studies have provided well-characterized values for many of the parameters in our experimental system, leaving only the values of the MWC parameters (K_A, K_I, and Δε_AI) to be determined. We note that while previous studies have obtained values for K_A, K_I, and L = e^{− βΔε_AI} (O’Gorman et al. 1980; Daber, Sochor, and Lewis 2011), they were either based upon in vitro biochemical experiments or in vivo conditions involving poorly characterized transcription factor copy numbers and gene copy numbers. These differences relative to our experimental conditions and fitting techniques led us to believe that it was important to perform our own analysis of these parameters. After inferring these three MWC parameters (see the supplemental Chapter 6 for details regarding the inference of Δε_AI, which was fitted separately from K_A and K_I), we were able to predict the input/output response of the system under a broad range of experimental conditions. For example, this framework can predict the response of the system at different repressor copy numbers R, repressor-operator affinities Δε_RA, inducer concentrations c, and gene copy numbers.

Results

Experimental Design

     We test our model by predicting the induction profiles for an array of strains that could be made using previously characterized repressor copy numbers and DNA binding energies. Our approach contrasts with previous studies that have parameterized induction curves of simple repression motifs, as these have relied on expression systems where proteins are expressed from plasmids, resulting in highly variable and unconstrained copy numbers (Murphy, Balázsi, and Collins 2007; Daber, Sharp, and Lewis 2009; Daber, Sochor, and Lewis 2011; Sochor 2014). Instead, our approach relies on a foundation of previous work as depicted in Fig. 1 (C). This includes work from our laboratory that used E. coli constructs based on components of the lac system to demonstrate how the Lac repressor (LacI) copy number R and operator binding energy Δε_RA affect gene expression in the absence of inducer (Garcia et al. 2011a). Rydenfelt et al. (2014) extended the theory used in that work to the case of multiple promoters competing for a given transcription factor, which was validated experimentally by Brewster et al. (2014), who modified this system to consider expression from multiple-copy plasmids as well as the presence of competing repressor binding sites.

    The present study extends this body of work by introducing three additional biophysical parameters – Δε_AI, K_A, and K_I – which capture the allosteric nature of the transcription factor and complement the results shown by Garcia et al. (2011a) and Brewster et al. (2014). Although the current work focuses on systems with a single site of repression, in the Materials & Methods, we utilize data from Brewster et al. (2014), in which multiple sites of repression are explored, to characterize the allosteric free energy difference Δε_AI between the repressor’s active and inactive states. This additional data set is critical because multiple degenerate sets of parameters can characterize an induction curve equally well, with the Δε_AI parameter compensated by the inducer dissociation constants K_A and K_I (see supplemental Chapter 6). After fixing Δε_AI as described in the Materials & Methods, we can use data from single-site simple repression systems to determine the values of K_A and K_I.

    We determine the values of K_A and K_I by fitting to a single induction profile using Bayesian inferential methods (Sivia and Skilling 2006). We then use Eq. 5 to predict gene expression for any concentration of inducer, repressor copy number, and DNA binding energy, and compare these predictions against experimental measurements. To obtain induction profiles for a set of strains with varying repressor copy numbers, we used modified lacI ribosomal binding sites from Garcia et al. (2011a) to generate strains with repressor copy number of R = 22 ± 4, 60 ± 20, 124 ± 30, 260 ± 40, 1220 ± 160, and 1740 ± 340 per cell on average, where the error denotes standard deviation of at least three replicates as measured by Garcia et al. (2011a). We note that R refers to the number of repressor dimers in the cell, which is twice the number of repressor tetramers reported by Garcia et al. (2011a); since both heads of the repressor are assumed to always be either specifically or non-specifically bound to the genome, the two repressor dimers in each LacI tetramer can be considered independently. Gene expression was measured using a Yellow Fluorescent Protein (YFP) gene, driven by a lacUV5 promoter. Each of the six repressor copy number variants were paired with the native O1, O2, or O3 lac operator (Oehler et al. 1994) placed at the YFP transcription start site, thereby generating eighteen unique strains. The repressor-operator binding energies (O1 Δε_RA =  − 15.3 ± 0.2 k_BT, O2 Δε_RA =  − 13.9 k_BT ± 0.2, and O3 Δε_RA =  − 9.7 ± 0.1 k_BT) were previously inferred by measuring the fold-change of the lac system at different repressor copy numbers, where the error arises from model fitting (Garcia et al. 2011a). Additionally, we were able to obtain the value Δε_AI = 4.5 k_BT by fitting to previous data as discussed in the Materials & Methods. We measure fold-change over a range of known IPTG concentrations c, using n = 2 inducer binding sites per LacI dimer and approximating the number of non-specific binding sites as the length in base-pairs of the E. coli genome, N_NS = 4.6 × 10⁶.

    Our experimental pipeline for determining fold-change using flow cytometry is shown in Fig. 3. Briefly, cells were grown to exponential phase, in which gene expression reaches steady state (???), under concentrations of the inducer IPTG ranging between 0 and 5000 μM. We measure YFP fluorescence using flow cytometry and automatically gate the data to include only single-cell measurements (see Materials & Methods). To validate the use of flow cytometry, we also measured the fold-change of a subset of strains using the established method of single-cell microscopy (see supplemental Chapter 6). We found that the fold-change measurements obtained from microscopy were indistinguishable from that of flow-cytometry and yielded values for the inducer binding constants K_A and K_I that were within error.

Determination of the in vivo MWC Parameters

     The three parameters that we tune experimentally are shown in Fig. 4 (A), leaving the three allosteric parameters (Δε_AI, K_A, and K_I) to be determined by fitting. We used previous LacI fold-change data (Brewster et al. 2014) to infer that Δε_AI = 4.5 k_BT (see Materials & Methods). Rather than fitting K_A and K_I to our entire data set of eighteen unique constructs, we performed Bayesian parameter estimation on data from a single strain with R = 260 and an O2 operator (Δε_RA =  − 13.9 k_BT (Garcia et al. 2011a)) shown in Fig. 4(D, white-faced points). Using Markov Chain Monte Carlo, we determine the most likely parameter values to be K_A = 139_− 22^+ 29 μM and K_I = 0.53_− 0.04^+ 0.04 μM, which are the modes of their respective distributions, where the superscripts and subscripts represent the upper and lower bounds of the 95^th percentile of the parameter value distributions (see Fig. 4 (B)). Unfortunately, we are not able to make a meaningful value-for-value comparison of our parameters to those of earlier studies (Daber, Sharp, and Lewis 2009; Daber, Sochor, and Lewis 2011) because of uncertainties in both gene copy number and transcription factor copy numbers in these studies, as illustrated in supplemental Chapter 6. We then predicted the fold-change for the remaining seventeen strains with no further fitting (see Fig. 4 (C – E)) together with the specific phenotypic properties described in Fig. 1(B) and discussed in detail below (see (Fig. 4 (F – J)). The shaded regions in Fig. 4 (C – E) denote the 95% credible regions. Factors determining the width of the credible regions are explored in the supplemental Chapter 6.

     We stress that the entire suite of predictions in Fig. 4 (C – J) is based upon the induction profile of a single strain. Our ability to make such a broad range of predictions stems from the fact that our parameters of interest – such as the repressor copy number and DNA binding energy – appear as distinct physical parameters within our model. While the single data set in Fig. 4 could also be fit using a Hill function, such an analysis would be unable to predict any of the other curves in the figure. Phenomenological expressions such as the Hill function can describe data, but lack predictive power and are thus unable to build our intuition, help us design de novo input-output functions, or guide future experiments (Kuhlman et al. 2007; Murphy, Balázsi, and Collins 2007).

Comparison of Experimental Measurements with Theoretical Predictions

    We tested the predictions shown in Fig. 4 by measuring fold-change induction profiles in strains with a broad range of repressor copy numbers and repressor binding energies as characterized in Garcia et al. (2011a). With a few notable exceptions, the results shown in Fig. 5 demonstrate agreement between theory and experiment. We note that there was an apparently systematic shift in the O3 Δε_RA =  − 9.7 k_BT strains Fig. 5 (C) and all of the R = 1220 and R = 1740 strains. This may be partially due to imprecise previous determinations of their Δε_RA and R values. By performing a global fit where we infer all parameters including the repressor copy number R and the binding energy Δε_RA, we found better agreement for these strains, although a discrepancy in the steepness of the response for all O3 strains remains (see supplemental Chapter 6). We considered a number of hypotheses to explain these discrepancies such as including other states (e.g. non-negligible binding of the inactive repressor), relaxing the weak promoter approximation, and accounting for variations in gene and repressor copy number throughout the cell cycle, but none explained the observed discrepancies. As an additional test of our model, we considered strains using the synthetic Oid operator which exhibits an especially strong binding energy of Δε_RA =  − 17 k_BT (Garcia et al. 2011a). The global fit agrees well with the Oid microscopy data, though it asserts a stronger Oid binding energy of Δε_RA =  − 17.7 k_BT (see supplemental Chapter 6).

    To ensure that the agreement between our predictions and data is not an accident of the strain we used to perform our fitting, we also inferred K_A and K_I from each of the other strains. As discussed in supplemental Chapter 6 and Fig. 4, the inferred values of K_A and K_I depend minimally upon which strain is chosen, indicating that these parameter values are highly robust. We also performed a global fit using the data from all eighteen strains in which we fitted for the inducer dissociation constants K_A and K_I, the repressor copy number R, and the repressor DNA binding energy Δε_RA (see supplemental Chapter 6). The resulting parameter values were nearly identical to those fitted from any single strain. For the remainder of this chapter, we continue using parameters inferred from the strain with R = 260 repressors and an O2 operator.

Predicting the Phenotypic Traits of the Induction Response

     The properties shown in Fig. 1 (i.e. the leakiness, saturation, dynamic range, [EC₅₀], and effective Hill coefficient) are of significant interest to synthetic biology. For example, synthetic biology is often focused on generating large responses (i.e. a large dynamic range) or finding a strong binding partner (i.e. a small [EC₅₀]) (Brophy and Voigt 2014; Shis et al. 2014). While these properties are all individually informative and when taken together, they capture the essential features of the induction response. We reiterate that a Hill function approach cannot predict these features a priori and furthermore requires fitting each curve individually. The MWC model, on the other hand, enables us to quantify how each trait depends upon a single set of physical parameters as shown by Fig. 4 (F – J).

    We define these five phenotypic traits using expressions derived from the model presented in Eq. 5. These results build upon extensive work by Martins and Swain (2011), who computed many such properties for ligand-receptor binding within the MWC model. We begin by analyzing the leakiness, which is the minimum fold-change observed in the absence of ligand, given by
$$ \text{leakiness} = \text{fold-change}(c=0) \\ = \left( 1+\frac{1}{1+e^{-\beta \Delta \varepsilon_{AI} }}\frac{R}{N_{NS}}e^{-\beta \Delta\varepsilon_{RA}} \right)^{-1}, \qquad(6)$$
and the saturation, which is the maximum fold change observed in the presence of saturating ligand,
$$ \begin{aligned} \text{saturation} &= \text{fold-change}(c \to \infty) \\ &= \left( 1+\frac{1}{1+e^{-\beta \Delta \varepsilon_{AI} } \left(\frac{K_A}{K_I}\right)^n }\frac{R}{N_{NS}}e^{-\beta \Delta\varepsilon_{RA} } \right)^{-1}. \end{aligned} \qquad(7)$$
Systems that minimize leakiness repress strongly in the absence of effector while systems that maximize saturation have high expression in the presence of effector. Together, these two properties determine the dynamic range of a system’s response, which is given by the difference
dynamic range = saturation − leakiness.   (8)

     These three properties are shown in Fig. 4 (F-H). We discuss these properties in greater detail in supplemental Chapter 6. Fig. 6 shows that the measurements of these three properties, derived from the fold-change data in the absence of IPTG and the presence of saturating IPTG, closely match the predictions for all three operators.

     Two additional properties of induction profiles are the [EC₅₀] and effective Hill coefficient, which determine the range of inducer concentration in which the system’s output goes from its minimum to maximum value. The [EC₅₀] denotes the inducer concentration required to generate a system response halfway between its minimum and maximum value,
$$ \text{fold-change}(c = [EC_{50}]) = \frac{\text{leakiness} + \text{saturation}}{2}. \qquad(9)$$

The effective Hill coefficient h, which quantifies the steepness of the curve at the [EC₅₀] , is given by
$$ h = \left( 2 \frac{d}{d \log c} \left[ \log \left( \frac{ \text{fold-change}(c) - \text{leakiness}}{\text{dynamic range}} \right) \right] \right)_{c = [EC_{50}]}. \qquad(10)$$

      Fig. 6 (D–E) shows how the [EC₅₀] and effective Hill coefficient depend on the repressor copy number. In supplemental Chapter 6, we discuss the analytic forms of these two properties as well as their dependence on the repressor-DNA binding energy. Fig. 6 (D-E) shows the estimated values of the [EC₅₀] and the effective Hill coefficient overlaid on the theoretical predictions. Both properties were obtained by fitting to each individual titration curve and computing the [EC₅₀] and effective Hill coefficient. We find that the predictions made with the single strain closely match those made for each of the strains with O1 and O2 operators, but the predictions for the O3 operator are markedly off. In the supplemental Chapter 6, we show that the large, asymmetric error bars for the O3 R = 22 strain arise from its nearly flat response, where the lack of dynamic range makes it impossible to determine the value of the inducer dissociation constants K_A and K_I, as can be seen in the uncertainty of both the [EC₅₀] and effective Hill coefficient. Discrepancies between theory and data for O3 are improved, but not fully resolved, by performing a global fit or fitting the MWC model individually to each curve (see supplemental Chapter 6). It remains an open question how to account for discrepancies in O3, in particular regarding the significant mismatch between the predicted and fitted effective Hill coefficients.

Data Collapse of Induction Profiles

     Our primary interest heretofore was to determine the system response at a specific inducer concentration, repressor copy number, and repressor-DNA binding energy. However, the cell does not necessarily “care about” the precise number of repressors in the system or the binding energy of an individual operator. The relevant quantity for cellular function is the fold-change enacted by the regulatory system. This raises the question: given a specific value of the fold-change, what combination of parameters will give rise to this desired response? In other words, what trade-offs between the parameters of the system will give rise to the same mean cellular output? These are key questions both for understanding how the system is governed and, as will become evident in the following chapters of this dissertation, can provide insight as to what parameters may be changing in response to a physiological or environmental perturbation. To address these questions, we follow the data collapse strategy used in a number of previous studies (Sourjik and Berg 2002; Keymer et al. 2006; Swem et al. 2008).

    The equilibrium states and statistical weights outlined in Fig. 2 (A) can be further coarse grained into two possible states – one state being where the promoter is occupied by the repressor and another being where the promoter is not occupied by the repressor (Fig. 7 (A)). As the transcriptionally active state and the states in which the repressor is bound are mutually exclusive, we can compute the probability of the repressor not being bound p_¬r to the promoter as
$$ p_{\neg r} = \frac{\neg r}{r + \neg r}. \qquad(11)$$
We can now take a similar approach as in Eq. 2 and define the fold-change as the probability of the repressor not being bound when repressor is expressed p_¬r(R > 0) relative to the probability when no repressor is expressed p_¬r(R = 0). As the later term is equal to 1, the fold-change in gene expression is directly equivalent to p_¬r expressed in Eq. 11. This form can be algebraically manipulated to the form
$$ \text{fold-change} = \frac{1}{1 + \frac{r}{\neg r}} = \frac{1}{1 + e^{-\beta F}} \qquad(12)$$
where F can be interpreted as the difference in free energy between the repressor bound and repressor not bound states,
F = k_BT[log¬r−logr].   (13)

     As Fig. 2 provides mathematical forms for r and ¬r, F can be directly computed as
$$ F = \frac{\Delta\varepsilon_{RA}}{k_BT} - \log \frac{\left(1+\frac{c}{K_A}\right)^n}{\left(1+\frac{c}{K_A}\right)^n+e^{-\beta \Delta\varepsilon_{AI} }\left(1+\frac{c}{K_I}\right)^n} - \log \frac{R}{N_{NS}}. \qquad(14)$$

     The first term in F denotes the repressor-operator binding energy, the second the contribution from the inducer concentration, and the last the effect of the repressor copy number. We note that elsewhere, this free energy has been dubbed the Bohr parameter since such families of curves are analogous to the shifts in hemoglobin binding curves at different pHs, known as the Bohr effect (Mirny 2010; Phillips 2015; Einav, Mazutis, and Phillips 2016).

     Instead of analyzing each induction curve individually, the free energy provides a natural means to simultaneously characterize the diversity in our eighteen induction profiles. Fig. 8 (A) demonstrates how the various induction curves from Fig. 4 (C-E) all collapse onto a single master curve, where points from every induction profile that yield the same fold-change are mapped onto the same free energy. Fig. 8 (B) reveals complete data collapse for the 216 data points in Fig. 5 (A – C), demonstrating the close match between the theoretical predictions and experimental measurements across all eighteen strains.

    There are many different combinations of parameter values that can result in the same free energy as defined in Eq. 14. For example, suppose a system originally has a fold-change of 0.2 at a specific inducer concentration, and then operator mutations increase the Δε_RA binding energy (Garcia et al. 2011a). While this serves to initially increase both the free energy and the fold-change, a subsequent increase in the repressor copy number could bring the cell back to the original fold-change level. Such trade-offs hint that there need not be a single set of parameters that evoke a specific cellular response, but rather that the cell explores a large but degenerate space of parameters with multiple, equally valid paths.

Discussion

     Since the early work by Monod, Wyman, and Changeux (Monod, Changeux, and Jacob 1963; Monod, Wyman, and Changeux 1965), an array of biological phenomena has been tied to the existence of macromolecules that switch between inactive and active states. Examples can be found in a wide variety of cellular processes, including ligand-gated ion channels (Auerbach 2012), enzymatic reactions (Velyvis et al. 2007; Einav, Mazutis, and Phillips 2016), chemotaxis (Keymer et al. 2006), quorum sensing (Swem et al. 2008), G-protein coupled receptors (Canals et al. 2012), physiologically important proteins (Milo et al. 2007; Levantino et al. 2012), and beyond. One of the most ubiquitous examples of allostery is in the context of gene expression, where an array of molecular players bind to transcription factors to influence their ability to regulate gene activity (Li et al. 2014). A number of studies have focused on developing a quantitative understanding of allosteric regulatory systems. The work of Martins and Swain (2011) and Marzen, Garcia, and Phillips (2013) analytically derives fundamental properties of the MWC model, including the leakiness and dynamic range described in this work, noting the inherent trade-offs in these properties when tuning the model’s parameters. Work in the Church and Voigt labs, among others, has expanded on the availability of allosteric circuits for synthetic biology (Lutz and Bujard 1997; Moon et al. 2012; Rogers et al. 2015; Rohlhill, Sandoval, and Papoutsakis 2017). Somewhat recently, Daber, Sharp, and Lewis (2009) theoretically explored the induction of simple repression within the MWC model and experimentally measured how mutations alter the induction profiles of transcription factors (Daber, Sochor, and Lewis 2011). Vilar and Saiz (2013) analyzed a variety of interactions in inducible lac-based systems including the effects of oligomerization and DNA folding on transcription factor induction. Other work has attempted to use the lac system to reconcile in vitro and in vivo measurements (Tungtur et al. 2011).

     Although this body of work has done much to improve our understanding of allosteric transcription factors, there have been few attempts to explicitly connect quantitative models to experiments. Here, we generate a predictive model of allosteric transcriptional regulation, and then test the model against a thorough set of experiments using well-characterized regulatory components. Specifically, we used the MWC model to build upon a well-established thermodynamic model of transcriptional regulation (Bintu, Buchler, Garcia, Gerland, Hwa, Kondev, Kuhlman, et al. 2005; Garcia et al. 2011a), allowing us to compose the model from a minimal set of biologically meaningful parameters. This model combines both theoretical and experimental insights; for example, rather than considering gene expression directly, we analyze the fold-change in expression, where the weak promoter approximation circumvents uncertainty in the RNAP copy number. The resulting model depended upon experimentally accessible parameters, namely, the repressor copy number, the repressor-DNA binding energy, and the concentration of inducer. We tested these predictions on a range of strains whose repressor copy number spanned two orders of magnitude and whose DNA binding affinity spanned 6 k_BT. We argue that one would not be able to generate such a wide array of predictions by using a Hill function, which abstracts away the biophysical meaning of the parameters into phenomenological parameters (Forsén and Linse 1995). Furthermore, our model reveals systematic relationships between behaviors that previously were only determined empirically.

     One such property is the dynamic range, which is of considerable interest when designing or characterizing a genetic circuit, is revealed to have an interesting property: although changing the value of Δε_RA causes the dynamic range curves to shift to the right or left, each curve has the same shape and, in particular, the same maximum value. This means that strains with strong or weak binding energies can attain the same dynamic range when the value of R is tuned to compensate for the binding energy. This feature is not immediately apparent from the IPTG induction curves, which show very low dynamic ranges for several of the O1 and O3 strains. Without the benefit of models that can predict such phenotypic traits, efforts to engineer genetic circuits with allosteric transcription factors must rely on trial and error to achieve specific responses (Rogers et al. 2015; Rohlhill, Sandoval, and Papoutsakis 2017). Other calculable properties, such as leakiness, saturation, [EC₅₀], and the effective Hill coefficient, agree well with experimental measurement. One exception is the titration profile of the weakest operator, O3. While performing a global fit for all model parameters marginally improves the prediction of all properties for O3 (see supplemental Chapter 6), a noticeable difference remains when inferring the effective Hill coefficient or the [EC₅₀]. We further tried including additional states (such as allowing the inactive repressor to bind to the operator), relaxing the weak promoter approximation, accounting for changes in gene and repressor copy number throughout the cell cycle (Jones, Brewster, and Phillips 2014), and refitting the original binding energies from Garcia et al. (2011a), but such generalizations were unable to account for the O3 data. It remains an open question as to how the discrepancy between the theory and measurements for O3 can be reconciled.

     Despite the diversity observed in the induction profiles of each of our strains, our data are unified by their reliance on fundamental biophysical parameters. In particular, we have shown that our model for fold-change can be rewritten in terms of the free energy, which encompasses all of the physical parameters of the system. This has proven to be an illuminating technique in a number of studies of allosteric proteins (Sourjik and Berg 2002; Keymer et al. 2006; Swem et al. 2008). Although it is experimentally straightforward to observe system responses to changes in effector concentration c, framing the input-output function in terms of c can give the misleading impression that changes in system parameters lead to fundamentally altered system responses. Alternatively, if one can find the “natural variable” that enables the output to collapse onto a single curve, it becomes clear that the system’s output is not governed by individual system parameters, but rather the contributions of multiple parameters that define the natural variable. When our fold-change data are plotted against the respective free energies for each construct, they collapse cleanly onto a single curve (see Fig. 8). This enables us to analyze how parameters can compensate each other. For example, rather than viewing strong repression as a consequence of low IPTG concentration c or high repressor copy number R, we can now observe that strong repression is achieved when the free energy F(c)≤  − 5k_BT, a condition which can be reached in a number of ways.

     While our experiments validated the theoretical predictions in the case of simple repression, we expect the framework presented here to apply much more generally to different biological instances of allosteric regulation. For example, we can use this model to study more complex systems such as when transcription factors interact with multiple operators (Bintu, Buchler, Garcia, Gerland, Hwa, Kondev, Kuhlman, et al. 2005). We can further explore different regulatory configurations such as corepression, activation, and coactivation, each of which are found in E. coli. This work can also serve as a springboard to characterize not just the mean, but the full gene expression distribution, and thus quantify the impact of noise on the system (Eldar and Elowitz 2010). Another extension of this approach would be to theoretically predict and experimentally verify whether the repressor-inducer dissociation constants K_A and K_I or the energy difference Δε_AI between the allosteric states can be tuned by making single amino acid substitutions in the transcription factor (Daber, Sharp, and Lewis 2009; Phillips 2015). Finally, we expect that the kind of rigorous quantitative description of the allosteric phenomenon provided here will make it possible to construct biophysical models of fitness for allosteric proteins similar to those already invoked to explore the fitness effects of transcription factor binding site strengths and protein stability (Gerland, Moroz, and Hwa 2002; Berg, Willmann, and Lässig 2004; Zeldovich and Shakhnovich 2008). In total, these results show that a thermodynamic formulation of the MWC model supersedes phenomenological fitting functions for understanding transcriptional regulation by allosteric proteins.

Materials & Methods

Bacterial Strains and DNA Constructs

All strains used in these experiments were derived from E. coli K12 MG1655 with the lac operon removed, adapted from those created and described in Garcia et al. (2011a). Briefly, the operator variants and YFP reporter gene were cloned into a pZS25 background which contains a lacUV5 promoter that drives expression as is shown schematically in Fig. 2. These constructs carried a kanamycin resistance gene and were integrated into the galK locus of the chromosome using λ Red recombineering (Sharan et al. 2009). The lacI gene was constitutively expressed via a P_LtetO-1 promoter (Lutz and Bujard 1997), with ribosomal binding site mutations made to vary the LacI copy number as described in Salis, Mirsky, and Voigt (2009) using site-directed mutagenesis (Quickchange II; Stratagene), with further details in Garcia et al. (2011a). These lacI constructs carried a chloramphenicol resistance gene and were integrated into the ybcN locus of the chromosome. Final strain construction was achieved by performing repeated P1 transduction (Thomason, Costantino, and Court 2007) of the different operator and lacI constructs to generate each combination used in this work. Integration was confirmed by PCR amplification of the replaced chromosomal region and by sequencing. Primers and final strain genotypes are listed in supplemental Chapter 6.

     It is important to note that the rest of the lac operon (lacZYA) was never expressed. The LacY protein is a transmembrane protein which actively transports lactose as well as IPTG into the cell. As LacY was never produced in our strains, we assume that the extracellular and intracellular IPTG concentration was approximately equal due to diffusion across the membrane into the cell as is suggested by previous work (Fernández-Castané et al. 2012).

     To make this theory applicable to transcription factors with any number of DNA binding domains, we used a different definition for repressor copy number than has been used previously. We define the LacI copy number as the average number of repressor dimers per cell, whereas in Garcia et al. (2011a), the copy number is defined as the average number of repressor tetramers in each cell. To motivate this decision, we consider the fact that the LacI repressor molecule exists as a tetramer in E. coli (Lewis et al. 1996) in which a single DNA binding domain is formed from dimerization of LacI proteins, so that wild-type LacI might be described as dimer of dimers. Since each dimer is allosterically independent (i.e. either dimer can be allosterically active or inactive, independent of the configuration of the other dimer) (Daber, Sharp, and Lewis 2009), a single LacI tetramer can be treated as two functional repressors. Therefore, we have simply multiplied the number of repressors reported in Garcia et al. (2011a) by a factor of two.

     A subset of strains in these experiments were measured using fluorescence microscopy for validation of the flow cytometry data and results. To aid in the high-fidelity segmentation of individual cells, the strains were modified to constitutively express an mCherry fluorophore. This reporter was cloned into a pZS4*1 backbone (Lutz and Bujard 1997) in which mCherry is driven by the lacUV5 promoter. All microscopy and flow cytometry experiments were performed using these strains.

Growth Conditions for Flow Cytometry Measurements

All measurements were performed with E. coli cells grown to mid-exponential phase in standard M9 minimal media (M9 5X Salts, Sigma-Aldrich M6030; 2 mM magnesium sulfate, Mallinckrodt Chemicals 6066-04; 100 μM calcium chloride, Fisher Chemicals C79-500) supplemented with 0.5% (w/v) glucose. Briefly, 500 μL cultures of E. coli were inoculated into Lysogeny Broth (LB Miller Powder, BD Medical) from a 50% glycerol frozen stock (-80^∘C) and were grown overnight in a 2 mL 96-deep-well plate sealed with a breathable nylon cover (Lab Pak - Nitex Nylon, Sefar America Inc. Cat. No. 241205) with rapid agitation for proper aeration. After approximately 12 to 15 hours, the cultures had reached saturation and were diluted 1000-fold into a second 2 mL 96-deep-well plate where each well contained 500 μL of M9 minimal media supplemented with 0.5% w/v glucose (anhydrous D-Glucose, Macron Chemicals) and the appropriate concentration of IPTG (Isopropyl β-D-1 thiogalactopyranoside Dioxane Free, Research Products International). These were sealed with a breathable cover and were allowed to grow for approximately eight hours. Cells were then diluted ten-fold into a round-bottom 96-well plate (Corning Cat. No. 3365) containing 90 μL of M9 minimal media supplemented with 0.5% w/v glucose along with the corresponding IPTG concentrations. For each IPTG concentration, a stock of 100-fold concentrated IPTG in double distilled water was prepared and partitioned into 100 μL aliquots. The same parent stock was used for all experiments described in this work.

Flow Cytometry

All fold-change measurements were collected on a Miltenyi Biotec MACSquant Analyzer 10 Flow Cytometer graciously provided by the Pamela Björkman lab at Caltech. Detailed information regarding the voltage settings of the photo-multiplier detectors can be found in the supplemental Chapter 6.

     Prior to each day’s experiments, the analyzer was calibrated using MACSQuant Calibration Beads (Cat. No. 130-093-607) such that day-to-day experiments would be comparable. All YFP fluorescence measurements were collected via 488 nm laser excitation coupled with a 525/50 nm emission filter. Unless otherwise specified, all measurements were taken over the course of two to three hours using automated sampling from a 96-well plate kept at approximately $4^\circ \, \hbox{-} \, 10^\circ$C on a MACS Chill 96 Rack (Cat. No. 130-094-459). Cells were diluted to a final concentration of approximately 4 × 10⁴ cells per μL which corresponded to a flow rate of 2,000-6,000 measurements per second, and acquisition for each well was halted after 100,000 events were detected. Once completed, the data were extracted and immediately processed using the following methods.

Unsupervised Gating of Flow Cytometry Data

Flow cytometry data will frequently include a number of spurious events or other undesirable data points such as cell doublets and debris. The process of restricting the collected data set to those data determined to be “real” is commonly referred to as gating. These gates are typically drawn manually and restrict the data set to those points which display a high degree of linear correlation between their forward-scatter (FSC) and side-scatter (SSC). The development of unbiased and unsupervised methods of drawing these gates is an active area of research (Lo, Brinkman, and Gottardo 2008; Aghaeepour et al. 2013). For our purposes, we assume that the fluorescence level of the population should be log-normally distributed about some mean value. With this assumption in place, we developed a method that allows us to restrict the data used to compute the mean fluorescence intensity of the population to the smallest two-dimensional region of the log (FSC) vs. log (SSC) space in which 40% of the data is found. This was performed by fitting a bivariate Gaussian distribution and restricting the data used for calculation to those that reside within the 40th percentile. This procedure is described in more detail in the supplemental Chapter 6.

Experimental Determination of Fold-Change

     For each strain and IPTG concentration, the fold-change in gene expression was calculated by taking the ratio of the population mean YFP expression in the presence of LacI repressor to that of the population mean in the absence of LacI repressor. However, the measured fluorescence intensity of each cell also includes the autofluorescence contributed by the weak excitation of the myriad protein and small molecules within the cell. To correct for this background, we computed the fold change as
$$ \text{fold-change} = \frac{\langle I_{R > 0} \rangle - \langle I_\text{auto}\rangle}{\langle I_{R = 0} \rangle - \langle I_\text{auto}\rangle}, \qquad(15)$$
where ⟨I_R > 0⟩ is the average cell YFP intensity in the presence of repressor, ⟨I_R = 0⟩ is the average cell YFP intensity in the absence of repressor, and ⟨I_auto⟩ is the average cell autofluorescence intensity, as measured from cells that lack the lac-YFP construct.

Bayesian Parameter Estimation

In this work, we determine the most likely parameter values for the inducer dissociation constants K_A and K_I of the active and inactive state, respectively, using Bayesian methods. We compute the probability distribution of the value of each parameter given the data D, which by Bayes’ theorem is given by
$$P(K_A, K_I \,\vert\, D) = \frac{P(D \,\vert\, K_A, K_I)P(K_A, K_I)}{P(D)}, \qquad(16)$$
where D is all the data composed of independent variables (repressor copy number R, repressor-DNA binding energy Δε_RA, and inducer concentration c) and one dependent variable (experimental fold-change). P(D ∣ K_A, K_I) is the likelihood of having observed the data given the parameter values for the dissociation constants, P(K_A, K_I) contains all the prior information on these parameters, and P(D) serves as a normalization constant, which we can ignore in our parameter estimation. Eq. 5 assumes a deterministic relationship between the parameters and the data, so in order to construct a probabilistic relationship as required by Eq. 16, we assume that the experimental fold-change for the i^th datum given the parameters is of the form
$$ \text{fold-change}_\text{exp}^{(i)} = \left( 1 + \frac{\left(1 + \frac{c^{(i)}}{K_A}\right)^2}{\left( 1 + \frac{c^{(i)}}{K_A}\right)^2 + e^{-\beta \Delta \varepsilon_{AI}} \left(1 + \frac{c^{(i)}}{K_I} \right)^2} \frac{R^{(i)}}{N_{NS}} e^{-\beta \Delta \varepsilon_{RA}^{(i)}}\right)^{-1} + \epsilon^{(i)}, \qquad(17)$$
where ϵ⁽ⁱ⁾ represents the departure from the deterministic theoretical prediction for the i^th data point. If we assume that these ϵ⁽ⁱ⁾ errors are normally distributed with mean zero and standard deviation σ, the likelihood of the data given the parameters is of the form
$$ \begin{aligned} P(D \vert K_A, K_I, \sigma) &= \frac{1}{(2\pi\sigma^2)^{\frac{n}{2}}}\times\\ &\prod\limits_{i=1}^n \exp \left[-\frac{(\text{fold-change}^{(i)}_\text{exp} - \text{fc}^\text{(theo)}(K_A, K_I, R^{(i)}, \Delta\varepsilon_{RA}^{(i)}, c^{(i)}))^2}{2\sigma^2}\right], \end{aligned} \qquad(18)$$
where fold-change_exp⁽ⁱ⁾ is the experimental fold-change and fc^(theo)( ⋯) is the theoretical prediction. The product $\prod\limits_{i=1}^n$ captures the assumption that the n data points are independent. Note that the likelihood and prior terms now include the extra unknown parameter σ. In applying Eq. 18, a choice of K_A and K_I that provides better agreement between theoretical fold-change predictions and experimental measurements will result in a more probable likelihood.

     Both mathematically and numerically, it is convenient to define $\tilde{k}_A = -\log \frac{K_A}{1\,\mu\text{M}}$ and $\tilde{k}_I = -\log \frac{K_I}{1\,\mu\text{M}}$ and fit for these parameters on a log scale. Dissociation constants are scale invariant, so that a change from 10 μM to 1 μM leads to an equivalent increase in affinity as a change from 1 μM to 0.1 μM. With these definitions, we assume for the prior P(k̃_A, k̃_I, σ) that all three parameters are independent. In addition, we assume a uniform distribution for k̃_A and k̃_I and a Jeffreys prior for the scale parameter σ. This yields the complete prior
$$P(\tilde{k}_A, \tilde{k}_I, \sigma) \equiv \frac{1}{(\tilde{k}_A^{\max} - \tilde{k}_A^{\min})} \frac{1}{(\tilde{k}_I^{\max} - \tilde{k}_I^{\min})}\frac{1}{\sigma}. \qquad(19)$$

These priors are maximally uninformative meaning that they imply no prior knowledge of the parameter values. We defined the k̃_A and k̃_A ranges uniform on the range of  − 7 to 7, although we note that this particular choice does not affect the outcome provided the chosen range is sufficiently wide.

     Putting all these terms together we can now sample from P(k̃_A, k̃_I, σ ∣ D) using Markov chain Monte Carlo to compute the most likely parameter as well as the error bars (given by the 95% credible region) for K_A and K_I.

Data Curation

     All of the data used in this work as well as all relevant code can be found at the website associated with the publication mentioned at the beginning of this chapter. Data were collected, stored, and preserved using the Git version control software in combination with off-site storage and hosting website GitHub. Code used to generate all figures and complete all processing step as and analyses are available on the GitHub repository. Many analysis files are stored as instructive Jupyter Notebooks. The scientific community is invited to fork our repositories and open constructive issues on the GitHub repository.

References