The Energetics of Molecular Adaptation in Transcriptional Regulation

TL;DR

This study introduces a framework linking mutations in an allosteric transcriptional repressor to changes in free energy, enabling prediction of gene expression responses and epistatic interactions through a parameter-based model.

Contribution

The paper presents a novel energetic framework that quantitatively connects mutations to functional responses in transcriptional regulation, validated with experimental data on LacI mutants.

Findings

Mutations can be categorized by their energetic effects with characteristic curves.

Single mutations alter specific model parameters, such as DNA affinity or allosteric response.

Double mutant induction profiles can be accurately predicted from single mutant data.

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 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 a subset of 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.

Tables

Table 1. Table 1: Inferred values of K A subscript 𝐾 𝐴 K_{A} , K I subscript 𝐾 𝐼 K_{I} , and Δ ​ ε A ​ I Δ subscript 𝜀 𝐴 𝐼 \Delta\varepsilon_{AI} for inducer binding mutants

Mutant	$K_A$	$K_I$	$\Delta {\epsilon }_{AI}$ [$k_{B}T$]	Reference
WT	${139}_{-22}^{+29}\mu$M	${0.53}_{-0.04}^{+0.04}\mu$M	4.5	(10)
F164T	${165}_{-65}^{+90}\mu$M	$3_{-3}^{+6}\mu$M	$1_{-2}^{+5}$	This study
Q294V	${650}_{-250}^{+450}\mu$M	$8_{-8}^{+8}\mu$M	$3_{-3}^{+6}$	This study
Q294K	$>1$ mM	${310}_{-60}^{+70}\mu$M	$-{3.11}_{-0.07}^{+0.07}$	This study
Q294R	$9_{-9}^{+20}\mu$M	$8_{-8}^{+20}\mu$M	$-{2.35}_{-0.09}^{+0.01}$	This study

Equations

$$\text{fold-change}={\left(1+\frac{R_A}{N_{NS}}{e}^{-\beta \Delta {\epsilon }_{RA}}\right)}^{-1},$$

$$p_{\text{act}}(c)=\frac{{\left(1+\frac{c}{K_A}\right)}^n}{{\left(1+\frac{c}{K_A}\right)}^n+e^{-\beta \Delta {\epsilon }_{AI}}{\left(1+\frac{c}{K_I}\right)}^n}.$$

$$\text{fold-change}={\left(1+e^{-\beta F}\right)}^{-1},$$

$$F=-k_{B}T\log p_{\text{act}}(c)-k_{B}T\log \left(\frac{R}{N_{NS}}\right)+\Delta {\epsilon }_{RA}.$$

$$\Delta F=-k_{B}T\log \left(\frac{p_{\text{act}}^{(\mathrm{mut})}(c)}{p_{\text{act}}^{(\mathrm{wt})}(c)}\right)$$

$$+(\Delta {\epsilon }_{RA}^{(\mathrm{mut})}-\Delta {\epsilon }_{RA}^{(\mathrm{wt})}).$$

$$\underset{c\to \infty }{\lim }\log \left(\frac{p_{\text{act}}^{(\mathrm{mut})}}{p_{\text{act}}^{(\mathrm{wt})}}\right)\approx \log \left[\frac{1+e^{-\beta \Delta {\epsilon }_{AI}^{(\mathrm{wt})}}{\left(\frac{K_A^{(\mathrm{wt})}}{K_I^{(\mathrm{wt})}}\right)}^n}{1+e^{-\beta \Delta {\epsilon }_{AI}^{(\mathrm{wt})}}{\left(\frac{K_A^{(\mathrm{mut})}}{K_I^{(\mathrm{mut})}}\right)}^n}\right],$$

$$g(\theta |y)\propto f(y|\theta )g(\theta ),$$

$$f(y|\theta )=\frac{1}{(2\pi {\sigma }^2{)}^{N/2}}\prod _i^N\exp \left[-\frac{(y_i-\mu (\theta ){)}^2}{2{\sigma }^2}\right],$$

$$g(\sigma )=\frac{1}{\varphi }\sqrt{\frac{2}{\pi }}\exp \left[-\frac{x}{2{\varphi }^2}\right];x\ge 0,$$

$$g(\mu ,\sigma |y)\propto g(\mu )g(\sigma )\frac{1}{(2\pi {\sigma }^2{)}^{N/2}}\prod _i^N\exp \left[\frac{-(y_i-\mu {)}^2}{2{\sigma }^2}\right],$$

$$F=-\log \left({\mu }^{-1}-1\right)$$

Full text

\significancestatement

We present a biophysical model of allosteric transcriptional regulation that directly links the location of a mutation within a repressor to the biophysical parameters that describe its behavior. We explore the phenotypic space of a repressor with mutations in either the inducer binding or DNA binding domains. Using the LacI repressor in E. coli, we make sharp, falsifiable predictions and use this framework to generate a null hypothesis for how double mutants behave given knowledge of the single mutants. Linking mutations to the parameters which govern the system allows for quantitative predictions of how the free energy of the system changes as a result, permitting coarse graining of high-dimensional data into a single-parameter description of the mutational consequences.

\dropcap

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 variables (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13). These models quantitatively describe numerous properties of input-output functions, such as the leakiness, saturation, dynamic range, steepness of response, and the [$EC_{50}$] – 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 (10). 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 result in 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 first present a theoretical framework for understanding how mutations in the repressor affect different parameters and alter the free energy of the system. The multi-dimensional parameter space of the aforementioned thermodynamic models is highly degenerate with multiple combinations of parameter values yielding the same phenotypic response. This degeneracy can be subsumed into the free energy of the system, transforming the input-output function into a one-dimensional description with the form of a Fermi function (14, 15). 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.

Results

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 (4, 6, 10). While we direct the reader to the SI text for a complete derivation, 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

[TABLE]

where $R_{A}$ is the number of active repressors per cell, $N_{NS}$ is the number of non-specific binding sites for the repressor, $\Delta\varepsilon_{RA}$ is the binding energy of the repressor to its specific binding site relative to the non-specific background, and $\beta$ 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_{\text{act}}(c)$ which captures the allosteric nature of the repressor and encodes the probability a repressor is in the active (repressive) state rather than the inactive state for a given inducer concentration $c$, namely,

[TABLE]

Here, $K_{A}$ and $K_{I}$ are the dissociation constants of the inducer to the active and inactive repressor, $\Delta\varepsilon_{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 (1) as $R_{A}=p_{\text{act}}(c)R$.

A key feature of (1) and (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 (16)],

[TABLE]

where

[TABLE]

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 (10) collapse onto a single master curve. The tight distribution about this curve reveals that 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 (3) and (4) can be used as a tool to identify mechanistic changes in the regulatory architecture (14). 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 (14, 15).

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 (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^{\mathrm{(mut)}}$ differs from that of the wild-type $F^{\mathrm{(wt)}}$, we define $\Delta F=F^{\mathrm{(mut)}}-F^{\mathrm{(wt)}}$, which has the form

[TABLE]

$\Delta 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 (1) and (2) into three distinct quantities, each with different sensitivities to parametric changes. By examining how a mutation changes the free energy changes 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 paper website which makes exploration of parameter space a simpler task.

The first term in (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 $\Delta\varepsilon_{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 changes to the inducer binding constants $K_{A}$ and $K_{I}$ alone alter the induction profiles and resulting 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_{\text{act}}(c)$ given by (2), meaning that $\Delta\varepsilon_{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 $\Delta F$ should be $0\,k_{B}T$ when $c=0$, illustrated by the overlapping red, purple, and grey curves in the right-hand plot of Fig. 2(A). However, if $\Delta\varepsilon_{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 $\Delta F$ when $c=0$. This is illustrated in Fig. 2 (B) where $\Delta\varepsilon_{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 $\Delta 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,

[TABLE]

which illustrates that $\Delta F$ returns to zero at saturating inducer concentration when $K_{A}/K_{I}$ is the same for both the mutant and wild-type repressors, so long as $\Delta\varepsilon_{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.

The second term in (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 magnitude of the change in free energy scales logarithmically with changing repressor copy number, a mutation which increases expression from 1 to 10 repressors per cell is more impactful from an energetic standpoint ($k_{B}T\log(10)\approx 2.3\,k_{B}T$) than an increase from 90 to 100 ($k_{B}T\log(100/90)\approx 0.1\,k_{B}T$). 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 (17).

The third and final term in (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 ($\Delta\varepsilon_{RA}^{\mathrm{(wt)}}>\Delta\varepsilon_{RA}^{\mathrm{(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_{B}T$, the repressor copy number must change by a factor of $\approx 2.3$ whereas the DNA binding energy must change by $1\,k_{B}T$.

The unique behavior of each quantity in (5) and its sensitivity with respect to the parameters makes $\Delta F$ useful as a diagnostic tool to classify mutations. Given a set of fold-change measurements, a simple rearrangement of (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.

DNA Binding Domain Mutations

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 (18, 19, 20, 21). 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 (5). However, experimental variables such as the repressor copy number or the number of specific binding sites were not known, making precise calculation of $\Delta F$ as presented here not tractable. Using this dataset as a guide, we chose a subset of the mutations and inserted them into our experimental strains of E. coli where these parameters are known and tightly controlled (4, 10).

We made three amino acid substitutions (Y20I, Q21A, and Q21M) that are critical for the DNA-repressor interaction. These mutations were introduced into the lacI sequence used in Garcia and Phillips 2011 (4) 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 ($\Delta\varepsilon_{RA}=-13.9\,k_{B}T$). 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 (4) for the wild-type sequence paired with the same ribosomal binding site sequences. In characterizing the effects of these DNA binding mutations, we take the repressor copy number to be unchanged. Any error introduced by this mutation 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 altered. This hypothesis appears to contradict the core principle of allostery in that ligand binding in one domain influences binding in another, suggesting that changing any 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 a Bayesian approach, the details of which are discussed in the Materials and Methods as well as the SI text. 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 Q21A 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 (10), illustrating that our model is imperfect in characterizing weakly repressing architectures. Including other parameters in the fit (such as $\Delta\varepsilon_{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 SI text).

Mutations Y20I and Q21A both weaken the affinity of the repressor to the DNA relative to the wild type strain ($-9.9^{+0.1}_{-0.1}\,k_{B}T$ and $-11.0^{+0.1}_{-0.1}\,k_{B}T$, 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_{B}T$. The mutation Q21M increases the strength of the DNA-repressor interaction relative to the wild-type repressor with a binding energy of $-15.43^{+0.07}_{-0.06}\,k_{B}T$, comparable to the affinity of the wild-type repressor to the native O1 operator sequence ($-15.3k_{B}T$). 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 (22, 4).

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 $\Delta 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 and Methods as well as the SI text 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 $\Delta F$. We direct the reader to the SI text for a detailed discussion of this systematic error. With a measure of $F^{\mathrm{(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 $\Delta F$ solely to changes in $\Delta\varepsilon_{RA}$.

The change in free energy for each mutant is shown in Fig. 3(C). It can be seen that the $\Delta F$ for each mutant is constant as a function of the inducer concentration and is concordant with the prediction generated from fitting $\Delta\varepsilon_{RA}$ to a single repressor copy number [red 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 $\Delta 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 $\Delta\varepsilon_{RA}$, one would expect a shift away from the prediction that scales logarithmically with the change in repressor copy number. However, as the $\Delta 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 $\Delta\varepsilon_{RA}$ is the only parameter modified by the DNA mutants examined.

Inducer Binding Domain Mutations

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 $\Delta\varepsilon_{AI}$ is harder to restrict to a single spatial region of the protein. As $K_{A}$, $K_{I}$, and $\Delta\varepsilon_{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 (F164T, Q294V, Q294R, and Q294K) that have been shown previously to alter binding to multiple allosteric effectors (5). 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 (4), 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 $\Delta\varepsilon_{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, Q294R and Q294K appear to have the most significant impact, with Q294R abolishing the characteristic sigmoidal titration curve entirely. It is notable that both Q294R and Q294K 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 Q294K and Q294R.

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 Q294R and Q294K, and undershoot the observed profiles for F164T and Q294V, especially when paired with the weak operator sequence O3. The change in the leakiness for Q294R and Q294K 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 F164T and Q294V, but not to the same degree as Q294K and Q294R. The disagreement is most evident for the weakest operator O3 [green lines in Fig. 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 (10, 22).

Including $\Delta\varepsilon_{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 (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 Tab. LABEL:tab:ind_params. As evident in the data alone, Q294R abrogates inducibility outright ($K_{A}\approx K_{I}$). For Q294K, the active state of the repressor can no longer bind inducer whereas the inactive state binds with weak affinity. The remaining two mutants, Q294V and F164T, both show diminished binding affinity of the inducer to both the active and inactive states of the repressor relative to the wild-type.

Given the collection of fold-change measurements, we computed the $\Delta 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 $\Delta\varepsilon_{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 $\Delta\varepsilon_{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). Q294K, F164T, and Q294V 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 $\Delta F$ at $c=0$ for Q294R and Q294K coupled with an inducer concentration dependence is a telling sign that $\Delta\varepsilon_{AI}$ must be significantly modified. This shift in $\Delta F$ is positive in all cases, indicating that $\Delta\varepsilon_{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 $\Delta\varepsilon_{AI}$ (Tab. LABEL:tab:ind_params) reveal both mutations Q294R and Q294K make the inactive state more favorable than the active state. Thus, for these two mutations, only $\approx 10\%$ of the repressors are active in the absence of inducer, whereas the basal active fraction is $\approx 99\%$ for the wild-type repressor (10).

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 individual mutation, 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 $\Delta F$ defined by (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 (F164T, Q294V, Q294K) and DNA binding mutation (Y20I, Q21A, Q21M). 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 manifest.

The precise agreement between the predictions and measurements for Q294K paired with either Q21A or Q21M is striking as Q294K drastically changed $\Delta\varepsilon_{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 Y20I-Q294K exhibits fold-change of approximately $1$ across all IPTG concentrations [Fig. 5(A)], these mutations in tandem make repression so weak it is beyond the limits which are detectable by our experiments. As a consequence, we are unable to estimate $\Delta 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 $\Delta F$ is likely accurate. One would be able to infer the $\Delta 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 (23, 24, 25). 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 (26, 27, 28). 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 $\Delta 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 $\Delta F=0$ at $c=0$. Deviations from this condition imply that more than the inducer binding constants must have changed. If this shift in $\Delta F$ off of [math] at $c=0$ is not constant across all inducer concentrations, we can surmise that the energy difference between the allosteric states $\Delta\varepsilon_{AI}$ must also be modified. We again see this effect for all of our inducer mutants. By examining the inferred $\Delta F$, we can immediately say that in addition to $K_{A}$ and $K_{I}$, $\Delta\varepsilon_{AI}$ must decrease relative to the wild-type value as $\Delta 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, Q294R and Q294K, make the inactive state of the repressor more stable than the active state, which is not the case for the wild-type repressor (10).

Our formulation of $\Delta F$ indicates that shifts away from [math] 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, $\Delta F$ is the same irrespective of the inducer concentration. Measurements of $\Delta F$ for these mutants with repressor copy numbers across three orders of magnitude yield approximately the same value, revealing that $\Delta\varepsilon_{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 $\Delta 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 $\Delta 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 $\Delta F$ is that a simple summation of the energetic contribution of each mutation should be sufficient to predict the double mutants (so 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 $\Delta 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 [math] when $c=0$, one could surmise that the modified parameter must not be $\Delta\varepsilon_{AI}$ nor $\Delta\varepsilon_{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 (4, 6, 10) along with all of the measurements presented in this work collapse onto the master curve defined by (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 wide range of regulatory scenarios. 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.

\matmethods

Bacterial Strains and DNA Constructs

All wild-type strains from which the mutants were derived were generated in previous work from the Phillips group (4, 10). Briefly, mutations were first introduced into the lacI gene of our pZS3*1-lacI plasmid (4) 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 SI text.

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 $\lambda$ Red recombineering (29) as described in Reference (4).

The mutants studied in this work were chosen from data reported in (5). 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 (30). For this reason, all mutants given here are with respect to our sequence and their positions are shifted by three to those studied in (5).

Flow Cytometry

All fold-change measurements were performed on a MACSQuant flow cytometer as described in Razo-Mejia et al. (10). Briefly, saturated overnight cultures 500 $\mu$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 $\mu$L of M9 minimal media supplemented with 0.5% w/v glucose (anhydrous D-Glucose, Macron Chemicals) and the appropriate concentration of IPTG (Isopropyl $\beta$-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${}_{\text{600nm}}\approx 0.3$. Cells were then diluted ten-fold into a round-bottom 96-well plate (Corning Cat. No. 3365) containing 90 $\mu$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 Reference (10).

The data was 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_{10}$ forward-scattering (FSC) and the $\log_{10}$ 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 Reference (10) 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 the SI text, 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 $\theta$ which represents $\Delta\varepsilon_{RA}$, $K_{A}$, $K_{I}$, and/or $\Delta\varepsilon_{AI}$ depending on the specific LacI mutant.

As prescribed by Bayes’ theorem, we are interested in the posterior probability distribution

[TABLE]

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 $\theta$ is captured by $f(y\,|\,\theta)$. All prior information we have about the possible values of $\theta$ are described by $g(\theta)$.

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 $\mu$ dictated by (1) with a variance $\sigma^{2}$,

[TABLE]

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 $\sigma$. As $\sigma$ must be positive and greater than zero, we define as a prior distribution a half-normal distribution with a standard deviation $\phi$,

[TABLE]

where $x$ is a given range of values for $\sigma$. A standard deviation of $\phi=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 [math] and $1.0$, and given our knowledge of the sensitivity of the experiment, it is reasonable to assume that the error will be closer to [math] than to $1.0$. Further justification of this choice of prior through simulation based methods are given in the SI text. The prior distribution for $\theta$ 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 SI text.

All statistical modeling and parameter inference was performed using Markov chain Monte Carlo (MCMC). Specifically, Hamiltonian Monte Carlo sampling was used as sis implemented in the Stan probabilistic programming language (31). 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 paper website.

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 $\mu$ and standard deviation $\sigma$. Using Bayes’ theorem, we can write the posterior distribution as

[TABLE]

where $y$ is a collection of fold-change measurements. The prior distribution for $\mu$ was chosen to be uniform between 0 and 1 while the prior on $\sigma$ was chosen to be half normal, as written in (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 $\mu$, the free energy was calculated as

[TABLE]

which is simply the rearrangement of (3). Using simulated data, we determined that when $\mu<\sigma$ or $(1-\mu)<\sigma$, 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-change 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 $\Delta F$. We further discuss details of this limitation in the SI text.

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.

\showmatmethods

\acknow

We thank Pamela Björkman, Rachel Galimidi, and Priyanthi Gnanapragasam for access and training for the use of the Miltenyi Biotec MACSQuant flow cytometer. The experimental efforts first took place at the Physiology summer course at the Marine Biological Laboratory in Woods Hole, MA, operated by the University of Chicago. We thank Ambika Nadkarni and Damian Dudka for their work on the project during the course. We also thank Suzannah Beeler, Justin Bois, Robert Brewster, Soichi Hirokawa, Heun Jin Lee, and Muir Morrison for thoughtful advice and discussion. This work was supported by La Fondation Pierre-Gilles de Gennes, the Rosen Center at Caltech, the NIH DP1 OD0002179 (Director’s Pioneer Award), R01 GM085286, and 1R35 GM118043 (MIRA). Nathan M. Belliveau was supported by a Howard Hughes Medical Institute International Student Research fellowship.

\showacknow