GlnK facilitates the dynamic regulation of bacterial nitrogen assimilation

TL;DR

This study combines experiments and modeling to elucidate how GlnK and GlnB proteins regulate nitrogen assimilation in E. coli, revealing GlnK's role as a nitrogen shock buffer during starvation.

Contribution

It introduces a biochemical model of GlnK and GlnB interactions that predicts enzyme responses under dynamic nitrogen conditions, highlighting GlnK's regulatory role.

Findings

GlnK is regulated under nitrogen-rich conditions but expressed during starvation.

GlnK acts as a buffer of nitrogen shock after starvation.

The model accurately predicts responses of GS, GlnB, and GlnK in various conditions.

Abstract

Ammonium assimilation in E. coli is regulated by two paralogous proteins (GlnB and GlnK), which orchestrate interactions with regulators of gene expression, transport proteins and metabolic pathways. Yet how they conjointly modulate the activity of glutamine synthetase (GS), the key enzyme for nitrogen assimilation, is poorly understood. We combine experiments and theory to study the dynamic roles of GlnB and GlnK during nitrogen starvation and upshift. We measure time-resolved in vivo concentrations of metabolites, total and post-translationally modified proteins, and develop a concise biochemical model of GlnB and GlnK that incorporates competition for active and allosteric sites, as well as functional sequestration of GlnK. The model predicts the responses of GS, GlnB and GlnK under time-varying external ammonium level in the wild type and two genetic knock-outs. Our results show that GlnK is tightly regulated under nitrogen-rich conditions, yet it is expressed during ammonium run-out and starvation. This suggests a role for GlnK as a buffer of nitrogen shock after starvation, and provides a further functional link between nitrogen and carbon metabolisms.

Tables

Table 1. Table 1: Model parameters

	UT/UR			AT/AR
		Literature	Fitted		Literature	Fitted
Maximum enzyme rates	$V_{\text{UT}}$	-	119.7 $\mu$M/min	$V_{\text{AT}}$	-	976.9 $\mu$M/min
	$V_{\text{UR}}$	-	14.7 $\mu$M/min	$V_{\text{AR}}$	-	473.4 $\mu$M/min
Michaelis-Menten constants	$K_{\text{m,GlnB}}$	3 $\mu$M (ninfa3)	-	$K_{\text{m,GS}}$	-	133.8 $\mu$M
	$K_{\text{m,GlnB-U}}$	2 $\mu$M (ninfa3)	-	$K_{\text{m,GS-A}}$	-	56.1 $\mu$M
	$K_{\text{m,GlnK}}$	-	9.9 $\mu$M
	$K_{\text{m,GlnK-U}}$	-	59.8 $\mu$M
Dissociation constants	$K_{\text{GLN,1}}$	90 $\mu$M (ninfa3)	-	$\overline{K_{\text{GlnB}}}$	15.0 $\mu$M (ninfa9)	-
	$K_{\alpha \text{-KG,1}}$	-	3.0 mM	$\overline{K_{\text{GlnB-U}}}$	0.25 $\mu$M (ninfa9)	-
				$\overline{K_{\text{GlnK}}}$	-	8.3 $\mu$M
				$\overline{K_{\text{GlnK-U}}}$	-	9.7 $\mu$M
				$K_{\text{GLN,2}}$	-	1.4 mM
				$K_{\alpha \text{-KG,2}}$	-	0.89 mM
Hill coefficients	$n_{\text{UT}}$	3 (hwa)	-	$n_{\text{AT}}$	-	2
Fractional sequestration	$x$	-	0.4

Table 2. Table S1: Information criteria for model selection of uridylylation model with and without sequestration using the WT and Δ Δ \Delta glnB data ( n = 24 𝑛 24 n=24 points)

Model	Number of fitted parameters ($k$)	RSS	AICc	BIC
without sequestration	4	60.3	108.5	111.1
with sequestration	6	21.7	90.8	92.9

Table 3. Table S2: MRM-MS signals of GlnK unlabelled/labelled signature peptides

Protein	Peptide	Internal standard (Y/N)	Q1	Q3	Retention	Collision energy
GlnK-1	GAEYSVNFLPK-1a	N	612.8	804.5	35.8	30
GlnK-1	GAEYSVNFLPK-1b	N	612.8	967.5	35.8	30
GlnK-1	GAEYSVNFLPK-1c	N	612.8	244.2	35.8	30
GlnK-1	GAEYSVNFLPK-1a-is	Y	616.8	812.5	35.8	30
GlnK-1	GAEYSVNFLPK-1b-is	Y	616.8	975.5	35.8	30
GlnK-1	GAEYSVNFLPK-1c-is	Y	616.8	252.2	35.8	30
GlnK(u)-1	GAEY(U)SVNFLPK-1a	N	765.9	244.2	33.2	30
GlnK(u)-1	GAEY(U)SVNFLPK-1b	N	765.9	804.5	33.2	30
GlnK(u)-1	GAEY(U)SVNFLPK-1a-is	Y	769.9	252.2	33.2	30
GlnK(u)-1	GAEY(U)SVNFLPK-1b-is	Y	769.9	812.5	33.2	30
GlnK-2	IFVAELQR-2a	N	488.3	715.4	31.8	30
GlnK-2	IFVAELQR-2b	N	488.3	616.3	31.8	30
GlnK-2	IFVAELQR-2c	N	488.3	862.5	31.8	30
GlnK-2	IFVAELQR-2a-is	Y	493.3	725.4	31.8	30
GlnK-2	IFVAELQR-2b-is	Y	493.3	626.3	31.8	30
GlnK-2	IFVAELQR-2c-is	Y	493.3	872.5	31.8	30

Equations

$$\frac{\mathrm{d}[\text{GlnB-U}]}{\mathrm{d}t}$$

$$\frac{\mathrm{d}[\text{GlnK-U}]}{\mathrm{d}t}$$

$$\frac{\mathrm{d}[\text{GS-A}]}{\mathrm{d}t}$$

$$v_{\text{UT,GlnB}}=\frac{V_{\text{UT}}\frac{[\text{GlnB}]}{K_{\text{m,GlnB}}}}{\left(1+\frac{[\text{GlnB}]}{K_{\text{m,GlnB}}}+\frac{[\text{GlnK}]}{K_{\text{m,GlnK}}}\right)\left(1+\frac{[\text{GLN}]}{K_{\text{GLN,1}}}\right)}$$

$$v_{\text{UR,GlnB-U}}=\frac{V_{\text{UR}}\frac{[\text{GlnB-U}]}{K_{\text{m,GlnB-U}}}}{\left(1+\frac{[\text{GlnB-U}]}{K_{\text{m,GlnB-U}}}+\frac{[\text{GlnK-U}]}{K_{\text{m,GlnK-U}}}\right)\left(1+\frac{K_{\text{GLN,1}}}{[\text{GLN}]}\right)}$$

$$v_{\text{UT,GlnK}}=\frac{V_{\text{UT}}\frac{[\text{GlnK}]}{K_{\text{m,GlnK}}}}{\left(1+\frac{[\text{GlnB}]}{K_{\text{m,GlnB}}}+\frac{[\text{GlnK}]}{K_{\text{m,GlnK}}}\right)\left(1+\frac{[\text{GLN}]}{K_{\text{GLN,1}}}\right)}$$

$$v_{\text{UR,GlnK-U}}=\frac{V_{\text{UR}}\frac{[\text{Glnk-U}]}{K_{\text{m,Glnk-U}}}}{\left(1+\frac{[\text{GlnB-U}]}{K_{\text{m,GlnB-U}}}+\frac{[\text{GlnK-U}]}{K_{\text{m,GlnK-U}}}\right)\left(1+\frac{K_{\text{GLN,1}}}{[\text{GLN}]}\right)}$$

$$[\text{GlnB}]=[{\text{GlnB}}_{\text{total}}](t)-[\text{GlnB-U}]$$

$$[\text{GlnK}]=[{\text{GlnK}}_{\text{total}}]\left(x+(1-x)\frac{[\alpha \text{-KG}{]}^{n_{\text{UT}}}}{K_{\text{α-KG,1}}^{n_{\text{UT}}}+[\alpha \text{-KG}{]}^{n_{\text{UT}}}}\right)$$

$$K_{\text{GlnB}}=\overline{K_{\text{GlnB}}}\left(1+\frac{[\alpha \text{-KG}]}{K_{\text{α-KG,2}}}\right)$$

$$K_{\text{GlnB-U}}=\overline{K_{\text{GlnB-U}}}\left(1+\frac{K_{\text{α-KG,2}}}{[\alpha \text{-KG}]}\right).$$

$$v_{\text{AT}}$$

$$+\frac{1}{\left(1+\frac{K_{\text{GlnK}}}{[\text{GlnK}]}+\frac{[\text{GlnB}]}{K_{\text{GlnB}}}\frac{K_{\text{GlnK}}}{[\text{GlnK}]}\right)}]\frac{1}{1+\frac{K_{\text{GLN,2}}}{[\text{GLN}]}}$$

$$v_{\text{AR}}$$

$$+\frac{1}{\left(1+\frac{K_{\text{GlnK-U}}}{[\text{GlnK-U}]}+\frac{[\text{GlnB-U}]}{K_{\text{GlnB-U}}}\frac{K_{\text{GlnK-U}}}{[\text{GlnK-U}]}\right)}]\frac{1}{1+\frac{[\text{GLN}]}{K_{\text{GLN,2}}}}$$

$$J(p)=\sum _i\sum _j\Vert \frac{M^{\text{mod}}}{M^{\text{total}}}-\frac{F^{\text{mod}}(p)}{M^{\text{total}}}\Vert \left(\frac{M^{\text{total}}}{M^{\text{error}}}\right),$$

$$\mathrm{GlnB}+{\mathrm{E}}_1$$

$$\mathrm{GlnK}+{\mathrm{E}}_1$$

$$GlnB-U+E_2$$

$$GlnK-U+{\mathrm{E}}_2$$

$$v_{\text{UT,GlnB}}$$

$$v_{\text{UT,GlnK}}$$

$$v_{\text{UR,GlnB-U}}$$

$$v_{\text{UR,GlnK-U}}$$

$$v_{\text{UT,GlnB}}$$

$$v_{\text{UT,GlnK}}$$

$$v_{\text{UR,GlnB-U}}$$

$$v_{\text{UR,GlnK-U}}$$

$$(e_0-c_1-c_2-c_3-c_4-c_5)s$$

$$(e_0-c_1-c_2-c_3-c_4-c_5)a$$

$$(e_0-c_1-c_2-c_3-c_4-c_5)b$$

$$c_{3}s$$

$$c_{4}s$$

$$c_{1}a$$

$$c_{1}b$$

$$\left(e_0-\frac{K_2}{a}{c}_2-c_2-\frac{K_1}{s}{c}_2-\frac{K_1}{s}\frac{b}{K_3}\frac{K_2}{a}{c}_2-\frac{b}{K_3}\frac{K_2}{a}{c}_2\right)s$$

$$\left(e_0-\frac{K_3}{b}{c}_5-\frac{a}{K_2}\frac{K_3}{b}{c}_5-\frac{K_1}{s}\frac{a}{K_2}\frac{K_3}{b}{c}_5-\frac{K_1}{s}{s}_5-c_5\right)a$$

$$c_2$$

$$c_5$$

$$v=k_1{c}_2+k_1^{\prime }{c}_5=\frac{1}{\left(1+\frac{K_1}{s}\right)}\left[\frac{k_1{e}_0}{(1+\frac{K_2}{a}+\frac{b}{K_3}\frac{K_2}{a})}+\frac{k_1^{\prime }{e}_0}{(1+\frac{K_3}{b}+\frac{a}{K_2}\frac{K_3}{b})}\right].$$

$$\left(e_0-c_1-\frac{b}{K_3}{c}_1-\frac{K_1}{s}\frac{b}{K_3}{c}_1-\frac{K_1}{s}\frac{a}{K_2}{c}_1-\frac{a}{K_2}{c}_1\right)s$$

$$v=k_1{c}_1=\frac{e_0}{\left(\frac{K_1}{s}+1\right)\left(1+\frac{a}{K_2}+\frac{b}{K_3}\right)}.$$

$$E_{\mathcal{D}}(\theta )=\underset{\theta }{\min }\left\{\sum _{j=1}^m\Vert (\mathbf{X}(t_i;\theta )-\tilde{\mathbf{X}}(t_i))/S{E}_i\Vert \right\},$$

Full text

GlnK facilitates the dynamic regulation of bacterial nitrogen assimilation

Adam Gosztolai${}^{\text{}1}$ , Jörg Schumacher${}^{\text{}2}$ , Volker Behrends${}^{\text{}3,4}$ , Jacob G Bundy${}^{\text{}4}$ ,

Franziska Heydenreich${}^{\text{}2,}\thanks{Current address: Laboratory of Biomolecular Research, Paul Scherrer Institute, 5232, Villigen, Switzerland and Department of Biology, ETH Zurich, 8093, Zurich, Switzerland}$ , Mark H Bennett${}^{\text{}2}$ , Martin Buck${}^{\text{}2}$ , Mauricio Barahona${}^{\text{}1,}\thanks{Correspondence: m.barahona@imperial.ac.uk.}$ Current address: Laboratory of Biomolecular Research, Paul Scherrer Institute, 5232, Villigen, Switzerland and Department of Biology, ETH Zurich, 8093, Zurich, SwitzerlandCorrespondence: m.barahona@imperial.ac.uk.

1. Department of Mathematics, Imperial College London, SW7 2AZ, London, United Kingdom.

2. Department of Life Sciences, Imperial College London, SW7 2AZ, London, United Kingdom.

3. Department of Life Sciences, University of Roehampton, SW15 5PU, London, United Kingdom.

4. Department of Surgery and Cancer, Imperial College London, SW7 2AZ, London, United Kingdom

Abstract

Ammonium assimilation in E. coli is regulated by two paralogous proteins (GlnB and GlnK), which orchestrate interactions with regulators of gene expression, transport proteins and metabolic pathways. Yet how they conjointly modulate the activity of glutamine synthetase (GS), the key enzyme for nitrogen assimilation, is poorly understood. We combine experiments and theory to study the dynamic roles of GlnB and GlnK during nitrogen starvation and upshift. We measure time-resolved in vivo concentrations of metabolites, total and post-translationally modified proteins, and develop a concise biochemical model of GlnB and GlnK that incorporates competition for active and allosteric sites, as well as functional sequestration of GlnK. The model predicts the responses of GS, GlnB and GlnK under time-varying external ammonium level in the wild type and two genetic knock-outs. Our results show that GlnK is tightly regulated under nitrogen-rich conditions, yet it is expressed during ammonium run-out and starvation. This suggests a role for GlnK as a buffer of nitrogen shock after starvation, and provides a further functional link between nitrogen and carbon metabolisms.

Introduction

To adapt to the highly variable nutrient conditions in their natural habitat, microorganisms have evolved a complex intracellular circuitry coupling signal transduction, membrane transport, gene expression and metabolism. A widely conserved example is the ammonium assimilatory system in Escherichia coli, which coordinates the uptake of nitrogen with other pathways involved in carbon assimilation and maintenance of cellular energy status (reitzer, ninfabook, heeswijk). External ammonium (E. coli’s preferred nitrogen source) is assimilated into glutamate and glutamine in two steps: first, ammonium and intracellular glutamate (GLU) are converted into glutamine (GLN) in a reaction catalysed by the enzyme glutamine synthetase (GS); second, glutamine is combined with $\alpha$-ketogluterate ($\alpha$-KG), a product of carbon metabolism, to yield two molecules of glutamate (for a net production of one GLU molecule) in a reaction catalysed by the enzyme glutamate synthase (GOGAT) (Fig. 1).

Although glutamine and glutamate are both key intermediates towards the production of amino acids and nucleotides, their wider roles in cellular metabolism are different. Whereas glutamate is the nitrogen source for most ($88\%$) reactions in the cell (reitzer), glutamine constitutes the primary signal controlling the activity of enzymes in the ammonium assimilation system. As a result, their regulatory requirements are different: homeostatic regulation of glutamate is essential to maintain growth rate, whereas glutamine levels must adjust rapidly to reflect cellular nitrogen status while at the same time remaining within physiologically acceptable bounds. In particular, when nutrient-rich conditions are suddenly restored after a period of starvation, the cell must be able to avoid glutamine shock (kustu, wingreen).

The modulation of GS activity is the primary means by which E. coli achieves these regulatory targets in the face of changing nitrogen levels. On one hand, as extracellular ammonium is used up during growth, the activity and abundance of GS is raised so as to maintain the growth rate as long as possible. On the other hand, when favourable conditions are suddenly restored, GS can be ‘turned off’. This nonlinear response of GS is controlled by the bifunctional enzyme adenylyl-transferase/adenylyl-removase (AT/AR), which can reversibly adenylylate (inactivate) or deadenylylate (activate) independently each of the 12 monomers that conform GS, thus modulating GS activity in a graded manner (stadtman). The activity of AT/AR is in turn regulated allosterically by glutamine and by two highly similar signal transduction proteins, GlnB and GlnK (encoded, respectively, by the glnB and glnK genes) (ninfa9, ninfabook). Importantly, GlnB and GlnK share a source of upstream regulation, themselves being reversibly uridylylated by a second bifunctional enzyme, uridilyl-transferase/uridylyl-removase (UT/UR), which enables GlnB and GlnK to control the enzyme AT/AR (ninfa1). Specifically, under glutamine limitation GlnB and GlnK are mostly uridylylated, leading to increased AR activity and inhibition of AT activity (Fig. 1). In addition to regulating AT/AR, GlnB and GlnK also regulate the expression of GS and bind $\alpha$-KG, communicating with carbon metabolism and signal transduction pathways (ninfabook).

Due to the apparent redundancy of having two signalling proteins (GlnB and GlnK) relaying the glutamine state of the cell, previous models (Bruggeman, doucette, ninfa3, straube) have ignored GlnK regulation of AT/AR, it being less potent than GlnB. This simplification is supported by mutational studies of the receptor interaction domains of GlnB and GlnK (ninfa10) and by in vitro studies (ninfa7). From this viewpoint, the role of GlnK is mostly circumscribed to regulating the membrane-bound channel protein AmtB, whose gene (amtB) lies in the same operon as glnK, causing glnK and amtB expression to be directly correlated (merrick3, merrick1, hwa).

However, various experiments suggest that GlnK also acts as an effective regulator of AT/AR. For instance, GlnB deficient strains (i.e., $\Delta$glnB mutants) exhibit unimpaired growth (ninfa8) with regulated GS (de)adenylylation (Heeswijk2), suggesting that GlnK can substitute for GlnB. It has also been found that GlnK does interact with AT/AR in vitro (ninfa7), albeit relatively weakly, and binds $\alpha$-KG, ATP and ADP similarly to GlnB (xu, ninfa13). Importantly, glnK expression is induced during nitrogen limitation, whereas glnB is expressed constitutively (Heeswijk2). Taken together, these observations raise questions about the in vivo regulation of AT/AR and about a complementary role for GlnK in nitrogen assimilation to provide not only added redundancy but also a flexible and rapid response to the consequences of variations in external ammonium levels.

To establish the relevance of GlnK in nitrogen regulation during external nitrogen variations, we combined experiment with modelling. We measured the growth rate (Fig. 2) and the in vivo temporal response of relevant metabolites, proteins and post-translational modification (PTM) protein states during ammonium run-out, starvation and a subsequent ammonium shock (Fig. 3). Measurements of the wild type (WT) strain were complemented with experiments with two genetic knock-outs: $\Delta$glnB (no GlnB) and $\Delta$glnK (no GlnK), as shown in Fig. 5.

Our experiments show that growth rate and ammonium uptake are not impaired by the removal of either GlnB or GlnK (Fig. 2) suggesting that GlnB and GlnK conform two complementary branches of nitrogen signal transduction regulating AT/AR. We also find that GlnK plays a limited role in regulating GS expression and adenylylation during ammonium run-out and starvation. However, GlnK levels rise steeply during run-out, exceeding those of GlnB four-fold during starvation (Fig. 3). This suggests a functional role of GlnK, post-starvation, until normal ammonium conditions are restored.

To confirm this hypothesis, we developed a concise, mechanistic model of GS regulation (Fig. 1). We show that competition between GlnB and GlnK for UT/UR active sites and AT/AR allosteric sites, as well as functional sequestration of GlnK are necessary to explain the in vivo uridylylation and adenylylation dynamics in the WT and to predict those in the two mutants. Our model predicts the dynamic response of GS adenylylation levels in the mutants, indicating that GlnK can substitute for GlnB. Moreover, we find that GlnK is a potent regulator of AT/AR post-starvation, when GlnK levels are high relative to GlnB. Our results agree with the role of GlnK as a regulator of the ammonium transporter AmtB upon nitrogen upshift (merrick1, merrick3), and with the decreased viability of glnK mutants following extended ammonium starvation (blauwkamp). Hence our work suggests a role of GlnK as an anticipatory buffer molecule to help avoid glutamine shock: GlnK is expressed before and during starvation in order to temporarily alleviate excess ammonium import and assimilation when ammonium becomes available again.

Materials and methods

Strains and growth conditions

Experiments were performed with E. coli K12 (NCM3722) (soupene). NCM3722$\Delta$glnB in-frame deletions were obtained through P1 phage transduction (miller), using strain JW2537 with the glnK knock-out from the Keio collection as the donor strain (baba). NCM3722$\Delta$glnB was verified by genomic locus sequencing. Pre-cultures were grown in Gutnick minimal media, supplemented with 10 mM NH4Cl, 0.4% (w/v) glucose and Ho-Le trace elements. Main cultures were inoculated in same media but with limiting 3 mM NH4Cl, resulting in NH4Cl starvation at mid log phase. Ammonium depletion and OD600 were determined at defined time points during the course of experiments as described by (Schumacher). Following one generation time (40 min) of growth arrest, NH4Cl was added to obtain a concentration of 3 mM (upshift).

Metabolite and protein measurements

Glutamine and $\alpha$-ketoglutarate were quantified using liquid chromatography-mass spectrometry (LC-MS) and nuclear magnetic resonance (NMR), as described previously (Schumacher, komorowski). Total abundance and uridylylation level of GlnK were determined using multiple reaction monitoring mass-spectrometry (MRM-MS) (picotti) with absolute standard protein quantification (PSAQ) (picard, Schumacher). Briefly, isotopically labelled GlnK protein standard was purified following IPTG induced over-expression in E. coli harbouring JW0440 from the ASKA(-) collection (glnK) (kitagawa), grown in the presence of labelled L-Arginine (13C6, 15N4) and L-Lysine (13C6, 15N2). Ni-affinity purified GlnK standard purity was judged visually by SYPRO Ruby Protein Gel Stain (Invitrogen) stained SDS-PAGE to be $>$90% and GlnK standard concentration determined. Isotopic labelling efficiency was determined by MRM-MS to be 100%. Total GlnK abundances were determined by the ratio of MRM-MS signals of GlnK unlabelled/labelled signature peptides, excluding the uridylylation site peptide (Table S1 in the Supporting Material) and intracellular concentrations calculated (Schumacher). GlnK uridylylation levels were directly measured by the abundance of the uridylylated GlnK peptide GAEYS(UMP)VNFLPK compared with the non-uridylylated peptide GAEYSVNFLPK. The GlnK protein concentrations derived from non-uridylylated peptides correspond well with GlnK protein concentrations derived independently from the sum of uridylylated and non-uridylylated peptide GAEYSVNFLPK (Fig. S8 in the Supporting Material). A similar procedure is used for quantifying GlnB, GS and their PTMs (komorowski).

Description of the model

The schematic of the system regulating nitrogen assimilation in E. coli is shown in Fig. 1. We use three coupled non-linear ODEs to describe the dynamics of the monomeric concentrations of uridylylated GlnB, uridylylated GlnK and adenylylated GS:

[TABLE]

where the right hand side contains the balance of fluxes corresponding to (de)uridylylation and (de)adenylylation of the corresponding proteins. These fluxes contain several kinetic parameters and depend non-linearly on the concentrations of substrates, products and allosteric effectors. We now describe these terms in detail.

Equations for UT/UR:

The UT and UR activities in the bifunctional enzyme UT/UR are treated as independent, unidirectional reactions obeying Michaelis-Menten kinetics. Due to the specificity of UT/UR active sites (Zhang), GlnB and GlnK are modelled as competing substrates. Our model considers only the experimentally measured monomeric concentration of the uridylylated species, and disregard the multimeric nature of GlnB and GlnK, since (de)uridylylation of GlnB monomers is non-cooperative (atkinson). Glutamine (GLN) acts as an allosteric activator to UR and an inhibitor to UT activities (ninfa1). These can be encapsulated in the following rate equations (see Section S1 in the Supporting Material for derivations):

[TABLE]

Here $V_{*}$ and $K_{\text{m,*}}$ are the maximal rates and Michaelis constants of the (de)uridylylation reactions, and $K_{\text{GLN,1}}$ is the binding affinity of GLN to the allosteric site of UT/UR. The concentration of GlnB (and similarly GlnK) is obtained from the conservation relation:

[TABLE]

where the total concentration is a time-varying experimental measurement (Fig. 3).

Sequestration of GlnK:

We assume that GlnK-AmtB complex formation regulates the amount of GlnK accessible to other reactions, and is controlled by $\alpha$-KG following a Hill-type form (hwa):

[TABLE]

where $[\text{GlnK}_{\text{total}}]$ and $[\text{GlnK}]$ are the total and accessible concentrations of GlnK, respectively; K_{\text{\alpha-KG,1}} is the dissociation constant for $\alpha$-KG; $n_{\text{UT}}$ is the Hill-coefficient; and $x$ is the minimum concentration of accessible GlnK, which indicates the ratio of AmtB to GlnK. If [\alpha\text{-KG}]\ll K_{\text{\alpha-KG,1}}, then $[\text{GlnK}]=x[\text{GlnK}_{\text{total}}]$, and the maximal amount of GlnK is sequestered; if [\alpha\text{-KG}]\gg K_{\text{\alpha-KG,1}}, then $[\text{GlnK}]\approx x[\text{GlnK}_{\text{total}}]$, so no GlnK is sequestered.

Equations for AT/AR:

We split the activity of AT/AR into independent AT and AR reactions, both obeying Hill-type kinetics. We consider that GlnB-U/GlnK-U and GlnB/GlnK are pairwise competing for allosteric sites, while glutamine can bind independently to a third allosteric site (jaggi, ninfa9, ninfa13). The binding of GlnB and $\alpha$-KG are synergistic, and that of GlnB-U and $\alpha$-KG are antagonistic (ninfa12) leading to the following activation parameters:

[TABLE]

Similar expressions are used for GlnK and GlnK-U. As above, we only describe GS-A monomers (without explicit GS-A complexes) (mutalik, Bruggeman, ninfa16). We neglect retroactivity from GlnB and GlnK on $\alpha$-KG, since the latter is two orders of magnitude more abundant, and from AT/AR on GlnB and GlnK, since GlnK level is expected to be significantly higher (ninfa9). These assumptions result in the following rate equations (see Section S2 in Supporting Material for a detailed derivation):

[TABLE]

where $V_{*}$ and $K_{m,*}$ are the maximal enzyme rate and Michaelis constants of (de)adenylylation, $K_{*}$ are the binding affinities of GlnB(-U), GlnK(-U) and GLN to the allosteric sites of AT/AR and $n_{\text{AT}}$ is the Hill-coefficient. A summary of these interactions is shown in Figure 4.

Simulations were performed in Matlab. An SBML format of the model is included in the Supporting Material where we refer the reader for further details.

Model parameters: fitting and sensitivity analysis

The parameters were taken from the literature whenever possible (6 out of 21) and the rest (15 out of 21) were fitted to the 36 time points of the WT data (Table 1).

To assess the robustness of the model and to evaluate the relevance of the parameters to be fitted, we performed a global sensitivity analysis using the eFAST algorithm (saltelli). This technique considers uncertainty in combinations of parameters and quantifies the deviation they introduce over the whole time series. See Section S3 in the the Supporting Material for details of the method and its results.

To fit the parameters of the model to the measured time series, we used the Squeeze-and-Breathe evolutionary Monte Carlo optimisation algorithm (beguerisse). This algorithm minimises a cost function $J$, which is the sum-of-squares deviation of the fractional PTM levels from the measurements, weighted by the relative errors. Mathematically,

[TABLE]

where $p$ denotes the vector of parameters; $M$ and $F$ denote measured and fitted, respectively; and the superscripts $mod$, $total$ and $error$ denote total, modified and error, respectively. The parameter values are presented in Table 1, rounded to the nearest decimal found in the confidence interval given by the distribution of optimised parameters obtained with the Squeeze-and-Breathe method (see Section S4 in the Supporting Material for details).

Results

In vivo measurements of E. coli temporal responses to nitrogen starvation and upshift

We measured concentrations of relevant metabolites ($\alpha$-KG and GLN) and active/inactive PTM states of the proteins GlnB, GlnK and GS in response to dynamic variations in external ammonium (see Materials and Methods). We use this novel experimental approach to measure PTM states (Schumacher, komorowski) in order to link metabolic status and regulatory mechanisms.

Three E. coli cultures (a WT strain and two isogenic strains with deleted glnK or glnB) were inoculated into minimal media with a (limiting) concentration of 3 mM ammonium (see Materials and Methods). The growing cultures depleted the external nitrogen levels during the run-out period eventually leading to the starvation period measured from the point when growth stopped and extending for 40 minutes (Fig. 2). The starvation period was followed by a sudden 3 mM ammonium upshift to investigate the shock response (0 min to 40 min). We sampled the cultures at specific external ammonium concentrations during run-out and at specific time points during starvation and upshift, and measured the metabolite and protein concentrations. The corresponding time courses are plotted on Fig. 3.

In vivo regulation of GS by GlnK

First, we investigated whether our experiments provide evidence for GlnK controlling the level of GS expression. We found no discernible difference in total GS level between $\Delta$glnK and WT (Fig. 3), indicating that GlnK is redundant to regulate GS expression in the WT. In addition, we found that GlnK and GS were two-fold over-expressed in $\Delta$glnB (Fig. 3). This observation provides an in vivo quantification of the reduced regulatory effect of GlnK on GS expression, which has been previously evidenced in genetic studies (ninfa7).

Next, we asked whether GlnK plays a relevant role as a regulator of GS adenylylation. The response of the $\Delta$glnB strain to external ammonium variations was consistent with GlnK substituting for GlnB in its role to control AT/AR (ninfa8). Although the growth of $\Delta$glnB was slower than the WT during ammonium run-out, it was comparable to $\Delta$glnK (Fig. 2), indicating that, relatively speaking, ammonium uptake was not limiting. We further found that GS was abruptly adenylylated upon nitrogen upshift in all strains (Figs.  5c,e,g), including the $\Delta$glnB mutant, despite higher GS levels in the latter (Fig. 3). This suggests that, after starvation, GlnK is a potent activator of GS adenylylation. In contrast, GS was not fully deadenylylated in $\Delta$glnB (Fig.  5e and Fig. S6a in the Supporting Material) under ammonium rich conditions (unlike in WT and in $\Delta$glnK, see Fig. 5c,e), showing that GlnK might be less efficient in activating the AR activity. Indeed, GlnK was almost fully uridylylated ($>$95%) during starvation (see Fig. S6c in the Supporting Material), consistent with the dominant GlnK species under starvation being the fully uridylylated GlnK trimer, GlnK3-U3, which is known to be a poor activator of deadenylylation in vitro (Heeswijk4).

Unlike in $\Delta$glnB, we observed no evidence of GlnK regulated GS adenylylation in the WT during ammonium run-out and starvation, the GS-A levels being similar between the WT and $\Delta$glnK (Fig. 5c,g). Despite this purported redundancy, the expression and uridylylation of GlnK are already induced when ammonium levels begin to decrease (Heeswijk2) (Fig. 3 and 5b). Since GlnB expression is constitutive, we find that during starvation the concentration of GlnK (12.5 $\mu$M) largely exceeds that of GlnB (3.3 $\mu$M) (Fig. 3) and similarly GlnK-U (11 $\mu$M) largely exceeds GlnB-U (2.9 $\mu$M) (Fig. 5a,b). Thus intracellular glutamine sensing by UT/UR operates before ammonium concentration becomes suboptimal for growth.

Taken together, these findings suggest a role for GlnK to prepare for a buffered return to subsequent normal nitrogen conditions after starvation, by regulating the inactivation of GS. Our measurements show that the $\Delta$glnK strain grows at a higher rate after ammonium upshift (Fig.  2). We hypothesise that this difference is due to additional ammonium being imported via the membrane-bound ammonium channel AmtB, which accumulates during nitrogen starvation and whose activity becomes unimpeded at low levels of GlnK.

To understand further the interaction between GlnB and GlnK during nitrogen run-out, starvation and upshift, we constructed a mathematical model motivated by our in vivo dynamic measurements of protein (total and PTM) and metabolite concentrations. To date, the use of directly determined PTM data in models of regulation of GlnB and GlnK has been lacking, yet such modification directly determines their functionality. We now present an inclusive modelling approach to link GlnB and GlnK activity via PTMs to metabolic status and control of nitrogen assimilation processes.

Mathematical modelling of the enzymatic network of nitrogen signalling and regulation

In view of the experimental data, we built an ordinary differential equation (ODE) model to describe the dynamics of PTM states, i.e., the (de)uridylylation of GlnB and GlnK and the (de)adenylylation of GS in response to temporal inputs (drivers). An ODE model was adopted due to the large copy numbers observed—the smallest being GlnB with an average of $\sim 10^{3}$ molecules per cell. We based our equations on the enzyme mechanisms reported in the literature with the aim of a parsimonious model, only including terms with an observable effect in our experiments. The resulting model remains mechanistic with a moderate number of parameters, facilitating direct comparison to experiments and parameter fitting.

The terms of the ODE model are derived from the biochemical details of competition for the active and allosteric sites, which is inherent to the uridylylation and adenylylation reactions, and from the formation of complexes between GlnK and AmtB. In addition, we use five measured input drives to isolate the responses pertaining to transcription and other parts of the E. coli metabolism. These include concentrations of metabolites GLN and $\alpha$-KG, which act as enzyme effectors and reflect the metabolic state of the cell, and the total concentrations of GlnB, GlnK and GS, which reflect the relevant transcriptional responses of the cell (Fig. 3). We further assume that key cofactors (e.g., ATP, NADPH) remain at homeostasis under our experimental conditions (doucette), so that their concentrations can be approximated to be constant and can thus be absorbed into the kinetic parameters of the equations.

The model parameters include: binding affinities ($K_{\text{d}}$) of various effectors to active and allosteric sites of the enzymes; Michaelis constants ($K_{\text{m}}$) and maximum enzyme rates ($V_{\text{max}}$); and parameters representing AmtB-GlnK complex formation (Table 1). Whenever possible we used parameters from the biochemical literature, and we fitted the rest to the WT experimental time-courses using the Squeeze-and-Breathe algorithm (beguerisse). To confirm that the parameters capture meaningfully the experimental observations, we also performed a global sensitivity analysis over the time series (saltelli). For details, see Materials and Methods and the Supporting Material.

To make predictions for the $\Delta$glnB and $\Delta$glnK mutants, we set the corresponding reaction rates to zero, fixing all other parameters at values fitted to the WT time-course data. We then generate predicted time series for the PTM protein states under the experimentally measured drives. These predictions are compared to the independently obtained measurements.

Modelling the mechanism of GlnB and GlnK uridylylation

The bifunctional enzyme UT/UR has two independent active sites catalysing uridylylation (UT) and deuridylylation (UR), which are highly specific for GlnB and GlnB-U, respectively (Zhang); it also has an allosteric site for glutamine (Fig. 4a). Due to the similarity between GlnB and GlnK, we assume similar specificity of UT/UR active sites to GlnK/GlnK-U, yet with different binding affinities. Hence we model GlnB and GlnK uridylylation as two monocycles, GlnB-UT/UR and GlnK-UT/UR, which operate in parallel but not independently, since the total amount of enzyme is distributed between them at any given time. We further consider that the binding of glutamine inhibits the uridylylation of GlnB/GlnK and activates the deuridylylation of GlnB-U/GlnK-U (ninfa1) (Eq. 2, Materials and Methods).

As shown in Fig. 5a, our model captures accurately the excitatory and inhibitory effects in the measured GlnB uridylylation under different conditions in the WT without the need for modelling ternary complexes explicitly as in the more complex model in Ref. (straube) (Fig. S5 in the Supporting Material). The parameters fitted to the WT data show that the affinities of GlnK and GlnK-U to AT/AR ($K_{\text{m}}=9.9\mu$M and 59.8$\mu$M, respectively) are significantly lower than those of GlnB and GlnB-U ($K_{\text{m}}=3\mu$M and 2$\mu$M, respectively), which are known experimentally (ninfa3). Hence GlnK only affects nitrogen assimilation when its concentration is several times higher than GlnB, in agreement with previous in vitro findings (ninfa7).

Although this model allows us to fit well the overall WT time measurements, it cannot capture the GlnK-U level in response to the ammonium upshift after starvation (Fig. 5a). Moreover, the model over-predicts the level of GlnK-U in the $\Delta$glnB knock-out (Fig. 5b). In the next section we discuss an additional mechanism, which is necessary to make accurate predictions of the observed GlnK-U levels.

Sequestration of GlnK is necessary to predict uridylylation levels

To accurately describe GlnK uridylylation dynamics, we take into account the interaction of GlnK with AmtB. The membrane-bound channel protein AmtB facilitates ammonium active import to the cell. AmtB channels are blocked when AmtB forms a complex with trimers of deuridylylated GlnK above 20-50 $\mu$M external ammonium concentrations (merrick1). Hence AmtB-facilitated transport is negligible in our experiments before complete run-out. However, the formation of complexes between AmtB and GlnK does affect the dynamical responses, since it provides a mechanism for the sequestration of GlnK and GlnK-U at the membrane (merrick4). Note that while fully uridylylated trimers cannot bind AmtB (merrick3), our measurements show that the uridylylation level of GlnK is around 30% after starvation (see Fig. S6c in the Supporting Material) suggesting that the GlnK trimers will be mostly mono-uridylylated, i.e. one uridylylated GlnK per GlnK trimer, on average. In our model, we assume that both fully deuridylylated GlnK trimers and partially uridylylated GlnK trimers may be functionally sequestered to the membrane. However, unlike the fully deuridylyated GlnK trimers, partially uridylylated GlnK trimers do not block AmtB channel activity.

We model AmtB-GlnK complex formation by allowing a proportion of the total GlnK to be sequestered with a Hill-function dependence (Eq. 4, Materials and Methods) on the concentration of $\alpha$-KG (merrick1). Note that, since high levels of $\alpha$-KG level promote the dissociation of the AmtB-GlnK complex, sequestration only occurs during external ammonium abundance, i.e., when the $\alpha$-KG levels are low (Fig. 3).

The model with sequestration fitted to the WT data is able to capture the GlnK-U levels after nitrogen upshift (Fig. 5a). To justify the necessity of the two additional parameters ($K_{\alpha\text{-KG},1}$, $x$) involved in GlnK sequestration, we use the fitted model to predict the uridylylation states of GlnK and GlnB in the $\Delta$glnK and $\Delta$glnB knock-outs, and compare them to the measured time-courses (Fig. 5b,c). This approach has been used before to develop systems biology models (link, doucette2). Our predictions agree with the measurements: without sequestration, the model overpredicts the levels of GlnK-U before starvation (time before -40 min) and after nitrogen upshift (time after 0 min) in the $\Delta$glnB mutant (Fig. 5d). To quantify this, we compute the error of both models (with and without sequestration) against the data (7), and rank both models according to two information theoretic criteria (AICc, BIC), which penalise the additional parameter complexity against any improvement in the error of the model. The extended model with sequestration is clearly selected by both information criteria (see Table S2 in the Supporting Material) lending further basis for the need for a GlnK sequestration term during nitrogen excess.

The higher growth rate observed in the $\Delta$glnK strain after starvation (Fig. 2) is also consistent with GlnK sequestration. Since amtB is still present in this strain, AmtB is expected to accumulate during nitrogen starvation, which provides unrestricted ammonium import and higher growth rate in the short term.

GlnK can regulate AT/AR in the absence of GlnB

Our results so far indicate that GlnK not only competes with GlnB as a substrate or UT/UR, but also provides a link between active membrane transport and carbon metabolism through the interaction with $\alpha$-KG. To further test this mechanism, we asked: is GlnK sufficient on its own to regulate the adenylylation enzyme AT/AR in the absence of GlnB? If so, one would expect that the predicted GlnK-U levels in the $\Delta$glnB strain should account correctly for GS-A levels and, hence, for the activity of GS in nitrogen assimilation.

To test this hypothesis, we extended our model to include GS adenylylation. Similarly to UT/UR, the two antagonistic activities of AT/AR are modelled as independent uni-directional reactions, justified by structural studies showing the two active sites to be well separated (jaggi), making interaction between the AT and AR domains less likely (Fig. 4b).

We model different facets of AT/AR regulation as found in the literature (Fig. 4b). Firstly, AT and AR activities are inhibited by the products, GS-A and GS, respectively (ninfa9). Secondly, in the absence of allosteric effectors AT/AR has no observable activity (ninfa9). From this state, the rates of AT and AR actitivites are regulated by five allosteric effectors, which no not affect the binding rates of substrates (ninfa13). In particular, GLN, GlnB and GlnB-U may bind to three distinct sites (jaggi, ninfa9, ninfa13). As above, we assume that GlnK(-U) has the same specificity as GlnB(-U) but different affinity, so that GlnB and GlnK compete for one allosteric site, and GlnB-U and GlnK-U compete for the other. Finally, to account for intramolecular signalling (jaggi), we allow the activation parameters of GlnB, GlnB-U, GlnK and GlnK-U to depend on $\alpha$-KG (Eqs. 5, Materials and Methods). We encapsulate these terms into the reaction schemes (Eq. 6), by generalising the classical allosteric activation model (cornish) for two competing effectors (see Section S2 in the Supporting Material for details).

The additional equation for GS-A was fitted to the WT data and captures well the observed variability (Fig. 5a). To avoid overfitting, we keep fixed all parameters related to UT/UR and sequestration, and fix the insensitive activation parameters $\overline{K_{\text{GlnB}}}$ and $\overline{K_{\text{GlnB-U}}}$ to the values reported by (ninfa9) (Table 1, Materials and Methods and Section S3 in the Supporting Material). In contrast, the activation parameters $\overline{K_{\text{GlnK}}}$ and $\overline{K_{\text{GlnK-U}}}$ are relatively insensitive in the WT, yet essential to predict GS-A levels in the $\Delta$glnB strain. The fitted activation constant $K_{\text{GLN,2}}$ for AT/AR is in the range reported in Ref. (ninfa9) and the Hill coefficient $n_{\text{AT}}$ is found to be as reported in Ref. (mutalik). Note also that the activation constants of GlnK and GlnK-U for AT/AR are much higher than those of GlnB and GlnB-U, consistent with their observed lower affinity (ninfa7).

To assess the model fitted to the WT data, we predict the level of GS-A for the two knock-outs in response to the measured drives. The predicted responses agree with the measured GS-A levels in both strains (Fig. 5b,c). The GlnB-deficient strain ($\Delta$glnB) is able to regulate GS activity in response to external changes in ammonium. The prediction for the $\Delta$glnB strain underestimates the GS-A level before starvation—we consider the possible source for this discrepancy in the Discussion. Given that GS regulation depends on the fine balance between GlnB/GlnB-U and GlnK/GlnK-U together with the level of GLN and $\alpha$-KG, our model demonstrates that GlnK can regulate GS under varying ammonium conditions.

Discussion

We have investigated the role of GlnK in the dynamic regulation of ammonium assimilation from a systems perspective, and have argued that paradoxical observations about GlnK stem from the subtle interplay of regulatory mechanisms on two levels: transcriptional and post-translational. To decouple both levels of regulation, we collected concurrent in vivo time-courses of metabolite, total protein and PTM protein concentrations in response to time-varying external ammonium levels. We analysed the data in conjunction with a mechanistic model derived from biochemical principles. The combined analysis allowed us to fit the WT dynamics of PTMs, and to used the fitted model to predict the dynamic responses of experimental data from genetic knock-outs, shedding light on experimental ambiguities about the regulatory role of GlnK.

We found that competition between GlnB and GlnK for the bifunctional enzyme UT/UR together with GlnK sequestration by AmtB describe the uridylylation responses in the WT and can predict the response of two mutants under dynamic ammonium perturbations. AmtB-GlnK complex formation not only blocks ammonium import (merrick1), but also rapidly sequesters GlnK, whereas the less abundant GlnK limits ammonium assimilation by facilitating more rapid deuridylylation of GlnK-U (Fig. 5) and thus higher GS-A levels. Thus GlnK dampens the assimilatory response to ammonium upshift, perhaps to avoid excess ammonium uptake causing the depletion of the glutamate pool (wingreen, kustu). The fact that sequestration is regulated by $\alpha$-KG (a product of the carbon metabolism), suggests that nitrogen shock responses could also help to avoid untoward impacts on carbon status.

Our experiments showed that during ammonium rich conditions, GlnK is present at low level and does not influence GS production and adenylylation, which is instead regulated by the more potent GlnB. However, due to the induction of glnK during ammonium run-out, GlnK outnumbers GlnB four-fold during starvation (Fig. 3). To examine the capability of GlnK for AT/AR regulation after ammonium upshift, we developed a mathematical model of GS adenylylation based on allosteric competition between GlnB and GlnK for AT/AR. The model predicts accurately the abrupt increase in GS-A levels after ammonium upshift in both the $\Delta$glnB and $\Delta$glnK strains (Fig. 5) showing that both GlnB and GlnK are independently sufficient for effective AT/AR regulation. Our model slightly underestimates the GS-A levels before starvation in the $\Delta$glnB strain. Since the model was fitted to the WT, where both GlnB and GlnK are present, this mismatch may indicate unaccounted-for interactions which are disrupted when deleting the glnB gene, e.g., GlnB and GlnK forming heterotrimers (Forchhammer, Heeswijk4, Jiang:2009).

To check the relevance of GlnK regulation of AT/AR, we used our model to simulate two ‘in silico’ cellular scenarios. Firstly, we considered our model for a ‘computational strain’ where GlnK does not interact with AT/AR (by setting $\overline{K_{\text{GlnK}}}\gg[\text{GlnK}]$ and $\overline{K_{\text{GlnK-U}}}\gg[\text{GlnK-U}]$, such that the terms involving GlnK and GlnK-U drop out in Eq. 6), and hence AT/AR is regulated only by GlnB. In this scenario where the regulatory role of GlnK is removed, the predicted GS-A levels in response to ammonium upshift after starvation are reduced (Fig. 6a), thus leading to an increase in GS activity and larger ammonium assimilation. This is consistent with the proposed role of GlnK as a buffer to ammonium shock, through the adenylylation (inactivation) of GS. Noting that the level of GS-A does not change during run-out and starvation, we also used our model to predict the response of another ‘computational strain’ (Fig. 6) where the expression of glnK is constitutive and unregulated (by fixing $[\text{GlnK}_{\text{total}}]$ as a constant). In this scenario, if GlnK abundance is kept fixed at the late-starved level (26 $\mu$M), GS-A levels remain elevated during run-out and starvation. This simulation supports the notion that GlnK must be kept down-regulated during nitrogen abundance, since early glnK expression would lead to reduced ammonium assimilation due to increased GS-A levels. During ammonium run-out and starvation, amtB (and conjointly glnK, being on the same operon) is expressed so as to increase ammonium uptake via AmtB channels. The resulting higher levels of GlnK induced by starvation serve also to control potential overshoots in ammonium levels following starvation not only by blocking AmtB channels, but also by directly reducing GS activity through adenylylation.

Beyond nitrogen regulation, our work exemplifies the need for in vivo experimental time series data in developing predictive models. Time-course inputs together with classic genetic knock-outs provide tools to gain insight into the dynamic aspects of regulation in signalling pathways. Although large-scale biochemical models can be powerful, their applicability is many times hampered due to the large number of unknown parameters they contain, which need to be fitted to scarce data. Here we derived a small set of mathematical modules commensurate with the observed dynamics, and justified their relevance by a global parameter sensitivity analysis. We used a recently developed parameter fitting algorithm (beguerisse), which is especially appropriate for time-courses, to fit the WT data, and we then used the fitted model to predict out-of-sample dynamical observations from the genetic knock-outs.

Although we have focussed here on E. coli, with the dual GlnB/GlnK proteins, we expect nitrogen assimilation and transport to be more strongly coupled in bacteria containing only one GlnB-like protein (e.g. cyanobacteria, azotobacter), enabling the tight coordination of carbon and nitrogen metabolisms. Similarly, patterns of sequestration may be a biophysical regulatory mechanism operating widely in other contexts, and understanding its spatiotemporal dynamics in the cell would help link gene and metabolic regulation responses to biophysical feedback to environmental inputs (alley2008biophysical).

Author contributions

AG, MBa, MBu and JS designed the study. AG performed the analysis under the supervision of MBa and MBu. JS designed and coordinated experiments performed by VB, JGB, FH, and MBe. AG, MBa, JS and MBu wrote the manuscript, with inputs from other authors.

Acknowledgements

We thank Mariano Beguerisse-Diaz for his intellectual contributions about parameter fitting. This work was funded by a BBSRC LoLa grant (BB/G020434/1). AG acknowledges funding through a PhD studentship under the BBSRC DTP at Imperial College London (BB/M011178/1). M. Barahona acknowledges funding from the EPSRC (EP/N014529/1 and EP/I032223/1).

Data availability

Data supporting this study is available at figshare with DOI: 10.6084/m9.figshare.4880003.

Conflict of interest

The authors declare that they have no conflict of interest.

Supporting Material:

GlnK facilitates the dynamic regulation of bacterial nitrogen assimilation

S1 Kinetic mechanism for an enzyme with multiple alternative substrates

Here we present the detailed mathematical formulation of the kinetics of an enzyme with two alternative substrates as depicted on Fig. 4a and described by the following reactions:

[TABLE]

where (S1a) and (S1b) represent uridylylation, whereas (S1c) and (S1d) represent deuridylylation. Here $k_{1}$, $k_{2}$, $k_{3}$ and $k_{4}$ are the catalytic rate constants, which are functions of the glutamine concentration. Furthermore, ${\operatorname{E_{1}}}$ and ${\operatorname{E_{2}}}$ denote the concentration of the enzyme UT/UR with UT and UR activities respectively. To derive the rate equations corresponding to reactions (S1), we assume that complex formation between the substrates and the enzyme occurs much faster than product formation, which is the standard assumption in Michaelis-Menten kinetics (cornish). Furthermore, since the enzyme active sites are shared between the two substrates, we may treat either substrate as the competitive inhibitor of the other. Using techniques described in section 6.2.1 of (cornish) we obtain the following system of four equations:

[TABLE]

where $K_{m,*}=(d_{i}+k_{i})/a_{i}$, $i\in\{1,2,3,4\}$ are the usual Michaelis constants, and $V_{\text{UT,*}}=k_{p}[E_{1}]$, $p\in\{1,2\}$ and $V_{UR,*}=k_{q}[E_{2}]$, $q\in\{3,4\}$ are the maximum enzyme velocities.

In addition to the competing substrates, glutamine (GLN) may also bind to the allosteric side of UT/UR, independently of the substrate binding at a rate $K_{\text{GLN}}=d_{5}/a_{5}$, where $a_{5}$ and $d_{5}$ are the association and dissociation constants of GLN. When bound, GLN activates the UR activity and inhibits the UT activity (ninfa1). We may describe this using the noncompetitive activation and inhibition model (see Section 6.2.3 and Table 6.1 of (cornish)), depicted in Fig. S1), which result in an additional term modifying the $V_{\text{m,*}}$’s of eqs. (S2):

[TABLE]

S2 Kinetic mechanism for an enzyme with competing allosteric effectors

The mathematical description of simple inhibition and activation can be found in standard textbooks on reaction (cornish). However, enzymes with multiple competing allosteric effectors such as AT/AR are rare, therefore we need to give a derivation of such a system. Fig. 4b illustrates graphically an enzyme which can bind two allosteric effectors (inhibitors or activators) and a substrate molecule simultaneously, but their presence do not influence each others’ binding. As in Michaelis-Menten kinetics, we assume that the reactions involving complex formation between the enzyme, the substrates and the effectors are in equilibrium independently and at a much shorter timescale than that at which the reaction takes place.

We begin with the case of multiple competing allosteric activators, which can be depicted by the reaction diagram on the left of Fig. S2. Let $E$ denote a general enzyme catalysing the conversion of a substrate $S$ to a product $P$. Assume that $E$ is in a neutral state unless one of the activators $A$ or $B$ (having different affinities) bind to the common allosteric site. When either effector is bound, product formation may take place, however the reaction rate depends on whether $A$ or $B$ is bound. To express the reaction rate mathematically as a function of effector concentrations note that we can define five different complexes whose concentrations are denoted as: $c_{1}=[SE]$, $c_{2}=[SEA]$, $c_{3}=[EA]$, $c_{4}=[EB]$, $c_{5}=[SEB]$. Also, let $s=[S]$, $a=[A]$, $b=[B]$. Then by the law of mass action it follows that:

[TABLE]

where $e_{0}=e+c_{1}+c_{2}+c_{3}+c_{4}+c_{5}$ is the total concentration of the enzyme and $K_{i}=d_{i}/a_{i}$. This is a linear system of equations with eight variables. There are seven equations, but two are linear combinations of the other five (system has rank 5) so we can express the concentration of the complexes as functions of three, which we choose to be $s$, $a$ and $b$. Substituting (S4d), (S4e), (S4f) and (S4g) into (S4a) we obtain:

[TABLE]

Rearranging, it follows that:

[TABLE]

Hence the catalytic rate of the enzyme is:

[TABLE]

Hence comparing with the case of a single activator (see Cornish-Bowden 2013, pg. 152-157) we obtain that there the effect of the two activators are almost additive except for an additional mixed term in the denominator.

We proceed similarly to obtain a relationship for the case of multiple competing allosteric inhibitors (see right of Fig. S2). Using rapid equilibrium assumption we obtain a set of equations identical to system (S4). Here we need to solve for complex $c_{1}$, so we write:

[TABLE]

from which it follows that the enzyme velocity is:

[TABLE]

So similarly to the case of single inhibitor we obtain that there is a reduction of the maximum velocity of the reaction leaving the $K_{1}$ unchanged.

S3 Global sensitivity analysis of the model

To assess the significance of the model parameters we conducted a sensitivity analysis. The objective of the latter is to quantify the change in the system output to a small perturbation in one or a combination of parameters, while the others are kept constant. Since the outputs are time series, the distance between two time series if of interest, therefore conventional (local) sensitivity analysis techniques, which are only concerned with variations in the steady state are not adequate. Instead, here we use the eFAST method (saltelli) to assess parameter sensitivities, which is a global sensitivity analysis technique, and is available in the Systems Biology Toolbox 2 for Matlab (schmidt). In essence, global sensitivity analyses are concerned with sampling the parameter space - usually randomly, according to some optimal strategy - in the vicinity of the fitted parameters and estimating the resulting variation in the output. The latter is usually measured by the pointwise difference between the nominal and perturbed output time series integrated over the time range of interest. The eFAST method in particular samples the parameter space along such a trajectory that allows the model output to be expressed as a Fourier series. Then the variance in the output can be decomposed into a sum of terms involving Fourier coefficients, which can be estimated using Monte Carlo techniques. In our analysis we used a 10% perturbation from the nominal parameters and $10^{5}$ samples.

Fig. S3 shows the relative sensitivities of the model parameters. Since GlnB/GlnB-U and GlnK/GlnK-U are important effectors of the AT/AR enzyme, GS-A levels will depend on the parameters that are specific to uridylylation. However, GS-A levels do not feed back into the model since we used glutamine concentration as a driver. Hence the parameters specific to adenylylation do not affect GlnB-U and GlnK-U states. Hence we grouped the parameters according to whether they are involved in uridylylation or adenylylation. As expected, the most sensitive parameters are the maximum enzyme velocities, Michaelis constants and the Hill-coefficient of the AT/AR reactions ($n_{\text{AT}}$), which has a large effect on setting the basal GS-A levels. Another parameter of high sensitivity is the percentage sequestration, $x$, which confirms the importance of modelling AmtB-GlnK complex formation. The parameters of lowest relative sensitivity are the activation constants $K_{\text{GlnB}}$, $K_{\text{GlnB-U}}$, hence we fixed these from the literature (see Table I in the paper). Although $K_{\text{GlnK-U}}$ was also found to be relatively insensitive, which is commensurate with the wide sample distribution in Fig. S4, we could not find parameters in the literature measured under similar conditions to our experiments. We leave these at their fitted values, because they are highly influential for our predictions.

S4 Evolutionary Monte Carlo optimisation algorithm for parameter fitting

The parameters were fitted using the ’Squeeze and Breathe’ evolutionary optimisation method (beguerisse). Let $\mathbf{X}(t)=\{x_{1}(t),\dots,x_{d}(t)\}$ denote the state of the system with $d$ variables at time $t$, which in our case are the concentrations of GlnB, GlnK and GS. The evolution of these variables is described by a set of ODEs, $\dot{\mathbf{X}}=f(\mathbf{X},t;\mathbf{\theta})$, where $\theta=\{\theta_{1},\dots\theta_{n}\}$ is a set of $n$ parameters. In essence, the objective of the parameter optimisation is to find parameters $\theta$ such that the distance between the solution $\mathbf{X}$ and an experimental dataset $\mathcal{D}=\{\widetilde{\mathbf{X}}(t_{i})|i=1,\dots,m\}$ of $m$ observations is minimised. To account for measurement error we define the following cost function, weighted by the standard error of the measurements:

[TABLE]

where $\|.\|$ is the Euclidean norm. Since the cost function depends on all $n$ parameters, its value lies in an $n$ dimensional space. The optimal parameters will be the coordinates corresponding the global minimum of $E_{\mathcal{D}}(\theta)$ over all possible parameters. The cost function $E_{\mathcal{D}}(\theta)$ defines a very rough landscape and as a result optimisation methods may get stuck in a local minimum. In fact, an algorithm which guarantees to find the best parameters does not exist. Therefore the objective of the parameter optimisation is to explore a large portion of this space to get as close to the global minimum as possible. The Squeeze and Breathe algorithm achieves this by first running local optimisation around random samples in the parameter space. These are then ranked according to optimality and culled keeping only the best few. The culled set is then used to obtain a posterior sampling distribution. The process is repeated until the difference between subsequent posterior distributions is small.

The fitted parameters are shown in Table I in the main paper and the histograms from the parameter sampling are shown in Fig. S4. The histograms show the number of times the parameter fitting algorithm converged to the particular parameter values, whereas the red asterisk show the parameter value with the minimum cost function. A narrow distribution around a fitted parameter reflects a well defined minimum in the explored parameter range, indicating that the fitted model has higher sensitivity to the variations of these parameters. On the other hand a wide distribution shows that the cost function landscape is shallow or contains many local minima. This indicates that the model might be less sensitive to the corresponding parameters. Fig. S4 show that most parameter have a narrow spread. An exception is $K_{\text{GlnK-U}}$, which show high variation, indicating low sensitivity of the model to this parameter (confirmed by sensitivity analysis). In the main paper we argue, that the low sensitivity of $K_{\text{GlnK-U}}$ in the WT model has biological origins, since the GlnK protein is less potent compared to GlnB. However, we do not $K_{\text{GlnK-U}}$ to a value found in the literature, since this parameter will be essential for accurate prediction of the GS adenylylation levels in the mutant strain containing no GlnB. The parameters to be fitted were chosen also based on the information provided by a global sensitivity analysis of the model presented in Section S3.

S5 Model selection: uridylylation reactions with and without sequestration

We compared the uridylylation models with and without GlnK sequestration (Eq. 4 in the main text) using different model selection criteria. All criteria unanimously selected the model with sequestration.

Our model is fitted to the WT time-series (training dataset). The fitted model is then used to predict the dynamical responses for the two mutants (test datasets) without using the measurements. The model with sequestration provides an improved fit both for the training dataset and an improved prediction of the test datasets, yet it contains 6 parameters (instead of 4 parameters for the model without sequestration).

To make the comparison precise, we obtain the statistical significance of the models with and without sequestration based on two information criteria:

• •

Akaike Information Criterion corrected for small datasets:

[TABLE]

• •

Bayesian Information Criterion:

[TABLE]

Here $n$ is the number of data points (24 using the WT and the $\Delta$glnB datasets); $k$ is the number of fitted parameters (4 without sequestration and 6 with sequestration); and RSS is the residual sum-of-squares deviation of the model from the data (Eq. 7 ).

Both criteria compare models based on their goodness-of-fit and penalise the number of parameters. The AICc criterion is theoretically more appropriate (since it compares the model to the true model), but has a bias toward models with higher complexity. As a further confirmation, we also computed the BIC criterion, which puts a higher penalty on the number of model parameters, and would thus favour the lower complexity model.

Table S1 shows that the model with sequestration is selected according to both the AICc and BIC criteria, as shown by lower values of both criteria for the model with sequestration.

S6 Comparison to the Straube model of GlnB uridylylation

The enzyme UT/UR has two active sites (Zhang), catalysing uridylylation (UT) of GlnB/GlnK and deuridylylation (UR) of GlnB-U and GlnK-U. Having two active sites both with two distinct substrates means there are four different ways in which ternary complexes can form. To test whether these contribute significantly to product formation in the UT and UR reactions we compared the Straube model (straube), which describes GlnB uridylylation and includes the ternary complex between UT/UR, GlnB and GlnB-U, and another model, where the contribution of the ternary complex is ignored. The latter is equivalent to the classical Goldbeter-Koshland model (gk) that assumes Michaelis-Menten kinetics with allosteric inhibition and activation by glutamine for the UR and UT reactions respectively. For both models we used literature parameters from the in vitro reconstituted GlnB-UT/UR system (ninfa3) except for the $V_{max}$, which had to be fitted. To make the best comparison we used the fact that $V_{\text{max}}=ke$, where $k$ is the catalytic rate and $e$ is the concentration of the enzyme, fixed the ratio between $V_{UT}/V_{UR}=k_{UT}/k_{UR}$ using literature values of $k_{UT}$ and $k_{UR}$ (ninfa3) and fitted only one parameter, corresponding to $e$. As Fig. S5 shows, both models produce indistinguishable results under most conditions except 30s after upshift where the Straube model, being a steady state model, performs worse due to strong transient dynamics. This suggests that the contribution of the ternary complexes are not likely to be significant.

Despite performing well on GlnB uridylylation using literature parameters, the Straube model could only produce an adequate fit to the GlnK data in the $\Delta$glnB strain when the basal UT activity (which is one of its parameters) was as high as 50%. We found this after an extensive parameter search using a non-linear least squares fitting procedure. This is much higher than the 1% basal UT activity reported in a previous in vitro study with GlnB and UT/UR (ninfa3). Furthermore, due to the number of complexes this model cannot be extended to treat GlnB and GlnK at the same time.

The high basal UT activity for GlnK in the Straube model highlights the asymmetry in the system with respect to glutamine. In other words, glutamine binding to UT/UR activates the UR activity and inhibits the UT activity by different relative amounts. We could account for this in our model by defining two different $K_{d}$ for the activation of UR and inhibition of UT reactions by glutamine.

S7 Additional Supplementary Figures and Tables