Unknown Knowns, Known Unknowns, and Unforeseen Consequences: Using Free Energy Shifts To Predict Mutant Phenotypes

Published as …

---

A version of this chapter originally appeared as Chure, G, Razo-Mejia, M., Belliveau, N.M., Kaczmarek, Zofii A., Einav, T., Barnes, Stephanie L., Lewis, M., and Phillips, R. (2019). Predictive shifts in free energy couple mutations to their phenotypic consequences. Proceedings of the National Academies of Sciences 116(37) DOI: https://doi.org/10.1073/pnas.1907869116. G.C., M.R.M, N.M.B., Z.A.K., and S.L.B designed the experiments and collected and analyzed data. G.C. developed the theoretical treatment of free energy shifts. G.C., M.R.M, N.M.B., Z.A.K., T.E., S.L.B., and R.P. designed the research project. G.C. and R.P. wrote the paper. M.L. provided guidance and advice.

Abstract

Mutation is a critical mechanism by which evolution explores the functional landscape of proteins. Despite our ability to experimentally inflict mutations at will, it remains difficult to link sequence-level perturbations to systems-level responses. Here, we present a framework centered on measuring changes in the free energy of the system to link individual mutations in an allosteric transcriptional repressor to the parameters which govern its response. We find that the energetic effects of the mutations can be categorized into several classes which have characteristic curves as a function of the inducer concentration. We experimentally test these diagnostic predictions using the well-characterized LacI repressor of Escherichia coli, probing several mutations in the DNA binding and inducer binding domains. We find that the change in gene expression due to a point mutation can be captured by modifying only the model parameters that describe the respective domain of the wild-type protein. These parameters appear to be insulated, with mutations in the DNA binding domain altering only the DNA affinity and those in the inducer binding domain altering only the allosteric parameters. Changing these subsets of parameters tunes the free energy of the system in a way that is concordant with theoretical expectations. Finally, we show that the induction profiles and resulting free energies associated with pairwise double mutants can be predicted with quantitative accuracy given knowledge of the single mutants, providing an avenue for identifying and quantifying epistatic interactions.

Introduction

Thermodynamic treatments of transcriptional regulation have been fruitful in their ability to generate quantitative predictions of gene expression as a function of a minimal set of physically meaningful parameters (Ackers and Johnson 1982; Buchler, Gerland, and Hwa 2003; Vilar and Leibler 2003; Garcia and Phillips 2011; Daber, Sharp, and Lewis 2009; Brewster et al. 2014; Weinert et al. 2014; Rydenfelt et al. 2014; Razo-Mejia et al. 2014, 2018; Bintu, Buchler, Garcia, Gerland, Hwa, Kondev, and Phillips 2005; Bintu, Buchler, Garcia, Gerland, Hwa, Kondev, Kuhlman, et al. 2005; Kuhlman et al. 2007). These models quantitatively describe numerous properties of input-output functions, such as the leakiness, saturation, dynamic range, steepness of response, and the [EC₅₀] – the concentration of inducer at which the response is half maximal. The mathematical forms of these phenotypic properties are couched in terms of a minimal set of experimentally accessible variables, such as the inducer concentration, transcription factor copy number, and the DNA sequence of the binding site (see Chapter 2 and Razo-Mejia et al. (2018)). While the amino acid sequence of the transcription factor is another controllable variable, it is seldom implemented in quantitative terms considering mutations with subtle changes in chemistry frequently yield unpredictable physiological consequences. In this work, we examine how a series of mutations in either the DNA binding or inducer binding domains of a transcriptional repressor influence the values of the biophysical parameters which govern its regulatory behavior.

     We build upon the results presented in Chapter 2 of this thesis and present a theoretical framework for understanding how mutations in the amino acid sequence of the repressor affect different parameters and alter the free energy of the system. We find that the parameters capturing the allosteric nature of the repressor, the repressor copy number, and the DNA binding specificity contribute independently to the free energy of the system with different degrees of sensitivity. Furthermore, changes restricted to one of these three groups of parameters result in characteristic changes in the free energy relative to the wild-type repressor, providing falsifiable predictions of how different classes of mutations should behave.

     Next, we test these descriptions experimentally using the well-characterized transcriptional repressor of the lac operon LacI in E. coli regulating expression of a fluorescent reporter. We introduce a series of point mutations in either the inducer binding or DNA binding domain. We then measure the full induction profile of each mutant, determine the minimal set of parameters that are affected by the mutation, and predict how each mutation tunes the free energy at different inducer concentrations, repressor copy numbers, and DNA binding strengths. We find in general that mutations in the DNA binding domain only influence DNA binding strength, and that mutations within the inducer binding domain affect only the parameters which dictate the allosteric response. The degree to which these parameters are insulated is notable, as the very nature of allostery suggests that all parameters are intimately connected, thus enabling binding events at one domain to be “sensed” by another.

     With knowledge of how a collection of DNA binding and inducer binding single mutants behave, we predict the induction profiles and the free energy changes of pairwise double mutants with quantitative accuracy. We find that the energetic effects of each individual mutation are additive, indicating that epistatic interactions are absent between the mutations examined here. Our model provides a means for identifying and quantifying the extent of epistatic interactions in a more complex set of mutations, and can shed light on how the protein sequence and general regulatory architecture coevolve.

Theoretical Model

This work considers the inducible simple repression regulatory motif depicted in Fig. 1 (A) from a thermodynamic perspective which has been thoroughly dissected and tested experimentally (Garcia and Phillips 2011; Brewster et al. 2014; Razo-Mejia et al. 2018) and is described in depth in Chapter 2. The result of this extensive theory-experiment dialogue is a succinct input-output function schematized in Fig. 1 (B) that computes the fold-change in gene expression relative to an unregulated promoter. This function is of the form
$$ \text{fold-change} = \left(1 + {R_A \over N_{NS}}e^{-\beta\Delta\varepsilon_{RA}}\right)^{-1}, \qquad(1)$$
where R_A is the number of active repressors per cell, N_NS is the number of non-specific binding sites for the repressor, Δε_RA is the binding energy of the repressor to its specific binding site relative to the non-specific background, and β is defined as ${1 \over k_B T}$ where k_B is the Boltzmann constant and T is the temperature. While this theory requires knowledge of the number of active repressors, we often only know the total number R which is the sum total of active and inactive repressors. We can define a prefactor p_act(c) which captures the allosteric nature of the repressor and encodes the probability that a repressor is in the active (repressive) state rather than the inactive state for a given inducer concentration c, namely,
$$ p_\text{act}(c) = {\left(1 + {c \over K_A}\right)^n \over \left(1 + {c \over K_A}\right)^n + e^{-\beta\Delta\varepsilon_{AI}}\left(1 + {c \over K_I}\right)^n}. \qquad(2)$$

Here, K_A and K_I are the dissociation constants of the inducer to the active and inactive repressor, Δε_AI is the energetic difference between the repressor active and inactive states, and n is the number of allosteric binding sites per repressor molecule (n = 2 for LacI). With this in hand, we can define R_A in Eq. 1 as R_A = p_act(c)R.

     A key feature of Eq. 1 and Eq. 2 is that the diverse phenomenology of the gene expression induction profile can be collapsed onto a single master curve by rewriting the input-output function in terms of the free energy F also called the Bohr parameter (Phillips 2015),
$$ \text{fold-change} = \frac{1}{1 + e^{-\beta F}}, \qquad(3)$$
where
$$ F = -k_BT \log p_\text{act}(c) - k_BT\log\left({R \over N_{NS}}\right) + \Delta\varepsilon_{RA}. \qquad(4)$$
Hence, if different combinations of parameters yield the same free energy, they will give rise to the same fold-change in gene expression, enabling us to collapse multiple regulatory scenarios onto a single curve. This can be seen in Fig. 1 (C) where eighteen unique inducer titration profiles of a LacI simple repression architecture collected and analyzed in Razo-Mejia et al. (2018) collapse onto a single master curve. The tight distribution about this curve reveals that the fold-change across a variety of genetically distinct individuals can be adequately described by a small number of parameters. Beyond predicting the induction profiles of different strains, the method of data collapse inspired by Eq. 3 and Eq. 4 can be used as a tool to identify mechanistic changes in the regulatory architecture (Swem et al. 2008). Similar data collapse approaches have been used previously in such a manner and have proved vital for distinguishing between changes in parameter values and changes in the fundamental behavior of the system (Swem et al. 2008; Keymer et al. 2006).

     Assuming that a given mutation does not result in a non-functional protein, it is reasonable to say that any or all of the parameters in Eq. 1 can be affected by the mutation, changing the observed induction profile and therefore the free energy. To examine how the free energy of a mutant F^(mut) differs from that of the wild-type F^(wt), we define ΔF = F^(mut) − F^(wt), which has the form
$$ \begin{aligned} \Delta F = -k_BT\log\left({p_\text{act}^\mathrm{(mut)}(c) \over p_\text{act}^\mathrm{(wt)}(c)}\right) &- k_BT \log\left({R^\mathrm{(mut)}\over R^\mathrm{(wt)}}\right)\\ &+ (\Delta\varepsilon_{RA}^\mathrm{(mut)}- \Delta\varepsilon_{RA}^\mathrm{(wt)}). \end{aligned} \qquad(5)$$
ΔF describes how a mutation translates a point across the master curve shown in Fig. 1 (C). As we will show in the coming paragraphs (illustrated in Fig. 2), this formulation coarse grains the myriad parameters shown in Eq. 1 and Eq. 2 into three distinct quantities, each with different sensitivities to parametric changes. By examining how a mutation changes the ΔF as a function of the inducer concentration, one can draw conclusions as to which parameters have been modified based solely on the shape of the curve. To help the reader understand how various perturbations to the parameters tune the free energy, we have hosted an interactive figure on the dedicated website for the publication which makes exploration of parameter space a simpler task.

     The first term in Eq. 5 is the log ratio of the probability of a mutant repressor being active relative to the wild type at a given inducer concentration c. This quantity defines how changes to any of the allosteric parameters – such as inducer binding constants K_A and K_I or active/inactive state energetic difference Δε_AI – alter the free energy F, which can be interpreted as the free energy difference between the repressor bound and unbound states of the promoter. Fig. 2 (A) illustrates how perturbations to the inducer binding constants K_A and K_I alone alter the induction profiles and free energy as a function of the inducer concentration. In the limit where c = 0, the values of K_A and K_I do not factor into the calculation of p_act(c) given by Eq. 2 meaning that Δε_AI is the lone parameter setting the residual activity of the repressor. Thus, if only K_A and K_I are altered by a mutation, then ΔF should be 0 k_BT when c = 0, illustrated by the overlapping red, purple, and grey curves in the right-hand plot of Fig. 2 (A). However, if Δε_AI is influenced by the mutation (either alone or in conjunction with K_A and K_I), the leakiness will change, resulting in a non-zero ΔF when c = 0. This is illustrated in Fig. 2 (B) where Δε_AI is the only parameter affected by the mutation.

     It is important to note that for a mutation which perturbs only the inducer binding constants, the dependence of ΔF on the inducer concentration can be non-monotonic. While the precise values of K_A and K_I control the sensitivity of the repressor to inducer concentration, it is the ratio K_A/K_I that defines whether this non-monotonic behavior is observed. This can be seen more clearly when we consider the limit of saturating inducer concentration,
$$ \lim\limits_{c \rightarrow \infty} \log\left({p_\text{act}^\mathrm{(mut)}\over p_\text{act}^\mathrm{(wt)}}\right) \approx \log\left[{1 + e^{-\beta\Delta\varepsilon_{AI}^\mathrm{(wt)}} \left({K_A^\mathrm{(wt)}\over K_I^\mathrm{(wt)}}\right)^n \over 1 + e^{-\beta\Delta\varepsilon_{AI}^\mathrm{(wt)}} \left({K_A^\mathrm{(mut)}\over K_I^\mathrm{(mut)}}\right)^n}\right], \qquad(6)$$
which illustrates that ΔF returns to zero at saturating inducer concentration when the ratio K_A/K_I is the same for both the mutant and wild-type repressors, so long as Δε_AI is unperturbed. Non-monotonicity can only be achieved by changing K_A and K_I and therefore serves as a diagnostic for classifying mutational effects reliant solely on measuring the change in free energy. A rigorous proof of this non-monotonic behavior given changing K_A and K_I can be found in supplemental Chapter 7.

     The second term in Eq. 5 captures how changes in the repressor copy number contributes to changes in free energy. It is important to note that this contribution to the free energy change depends on the total number of repressors in the cell, not just those in the active state. This emphasizes that changes in the expression of the repressor are energetically divorced from changes to the allosteric nature of the repressor. As a consequence, the change in free energy is constant for all inducer concentrations, as is schematized in Fig. 2 (C). Because the magnitude of the free energy shift logarithmically proportional to the relative change in repressor copy number, a mutation which increases expression from 1 to 10 repressors per cell is more impactful from an energetic standpoint (k_BTlog (10) ≈ 2.3 k_BT) than an increase from 90 to 100 (k_BTlog (100/90) ≈ 0.1 k_BT). Appreciable changes in the free energy only arise when variations in the repressor copy number are larger than or comparable to an order of magnitude. Changes of this magnitude are certainly possible from a single point mutation, as it has been shown that even synonymous substitutions can drastically change translation efficiency (Frumkin et al. 2018).

     The third and final term in Eq. 5 is the difference in the DNA binding energy between the mutant and wild-type repressors. All else being equal, if the mutated state binds more tightly to the DNA than the wild type (Δε_RA^(wt) > Δε_RA^(mut)), the net change in the free energy is negative, indicating that the repressor bound states become more energetically favorable due to the mutation. Much like in the case of changing repressor copy number, this quantity is independent of inducer concentration and is therefore also constant (Fig. 2 (D)). However, the magnitude of the change in free energy is linear with DNA binding affinity while it is logarithmic with respect to changes in the repressor copy number. Thus, to change the free energy by 1 k_BT, the repressor copy number must change by a factor of  ≈ 2.3 whereas the DNA binding energy must change by 1 k_BT.

     The unique behavior of each quantity in Eq. 5 and its sensitivity with respect to the parameters makes ΔF useful as a diagnostic tool to classify mutations. Given a set of fold-change measurements, a simple rearrangement of Eq. 3 permits the direct calculation of the free energy, assuming that the underlying physics of the regulatory architecture has not changed. Thus, it becomes possible to experimentally test the general assertions made in Fig. 2.

Results

DNA Binding Domain Mutants

With this arsenal of analytic diagnostics, we can begin to explore the mutational space of the repressor and map these mutations to the biophysical parameters they control. As one of the most thoroughly studied transcription factors, LacI has been subjected to numerous crystallographic and mutational studies (Lewis et al. 1996; Daber, Sharp, and Lewis 2009; Daber, Sochor, and Lewis 2011). One such work generated a set of point mutations in the LacI repressor and examined the diversity of the phenotypic response to different allosteric effectors (Daber, Sochor, and Lewis 2011). However, several experimental variables were unknown, precluding precise calculation of ΔF as presented in the previous section. In Daber, Sochor, and Lewis (2011), the repressor variants and the fluorescence reporter were expressed from separate plasmids. As the copy numbers of these plasmids fluctuate in the population, both the population average repressor copy number and the number of regulated promoters were unknown. Both of these quantities have been shown previously to significantly alter the measured gene expression, and calculation of ΔF is dependent on knowledge of their values. While the approach presented in Daber, Sochor, and Lewis (2011) considers the Lac repressor as an MWC molecule, the copy numbers of the repressor and the reporter gene were swept into an effective parameter R/K_DNA, hindering our ability to distinguish between changes in repressor copy number or in DNA binding energy. To test our hypothesis of free energy differences resulting from various parameter perturbations, we used the data set in Daber, Sochor, and Lewis (2011) as a guide and chose a subset of the mutations to quantitatively dissect. To control copy number variation, the mutant repressors and the reporter gene were integrated into the E. coli chromosome where the copy numbers are known and tightly controlled (Razo-Mejia et al. 2018; Garcia and Phillips 2011). Furthermore, the mutations were paired with ribosomal binding sites where the level of translation of the wild-type repressor had been directly measured previously (Garcia and Phillips 2011).

     We made three amino acid substitutions (Y17I, Q18A, and Q18M) that are critical for the DNA-repressor interaction. These mutations were introduced into the lacI sequence used in Garcia and Phillips (2011) with four different ribosomal binding site sequences that were shown (via quantitative Western blotting) to tune the wild-type repressor copy number across three orders of magnitude. These mutant constructs were integrated into the E. coli chromosome harboring a Yellow Fluorescent Protein (YFP) reporter. The YFP promoter included the native O2 LacI operator sequence which the wild-type LacI repressor binds with high specificity (Δε_RA =  − 13.9 k_BT). The fold-change in gene expression for each mutant across twelve concentrations of IPTG was measured via flow cytometry. As we mutated only a single amino acid with the minimum number of base pair changes to the codons from the wild-type sequence, we find it unlikely that the repressor copy number was drastically altered from those reported in Garcia and Phillips (2011) for the wild-type sequence paired with the same ribosomal binding site sequence. In characterizing the effects of these DNA binding mutations, we take the repressor copy number to be unchanged. Any error introduced by this assumption should be manifest as a larger than predicted systematic shift in the free energy change when the repressor copy number is varied.

     A naïve hypothesis for the effect of a mutation in the DNA binding domain is that only the DNA binding energy is affected. This hypothesis appears to contradict the core principle of allostery in that ligand binding in one domain influences binding in another, suggesting that changing parameter modifies them all. The characteristic curves summarized in Fig. 2 give a means to discriminate between these two hypotheses by examining the change in the free energy. Using a single induction profile (white-faced points in Fig. 3), we estimated the DNA binding energy using Bayesian inferential methods, the details of which are thoroughly discussed in the Materials & Methods as well as in the supplemental Chapter 7. The shaded red region for each mutant in Fig. 3 represents the 95% credible region of this fit whereas all other shaded regions are 95% credible regions of the predictions for other repressor copy numbers. We find that redetermining only the DNA binding energy accurately captures the majority of the induction profiles, indicating that other parameters are unaffected. One exception is for the lowest repressor copy numbers (R = 60 and R = 124 per cell) of mutant Q18A at low concentrations of IPTG. However, we note that this disagreement is comparable to that observed for the wild-type repressor binding to the weakest operator in Razo-Mejia et al. (2018), illustrating that our model is imperfect in characterizing weakly repressing architectures. Including other parameters in the fit (such as Δε_AI) does not significantly improve the accuracy of the predictions. Furthermore, the magnitude of this disagreement also depends on the choice of the fitting strain (see supplemental Chapter 7).

     Mutations Y17I and Q18A both weaken the affinity of the repressor to the DNA relative to the wild type strain with binding energies of  − 9.9_− 0.1^+ 0.1 k_BT and  − 11.0_− 0.1^+ 0.1 k_BT, respectively. Here we report the median of the inferred posterior probability distribution with the superscripts and subscripts corresponding to the upper and lower bounds of the 95% credible region. These binding energies are comparable to that of the wild-type repressor affinity to the native LacI operator sequence O3, with a DNA binding energy of  − 9.7 k_BT. The mutation Q18M increases the strength of the DNA-repressor interaction relative to the wild-type repressor with a binding energy of  − 15.43_− 0.06^+ 0.07 k_BT, comparable to the affinity of the wild-type repressor to the native O1 operator sequence ( − 15.3 k_BT). It is notable that a single amino acid substitution of the repressor is capable of changing the strength of the DNA binding interaction well beyond that of many single base-pair mutations in the operator sequence (Barnes et al. 2019).

     Using the new DNA binding energies, we can collapse all measurements of fold-change as a function of the free energy as shown in Fig. 3 (B). This allows us to test the diagnostic power of the decomposition of the free energy described in Fig. 2. To compute the ΔF for each mutation, we inferred the observed mean free energy of the mutant strain for each inducer concentration and repressor copy number (see Materials & Methods as well as the supplemental Chapter 7 for a detailed explanation of the inference). We note that in the limit of extremely low or high fold-change, the inference of the free energy is either over- or under-estimated, respectively, introducing a systematic error. Thus, points which are close to these limits are omitted in the calculation of ΔF. We direct the reader to the supplemental Chapter 7 for a detailed discussion of this systematic error. With a measure of F^(mut) for each mutant at each repressor copy number, we compute the difference in free energy relative to the wild-type strain with the same repressor copy number and operator sequence, restricting all variability in ΔF solely to changes in Δε_RA.

     The change in free energy for each mutant is shown in Fig. 3 (C). It can be seen that the ΔF for each mutant is constant as a function of the inducer concentration and is concordant with the prediction generated from fitting Δε_RA to a single repressor copy number (orange lines Fig. 3 (C)]) This is in line with the predictions outlined in Fig. 2 (C) and (D), indicating that the allosteric parameters are “insulated,” meaning they are not affected by the DNA binding domain mutations. As the ΔF for all repressor copy numbers collapses onto the prediction, we can say that the expression of the repressor itself is the same or comparable with that of the wild type. If the repressor copy number were perturbed in addition to Δε_RA, one would expect a shift away from the prediction that scales logarithmically with the change in repressor copy number. However, as the ΔF is approximately the same for each repressor copy number, it can be surmised that the mutation does not significantly change the expression or folding efficiency of the repressor itself. These results allow us to state that the DNA binding energy Δε_RA is the only parameter modified by the DNA mutants examined.

Inducer Binding Domain Mutants

Much as in the case of the DNA binding mutants, we cannot safely assume a priori that a given mutation in the inducer binding domain affects only the inducer binding constants K_A and K_I. While it is easy to associate the inducer binding constants with the inducer binding domain, the critical parameter in our allosteric model Δε_AI is harder to restrict to a single spatial region of the protein. As K_A, K_I, and Δε_AI are all parameters dictating the allosteric response, we consider two hypotheses in which inducer binding mutations alter either all three parameters or only K_A and K_I.

     We made four point mutations within the inducer binding domain of LacI (F161T, Q291V, Q291R, and Q291K) that have been shown previously to alter binding to multiple allosteric effectors (Daber, Sharp, and Lewis 2009). In contrast to the DNA binding domain mutants, we paired the inducer binding domain mutations with the three native LacI operator sequences (which have various affinities for the repressor) and a single ribosomal binding site sequence. This ribosomal binding site sequence, as reported in Garcia and Phillips (2011), expresses the wild-type LacI repressor to an average copy number of approximately 260 per cell. As the free energy differences resulting from point mutations in the DNA binding domain can be described solely by changes to Δε_RA, we continue under the assumption that the inducer binding domain mutations do not significantly alter the repressor copy number.

     The induction profiles for these four mutants are shown in Fig. 4 (A). Of the mutations chosen, Q291R and Q291K appear to have the most significant impact, with Q291R abolishing the characteristic sigmoidal titration curve entirely. It is notable that both Q291R and Q291K have elevated expression in the absence of inducer compared to the other two mutants paired with the same operator sequence. Panel (A) in Fig. 2 illustrates that if only K_A and K_I were being affected by the mutations, the fold-change should be identical for all mutants in the absence of inducer. This discrepancy in the observed leakiness immediately suggests that more than K_A and K_I are affected for Q291K and Q291R.

     Using a single induction profile for each mutant (shown in Fig. 4 as white-faced circles), we inferred the parameter combinations for both hypotheses and drew predictions for the induction profiles with other operator sequences. We find that the simplest hypothesis (in which only K_A and K_I are altered) does not permit accurate prediction of most induction profiles. These curves, shown as dotted lines in Fig. 4 (A), fail spectacularly in the case of Q291R and Q291K, and undershoot the observed profiles for F161T and Q291V, especially when paired with the weak operator sequence O3. The change in the leakiness for Q291R and Q291K is particularly evident as the expression at c = 0 should be identical to the wild-type repressor under this hypothesis. Altering only K_A and K_I is not sufficient to accurately predict the induction profiles for F161T and Q291V, but not to the same degree as Q291K and Q291R. The disagreement is most evident for the weakest operator O3 green lines in 4 (A), though we have discussed previously that the induction profiles for weak operators are difficult to accurately describe and can result in comparable disagreement for the wild-type repressor (Razo-Mejia et al. 2018).

     Including Δε_AI as a perturbed parameter in addition to K_A and K_I improves the predicted profiles for all four mutants. By fitting these three parameters to a single strain, we are able to accurately predict the induction profiles of other operators as seen by the shaded lines in Fig. 4 (A). With these modified parameters, all experimental measurements collapse as a function of their free energy as prescribed by Eq. 3 (Fig. 4 (B)). All four mutations significantly diminish the binding affinity of both states of the repressor to the inducer, as seen by the estimated parameter values reported in Table 3.1. As evident in the data alone, Q291R abrogates inducibility outright (K_A ≈ K_I). For Q291K, the active state of the repressor can no longer bind inducer whereas the inactive state binds with weak affinity. The remaining two mutants, Q291V and F161T, both show diminished binding affinity of the inducer to both the active and inactive states of the repressor relative to the wild-type.

Inferred values of K_A, K_I, and Δε_AI for inducer binding mutants

Mutant	KA	KI	ΔεAI [kBT]	Reference
WT	139 − 22 + 29 μM	0.53 − 0.04 + 0.04 μM	4.5	Razo-Mejia et al. (2018)
F161T	165 − 65 + 90 μM	3 − 3 + 6 μM	1 − 2 + 5	This study
Q291V	650 − 250 + 450 μM	8 − 8 + 8 μM	3 − 3 + 6	This study
Q291K	> 1 mM	310 − 60 + 70 μM	− 3.11 − 0.07 + 0.07	This study
Q291R	9 − 9 + 20 μM	8 − 8 + 20 μM	− 2.35 − 0.09 + 0.01	This study

     Given the collection of fold-change measurements, we computed the ΔF relative to the wild-type strain with the same operator and repressor copy number. This leaves differences in p_act(c) as the sole contributor to the free energy difference, assuming our hypothesis that K_A, K_I, and Δε_AI are the only perturbed parameters is correct. The change in free energy can be seen in Fig. 4 (C). For all mutants, the free energy difference inferred from the observed fold-change measurements falls within error of the predictions generated under the hypothesis that K_A, K_I, and Δε_AI are all affected by the mutation (shaded curves in Fig. 4 (C)). The profile of the free energy change exhibits some of the rich phenomenology illustrated in Fig. 2 (A) and (B). Q291K, F161T, and Q291V exhibit a non-monotonic dependence on the inducer concentration, a feature that can only appear when K_A and K_I are altered. The non-zero ΔF at c = 0 for Q291R and Q291K coupled with an inducer concentration dependence is a telling sign that Δε_AI must be significantly modified. This shift in ΔF is positive in all cases, indicating that Δε_AI must have decreased, and that the inactive state has become more energetically favorable for these mutants than for the wild-type protein. Indeed the estimates for Δε_AI (Table 3.1) reveal both mutations Q291R and Q291K make the inactive state more favorable than the active state. Thus, for these two mutations, only  ≈ 10% of the repressors are active in the absence of inducer, whereas the basal active fraction is  ≈ 99% for the wild-type repressor (Razo-Mejia et al. 2018).

     We note that the parameter values reported here disagree with those reported in Daber, Sochor, and Lewis (2011). This disagreement stems from different assumptions regarding the residual activity of the repressor in the absence of inducer and the parametric degeneracy of the MWC model without a concrete independent measure of Δε_AI. A detailed discussion of the difference in parameter values between our previous work (Garcia and Phillips 2011; Razo-Mejia et al. 2018), that of Daber, Sochor, and Lewis (2011), and those of other seminal works can be found in the supplemental Chapter 7.

     Taken together, these parametric changes diminish the response of the regulatory architecture as a whole to changing inducer concentrations. They furthermore reveal that the parameters which govern the allosteric response are interdependent and no single parameter is insulated from the others. However, as only the allosteric parameters are changed, one can say that the allosteric parameters as a whole are insulated from the other components which define the regulatory response, such as repressor copy number and DNA binding affinity.

Predicting Effects of Pairwise Double Mutations

Given full knowledge of each mutation individually, we can draw predictions of the behavior of the pairwise double mutants with no free parameters based on the simplest null hypothesis of no epistasis. The formalism of ΔF defined by Eq. 5 explicitly states that the contribution to the free energy of the system from the difference in DNA binding energy and the allosteric parameters are strictly additive. Thus, deviations from the predicted change in free energy would suggest epistatic interactions between the two mutations.

      To test this additive model, we constructed nine double mutant strains, each having a unique inducer binding (F161T, Q291V, Q291K) and DNA binding mutation (Y17I, Q18A, Q18M). To make predictions with an appropriate representation of the uncertainty, we computed a large array of induction profiles given random draws from the posterior distribution for the DNA binding energy (determined from the single DNA binding mutants) as well as from the joint posterior for the allosteric parameters (determined from the single inducer binding mutants). These predictions, shown in Fig. 5 (A) and (B) as shaded blue curves, capture all experimental measurements of the fold-change (Fig. 5 (A)) and the inferred difference in free energy (Fig. 5 (B)). The latter indicates that there are no epistatic interactions between the mutations queried in this work, though if there were, systematic deviations from these predictions would shed light on how the epistasis is manifested.

     The precise agreement between the predictions and measurements for Q291K paired with either Q18A or Q18M is striking as Q291K drastically changed Δε_AI in addition to K_A and K_I. Our ability to predict the induction profile and free energy change underscores the extent to which the DNA binding energy and the allosteric parameters are insulated from one another. Despite this insulation, the repressor still functions as an allosteric molecule, emphasizing that the mutations we have inserted do not alter the pathway of communication between the two domains of the protein. As the double mutant Y17I-Q291K exhibits fold-change of approximately 1 across all IPTG concentrations (Fig. 5 (A)), these mutations in tandem make repression so weak that it is beyond the limits which are detectable by our experiments. As a consequence, we are unable to estimate ΔF nor experimentally verify the corresponding prediction (grey box in Fig. 5 (B)). However, as the predicted fold-change in gene expression is also approximately 1 for all c, we believe that the prediction shown for ΔF is likely accurate. One would be able to infer the ΔF to confirm these predictions using a more sensitive method for measuring the fold-change, such as single-cell microscopy or colorimetric assays.

Discussion

Allosteric regulation is often couched as “biological action at a distance.” Despite extensive knowledge of protein structure and function, it remains difficult to translate the coordinates of the atomic constituents of a protein to the precise parameter values which define the functional response, making each mutant its own intellectual adventure. Bioinformatic approaches to understanding the sequence-structure relationship have permitted us to examine how the residues of allosteric proteins evolve, revealing conserved regions which hint to their function. Co-evolving residues reveal sectors of conserved interactions which traverse the protein that act as the allosteric communication channel between domains (Süel et al. 2003; McLaughlin Jr et al. 2012; Reynolds, McLaughlin, and Ranganathan 2011). Elucidating these sectors has advanced our understanding of how distinct domains “talk” to one another and has permitted direct engineering of allosteric responses into non-allosteric enzymes (Poelwijk et al. 2011; Raman, White, and Ranganathan 2016). Even so, we are left without a quantitative understanding of how these admittedly complex networks set the energetic difference between active and inactive states or how a given mutation influences binding affinity. In this context, a biophysical model in which the various parameters are intimately connected to the molecular details can be of use and can lead to quantitative predictions of the interplay between amino-acid identity and system-level response.

     By considering how each parameter contributes to the observed change in free energy, we are able to tease out different classes of parameter perturbations which result in stereotyped responses to changing inducer concentration. These characteristic changes to the free energy can be used as a diagnostic tool to classify mutational effects. For example, we show in Fig. 2 that modulating the inducer binding constants K_A and K_I results in non-monotonic free energy changes that are dependent on the inducer concentration, a feature observed in the inducer binding mutants examined in this work. Simply looking at the inferred ΔF as a function of inducer concentration, which requires no fitting of the biophysical parameters, indicates that K_A and K_I must be modified considering those are the only parameters which can generate such a response.

     Another key observation is that a perturbation to only K_A and K_I requires that the ΔF = 0 k_BT at c = 0. Deviations from this condition imply that more than the inducer binding constants must have changed. If this shift in ΔF off of 0 k_BT at c = 0 is not constant across all inducer concentrations, we can surmise that the energy difference between the allosteric states Δε_AI must also be modified. We again see this effect for all of our inducer mutants. By examining the inferred ΔF, we can immediately say that in addition to K_A and K_I, Δε_AI must decrease relative to the wild-type value as ΔF > 0 at c = 0. When the allosteric parameters are fit to the induction profiles, we indeed see that this is the case, with all four mutations decreasing the energy gap between the active and inactive states. Two of these mutations, Q291R and Q291K, make the inactive state of the repressor more stable than the active state, which is not the case for the wild-type repressor (Razo-Mejia et al. 2018).

      Our formulation of ΔF indicates that shifts away from 0 k_BT that are independent of the inducer concentration can only arise from changes to the repressor copy number and/or DNA binding specificity, indicating that the allosteric parameters are untouched. We see that for three mutations in the DNA binding domain, ΔF is the same irrespective of the inducer concentration. Measurements of ΔF for these mutants with repressor copy numbers across three orders of magnitude yield approximately the same value, revealing that Δε_RA is the sole parameter altered via the mutations.

     We note that the conclusions stated above can be qualitatively drawn without resorting to fitting various parameters and measuring the goodness-of-fit. Rather, the distinct behavior of ΔF is sufficient to determine which parameters are changing. Here, these conclusions are quantitatively confirmed by fitting these parameters to the induction profile, which results in accurate predictions of the fold-change and ΔF for nearly every strain across different mutations, repressor copy numbers, and operator sequence, all at different inducer concentrations. With a collection of evidence as to what parameters are changing for single mutations, we put our model to the test and drew predictions of how double mutants would behave both in terms of the titration curve and free energy profile.

     A hypothesis that arises from our formulation of ΔF is that a simple summation of the energetic contribution of each mutation should be sufficient to predict the double mutants (as long as they are in separate domains). We find that such a calculation permits precise and accurate predictions of the double mutant phenotypes, indicating that there are no epistatic interactions between the mutations examined in this work. With an expectation of what the free energy differences should be, epistatic interactions could be understood by looking at how the measurements deviate from the prediction. For example, if epistatic interactions exist which appear as a systematic shift from the predicted ΔF independent of inducer concentration, one could conclude that DNA binding energy is not equal to that of the single mutation in the DNA binding domain alone. Similarly, systematic shifts that are dependent on the inducer concentration (i.e. not constant) indicate that the allosteric parameters must be influenced. If the expected difference in free energy is equal to 0 k_BT when c = 0, one could surmise that the modified parameter must not be Δε_AI nor Δε_RA as these would both result in a shift in leakiness, indicating that K_A and K_I are further modified.

     Ultimately, we present this work as a proof-of-principle for using biophysical models to investigate how mutations influence the response of allosteric systems. We emphasize that such a treatment allows one to boil down the complex phenotypic responses of these systems to a single-parameter description which is easily interpretable as a free energy. The general utility of this approach is illustrated in Fig. 6 where gene expression data from previous work along with all of the measurements presented in this work collapse onto the master curve defined by Eq. 3. While our model coarse grains many of the intricate details of transcriptional regulation into two states (one in which the repressor is bound to the promoter and one where it is not), it is sufficient to describe a swath of regulatory scenarios. As discussed in the supplemental Chapter 7, any architecture in which the transcription-factor bound and transcriptionally active states of the promoter can be separated into two distinct coarse-grained states can be subjected to such an analysis.

     Given enough parametric knowledge of the system, it becomes possible to examine how modifications to the parameters move the physiological response along this reduced one-dimensional parameter space. This approach offers a glimpse at how mutational effects can be described in terms of energy rather than Hill coefficients and arbitrary prefactors. While we have explored a very small region of sequence space in this work, coupling of this approach with high-throughput sequencing-based methods to query a library of mutations within the protein will shed light on the phenotypic landscape centered at the wild-type sequence. Furthermore, pairing libraries of protein and operator sequence mutants will provide insight as to how the protein and regulatory sequence coevolve, a topic rich with opportunity for a dialogue between theory and experiment.

Materials & Methods

Bacterial Strains and DNA Constructs

All wild-type strains from which the mutants were derived were generated in previous work from the Phillips group (Garcia and Phillips 2011; Razo-Mejia et al. 2018). Briefly, mutations were first introduced into the lacI gene of our pZS3*1-lacI plasmid (Garcia and Phillips 2011) using a combination of overhang PCR Gibson assembly as well as QuickChange mutagenesis (Agligent Technologies). The oligonucleotide sequences used to generate each mutant as well as the method are provided in the supplemental Chapter 7.

     For mutants generated through overhang PCR and Gibson assembly, oligonucleotide primers were purchased containing an overhang with the desired mutation and used to amplify the entire plasmid. Using the homology of the primer overhang, Gibson assembly was performed to circularize the DNA prior to electroporation into MG1655 E. coli cells. Integration of LacI mutants was performed with λ Red recombineering as described in Sharan et al. (2009) and Garcia and Phillips (2011).

      The mutants studied in this work were chosen from data reported in Daber, Sochor, and Lewis (2011). In selecting mutations, we looked for mutants which suggested moderate to strong deviations from the behavior of the wild-type repressor. We note that the variant of LacI used in this work has an additional three amino acids (Met-Val-Asn) added to the N-terminus than the canonical LacI sequence reported in Farabaugh (1978). To remain consistent with the field, we have identified the mutations with respect to their positions in the canonical sequence and those in Daber, Sochor, and Lewis (2011). However, their positions in the raw data files correspond to that of our LacI variant and is noted in the README files associated with the data.

Flow Cytometry

All fold-change measurements were performed on a MACSQuant flow cytometer as described in Razo-Mejia et al. (2018). Briefly, saturated overnight cultures 500 μL in volume were grown in deep-well 96-well plates covered with a breathable nylon cover (Lab Pak - Nitex Nylon, Sefar America, Cat. No. 241205). After approximately 12 to 15 hr, the cultures 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 8 hours until the OD_600nm ≈ 0.3. 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.

     The flow cytometer was calibrated prior to use with MACSQuant Calibration Beads (Cat. No. 130-093-607). During measurement, the cultures were held at approximately 4^∘ C by placing the 96-well plate on a MACSQuant ice block. All fluorescence measurements were made using a 488 nm excitation wavelength with a 525/50 nm emission filter. The photomultiplier tube voltage settings for the instrument are the same as those used in Razo-Mejia et al. (2018), and are listed in supplemental Chapter 6.

     The data were processed using an automatic unsupervised gating procedure based on the front and side-scattering values, where we fit a two-dimensional Gaussian function to the log₁₀ forward-scattering (FSC) and the log₁₀ side-scattering (SSC) data. Here we assume that the region with highest density of points in these two channels corresponds to single-cell measurements and consider data points that fall within 40% of the highest density region of the two-dimensional Gaussian function. We direct the reader to Razo-Mejia et al. (2018) and supplemental Chapter 6 for further detail and comparison of flow cytometry with single-cell microscopy.

Bayesian Parameter Estimation

We used a Bayesian definition of probability in the statistical analysis of all mutants in this work. In supplemental Chapter 7, we derive in detail the statistical models used for the various parameters as well as multiple diagnostic tests. Here, we give a generic description of our approach. To be succinct in notation, we consider a generic parameter θ which represents Δε_RA, K_A, K_I, and/or Δε_AI depending on the specific LacI mutant.

     As prescribed by Bayes’ theorem, we are interested in the posterior probability distribution
g(θ | y) ∝ f(y | θ)g(θ),   (7)
where we use g and f to represent probability densities over parameters and data, respectively, and y to represent a set of fold-change measurements. The likelihood of observing our dataset y given a value of θ is captured by f(y | θ). All prior information we have about the possible values of θ are described by g(θ).

     In all inferential models used in this work, we assumed that all experimental measurements at a given inducer concentration were normally distributed about a mean value μ dictated by Eq. 1 with a variance σ²,
$$ f(y\,\vert\, \theta) = {1 \over (2\pi\sigma^2)^{N/2}}\prod\limits_i^N \exp\left[-{(y_i - \mu(\theta))^2 \over 2\sigma^2}\right], \qquad(8)$$
where N is the number of measurements in the data set y.

     This choice of likelihood is justified as each individual measurement at a given inducer concentration is a biological replicate and independent of all other experiments. By using a Gaussian likelihood, we introduce another parameter σ. As σ must be positive and greater than zero, we define as a prior distribution a half-normal distribution with a standard deviation ϕ,
$$ g(\sigma) = { {1 \over \phi}\sqrt{2 \over \pi} }\exp\left[-{x \over 2\phi^2}\right]\,;\, x \geq 0, \qquad(9)$$
where x is a given range of values for σ. A standard deviation of ϕ = 0.1 was chosen given our knowledge of the scale of our measurement error from other experiments. As the absolute measurement of fold-change is restricted between 0 and 1, and given our knowledge of the sensitivity of the experiment, it is reasonable to assume that the error will be closer to 0 than to 1. Further justification of this choice of prior through simulation based methods are given in the supplemental Chapter 7. The prior distribution for θ is dependent on the parameter and its associated physical and physiological restrictions. Detailed discussion of our chosen prior distributions for each model can also be found in the supplemental Chapter 7.

     All statistical modeling and parameter inference was performed using Markov chain Monte Carlo (MCMC). Specifically, Hamiltonian Monte Carlo sampling was used as is implemented in the Stan probabilistic programming language (Carpenter et al. 2017). All statistical models saved as .stan models and can be accessed at the GitHub repository associated with this work (DOI: 10.5281/zenodo.2721798) or can be downloaded directly from the website associated with the publication this chapter is based on.

Inference of Free Energy From Fold-Change Data

While the fold-change in gene expression is restricted to be between 0 and 1, experimental noise can generate fold-change measurements beyond these bounds. To determine the free energy for a given set of fold-change measurements (for one unique strain at a single inducer concentration), we modeled the observed fold-change measurements as being drawn from a normal distribution with a mean μ and standard deviation σ. Using Bayes’ theorem, we can write the posterior distribution as
$$ g(\mu, \sigma\,\vert y) \propto g(\mu)g(\sigma){1\over(2\pi\sigma^2)^{N/2}}\prod\limits_i^N \exp\left[{-(y_i - \mu)^2 \over 2\sigma^2}\right] \qquad(10)$$
where y is a collection of fold-change measurements. The prior distribution for μ was chosen to be uniform between 0 and 1 while the prior on σ was chosen to be half normal, as written in Eq. 9. The posterior distribution was sampled independently for each set of fold-change measurements using MCMC. The .stan model for this inference is available on the paper website.

     For each MCMC sample of μ, the free energy was calculated as
F =  − log (μ^− 1−1)   (11)
which is simply the rearrangement of Eq. 3. Using simulated data, we determined that when μ < σ or (1 − μ) < σ, the mean fold-change in gene expression was over or underestimated for the lower and upper limit, respectively. This means that there are maximum and minimum levels of fold-changes that can be detected using flow cytometry which are set by the distribution of fold-change measurements resulting from various sources of day-to-day variation. This results in a systematic error in the calculation of the free energy, making proper inference beyond these limits difficult. This bounds the range in which we can confidently infer this quantity with flow cytometry. We hypothesize that more sensitive methods, such as single cell microscopy, colorimetric assays, or direct counting of mRNA transcripts via Fluorescence In Situ Hybridization (FISH) would improve the measurement of ΔF. We further discuss details of this limitation in the supplemental Chapter 7.

Data and Code Availability

All data was collected, stored, and preserved using the Git version control software. Code for data processing, analysis, and figure generation is available on the [GitHub repository] (https://www.github.com/rpgroup-pboc/mwc_mutants) or can be accessed via the paper website. Raw flow cytometry data is stored on the CaltechDATA data repository and can be accessed via DOI 10.22002/D1.1241.

References