Transient Accumulation of NO2 - and N2O during Denitrification Explained by Assuming Cell Diversification by Stochastic Transcription of Denitrification Genes

Denitrifying bacteria accumulate NO2−, NO, and N2O, the amounts depending on transcriptional regulation of core denitrification genes in response to O2-limiting conditions. The genes include nar, nir, nor and nosZ, encoding NO3−-, NO2−-, NO- and N2O reductase, respectively. We previously constructed a dynamic model to simulate growth and respiration in batch cultures of Paracoccus denitrificans. The observed denitrification kinetics were adequately simulated by assuming a stochastic initiation of nir-transcription in each cell with an extremely low probability (0.5% h-1), leading to product- and substrate-induced transcription of nir and nor, respectively, via NO. Thus, the model predicted cell diversification: after O2 depletion, only a small fraction was able to grow by reducing NO2−. Here we have extended the model to simulate batch cultivation with NO3−, i.e., NO2−, NO, N2O, and N2 kinetics, measured in a novel experiment including frequent measurements of NO2−. Pa. denitrificans reduced practically all NO3− to NO2− before initiating gas production. The NO2− production is adequately simulated by assuming stochastic nar-transcription, as that for nirS, but with a higher probability (0.035 h-1) and initiating at a higher O2 concentration. Our model assumes that all cells express nosZ, thus predicting that a majority of cells have only N2O-reductase (A), while a minority (B) has NO2−-, NO- and N2O-reductase. Population B has a higher cell-specific respiration rate than A because the latter can only use N2O produced by B. Thus, the ratio BA is low immediately after O2 depletion, but increases throughout the anoxic phase because B grows faster than A. As a result, the model predicts initially low but gradually increasing N2O concentration throughout the anoxic phase, as observed. The modelled cell diversification neatly explains the observed denitrification kinetics and transient intermediate accumulations. The result has major implications for understanding the relationship between genotype and phenotype in denitrification research.


Introduction
The dissimilative reduction of nitrate (NO À 3 ) to nitrite (NO À 2 ), nitric oxide (NO), nitrous oxide (N 2 O), and finally to N 2 (denitrification) is an indispensable process in the nitrogen cycle, returning N to the atmosphere as N 2 . However, denitrification significantly leaks the gaseous intermediates NO and N 2 O, both with serious consequences for the environment. N 2 O catalyses depletion of the stratospheric ozone [1] and causes global warming, contributing~10% to the anthropogenic climate forcing [2]. Data suggests that since the 1950s, the atmospheric N 2 O has been increasing, and before being photolysed in the stratosphere, the gas persists for an average~120 years in the troposphere [3].~70% of global N 2 O emissions are tentatively attributed to microbial nitrification and denitrification in soils [4], where denitrification, generally, is considered a more dominant source [5].
To mitigate N 2 O emissions, we need to understand the physiology of denitrifiers To devise robust strategies for mitigating global N 2 O emissions, a good understanding of its primary source is imperative, i.e., genetics, physiology, and regulatory biology of denitrifiers. Any knowledge of the environmental controllers of N 2 O is incomplete without understanding the causal relationships of such controllers at the physiological level [6].
The biogeochemical models developed for understanding the ecosystem controls of denitrification and N 2 O emissions treat the denitrifying community of soils and sediments as a single homogenous unit with certain characteristic responses to O 2 and NO À 3 concentrations [6,7]. Natural denitrifying communities, however, are mixtures of organisms with widely different denitrification regulatory phenotypes [8]. The regulatory response of such mixtures is not necessarily equal to the 'sum of its components' because there will be interactions, not the least, via the intermediates NO and NO À 2 . Hence, it is probably a mission impossible to predict the regulatory responses of complex communities based on their phenotypic composition.
Nevertheless, investigations of the regulation in model organisms like Pa. denitrificans provide us with essential concepts, enhancing our ability to understand the regulatory responses of mixed communities and to generate meaningful hypotheses. Thus, future biogeochemical models of N 2 O and NO emissions are expected to have more explicit simulations of the regulatory networks involved, and a first attempt has recently been published [9].
Simulating the cell diversification in response to impending anoxia to analyse its implications for NO À 2 , N 2 , and N 2 O kinetics Dynamic modelling has been used to a limited extent to analyse various denitrification phenotypes; for example, to analyse NO À 3 and NO À 2 reduction and gas-kinetic data for individual strains [10] and mixtures of selected phenotypes [11]; to model the consequence of competition for electrons between denitrification reductases [12,13]; to investigate the control of O 2 on denitrification enzymes and inhibition of cytochrome c oxidase by NO in Agrobacterium tumefaciens [14]; and to examine the effect of copper availability on N 2 O reduction in Paracoccus denitrificans [15]. In our previous model [16], we simulated O 2 and N 2 kinetics from batch incubations of Pa. denitrificans [8,17] to test if a postulated cell diversification, driven by stochastic initiation of nirS, could explain the N 2 production kinetics in NO À 2 -supplemented media. The available data also contained NO À 3 -supplemented treatments but NO À 3 and NO À 2 were not monitored, and the experiment provided no information about the N 2 O kinetics, except that the concentrations were extremely low (below the detection limit of the thermal conductivity detector used). Recently, a neat dataset was generated from batch incubations supplemented with NO À 3 , with frequent measurements of NO À 2 and a more sensitive detection of N 2 O by an electron capture detector [18]. That encouraged us to extend our previous model and simulate the cell diversification during transition from oxic to anoxic conditions, targeting the regulation of Nar and cNor/NosZ (N 2 O emissions) in Pa. denitrificans.

Regulatory network of denitrification in Paracoccus denitrificans
Pa. denitrificans is a facultative anaerobe capable of reducing NO À 3 all the way to N 2 : In response to impending anoxic conditions, the organism sustains respiratory metabolism by producing the membrane-bound cytoplasmic nitrate reductase (Nar), cytochrome cd 1 nitrite reductase (NirS), cytochrome c dependent nitric oxide reductase (cNor), and nitrous oxide reductase (NosZ). Transcription of the genes encoding these reductases (narG, nirS, norBC, and nosZ, respectively) are regulated by the FNR-type proteins FnrP, NarR, and NNR. FnrP contains a 4Fe-4S cluster for sensing O 2 , and NNR harbours a NO-sensing haem; NarR, however, is poorly characterised and is most likely a NO À 2 -sensor [19][20][21]. All these sensors remain inactive during aerobic growth conditions [19].
Transcription of denitrification genes in Pa. denitrificans. FnrP and NarR facilitate a product-induced transcription of the nar genes, and NNR facilitates a product-induced transcription of the nirS genes (Fig 1, see P 1 and P 2 ): Low oxygen concentration ([O 2 ]) activates the self-regulating FnrP, which induces nar transcription in coaction with NarR. The self-regulating NarR was previously assumed to be activated by either NO À 3 or NO À 2 [21], but a recent proteomics study indicates that NO À 2 is the activator [19]. Thus once a cell starts producing traces of NO À 2 , nar expression becomes autocatalytic. Transcription of nirS is induced by NNR, which is apparently inactivated by O 2 [22,23], but under anoxic/micro-oxic conditions, NNR is activated by NO. Thus, once traces of NO are produced, the expression of nirS also becomes autocatalytic [19,20]. In contrast, nor transcription is substrate (NO) induced via NNR while nosZ is equally induced by NNR or FnrP [24]. High concentrations of NO may constrain nar transcription by inactivating FnrP [25] and, like O 2 , render NosZ dysfunctional by inactivating the Cu Z subunit of the reductase [26], but these observations are ignored in our model because Pa. denitrificans restricts [NO] to very low levels.

Entrapment of cells in anoxia: The underlying hypothesis and modelling
Denitrification proteome, once produced in response to an anoxic spell, is likely to linger within the cells under subsequent oxic conditions, ready to be used if anoxia recurs. But the proteome will be diluted by aerobic growth because the transcription of denitrification genes is inactivated under oxic conditions [20]. Hence, a population growing through many generations under fully oxic conditions is expected to undertake de novo synthesis of denitrification enzymes when confronted with anoxia. Batch cultivations of such aerobically raised Pa. denitrificans provided indirect evidence for a novel claim that, in response to anoxia, only a small fraction of the incubated population is able to produce denitrification proteome [8,17,27,28]. Our dynamic modelling of Regulatory network of denitrification in Pa. denitrificans. The network is driven by four core enzyme-complexes: Nar (transmembrane nitrate reductase encoded by the narG gene), NirS (cytochrome cd 1 nitrite reductase encoded by nirS), cNor (NO reductase encoded by norBC), and NosZ (N 2 O reductase encoded by nosZ). When anoxia is imminent, the low [O 2 ] is sensed by FnrP, which in some interplay with NarR induces nar transcription. NarR is activated by NO À 2 ; thus once a cell starts producing traces of NO À 2 , nar expression becomes autocatalytic (see P 1 ). Transcription of nirS is induced by NNR, which is activated under anoxic/micro-oxic conditions by NO; thus once traces of NO are produced, the expression of nirS also becomes autocatalytic (see P 2 ) [20]. The activated P 2 will also induce nor and nosZ transcription via NNR. The transcription of nosZ, however, can also be induced equally and independently by FnrP [24]. Micromolar concentrations of NO may inactivate both FnrP [25] and NosZ [26]. These observations, however, are ignored for our modelling because Pa. denitrificans restricts NO to nanomolar levels. Bergaust et al.'s [17] NO À 2 -supplemented incubations corroborated this, suggesting that a probabilistic function (specific probability = 0.005 h -1 ) resulting in the recruitment of 3.8-16.1% of all cells to denitrification is adequate to explain the measured N 2 kinetics [16].
Our model was based on the hypothesis that the entrapment of a large fraction in anoxia is due to a low probability of initiating nirS transcription, which in response to O 2 depletion is possibly mediated through a minute pool of intact NNR, crosstalk with other factors (such as FnrP), unspecific reduction of NO À 2 to NO by Nar, and/or through non-biologically formed traces of NO found in a NO À 2 -supplemented medium. Regardless of the exact mechanism(s), once nirS transcription is initiated, the positive feedback via NO/NNR (Fig 1, see P 2 ) would allow the product of a single transcript of nirS to induce a subsequent burst of nirS transcription. The activated positive feedback will also help induce nor and nosZ transcription via NNR, rapidly transforming a cell into a full-fledged denitrifier. We further hypothesised that recruitment to denitrification will only be possible as long as a minimum of O 2 is available because, since Pa. denitrificans is non-fermentative, the synthesis of first molecules of NirS will depend on energy from aerobic respiration.
The above hypothesis was modelled by segregating the culture into two pools (subpopulations): one for the cells without (N D− ) and the other with denitrification enzymes (N D+ ). Initially, all cells were N D− , growing by consuming O 2 . As [O 2 ] fell below a certain threshold, N D− recruited to N D+ with a constant probability (h -1 ), assumed to be that of the nirS transcriptional activation, and the recruitment halted as O 2 was completely exhausted, assuming lack of energy (ATP) for enzyme synthesis.

Underlying assumptions and aims of the present modelling
The present model is an extension of that developed in Hassan et al. [16]. Here we have divided the respiring culture into four pools (Fig 2A): 1. Z − : cells without Nar, NirS, and cNor 2. Z Na : cells with Nar 3. Z NaNi : cells with Nar, NirS, and cNor

Z Ni : cells with NirS and cNor
All these subpopulations are assumed to scavenge O 2 (if present) and produce NosZ in response to impending anoxia. The latter because the nosZ genes are readily induced by the O 2 -sensor FnrP [24].
The Z − pool (Fig 2A) contains the inoculum that grows by aerobic respiration. As [O 2 ] falls below a critical threshold [empirically determined, 18], the cells within Z − are assumed to start synthesising Nar with a certain probability and populate the Z Na pool. The aim here is to investigate whether, like for nirS, the initiation of nar transcription (by some combined activity of FnrP and NarR) can also be explained as a probabilistic phenomenon, quickly differentiating a cell into a full-fledge NO À 3 scavenger through product (NO À 2 ) induced transcription via NarR (Fig 1, see P 1 ). If so, we were interested to estimate what fraction of the cells is required to adequately simulate the measured data (NO À 2 production), aiming at scrutinising the general assumption that all cells in batch cultures produce Nar in response to impending anoxia.
Next, when [O 2 ] is further depleted to another critical threshold [18], the Z − and Z Na cells are assumed to initiate nirS transcription with a low per hour probability and, thereby, populate the Z Ni and Z NaNi pools, respectively. As explained above for our previous model, NirS + cNor production is assumed to be a) coordinated because the transcription of both nirS and nor is induced by NO via the NO-sensor NNR (Fig 1), and b) stochastic because the initial transcription of nirS (paving the way for the autocatalytic expression of NirS and substrate-induced nor transcription) happens in the absence of NO or at too low [NO] to be sensed by NNR.
Synthesis of denitrification enzymes requires energy, which all the subpopulations can obtain by respiration only. Hence, the initiation of the autocatalytic expression of nar and nirS The squares represent state variables, the circles the rate of change of the state variables, the edges (thicker arrows) depict flows into or out of the state variables, the shaded ovals auxiliary variables, and the arrows portray mutual dependencies between the variables. All feedback relationships among the three model sectors could not be shown; however, for illustration the feedback relationships of one sub-population (Z − ) are shown (dashed arrows). Within each square (state variable), t 0 refers to the initial value. (i.e., recruitment to Z Na and Z NaNi /Z Ni , respectively, Fig 2A) depends on the availability of the relevant terminal e --acceptor(s) above a critical concentration to sustain a minimum of respiration. For Z − , the only relevant e --acceptors are O 2 and the traces of N 2 O produced by Z Ni and Z NaNi . The same applies For Z Na , but in addition, this subpopulation can also obtain energy by reducing NO À 3 , if present. In our previous model [16], we assumed that recruitment to denitrification was sustained by energy from O 2 -respiration only; not NO À 3 because we simulated NO À 2 -supplemented treatments, and not by N 2 O because we naively assumed that the pool of this e --acceptor was insignificant (N 2 O concentrations were below the detection limit of the system used for those experiments). However, the present model assumes that the recruitment from Z − to Z Na and Z − to Z Ni is sustained by both O 2 -and N 2 O-reduction, and the recruitment from Z Na to Z NaNi is sustained by O 2 -, N 2 O-and NO À 3 -reduction, when above a critical minimum (ve À min ). The default value for ve À min was set to an arbitrary low value (= 0.44% of the maximum e --flow rate to O 2 ), and we have investigated the consequences of increasing, decreasing, and setting ve À min = 0. The expressions of nar and nirS + nor (recruitments to Z Na and Z NaNi /Z Ni , respectively, Fig  2A) are modelled as instantaneous discrete-events in each cell, thus ignoring the time-lag from the initiation of gene transcription till the cell is fully equipped with the reductase(s) in question. That is because the lag observed between the emergence of denitrification gene transcripts and the subsequent gas products suggests that the synthesis of denitrification enzymes takes less than half an hour [17,18], which is negligible for our purposes here.
The main purpose of the present modelling is to investigate if a full-fledged model, including all four functional denitrification reductases, could adequately simulate the observed kinetics and stoichiometry of denitrification products [18]. These cultures reduced all available NO À 3 to NO À 2 prior to the onset of gas production and accumulated traces of N 2 O throughout the anoxic phase, as illustrated in S1 Fig In particular, we were interested to investigate the NO À 2 kinetics, controlled by narand nirS transcription, and to test if the peculiar N 2 O kinetics (low, but increasing concentrations throughout the anoxic phase) could be explained by our modelled cell diversification.

An overview of the modelled experiment
Batch incubation. Qu [18] incubated Pa. denitrificans (DSM-413) at 20°C using 50 mL Sistrom's [29] medium in 120 mL gas-tight vials. Either succinate or butyrate (5 mM) was used as the main carbon source, enough to secure consumption of all available e --acceptors. After distribution of the medium, each vial was loaded with a magnetic stirring bar, sterilised through autoclaving, supplemented with 2 mM KNO 3 , and was tightly sealed. To remove O 2 and N 2 from the headspace, the headspace air was evacuated and replaced by helium (He) through several cycles of evacuation and He-filling (He-washing). Some vials were supplemented with oxygen to reach 7 vol.% O 2 in headspace (treatment designated 7% O 2 ). The remaining vials received no O 2 (designated 0% O 2 , although there were traces of O 2 present despite the He washing). For each treatment (i.e., C source and initial O 2 ), there were three replicates, and each vial was inoculated with 2.2×10 8 aerobically grown cells.
NO À 2 and gas measurement. Gases (CO 2 , O 2 , NO, N 2 O, and N 2 ) were monitored by frequent sampling of the headspace, using an improved version of the robotised incubation system [30]. In short, the system draws gas samples from the headspace (peristaltic pumping) via the septum pierced by a needle, filling three loops used for injecting samples into the two GC columns and the chemiluminescence NO analyser. The sample drawn is replaced by He (reversing the peristaltic pump), thus securing~1 atm pressure. The primary improvements of the new system are a more sensitive detection of N 2 O (by an electron capture detector), lower sampling volumes (~1 mL), and lower leaks of O 2 and N 2 through the sampling system (4 nmol O 2 and 12 nmol N 2 per sampling, which is~20% of that for the old system).
To extract samples for measuring NO À 2 without tampering the original vials, identical (parallel) vials were prepared for each treatment. Using sterile syringes, samples of 0.1 mL were regularly drawn from the liquid-phase of the parallel vials and immediately analysed for NO À 2 . Results for one of the treatments are shown in S1 Fig, illustrating the complete reduction of NO À 3 to NO À 2 prior to the onset of significant N-gas production. In previous experiments [17], N 2 O concentrations were below the detection limit of the system, but thanks to the new system, the N 2 O kinetics were monitored with a reasonable precision.

The model
The model is constructed in Vensim DSS 6.2 Double Precision (Ventana Systems, inc. http:// vensim.com/) using techniques from the field of system dynamics [31].
Cell diversification and growth. The respiring population is divided into four subpopulations, according to their reductases (Fig 2A): 1) Z − : cells without Nar, NirS, and cNor; 2) Z Na : cells with Nar; 3) Z NaNi : cells with Nar, NirS, and cNor; and 4) Z Ni : cells with NirS and cNor. All the subpopulations are assumed to equally respire O 2, if present, and express nosZ in response to oxygen depletion [24]. Z − contains the inoculum (= 2.2×10 8 cells) that grows by aerobic respiration. As O 2 is depleted, the Z − cells populate the other pools by producing Nar and/or NirS + cNor.
The recruitment from Z − to Z Na (R Na , Fig 2A) takes place first: is a conditional specific probability (h -1 ) for any Z − cell to initiate nar transcription (quickly transforming it into a NO À 3 scavenger through autocatalytic gene expression, see (h -1 ) where r Na (h -1 ) is a constant specific probability for a cell to initiate nar transcription once O 2 concentration in the aqueous-phase ([O 2 ] aq , mol L -1 ) falls below a critical concentration ([O 2 ] na ), empirically determined as the [O 2 ] aq (= 4.75×10 −5 mol L -1 ) at the outset of NO À 2 accumulation in the medium [18]. The second condition for a cell to produce first molecules of Nar is a minimum of e --flow to an e --acceptor (ve À min , mol ecell -1 h -1 ), assumed to generate minimum ATP required for protein synthesis. ve À O 2 and ve À N 2 O (mol ecell -1 h -1 ) are the cell-specific velocities of e --flow to O 2 and N 2 O, respectively. The latter is weighed down by 0.5 because mole ATP per mole etransferred to NO À x /NO x is lower for denitrification than for aerobic respiration [17,20]. For a Z − cell, ve À NO À 2 and ve À NO are not considered here, since such a cell is assumed to have no NirS and cNor. The fraction of the cells that successfully produces Nar (F Na ) is calculated based on the integral of the recruitment (Eq 1): where t Na is the time-window available for the recruitment. In theory, t Na is the time-period when ½O 2 aq < ½O 2 na AND ðve À Since the e --flow to N 2 O started after all NO À 3 had been reduced to NO À 2 (S1 Fig), the recruitment based on ve À N 2 O would be inconsequential for the simulated (and measured) NO À 2 kinetics. Therefore, to calculate the functional F Na actually responsible for producing NO À 2 , we ignored the N 2 O-sustained recruitment, thus considering t Na to be the time when ½O 2 aq < ½O 2 na AND ve À O 2 > ve À min . Next, the cells within Z Na and Z − are recruited to Z NaNi and Z Ni (R NaNi and R Ni , respectively, Fig 2A), as they are assumed to stochastically initiate nirS transcription, paving the way for NO/NNR mediated autocatalytic expression of nirS + nor (Fig 1). In principle, the rates of both these recruitments are modelled as that of the recruitment from Z − to Z Na (Eqs 1 and 2): a) Both trigger as O 2 falls below another critical concentration ([O 2 ] ni ), low enough to activate NNR to induce nirS transcription; [O 2 ] ni (= 1.16×10 −5 mol L -1 ) is empirically determined as the O 2 concentration at the outset of NO accumulation [18]. b) Both continue as long as a minimum of e --flow to the relevant terminal e --acceptor is possible, sustaining the respiratory metabolism to generate ATP for protein synthesis: where r Ni is a constant specific probability (h -1 ) for the initiation of nirS transcription. ve À NO À 3 and ve À N 2 O are multiplied with 0.5 for the same reasons as described for Eq 2. The recruitment from Z − to Z Ni (R Ni , Fig 2A) is modelled as a product of Z − and a conditional specific probability, r Ni (O 2 ,N 2 O), which is different from Eq 5 only in that ve À NO À 3 is omitted, since Z − do not possess Nar: The fraction that successfully produced NirS + cNor (F Ni ) is calculated based on the integral of R NaNi and R Ni : where t NaNi is the duration of the recruitment from Z Na to Z NaNi , i.e., when ½O 2 aq < (Eqs 4 and 5), F Na is the fraction recruited to the pool of Nar positive cells (Z Na , Eq 3), and t Ni is the duration of the recruitment from Z − to Z Ni , i.e., when ½O 2 aq < ½O 2 ni AND ðve À (Eqs 6 and 7). Each of the populations will grow depending on the rates of e --flow to the various e --acceptors they are able to use: (cells h -1 ) where Ye À X (cells mol -1 eto X = O 2 or NO À x /NO x ) is the growth yield determined under the actual experimental conditions, and ve À X (mol ecell -1 h -1 ) is the cell-specific velocity of e --flow to X (O 2 or NO À x /NO x ), which depends on the concentration of the e --acceptor (see Eqs 17, 20 and 28). For NO À 3 and NO À 2 , a restricted velocity (ve À NO À x res ) is used so that when electrons flow to O 2 , NO À 3 , and NO À 2 simultaneously, the total ve − per cell does not exceed the maximum electrons that the TCA cycle (ve À maxTCA ) can deliver per hour (see Eqs 21 and 22). O 2 kinetics. O 2 is initially present in the headspace (O 2 g , mol, initialised according to the experiment, see Table 1) but is transported to the liquid-phase (O 2 aq ) due to its consumption therein (Fig 2B). The transport rate (Tr O 2 ) is modelled according to Molstad et al. [30]: In addition, the model simulates the changes in O 2 g due to sampling. The robotised incubation system used monitors gas concentrations by sampling the headspace, where each sampling alters the concentrations in a predictable manner: a fraction of O 2 g is removed and replaced by He (dilution), but the sampling also results in a marginal leakage of O 2 through the tubing and membranes in the injection system. The net change in O 2 g (ΔO 2(S) ) as a result of each sampling is calculated as: (mol h -1 ) where Z − , Z Na , Z NaNi , and Z Ni (cells) are all the sub-populations present (described above); thus, we assume that all cells have the same potential to consume O 2 . is modelled as a Michaelis-Menten function of oxygen concentration: (mol ecell -1 h -1 ) where ve À maxO 2 (mol ecell -1 h -1 ) is the maximum velocity of e --flow to O 2 per cell (determined under the actual experimental conditions), [O 2 ] aq (mol L -1 ) is the O 2 concentration in the liquid-phase, and K mO 2 (mol L -1 ) is the half-saturation constant for O 2 reduction.
Denitrification kinetics. The NO À 3 and NO À 2 pools (mol, Fig 2C) are initialised according to the experiment (Table 1; NO À 2 = 0). The kinetics of these nitrogen oxyanions (NO À x ) are modelled as: (mol h -1 ) where Rr NO À x (mol h -1 ) is the reduction rate, Z Na + Z NaNi (cells) is the total number of cells with Nar, Z NaNi + Z Ni (cells) is the total NirS active population, and v NO À x (mol cell -1 h -1 ) is the cell-specific velocity of NO À x consumption, obtained by the velocity of e --flow to NO À x 1 mol NO À 3 2 mol e À and 1 mol NO À 2 1 mol e À . The latter is modelled as a Michaelis-Menten function of NO À x concentration: (mol ecell -1 h -1 ) where ve À maxNO À x (mol ecell -1 h -1 ) is the maximum velocity of e --flow to NO À x per cell (determined under the actual experimental conditions), ½NO À x aq (mol L -1 ) is the NO À x concentration in the aqueous-phase, and K mNO À x (mol L -1 ) is the half-saturation constant for NO À x reduction. The velocity of NO À 3 and NO À 2 consumption had to be restricted (ve À NO À x res ) to ensure that when electrons flow to O 2 , NO À 3 , and NO À 2 simultaneously, the total ve − per cell does not exceed an estimated maximum delivery of electrons from the TCA cycle (ve À maxTCA ). In competition for electrons, O 2 is prioritised [20], followed by NO À 3 and NO À 2 , respectively [18]: (mol ecell -1 h -1 ) where ve À NO À 3 res is the realised e --flow to NO À 3 , limited either by available NO À 3 or the availability of electrons (due to competition with O 2 ); ve À NO À 2 res is the realised e --flow to NO À 2 . Such competition for electrons was not implemented for ve À NO and ve À N 2 O because at the onset of NO-, N 2 O-and N 2 production, the total velocity of e --flow to all available e --acceptors (as predicted by the enzyme kinetics alone) never exceeded ve À maxTCA . Gas consumption and production takes place in the aqueous phase, but the gases are transported between aqua and the headspace depending on their concentrations in the two phases. Each gas in aqua, X aq (molN, Fig 2C), is modelled as a function of production, consumption (not applicable to N 2 ), and the net transport, where N 2 O aq and N 2 aq are initialised with zero, and NO aq is initialised with a negligible 1×10 −25 mol to avoid division by zero (in Eq 28).
(molN h -1 ) where Rr NO x (molN h -1 ) is the relevant NO À x /NO x reduction rate, and Tr N X represents the gas transport rate between aqua and the headspace (Eq 29; N.B. Tr N X < 0 for the net transport from aqua to the headspace).
The reduction of NO to N 2 O (Rr NO ) and N 2 O to N 2 (Rr N 2 O ) is modelled likewise as a function of the number of relevant cells and the velocity of e --flow to NO and N 2 O (mol ecell -1 h -1 ), respectively: (mol cell -1 h -1 ) where ve À maxNO (mol ecell -1 h -1 ) is the empirically determined maximum velocity of e --flow to NO per cell, [NO] aq (mol L -1 ) is the NO concentration in the liquid-phase, and K 1NO & K 2NO (mol L -1 ) are the equilibrium dissociation constants for the cNor/NOand cNor/(NO) 2 complex, respectively.
The transport of NO, N 2 O, and N 2 between the liquid and the headspace (Eqs 23-25) is modelled as: where k t is the empirically determined coefficient for the transport of each gas between the headspace and the liquid, k H(N) (molN L -1 atm -1 ) is the solubility of NO, N 2 O, or N 2 in water at 20°C, P N (= [N] g ×R×T, atm) is the partial pressure of each gas in the headspace, and [N] aq (mol L -1 ) represents the concentration of each gas in the liquid-phase.
The amount of NO and N 2 O in the headspace (NO x g , molN, Fig 2C) is a function of transport (Eq 29) and the disturbance by gas sampling. The latter is simulated as discrete events at time-points given as input to the model (equivalent to the sampling times in the experiment): where ΔNO x(S) is the net change in the amount of NO x g (molN), D (dilution) is the fraction of each gas removed and replaced by equal amount of He, and t s (h) is the time taken to complete each sampling. For N 2 , the model ignores the sampling loss because the N 2 production data to be compared with the model output are already corrected for the sampling disturbance [30]. Thus, the model estimates somewhat higher N 2 concentrations than that experienced by the organisms, which is acceptable, since the concentration of N 2 is unlikely to have consequences for the metabolism.

Parameterisation
Most of the parameter values used in the model are well established in the literature (see Table 2); however, uncertain parameters include K mO 2 , K mN 2 O , ve À maxO 2 , and ve À min . The maximum cell-specific velocity of e --flow to NO À 3 1×10 −14 mol ecell -1 h -1 [18] ve À maxNO À
K mN 2 O . In vitro studies of NosZ from Pa. denitrificans estimate the values for K mN 2 O = 5 μM at 22°C and pH 7.1 [42] and 6.7 μM at 25°C and pH 7.1 [43]. When our model was simulated with K mN 2 O in this range, given our empirically estimated ve À maxN 2 O [24], the simulated N 2 O reached concentrations much higher than that measured (see Results/Discussion). A more adequate parameter value (= 0.6 μM) was found by optimising K mN 2 O in Vensim. The value is within the range determined for soil bacterial communities [44].
ve À maxO 2 (Eq 17) could be estimated using the empirically determined cell yield per mole of electrons to O 2 (Ye À O 2 , cells per mol e -) and the maximum specific growth rate (μ, h -1 ): We are confident about the yields for the two C-substrates used, but the empirically determined μ for the butyrate treatments is suspiciously low (= 0.067 h -1 ), providing ve À maxO 2 = 2.45×10 −15 mol ecell -1 h -1 . Simulations with this value grossly underestimated the rate of O 2 depletion compared to measured, which forced us to estimate the parameter value by optimisation: ve À maxO 2 = 4.42×10 −15 and 4.22×10 −15 mol ecell -1 h -1 for the succinate-and butyrate treatments, respectively. These values give μ = 0.22 and 0.12 h -1 , respectively: for the succinate treatments, the value is very close to that empirically determined (= 0.2 h -1 ); for the butyrate treatments, the value seems more realistic than 0.067 h -1 .
ve À min (Eqs 2, 5 and 7) is the per cell velocity of e --flow to O 2 (ve À O 2 ) assumed to generate minimum ATP required for synthesising the initial molecules of denitrification enzymes. Since we lack any empirical or other estimations for this parameter, it is arbitrarily assumed to be the ve À O 2 when [O 2 ] aq reaches 1 nM. At this concentration, ve À min is determined by the Michaelis-Menten equation ve À min ¼ ve À maxO 2 Â½O 2 aq ðK mO 2 þ½O 2 aq Þ , using ve À maxO 2 and K mO 2 given above. The values obtained for the succinate-and butyrate-supplemented treatments = 1.96×10 −17 and 1.87×10 −17 mol ecell -1 h -1 , respectively, which for both the cases is 0.44% of ve À maxO 2 . To investigate the impact of ve À min on the model behaviour (r Na and r Ni , Eqs 1, 2, 4, 5, 6 and 7), sensitivity analyses were performed by simulating the model with ve À min corresponding to [O 2 ] aq = 5×10 −9 , 5×10 −10 , and 0 mol L -1 (see Results/Discussion).

Results/Discussion
Low probabilistic initiation of nar transcription, resulting in the fraction of the population with Nar < 100% To test the assumption of a single homogeneous population with all cells producing Nar in response to O 2 depletion, we simulated the model with the specific probability for a Z − cell to initiate nar transcription (r Na ) = 4 h -1 . This resulted in 98% of the cells possessing Nar within an hour (see Eqs 1-3). Evidence suggests that less than half an hour is required to synthesise denitrification enzymes [17,18], but an hour's time is assumed here to allow margin for error. The results show that, for all the treatments, the simulated NO À 2 production (mol vial -1 ) grossly overestimates that measured (Fig 3).
To find a reasonable parameter value, we optimised r Na for the 0% O 2 treatments, so that the simulated NO À 2 production matches that measured. The results (Table 3) suggest that a low probabilistic initiation of nar transcription (average r Na = 0.035 h -1 ) is adequate to simulate the measured NO À 2 kinetics (Fig 3). In the Butyrate, 7% O 2 treatment (Fig 3B), the simulated NO À 2 starts earlier, but the rate of accumulation is similar to that measured.
Once O 2 falls below a certain threshold, the production of Nar is assumed to trigger with r Na = 0.035 h -1 and last until a minimum of respiration is sustained by the e --flow to O 2 and N 2 O (ve À O 2 and ve À N 2 O ), assumed to fulfil the ATP needs for Nar production (Eqs 1 and 2). But To test the assumption of a single homogeneous population with almost all cells expressing nar in response to O 2 depletion, we forced our model to achieve 98% Narpositive cells (Z Na ) within an hour by setting the specific-probability of initiating nar transcription (r Na ) = 4 h -1 . This resulted in grossly overestimated rates of NO À 2 accumulation for all treatments (grey curves). In contrast, we simulated the model with r Na = 0.035 h -1 obtained through optimisation, resulting in a reasonable agreement with measurements for all treatments (except for an apparent time frameshift for the Butyrate, 7% O 2 treatment). Cell Diversification in Response to Anoxia in Paracoccus denitrificans the production of Nar sustained by ve À N 2 O was inconsequential for simulating the measured NO À 2 production, since NO À 3 was already exhausted when N 2 O started accumulating (i.e., when ve À N 2 O > 0). For this reason, the fraction that produced Nar (F Na , Eq 3 and Table 4) is calculated as functional (= 0.23-0.43) and theoretical (= 0.56-0.81), where the first is the fraction actually responsible for NO À 2 production (sustained by ve À O 2 ), but the latter also incorporates the fraction that produced Nar after the exhaustion of NO À 3 (sustained by ve À O 2 as well as ve À N 2 O ). The rationale behind calculating the theoretical F Na is the empirical data indicating that Nar transcription is not turned off in response to NO À 3 depletion [18]. Although our model cannot test the theoretical F Na , but the functional F Na suggests that, contrary to the common assumption, the measured NO À 2 kinetics can be neatly explained by only 23-43.3% of the population producing Nar in response to O 2 depletion.

Very low probabilistic initiation of nirS transcription
When we optimised the specific probability of nirS transcriptional activation (r Ni , see Eqs 4, 5, 6 and 7) to fit the measured data, the average r Ni = 0.004 h -1 (Table 3) adequately simulated the measured NO À 2 depletion and N 2 accumulation (Fig 4). The recruitment to denitrification lasted for 19.5-47.3 h, i.e., the time when [O 2 ] was below a critical concentration and the velocity of e --flow to O 2 and the relevant NO À x /NO x remained above a critical minimum (Eqs 4, 5, 6 and 7). The resulting fraction recruited to denitrification (F Ni , see Eq 8 and Table 4) was 0.08-0.18, the bulk of which depended on the e --flow to NO À 3 and N 2 O (instead of aerobic respiration).
To test whether the measured data could be explained without the recruitment sustained by NO À 3 and N 2 O respiration, we also simulated the model with the recruitment as a function of Table 3. Specific-probability of nar and nirS transcriptional initiation (r Na and r Ni , respectively) estimated for each treatment by optimisation (best match between the simulated and measured data). Avg. = 0.035 Avg. = 0.004 *Treatment refers to the C-source, initial oxygen concentration in the headspace (measured as headspace-vol.%), and initial NO À 3 concentration in the medium (mM).
doi:10.1371/journal.pcbi.1004621.t003 Table 4. The fraction of the population with Nar (F Na ) and NirS (F Ni ) estimated based on the optimal specific-probability of nar and nirS transcriptional initiation (r Na and r Ni ), respectively.   . To test the model's sensitivity to this parameter, we estimated r Na and r Ni by optimisation for different values of ve À min , relative to the default value = 1.95×10 −17 mol ecell -1 h -1 . For all cases, the model was able to adequately simulate the measured N 2 kinetics by moderate adjustments of r Na and r Ni . Table 5 shows the average optimal values of r Na and r Ni , obtained by fitting the simulated N 2 kinetics to the data, for different values of ve À min . S3 Fig shows adequate simulations of the measured N 2 kinetics assuming ve À min = 0, with optimised r Na = 0.033 h -1 and r Ni = 0.0033 h -1 . Thus, although assuming ve À min > 0 appears logical, it is not necessary to explain the measured data.

N 2 O kinetics
To simulate the N 2 O kinetics, we used ve À maxN 2 O = 5.5×10 −15 mol ecell -1 h -1 , empirically determined under similar experimental conditions as simulated here [24], and adopted the literature values for K mN 2 O [= 5 and 7 μM 42,43, respectively]. But with K mN 2 O = 5 μM, the model predicted N 2 O accumulation~10-20 times higher than measured for the~0% and~2-3 times higher for the 7% O 2 treatments (Fig 5). This forced us to simulate the model with the parameter value estimated by optimisation, providing the average K mN 2 O = 0.6 μM.
The measured N 2 O shows a conspicuous increase throughout the entire active denitrification period, and this phenomenon is neatly captured by the model. The reason for this model prediction is that the number of N 2 O producing cells (Z NaNi + Z Ni , Fig 2A) is low to begin with compared to the number of N 2 O consuming cells (Z − + Z Na + Z NaNi + Z Ni ), but the fraction of N 2 O producers will increase during the anoxic phase for two reasons: one is the recruitment to Z NaNi & Z Ni , another is the fact that the model predicts approximately three times faster cellspecific growth rate for Z NaNi & Z Ni than for Z − & Z Na (ve À N 2 O is identical for all groups, while ve À NO À 2 and ve À NO are both zero for Z − & Z Na but for Z NaNi & Z Ni , it holds that ve À To illustrate this phenomenon, we ran the model, assuming that the Z − & Z Na cells had no N 2 O reductase, resulting in a) constant N 2 O concentration throughout the entire anoxic phase and b) much higher N 2 O concentrations than measured (Fig 5). The overestimation is a trivial result, easily avoidable by increasing ve À maxN 2 O or decreasing K mN 2 O moderately. However, the prediction of a constant N 2 O concentration is clearly in conflict with the experimental data, and no parameterisation could force the model to reproduce this phenomenon other than the differential expression of denitrification genes. Hence, although there is room for further refinements, our default assumption regarding differential expression of NirS and NosZ explains the observed N 2 O kinetics: 1) abrupt initial accumulation to very low levels due to recruitment of relatively small numbers to the N 2 O  Table 2, i.e., K mN 2 O = 0.6 μM (estimated through optimisation) and ve À maxN 2 O = 5.5×10 −15 mol ecell -1 h -1 [24]. In contrast, each inserted panel shows the simulated N 2 O assuming 1) N 2 O consumption only by the cells producing N 2 O (Z NaNi + Z Ni ), and 2) the literature value for K mN 2 O = 5 μM [42]. The results show that the default simulation best explains the measured N 2 O kinetics, assuming its production by a small fraction (Z NaNi + Z Ni ) and consumption by the entire population (Z − + Z Na + Z NaNi + Z Ni ). producing pools (Z NaNi & Z Ni ), and 2) increasing N 2 O concentration due to recruitment and faster cell-specific growth of Z NaNi & Z Ni than that of the cells only consuming N 2 O (Z − + Z Na ).
This modelling exercise sheds some light on the possible role of regulatory biology of denitrification in controlling N 2 O emissions from soils. If all cells in soils had the same regulatory phenotype as Pa. denitrificans, their emission of N 2 O would probably be miniscule, and soils could easily become strong net sinks for N 2 O because the majority of cells would be 'truncated denitrifiers' with only N 2 O reductase expressed. It remains to be tested, however, if the regulatory phenotype of Pa. denitrificans is a rare or a common phenomenon among full-fledged denitrifiers. We foresee that further exploration of denitrification phenotypes will unravel a plethora of response patterns.

Conclusion
Using dynamic modelling, we have demonstrated that the denitrification kinetics in Pa. denitrificans can be adequately explained by assuming low probabilistic transcriptional activation of the nar and nirS genes and a subsequent autocatalytic expression of the enzymes. Such autocatalytic gene expressions are common in prokaryotes, rendering a population heterogeneous because of the stochastic initiation of gene transcription, with a low probability [45]. For N 2 O kinetics, our hypothesis was that a) the gas is produced by a fraction of the incubated population that is able to initiate nirS transcription with a certain probability, leading to a coordinated expression of nirS + nor via NO [16], and b) N 2 O is consumed by the entire population because, in response to anoxia, nosZ is readily induced by FnrP [24]. Our model corroborated this hypothesis by reasonably simulating the N 2 O kinetics with the specific-probability of nirS transcriptional activation = 0.004 h -1 , resulting in 7.7-22.1% of the population producing NirS + cNor (hence N 2 O), but all cells producing NosZ (hence equally consuming N 2 O).
Supporting Information S1 Dynamic Model. The folder contains the dynamic model used in this study 'Hassan_et_ al_2015_Pa._denitrificans.mdl'. The model requires Vensim (Double Precision), which is available at http://vensim.com/free-download/. The zip folder also contains files with the empirical data; these files are automatically loaded into the model when it is run. (ZIP) S1 Fig. Pa. denitrificans gas and NO À 2 kinetics. Typical gas kinetics (O 2 , NO, N 2 O, N 2 ) and NO À 2 accumulation in Pa. denitrificans during the transition from aerobic respiration to denitrification; batch cultures, n = 3; 20°C; Sistrom's medium; 2 mM KNO 3 and 7 vol% initial O 2 in the headspace. All the available NO À 3 (100 μmol vial -1 ) was recovered as NO À 2 before the onset of N-gas production. In previous experiments [17], N 2 O concentrations were below the detection limit of the system, but thanks to a new system with electron capture detector, the N 2 O kinetics were monitored with reasonable precision. Adapted from [18]. (TIF) S2 Fig. Comparison of measured and simulated data assuming stochastic initiation of nirS transcription with aerobic respiration being the only energy source for producing NirS + cNor. In each panel, the measured NO À 2 depletion (sub-panel) and N 2 accumulation (main panel; n = 3-4) are compared with simulations. The simulations here are to be compared with the default simulations (Fig 4), which were run assuming that the coordinated NirS + cNor production (via nirS transcriptional activation) is sustained by the energy generated by O 2 as well as NO À 3 and/or N 2 O reduction. The default simulations provided an average specific-probability of nirS transcriptional activation (r Ni ) = 0.004 h -1 (Eqs 4, 5, 6 and 7) by optimisation, allowing 7.7-22.1% of the population to produce NirS + cNor (Eq 8) in 19.5-47.3 h. To match the measured data here, the average r Ni had to be raised to 0.012 h -1 , since the time available for the enzyme synthesis shrank (= 3.5-16 h) due to a rapid exhaustion of O 2 . Comparatively, the assumption that the ATP from NO À 3 and/or N 2 O reduction should help cells produce denitrification enzymes seems more plausible and provide better agreement with the measured data. (TIF) S3 Fig. Measured vs. simulated N 2 kinetics assuming ve À min = 0. The default simulations are carried out assuming that for a cell to produce first molecules of Nar and NirS, a minimum of e --flow to an available e --acceptor (ve À min , mol ecell -1 h -1 ) is necessary to generate a minimum of ATP required for protein synthesis (Eqs 1, 2, 4, 5, 6 and 7). Although assuming ve À min > 0 seems logical, the measured N 2 kinetics are adequately simulated here with ve À min = 0. This shows that the assumption is not necessary to explain the measured data. (TIF)

Author Contributions
Conceived and designed the experiments: LLB LRB. Performed the experiments: ZQ LLB. Analyzed the data: JH LRB. Contributed reagents/materials/analysis tools: JH. Wrote the paper: JH LRB.