Agent-Based Modeling of Oxygen-Responsive Transcription Factors in Escherichia coli

In the presence of oxygen (O2) the model bacterium Escherichia coli is able to conserve energy by aerobic respiration. Two major terminal oxidases are involved in this process - Cyo has a relatively low affinity for O2 but is able to pump protons and hence is energetically efficient; Cyd has a high affinity for O2 but does not pump protons. When E. coli encounters environments with different O2 availabilities, the expression of the genes encoding the alternative terminal oxidases, the cydAB and cyoABCDE operons, are regulated by two O2-responsive transcription factors, ArcA (an indirect O2 sensor) and FNR (a direct O2 sensor). It has been suggested that O2-consumption by the terminal oxidases located at the cytoplasmic membrane significantly affects the activities of ArcA and FNR in the bacterial nucleoid. In this study, an agent-based modeling approach has been taken to spatially simulate the uptake and consumption of O2 by E. coli and the consequent modulation of ArcA and FNR activities based on experimental data obtained from highly controlled chemostat cultures. The molecules of O2, transcription factors and terminal oxidases are treated as individual agents and their behaviors and interactions are imitated in a simulated 3-D E. coli cell. The model implies that there are two barriers that dampen the response of FNR to O2, i.e. consumption of O2 at the membrane by the terminal oxidases and reaction of O2 with cytoplasmic FNR. Analysis of FNR variants suggested that the monomer-dimer transition is the key step in FNR-mediated repression of gene expression.


Introduction
The bacterium Escherichia coli is a widely used model organism to study bacterial adaptation to environmental change. As an enteric bacterium, E. coli has to cope with an O 2 -starved niche in the host and an O 2 -rich environment when excreted. In order to exploit the energetic benefits that are conferred by aerobic respiration, E. coli has two major terminal oxidases: cytochrome bd-I (Cyd) and cytochrome bo9 (Cyo) that are encoded by the cydAB and cyoABCDE operons, respectively [1,2]. Cyd has a high affinity for O 2 and is induced at low O 2 concentrations (micro-aerobic conditions), whereas Cyo has a relatively low affinity for O 2 and is predominant at high O 2 concentrations (aerobic conditions) [3]. These two terminal oxidases contribute differentially to energy conservation because Cyo is a proton pump, whereas Cyd is not [1,2]; however, the very high affinity of Cyd for O 2 allows the bacterium to maintain aerobic respiration at nanomolar concentrations of O 2 , thereby maintaining aerobic respiratory activity rather than other, less favorable, metabolic modes [4][5][6].
The transcription factors, ArcA and FNR, regulate cydAB and cyoABCDE expression in response to O 2 supply [7]. FNR is an iron-sulfur protein that senses O 2 in the cytoplasm [8,9]. In the absence of O 2 the FNR iron-sulfur cluster is stable and the protein forms dimers that are competent for site-specific DNA-binding and regulation of gene expression [10]. The FNR iron-sulfur cluster reacts with O 2 in such a way that the DNA-binding dimeric form of FNR is converted into a non-DNA-binding monomeric species [10]. Under anaerobic conditions, FNR acts as a global regulator in E. coli [11][12][13], including the cydAB and cyoABCDE operons, which are repressed by FNR when the O 2 supply is restricted [7]. Under aerobic conditions, repression of cydAB and cyoABCDE is relieved and Cyd and Cyo proteins are synthesized [3]. In contrast, ArcA responds to O 2 availability indirectly via the membrane-bound sensor ArcB. In the absence of O 2 ArcB responds to changes in the redox state of the electron transport chain and the presence of fermentation products by autophosphorylating [14][15][16]. Phosphorylated ArcB is then able to transfer phosphate to the cytoplasmic ArcA regulator (ArcA,P), which then undergoes oligomerization to form a tetra-phosphorylated octomer that is capable of binding at multiple sites in the E. coli genome [17,18], including those in the promoter regions of cydAB and cyoABCDE to enhance synthesis of Cyd and inhibit production of Cyo [7,17]. Because the terminal oxidases (Cyd and Cyo) consume O 2 at the cell membrane, a feedback loop is formed that links the activities of the oxidases to the regulatory activities of ArcA and FNR (Figure 1). These features of the system -combining direct and indirect O 2 sensing with ArcA,P and FNR repression of cyoABCDE, and ArcA,P activation and FNR repression of cydAB -result in maximal Cyd production when the O 2 supply is limited (micro-aerobic conditions) and maximal Cyo content when O 2 is abundant (aerobic conditions) [3].
Although the cellular locations of the relevant genes (cydAB and cyoABCDE), the regulators (ArcBA and FNR) and the oxidases (Cyd and Cyo) are likely to be fundamentally important in the regulation of this system, the potential significance of this spatial organization has not been investigated. Therefore, a detailed agent-based model was developed to simulate the interaction between O 2 molecules and the electron transport chain components, Cyd and Cyo, and the regulators, FNR and ArcBA, to shed new light on individual events within local spatial regions that could prove to be important in regulating this core component of the E. coli respiratory process.

Results/Discussion
The starting point for this work was the suggestion that spatial effects play an important role in controlling the response of the E. coli transcription factor FNR to changes in O 2 availability [15]. Early work showed that when dissolved O 2 was detectable in the culture medium the activity of FNR decreased, exhibiting ,50% activity in the range 2-5 mM dissolved O 2 and that under these conditions the external concentration of O 2 was equivalent to that in the bacterial cytoplasm [19,20]. However, when the cultures were supplied with low concentrations of O 2 in the input gas, dissolved O 2 in the culture medium could not be measured conventionally, because the O 2 is consumed by the respiratory activity of the bacteria. Therefore, a physiological measure of O 2 availability (the aerobiosis scale) was adopted to investigate the effects of low O 2 concentrations on bacteria [21]. On the aerobiosis scale, the minimum O 2 input that results in undetectable excretion of the fermentation product acetate under carbonlimiting conditions is defined as 100% aerobiosis (100% AU). As the amount of O 2 supplied to cultures decreases the specific rate of acetate production (q acetate ) increases to a maximum under anaerobic conditions (0% AU). Between these limits (0-100% AU) lies the micro-aerobic range, defined by the linear decrease in q acetate as the O 2 transfer rate increases, i.e. there is an inverse correlation between q acetate and aerobiosis [21]. When the O 2 supply exceeds the minimum required to abolish acetate excretion the AU value extends beyond 100%, reaching 217% when the culture medium is O 2 saturated. This physiological measurement of O 2 availability is reliable for cultures grown at values as low as 4% AU [21]. In previous experiments using this approach steadystate chemostat cultures were established at fixed points on the aerobiosis scale and samples were taken for measurement of the numbers of Cyd and Cyo molecules per bacterial cell (Table 1) [3]. In addition, Western blotting showed that the concentration of FNR in the cell was constant at 3000 protomers per bacterium across the aerobiosis scale ( Table 1).
The agent-based model was constructed to simulate interactions between three groups of agents: O 2 molecules, terminal oxidases, (Cyd and Cyo) and regulators (FNR and ArcBA), which interact with O 2 directly or indirectly. As shown in Figure 2, the initial state of the model is defined in the file 0.xml, which contains information on agent properties such as Id, type, status, position etc. -the data in Table 1 was used to provide numbers of Cyd, Cyo, ArcA and FNR molecules. There is no information available on the abundance of ArcB in the cell and thus it was assumed that there are 1000 ArcB molecules per cell based on the ,10:1 ratio of response-regulator to sensor kinase of another E. coli twocomponent system, PhoBR [22]. At the beginning of the simulation, the 0.xml file is read to establish the 3-D E. coli cell. From the initial state, the agents start moving randomly within their 3-D activity space, and interact with each other according to a set of pre-defined interaction rules and interaction radii (Tables 2  and 3), leading to an emergent state.
The algorithm of supplying O 2 to the cell responds to given AU levels and automatically calculates the rate of O 2 molecules supplied to the cell. Whilst the model runs, agent information is updated and stored as a series of XML files for further analysis and visualization.

Modeling the regulatory response to O 2 availability
The dynamics of the system were investigated by running the simulation through two cycles of transitions from 0-217% AU. Figure 3a shows a top view of a 3-D E. coli cell at 0% AU (steadystate anaerobic conditions). Under these conditions, the FNR molecules are present as dimers, all ArcB molecules are phosphorylated and the ArcA is octameric. The DNA binding sites for ArcA (120 in the model) and FNR (350 in the model) in the nucleoid are fully occupied. The number of ArcA sites was chosen from the data reported by Liu and De Wulf [18]. The model must include a mechanism for ArcA,P to leave regulated promoters. Upon introduction of O 2 into anaerobic steady-state chemostat cultures ,5 min was required to inactivate ArcAmediated transcription [15]. In the agent-based model presented here, each iteration represents 0.2 sec. Therefore, assuming that ArcA,P leaving the 120 DNA sites is a first order process, then t K is ,45 sec, which is equivalent to ,0.3% ArcA,P leaving the DNA per iteration ( Table 3). The number of FNR binding sites was based on ChIP-seq and ChIP-Chip measurements, which detected ,220 FNR sites and a genome sequence analysis that predicted ,450 FNR sites; thus a mid-range value of 350 was chosen [23][24][25]. Interaction with O 2 causes FNR to dissociate from the DNA (Table 3). Under fully aerobic conditions (217% AU) the FNR dimers are disassembled to monomers, and the different forms of ArcA coexist (Figure 3b). The ArcA-and FNR-DNA binding sites in the nucleoid are mostly unoccupied due to the lower concentrations of FNR dimers and ArcA octamers.

Author Summary
The model bacterium Escherichia coli has a modular electron transport chain that allows it to successfully compete in environments with differing oxygen (O 2 ) availabilities. It has two well-characterized terminal oxidases, Cyd and Cyo. Cyd has a very high affinity for O 2 , whereas Cyo has a lower affinity, but is energetically more efficient. Expression of the genes encoding Cyd and Cyo is controlled by two O 2 -responsive regulators, ArcBA and FNR. However, it is not clear how O 2 molecules enter the E. coli cell and how the locations of the terminal oxidases and the regulators influence the system. An agent-based model is presented that simulates the interactions of O 2 with the regulators and the oxidases in an E. coli cell. The model suggests that O 2 consumption by the oxidases at the cytoplasmic membrane and by FNR in the cytoplasm protects FNR bound to DNA in the nucleoid from inactivation and that dimerization of FNR in response to O 2 depletion is the key step in FNR-mediated repression. Thus, the focus of the agent-based model on spatial events provides information and new insight, allowing the effects of dysregulation of system components to be explored by facile addition or removal of agents. Examination of the system as it transits from 0% to 217% AU showed that the DNA-bound, transcriptionally active FNR was initially protected from inactivation by consumption of O 2 at the cell membrane by the terminal oxidases and by reaction of O 2 with the iron-sulfur clusters of FNR dimers in the bacterial cytoplasm -the progress of this simulation is shown in Video S1.
This new insight into the buffering of the FNR response could serve a useful biological purpose by preventing pre-mature switching off of anaerobic genes when the bacteria are exposed to low concentration O 2 pulses in the environment.
In the various niches occupied by E. coli, the bacterium can experience the full range of O 2 concentrations from zero, in the   and tetramers combining to form ArcA octamers, their numbers dropped after initially increasing. The rate at which the ArcA octomers accumulated (ArcA activation) after O 2 withdrawal was slower than the rate of ArcA inactivation (Figures 4b and c). In this implementation of the modeled transition cycle, the numbers of The third column lists the number of molecules used to initiate the model. 1 See Table 4. 2 The ArcA numbers were reported by Rolfe et al. [3]. 3 The ArcB numbers are assumed based on the ,10:1 ratio of response-regulator to sensor-kinase for another E. coli two-component system PhoB-PhoR [22]. 4 Numbers are listed in Table 1  ArcA octamer bound to DNA is assigned a probability of 0.3% to leave the DNA in every iteration. This 'off rate' is required because ArcA,P dephosphorylation occurs by the action of ArcB at the cell membrane (see text). Cyo The interaction radii (nm) were defined and refined as described in the Methods. ArcA octamers in the cytoplasm and bound to DNA did not reach that observed in the initial state before the second cycle of O 2 supply began, indicating that a longer period is required to return to the fermentation state.
The numbers of FNR dimer bound to binding sites and free FNR dimer (cytoplasmic FNR dimer) decreased when O 2 was supplied to the system (Figures 4g-h), but the rate was slower than that for ArcA inactivation, consistent with O 2 consumption at the membrane, which can be sensed by ArcB to initiate inactivation of ArcA, but lowers the signal for inactivation of FNR. When O 2 was removed from the system (from iteration 5001) FNR was activated over a similar timeframe to ArcA (Figures 4b and g), which was again consistent with previous observations [15]. As with ArcA, free FNR dimers and FNR monomers did not fully return to their initial states after O 2 supply was withdrawn in the model, indicating that further iterations are required to reach steadystate (Figure 4h-i). These results clearly indicate that the model is self-adaptive to the changes in O 2 availability, and the reproducible responses prove the reliability and robustness of the model. The ArcBA system simulated in this model is based on a preliminary biological assumption, and the agent-based model presented here should prove a reliable and flexible platform for exploring the key components of the system and testing new experimental findings.

Model validation
In order to validate the model with biological measurements of FNR DNA-binding activity estimated using an FNR-dependent lacZ reporter, the ArcBA system agents were removed from the model by setting their agent numbers to zero. The ArcBA system is an indirect O 2 sensor and does not consume O 2 , hence the FNR system was not affected by withdrawing ArcBA from the model, but this simplification increased simulation speed.
The O 2 step length and other model parameters were estimated using the experimental data obtained at 31% AU. Using the estimated O 2 step length at 31% AU and defining the step length of O 2 molecule, S O2 , as 0 at 0% AU, a linear model, S O2~k |C O2 , was constructed to predict the step lengths of O 2 at other AU levels, where k = 2.1 and C O2 represents the O 2 concentration at different AU levels ( Table 4). The O 2 step lengths predicted by this model were used to validate the model at 85%, 115% and 217% AU, and the accuracy of the linear model was shown by the good correlation between the model and experimental data.
Profiles of five repetitive simulations in which the simplified model was used to predict the numbers of active FNR dimers in steady-state cultures of bacteria grown at different AU values are presented in Figure 5. At 31% AU, the model implied that FNRmediated gene expression is unaffected compared to an anaerobic culture (0% AU), i.e. the number of FNR binding sites occupied in the nucleoid remained unchanged (Figures 5a and e). Even at 85% AU, ,80% of the FNR-binding sites remained occupied (Figures 5b and f). It was only when the O 2 supply was equivalent to .115% AU that occupation of the FNR-binding sites in the nucleoid decreased (Figures 5 c, d, g and h). These outputs matched the FNR activities calculated from the measurements of an FNR-dependent reporter (Table 5) and thus demonstrate the abilities of the model to simulate the general behavior of FNR dimers in steady-state cultures of E. coli.
A second validation approach using two FNR variants that are compromised in their ability to undergo monomer-dimer transitions was adopted. The FNR variant FNR I151A can acquire an iron-sufur cluster in the absence of O 2 , but subsequent dimerization is impaired [26]. The FNR D154A variant can also acquire an iron-sulfur cluster under anaerobic conditions, but does not form monomers in the presence of O 2 [26]. To mimic the behavior of these two FNR variants the interaction radius for FNR dimer formation was changed in the model. Thus, the interaction distance for wild-type FNR monomers, which was initially set at 6 nm (r 3 , Table 3) was increased to 2000 nm for the FNR D154A variant, essentially fixing the protein as a dimer, or decreased to 2.5 nm for the FNR I151A variant, making this protein predominantly monomeric under anaerobic conditions. The results of simulations run under aerobic (217% aerobiosis) and anaerobic conditions (0% aerobiosis) suggested that under aerobic conditions wild-type FNR and FNR I151A should be unable to inhibit transcription from an FNR-repressed promoter (i.e. the output from the reporter system is 100%), whereas FNR D154A should retain ,50% activity (Table 6). Under anaerobic conditions, wild-type FNR was predicted to exhibit maximum repressive activity (i.e. 0% reporter output), whereas FNR I151A and FNR D154A mediated slightly enhanced repression compared to the simulated aerobic conditions (Table 6). To test the accuracy of these predictions, the ability of wild-type FNR, FNR I151A and FNR D154A to repress transcription of a synthetic FNR-regulated promoter (FFgalD4) under aerobic and anaerobic conditions was tested [27]. The choice of a synthetic FNR-repressed promoter was made to remove complications that might arise due to ironsulfur cluster incorporation influencing the protein-protein interactions between FNR and RNA polymerase; in the reporter system chosen FNR simply occludes the promoter of the reporter gene and as such DNA-binding by FNR controls promoter   activity. The experimental data obtained matched the general response of the FNR variants in the simulation, but not very precisely for FNR D154A, with the experimental data indicating more severe repression by FNR D154A under both aerobic and anaerobic conditions than predicted (Table 6). This suggested that the interaction radius (r 2 = 5 nm; Table 3), which controls the binding of FNR to its DNA target required adjustment to enhance DNA-binding of the FNR D154A variant. Therefore, the simulations were rerun after adjusting r 2 to 7 nm for all the FNR proteins considered here. The results of the simulations for both FNR variants now matched the experimental data well (Table 6). However, it was essential to ensure that the adjustment to r 2 did not significantly influence the model output for wild-type FNR. Therefore, simulations of the behaviour of wild-type FNR at 31, 85, 115 and 217% aerobiosis were repeated using the adjusted r 2 value of 7 nm. The model output was very similar to those obtained when r 2 was at the initial value of 5 nm (Table 7). These analyses imply that for FNR D154A, which is essentially fixed in a dimeric state, the rate of binding to the target DNA governs transcriptional repression, but for wild-type FNR the upstream monomer-dimer transition is the primary determinant controlling the output from the reporter.

Concluding remarks
The FNR switch has been the subject of several attempts to integrate extensive experimental data into coherent models that account for changes in FNR activity and target gene regulation in response to O 2 availability [15,[28][29][30][31]. These models have provided estimates of active and inactive FNR in E. coli cells exposed to different O 2 concentrations and the dynamic behavior of the FNR switch. The ability of FNR to switch rapidly between active and inactive forms is essential for it to fulfill its physiological role as a global regulator and the models are able to capture this dynamic behavior. Thus, it is thought that the 'futile' cycling of FNR between inactive and active forms under aerobic conditions has evolved to facilitate rapid activation of FNR upon withdrawal of O 2 and hence the physiological imperative for rapid activation has determined the structure of the FNR regulatory cycle [30,31]. However, it is less clear from these approaches how the system avoids undesirable switching between active and inactive states at low O 2 availabilities (micro-aerobic conditions, .0%-,100% AU). To achieve rapid FNR response times it has been suggested that minimizing the range of O 2 concentrations that constitute a micro-aerobic environment, from the viewpoint of FNR, is advantageous [31]. Unlike previous models of the FNR switch, The model simulation data was taken by averaging the model outputs at steady-state. For these results, the values shown are the averages and standard deviations. The last column shows experimentally determined transcription from the FNR-dependent FF-41.5-lacZ promoter; the transcriptional output at 0% AU was set to 100. 1 The total numbers of FNR dimers were calculated from the Western blots and FNR-dependent promoter activities. 2 The measured output from an FNR-dependent promoter was reported by [15,38]. doi:10.1371/journal.pcbi.1003595.t005 For the simulations, the standard deviations were calculated from 5 repeats for FNR and 13 repeats for FNR I151A and FNR D154A. The aerobic condition in simulation was modelled at AU level 217%, and the anaerobic condition was modelled at AU level 0%. r 3 is the interaction radius between the FNR dimer and its cognate DNAbinding site.
the agent-based model described here recognizes the importance of geometry and location in biology. This new approach reveals that spatial effects play a role in controlling the inactivation of FNR in low O 2 environments. Consumption of O 2 by terminal oxidases at the cytoplasmic membrane and reaction of O 2 with the iron-sulfur clusters of FNR in the cytoplasm present two barriers to inactivation of FNR bound to DNA in the nucleoid, thereby minimizing exposure of FNR to micro-aerobic conditions by maintaining an essentially anaerobic cytoplasm for AU values up to ,85%. It is suggested that this buffering of FNR response makes the regulatory system more robust by preventing large amplitude fluctuations in FNR activity when the bacteria are exposed to micro-aerobic conditions or experience environments in which they encounter short pulses of low O 2 concentrations. Furthermore, investigation of FNR variants with altered oligomerization properties suggested that the monomer-dimer transition, mediated by iron-sulfur cluster acquisition, is the primary regulatory step in FNR-mediated repression of gene expression. It is expected that the current model will act as a foundation for future investigations, e.g. predicting the effects of adding or removing a class of agent to identify the significant regulatory components of the system.

Measurement of the rate of O 2 supply
Knowledge of the rate of O 2 supply, R O2 , to the E. coli cells was required in order to simulate the response of the regulators of cydAB and cyoABCDE to different O 2 availabilities. Therefore, uninoculated chemostat vessels were used to measure dissolved O 2 concentrations, D O2 , as a function of the percentage O 2 in the input gas, P i , in the absence of bacteria. This allowed the rate at which O 2 dissolves in the culture medium to be calculated from the equation: D O2~RO2 |P i , yielding R O2 = 5.898 mmol/L/min. The number of O 2 molecules distributed to a single bacterial cell was then calculated from the following equation: Num O2~DO2 |N A |V cell |n (where, N A is the Avogadro constant (6.022610 23 ), V cell is the volume of E. coli cell (0.3925 mm 3 ) and as a constant for this equation, n (3.3610 29 ) includes the unit transformations, min to sec (60 21 ) and mmol to mol (10 26 ), and the time unit represented by an iteration (0.2 sec).

Control of agent mobility
In the model the individual agents (Cyd, Cyo, ArcB, ArcA, FNR and O 2 ) are able to move and interact within the confines of their respective locations in a 3-D-cylinder representing the E. coli cell. To control the velocity of agents, the maximal distances they can move in 3-D space during one iteration (step length) were predefined (Table 4). Thus, a step length is pre-defined in program header file (.h) and for each movement, this is multiplied by a randomly generated value within [0,1] to obtain a random moving distance, which in turn is directed towards a 3-D direction (movement vector) that was also randomly generated within defined spatial regions. An example is shown in Figure 6 to illustrate the movements of an O 2 molecule when it enters the cell.

Estimating interaction radii
Interactions between agents depend upon the biological rules governing their properties and being in close enough proximity to react. The interaction radius of an agent encapsulates the 3-D space within which reactions occur. As the interaction radii cannot be measured, they were first estimated on the basis of known biological properties. For the radii r 1…4 , r 12 and r 13 (Table 3), arbitrary values were initially set at 31% AU, and the model was then trained to match the experimental result for the number of FNR dimers at 31% AU ( Table 5). The modeled output of FNR dimer number at steady-state was compared with the experimental data, and the difference suggested re-adjustment of interaction radii. The adjusted radii were then tested against the FNR dimer numbers at 85%, 115% and 217% AU (Table 5) during model validation, and the results indicate that the interaction radii values are capable of describing the behavior of the system. The Figure 6. The 3-D movement of an O 2 molecule during five successive iterations. Within a pre-defined limit (step length) the agent moves a random distance per iteration. For O 2 this is also used to imitate the diffusion along the concentration gradient according to Fick's first law, in which the flux goes from regions of high concentration to regions of low concentration with a magnitude that is proportional to the concentration gradient (spatial derivative). The adjustment of step length of O 2 (Table 4) affects the spatial moving speed, which is simply represented as the greater the step length, the faster the movement. The angle range was defined for O 2 molecules that are outside the cell to enable them move towards cell. The O 2 molecules inside the cell and all other molecules move in any direction within their defined spatial regions. doi:10.1371/journal.pcbi.1003595.g006 interaction radii of Cyd and Cyo with O 2 reflect their relative affinities for O 2 (i.e. Cyd has a high O 2 affinity and thus reacts more readily, 7 nm interaction radius, than Cyo, which has a lower affinity for O 2 , 3 nm interaction radius). As, thus far, no accurate biological data is available for ArcBA system, the radii r 5…11 were arbitrarily defined and were refined by training the model to match current biological expectations.

Model description
The rod-shaped E. coli cell was modeled as a cylinder (500 nm62000 nm) [32] with the nucleoid represented as a sphere with a diameter of 250 nm at the centre of the cell. The experimentally-based parameters and locations of the agents in their initial state are listed in Table 2. As the number of ArcB molecules has not been determined experimentally, this value was arbitrarily assigned (see above). The interaction rules for the agents are shown in Table 3 (additional descriptions of an exemplar agent (O 2 ) and the rules for ArcBA and FNR are provided in Figure S1, Table S1 and Text S1). These rules, combined with the interaction radii, determine the final status of the system. The scale of the model is such that high performance computers are required to implement it, and the flexible agent-based supercomputing framework, FLAME (http://www.flame.ac.uk) acted as the framework to enable the simulation [33,34]. For more information on FLAME see Figure S2 and Text S2.

Measurement of the activities of FNR and FNR variants in vivo
Plasmids encoding the FNR variants were constructed by sitedirected mutagenesis (Quikchange, Agilent) of pGS196, which contains a 5.65 kb fragment of wild-type fnr ligated into pBR322 [35]. The three isogenic plasmids pGS196 (FNR), pGS2483 (FNR I151A) and pGS2405 (FNR D154A) were used to transform E. coli JRG4642 (an fnr lac mutant strain) containing a pRW50-based reporter plasmid carrying the lac-operon under the control of the FFgalD4 promoter [27].
b-Galactosidase assays were carried out as described previously on strains grown in LBK medium at pH 7.2 containing 20 mM glucose [36,37]. Cultures were grown either aerobically (25 ml culture in a 250 ml flask at 250 rpm agitation with 1:100 inoculation) or anaerobically (statically in a fully sealed 17 ml tube with 1:50 inoculation). Cultures (three biological replicates) were grown until mid-exponential phase (OD 600 = 0.35) before assaying for b-galactosidase activity. Figure S1 Stategraph for O 2 molecules. In order to describe the model clearly, every agent is given a formal description to illustrate its states, memory, functions, and relevant messages that it sends out or receives from other agents (see Table  S1). The stategraph for an oxygen agent is shown in the diagram. (TIFF) Figure S2 Building a FLAME simulation file. The agent definition (written in XMML) is parsed by a FLAME model parser, called xparser, which generates the simulation code. In the GCC environment, the code is compiled with the message board library, libmboard. The initial agent population settings are set in 0.xml file as the starting status of the model. (TIFF) Text S1 Additional description of interaction rules for the regulatory systems, ArcBA and FNR.