Dynamic Boolean modelling reveals the influence of energy supply on bacterial efflux pump expression.

Antimicrobial resistance (AMR) is a global health issue. One key factor contributing to AMR is the ability of bacteria to export drugs through efflux pumps, which relies on the ATP-dependent expression and interaction of several controlling genes. Recent studies have shown that significant cell-to-cell ATP variability exists within clonal bacterial populations, but the contribution of intrinsic cell-to-cell ATP heterogeneity is generally overlooked in understanding efflux pumps. Here, we consider how ATP variability influences gene regulatory networks controlling expression of efflux pump genes in two bacterial species. We develop and apply a generalizable Boolean modelling framework, developed to incorporate the dependence of gene expression dynamics on available cellular energy supply. Theoretical results show that differences in energy availability can cause pronounced downstream heterogeneity in efflux gene expression. Cells with higher energy availability have a superior response to stressors. Furthermore, in the absence of stress, model bacteria develop heterogeneous pulses of efflux pump gene expression which contribute to a sustained sub-population of cells with increased efflux expression activity, potentially conferring a continuous pool of intrinsically resistant bacteria. This modelling approach thus reveals an important source of heterogeneity in cell responses to antimicrobials and sheds light on potentially targetable aspects of efflux pump-related antimicrobial resistance.

Introduction

Antimicrobial resistance (AMR) is a huge global health challenge, contributing to over 700 000 deaths worldwide [1] and estimated to exceed cancer mortality by 2050 if neglected [2]. Protein complexes called efflux pumps are an important mechanism for intrinsic (and acquired) resistance to antibiotics [3-5]. Efflux pumps are capable of transporting a range of substances from the cell compartment to the extracellular medium, reducing intracellular accumulation of substances including antibiotics, and thus conferring resistance. Their central role in facilitating clinically relevant AMR makes efflux pumps an attractive drug target, with much research directed at efflux pump inhibitors as a novel therapeutic to overcome efflux-mediated resistance [6-10].

A current mystery in the activity of efflux pumps is the causes and effects of cell-to-cell heterogeneity in their expression. Recent data demonstrate efflux pump gene expression is not homogeneous throughout the population—instead, bacterial populations display significant cell-to-cell variability [11-13]. Heterogeneity in efflux activity has been suggested to have further downstream influence on AMR. One study, by Sánchez-Romero & Casadesús [13], reported a correlation between bacteria with high efflux pump activity in the absence of antibiotic compounds and increased resistance to antibiotics. Additionally, a study by El Meouche & Dunlop [11] has shown a sub-population of bacteria with higher efflux pump expression has lower expression of mismatch repair genes, facilitating mutations which can lead to permanent genetic changes bestowing resistance. However, the underlying dynamics and influences of efflux heterogeneity in cellular populations remain to be uncovered.

One potential source of heterogeneity in gene expression dynamics and interactions is cell-to-cell differences in available energy supply [14,15]. Heterogeneity in adenosine triphosphate (ATP) levels in bacteria (and other species; see §3) is increasingly being characterized experimentally [16,17]. The quantified range of intracellular ATP can cover almost an order of magnitude, varying between 0.32 and 2.76 mM in Escherichia coli [16], for example. Although not uncontroversial, the ubiquity of cell-to-cell variability in ATP across life suggests it is a genuine effect [18-20]. Following theoretical work showing the potentially profound influence of this variability on gene regulatory networks (GRNs) [14,15], we here investigate the hypothesis that cell-to-cell differences in available energy can contribute to cell heterogeneity in efflux pump expression, intrinsic resistance to antimicrobials, and response to different environmental stressors.

We previously applied an energy-dependent ordinary differential equation (ODE) framework to a simple GRN, describing the behaviour of many naturally occurring decision-making circuits [14]. We predicted that differences in cellular energy levels will cause differences in the dynamics and stable outcomes of cellular decision-making. However, while ODE models can in principle be used for GRNs of arbitrary size, the associated parameter space rapidly increases with network size, making it harder to investigate general principles. Here, we develop and use an alternative theoretical framework based on Boolean network models [21,22], simplifications that have been successfully applied to diverse biological systems [23-25] (although to date not often used to understand bacterial virulence and resistance mechanisms [26,27]). Despite this widespread and successful use, there has been little consideration of the fact that the processes in Boolean GRN models correspond biologically to ATP-dependent processes [21,26,28]. This ATP dependence is important, because it generally influences the dynamics of gene expression [29]. We therefore proceed by developing a simple but highly generalizable modification for including energy dependence in Boolean models of regulatory networks, and use this theoretical framework to explore the effects of ATP variability in experimentally derived GRN models of efflux pump expression in E. coli and Salmonella.

Results

Literature-based model construction accounting for energy availability

We first sought to construct a Boolean modelling framework describing efflux pump expression dynamics, supported by data and yielding verifiable predictions about biological behaviour. We considered AcrAB-TolC, an efflux pump associated with clinically important drug resistance in E. coli [30-32] and Salmonella [33-36]. TolC is constitutively expressed [37] and predicted to be integrated within at least seven additional efflux transport systems in E. coli [38], suggesting the dynamics of acrAB are more limiting on producing the complete efflux pump system within each cell. We performed a broad literature review (see the electronic supplementary material, S1) of efflux pump gene regulation in E. coli and Salmonella to define a core GRN of acrAB in each species (figure 1b,c). This review defined the structure of our model i.e. the nodes, regulatory interactions and edge weights in the corresponding wiring diagrams (electronic supplementary material, figure S1A,B).

Figure 1

Boolean network modelling of gene regulatory network architectures controlling efflux pump gene expression. (a) Schematic overview of the Boolean network framework. Regulatory architectures (i) are converted into wiring diagram form (ii). These wiring diagrams include ‘ghost nodes’ (G1) if necessary, providing a constant signal to ensure constitutively expressed genes are ‘on’ unless switched ‘off’ by a regulator. The considered network is repeatedly simulated over a set number of timesteps (iii), from some initial configuration, producing individual and averaged pictures of network behaviour through transition matrices and expression dynamics of each node per timestep (iv). (b,c) Coarse-grained E. coli (b) and Salmonella (c) AcrAB efflux pump gene regulatory networks, containing a collection of genetic regulators (blue circles) that interact with each other (directed edges) and govern the expression level of efflux genes acrAB. Full layouts with parameterizations are given in the electronic supplementary material, figure S1. Antibiotics, or other noxious stressor substances (lightning bolt), can reversibly bind to a gene product (transcription factor) removing its ability to bind to DNA and regulate its target (displayed through orange edges), with consequent downstream effects through the network. Removal of the stressor is integrated through a repressive interaction from acrAB to the stress node (dotted edges).

Escherichia coli contains a global activator of efflux genes acrAB, called marA (of the marRAB operon), which is repressed locally by the product of marR [39]. Similarly, in Salmonella, a global activator ramA is repressed locally by RamR, the product of ramR [40]. In both networks, repressor AcrR, the gene product of acrR, down-regulates its own expression and prevents over-expression of efflux genes acrAB [41] (see figure 1b,c for each complete GRN). acrAB encodes the proteins AcrA and AcrB [42], which can assemble together with TolC to form the AcrAB-TolC system [43,44] that provides the physical means to expel unwanted substances from the cell (such as toxic metabolic intermediates or bile salts [45,46]).

The expression of efflux pumps is inducible by the presence of environmental signals, or chemical compounds (such as antibiotics), which we will generally refer to as ‘stressors’. The inducibility of efflux genes enables bacteria to adapt to a wide range of environmental conditions. Noxious substances, such as antibiotics, can bind to GRN products MarR, RamR and AcrR, inhibiting their ability to bind to DNA, and initiating a response cascade that leads to the upregulation of efflux genes acrAB [47-52]. In our model the expression of acrAB corresponds to the capacity of the assembled efflux pump to export molecules.

To avoid the aforementioned parameterization issues in continuous models, we developed a general Boolean framework that captures the influence of energy availability. Our model represents genes as nodes in a network, where each gene can be ‘on’ (expressed) or ‘off’ (not expressed) at a given time. Edges describe activating and repressing interactions, and the state of a network changes over time with asynchronous update rules following these interactions (see Methods), capturing the stochastic nature of gene regulatory dynamics. Edge weights for our model were chosen through a manual tuning process, seeking to reproduce our ‘training’ observations given the imposed structure of the regulatory network from the literature review. In a Boolean network, edge weights encode logical interactions rather than specific physical rates or quantities, and therefore there is substantial flexibility in the particular values used (see §3).

To explicitly capture the fact that the fundamental processes (hence all gene interaction ‘arrows’) involved in gene expression depend on ATP, we modulate the rate with which regulatory interactions from node X to node Y are applied according to cellular energy levels (see Methods). This reflects the fact that, for example, cells with low ATP will have a lower rate of transcription elongation per timestep [53], slowing the dynamics of the biochemical intermediates involved in these regulatory interactions, and thus slowing the interactions themselves [14]. While fast fluctuations in ATP supply may exist within cells, as gene expression processes are dynamically slower, we assume the network is exposed to a time-averaged ATP level, and thus define three energy levels ‘high’, ‘intermediate’ and ‘low’, respectively, involving a characteristic mean rate of 1, 1/2 and 1/10 per simulation, corresponding to the aforementioned order of magnitude range observed in biology (see Methods).

Efflux network models reproduce experimentally observed behaviour in Escherichia coli and Salmonella

Having constructed our model based on a set of experimental observations, we next asked whether it made predictions that could be tested using independent experimental data. In the electronic supplementary material, we present a table assessing the predictions of our model, for our choice of edge weights, with respect to a multiscale set of experiments at both the population and single-cell level, where agreement is demonstrated across a wide variety of observations (electronic supplementary material, table S1). The overall behaviour of both systems—which we will describe in detail below—predicts a range of coarse- and fine-grained features observed in independent experiments, including pronounced gene expression heterogeneity at the single-cell level, prominent increases in gene expression when a short or long stress signal is provided, the expression of repressive genes in response to stress [41,54]), and the details of expression levels of different actors in the GRN in the presence and absence of stress.

Stress-free efflux network behaviour with a maximal cellular energy budget

To dissect the dynamics of these networks in detail, we first explored the stress-free population-level behaviour of each network. In the absence of stress, the population-wide mean behaviour in expression levels approaches a steady state (pre-stress dynamics in figure 2). We found at high energy availability, initial conditions had no influence on this population-level equilibrium of the E. coli and Salmonella networks (see the electronic supplementary material, S5). We therefore proceed by investigating network behaviour starting from a particular initial state—which models an ‘unstressed’ baseline condition—in the rest of the study and focus on the E. coli network before comparing the behaviour of the two networks.

Figure 2

Mean expression dynamics for E. coli from the unstressed initial condition in response to a short or long stress period. Time-series dynamics of mean expression (coloured line for low (red), intermediate (green) and high (blue) energy levels) for ((a), (c)) efflux genes acrAB and ((b), (d)) remaining regulatory components in E. coli; all panels include stressor behaviour. Behaviour is shown when the network is exposed to a (a)–(b) short (t = 150) or (c)–(d) long (t = 150 − 174) stress period, from the unstressed initial condition. The initial time region of each regulatory component prior to stress application (grey shaded region) displays the complete stress-free behaviour. Note: standard deviation (corresponding coloured shaded area) is included within efflux panels to represent the uncertainty, but the actual expression level cannot exceed the [0, 1] domain.

To explore the variability within the equilibrium population without energetic limitations, we calculated the expression level coefficient of variation (standard deviation divided by mean) for each genetic component in the E. coli network. In the absence of stress, expression level coefficient of variation (ELCV), corresponding to population variability in gene expression, ranged between 1.0 and 1.35 at equilibrium (see stress-free period in figure 2). This diversity is owing to population-level asynchrony in expression dynamics; in simulated individual bacteria, network components are expressed in heterogeneous pulses (pre-stress dynamics in figure 3). This diversity, particularly in acrAB, suggests that the network structure itself may support a degree of ‘bet-hedging’ in the absence of stress, where some members of a population are primed to respond to stressors at any given time, agreeing with experimentally observed behaviour in Salmonella enterica [13] (see §3).

Figure 3

Simulations predict E. coli expresses acrAB in heterogeneous pulses in the absence of stress. Qualitative plots displaying acrAB ‘on’/‘off’ (dark/light grey) dynamics for 20 randomly selected individual E. coli simulations, when exposed to a (a)–(c) short or (d)–(f) long stressor period. Panels display acrAB dynamics between timesteps t = 120 −250 at low, intermediate and high energy availability. The short stressor period is a single timestep pulse activated at t = 150 (blue line) and the long stressor period is activated at t = 150 (blue line) until t = 174.

Efflux network behaviour with a maximal cellular energy budget after exposure to a stressor

We next asked how exposure to a stressor affects network behaviour. Following a simulated stress-free period as in the previous section, we simulated exposure to a stressor for either a short (one timestep) or long (25 timesteps) stress period. In response to stress, the mean expression levels of all genes in the system change, with different rates and magnitude (figure 2). acrAB mean expression rises throughout stress exposure, thus becoming higher for the longer stress duration (figure 2a,c). After the long stress period, the high energy population displays homogeneous behaviour (active efflux expression), increasing the capacity of each cell to remove the stressor (see figure 3d to qualitatively observe this behaviour). The high energy population can then rapidly, and collectively, expel the noxious substance, promptly restoring baseline conditions in mean expression following the steps discussed previously.

When stimulated by a short stress (figure 2a), high variability is observed in efflux pump activity throughout the entire stress response, corresponding to diverse responses across the bacterial population. Individual simulation dynamics (figure 3a) suggest this behaviour originates from sub-populations of cells developing unsynchronized responses to stress, and a delayed development of the efflux response. The unsynchronized dynamics prolong the population-level stress response, potentially benefiting a clonal population by providing a hedging mechanism against future stress.

Once the simulated stressor is removed, repressive transcription factors resume downregulation of their GRN target, initially reducing the mean expression level of marRA and acrR, followed by a reduction in acrAB mean expression (figure 2). This process switches off the response cascade, allowing acrAB mean expression to return to the pre-stress population equilibrium (figure 2).

We next asked when the population of E. coli cells reached an equilibrium in acrAB expression after each perturbation—that is, when the statistics of efflux gene expression levels no longer changed with time in the population after the stress period. We found that the number of timesteps for acrAB to attain a post-stress population equilibrium, in the maximal energy case, was similar after a short or long stress (nine timesteps versus eight timesteps, respectively, figure 2a,c), and the equilibrium eventually attained was the same for both stress durations. These timings differ when energy is limited, as described in the following sections.

To understand these dynamics in more depth, we next asked how the different patterns of gene expression between which a network can transition vary when a cell is exposed to stress (reflecting the diversity of behaviours supported by the network). We define ‘accessible transitions’ to be the set of expression patterns that can be reached directly from a given network state. To explore these transition patterns, we analysed the transition matrices describing the dynamics of the E. coli regulatory network in the presence and absence of a stressor, visualized through heatmap images (figure 4; see figure 1a for illustration of the process to yield these outputs).

Figure 4

Escherichia coli transition dynamics depend on cellular energy and the presence of stress. Heatmap plots showing E. coli’s transition matrix. Each row contains the transition probability distribution from each of possible global states (rows) to each subsequent state (columns), in the absence (a)–(c) and presence (d)–(f) of a stressor, with increasing energy availability; here, owing to the coupling of marR and marA (as they form the marRA operon), we do not artificially implement unattainable network states, for example, ‘1011’. Each matrix element is coloured to indicate the probability p (using a logarithmic scale) that global state A transitions to global state B from timestep t = τ to t = τ + 1. Hierarchical clustering of global states is calculated using the destination possibilities (t = τ + 1) for each t = τ global state, with clustering employing a Euclidean distance metric to determine the distances between each element. Global state binary vectors are displayed in the order [marR, marA, acrR, acrAB].

In the absence of stress, all transitions are accessible (figure 4a). When a stress is applied, the accessible transitions for each network state are significantly restricted (compare figure 4a to d). The presence of a stressor removes the negative regulation within E. coli’s regulatory architecture (figure 1b), and consequently the ability to actively switch components ‘off’. The stressed network contains a steady-state attractor—a single state that when reached stays constant—corresponding to all components of the network being ‘on’. Under stress-free conditions no such single steady-state attractor exists (figure 4a). The presence of this attractor suggests that an E. coli cell under stress shifts its regulatory poise to allow the expression of each GRN component, allowing a response to the stressor and a faster dynamic response when the stress is removed.

Efflux network behaviour is modulated by available cellular energy

We next sought to understand how energy availability changes network behaviour in E. coli cells, in environments with and without a stressor. We first explored how the amount of energy available to a cell affects the accessible network transitions. Both with and without stress, the accessible transitions for each network state are the same across all energy levels, but the probabilities of these transitions vary dramatically with energy (figure 4). At lower energies, it is more likely that the system remains in the same state from one timestep to the next; increasing energy distributes the probability of each accessible transition more evenly among the total accessible transitions.

When a stressor is present, accessible transitions are restricted to subsets of the stress-free case at each energy level (figure 4a–c versus d–f). As energy increases, the likelihood of transitioning between supported network behaviours expands, transitioning more directly towards the steady-state attractor, ‘1111’. Therefore, increasing energy availability supports more rapid shifting of behaviour in the E. coli efflux network (figure 4), particularly towards a high efflux state.

Next, we asked whether the equilibrium mean expression of efflux pump genes changed at different energy levels, both in the presence and absence of stress. Perhaps counterintuitively, we found that increasing energy availability decreased the mean expression equilibrium of acrAB in simulated stress-free conditions (figure 2). In the presence of a stressor, the population response is more rapid and of higher magnitude as energy availability increases (figure 2). Simulated single-cell dynamics predict cell-to-cell variability in the time taken to transition to an active efflux expression state during stress exposure in each energy state (figure 3). The delayed response is more pronounced as energy decreases (figure 3d versus f). Following either stress period, the network returns to the pre-stress equilibrium more rapidly at higher energy levels (figure 2).

Energy levels influence heterogeneity in individual expression level dynamics

To investigate how energy level influences the detailed dynamics of single-cell pulses of efflux pump genes acrAB, we calculated the mean pulse length in acrAB expression at different energy levels (figure 5). Intuitively, from time-scaling effects, mean pulse length increased as energy decreased (figure 5a), agreeing with qualitative observations (figure 3a–c). The introduction of a short stress results in only a very small transient increase in mean pulse length, more pronounced at higher energy (figure 5). Pulse lengths during the long stress period are more dependent on energy level, with pulse statistics in high-energy cells changing more dramatically than in low-energy cells.

Figure 5

Energy availability and stress modulate E. coli acrAB pulse variability. Scatter plots displaying the (logarithmic) mean and coefficient of variation of acrAB pulse lengths (coloured dots for low (red), intermediate (green) and high (blue) energy levels) starting at each timestep in E. coli. Statistics are displayed between timesteps t = 50–250 to observe pulse behaviour before, during and after a (a)–(b) short or (c)–(d) long stress period. The short stressor period is a single timestep pulse activated at t = 150 and the long stressor period is activated at t = 150 until t = 174.

We next addressed our central hypothesis, whether energy availability could be a cause of the observed cell-to-cell variability in efflux pump expression. To explore this, we calculated the dynamics of the ELCV for E. coli acrAB at each energy level. As seen previously, substantial cell-to-cell variability exists in the system in the absence of a stressor, reflecting asynchrony and intrinsic noise in the system (figure 2). This variability is higher for low-energy cells (ELCV 1.34 at low energy, 1.01 at high energy). Exposure to stress reduces this variability, as the population synchronizes and responds in a similar way to the imposed stressor. Exposure to a longer stress period imposed a more substantial decrease in efflux variability than a shorter stress at all energy levels (figure 2a,c) and, as observed for the maximal energy case (§2.4), the ELCV converged to zero at intermediate to high energy (figure 2c). After the stress exposure, the original level of variability is recovered, allowing the same potential ‘hedging’ against further stress in the future (figure 2a) [55-58].

Next, we quantitatively investigated how energy level influences the variability of efflux pump gene expression pulses, by calculating the coefficient of variation of pulse lengths (PLCV) starting from each timestep in acrAB (figure 5). In the absence of stress, heterogeneity in pulse length is higher at lower energy (figure 5b,d), while high energy levels display relatively low PLCV <1, corresponding to lower variability when dynamics are faster and individual events less rare. As with expression level dynamics, when stimulated by a short stress period, PLCV decreases more dramatically at lower energy levels (figure 5b), as the population becomes more synchronized by the imposition of an external stress. For longer stress periods (figure 5d), differences between different energy levels are more dramatic. Pulse lengths remain highly variable for the low energy population, but PLCV decreases to zero at higher energy. Thus, higher energy availability allows more concerted and synchronous population-level behaviour in dealing with the stressor. Following stress exposure, PLCV behaviour returns to its pre-stress level for each energy level.

Escherichia coli and Salmonella efflux networks share functional similarity

We next compared the behaviour of the Salmonella network to identify similarities and differences between bacterial species. We first explored the capacity of the Salmonella network to transition between expression states. As with E. coli, the Salmonella network architecture supports the same accessible transitions (electronic supplementary material, figure S5). Clustering by the set of accessible transitions for each global state at a given energy level shows analogous features in both bacterial species, with transitions between states limited by energy (electronic supplementary material, figure S5).

Next, we asked how Salmonella’s acrAB efflux pump genes respond to a stressor, compared to the homologous efflux pump in E. coli, by analysing the population-level activity of acrAB in both species in response to a short and long stress period (figure 2; electronic supplementary material, figure S3 for E. coli and Salmonella respectively). Simulation results predict both species respond more dynamically to stress at higher energy, regardless of stress duration and, at each energy level, rapidly return to baseline conditions after the stress period (compare figure 2; electronic supplementary material, figure S3).

Salmonella displays several other analogous behaviours compared to E. coli. At maximal energy availability, the number of timesteps for Salmonella to return to the mean expression baseline after a short and long stress period is similar (electronic supplementary material, figure S3), as observed for E. coli (see §2.3). Other analogous behaviours include a zero acrAB ELCV at the maximum stress response to the long stress period (electronic supplementary material, figure S3(iii)), and the observation of heterogeneous pulses in network components (electronic supplementary material, figures S7–S8).

Simulations predict the stress-free acrAB ELCV is lower in Salmonella at low and intermediate energy, but greater at high energy, compared to E. coli (compare figure 2a; electronic supplementary material, figure S3(i)). This behaviour is a downstream consequence of the mean and variability in global activator expression in each species, which is modulated by the energy level, and the decoupled expression of Salmonella’s activator ramA and repressor ramR, compared with E. coli’s coupled marRA operon (figure 1b,c).

Together, this suggests the GRNs for both bacterial species in this study display behavioural similarities in the presence and absence of a stressor, but despite structural similarities in their architectures (figure 1) there are differences in the dynamics of efflux pump expression (figure 2; electronic supplementary material, figure S3).

Discussion

We have applied a Boolean modelling framework to coarse-grained E. coli and Salmonella efflux pump gene regulatory networks, to explore how efflux pump activity depends on the energy available to fuel the expression of efflux proteins and their regulators. We find quantitative support for our hypothesis that differences in available energy can contribute to cell heterogeneity in efflux pump expression, and dissect this heterogeneity in different mechanisms (see table 1 for a summary).

Table 1

Summary of energy availability effects on regulatory network behaviour. (The table displays a summary of our main results from applying the modified Boolean network framework to E. coli and Salmonella acrAB GRNs. The magnitude, under low, intermediate and high energy availability, of each key result is presented through a discretized scale: low (), medium () and high ().)

	low energy	intermediate energy	high energy
speed of network dynamics	∙	∙∙	∙∙∙
diversity of supported network behaviour in the absence of stress	∙∙∙	∙∙∙	∙∙∙
diversity of supported network behaviour in the presence of stress	∙	∙	∙
cell-to-cell variability in population-level acrAB expression in the absence of stress	∙∙∙	∙∙∙	∙∙∙
cell-to-cell variability in population-level acrAB expression in the presence of stress	∙∙	∙	∙
cell-to-cell variability in pulses of acrAB expression in the absence of stress	∙∙∙	∙∙	∙
cell-to-cell variability in pulses of acrAB expression in the presence of stress	∙∙∙	∙	∙

In the absence of cell-to-cell energy variability, our model suggests (supported by experiments) that substantial intrinsic variability exists in efflux gene expression. The imposition of a population-wide stressor reduces this intrinsic variability by partially synchronizing cell behaviour; stochastic effects after stress exposure recover the original variability. This intrinsic noise can potentially be acting to hedge the population against future stress, while maintaining cells that invest less resource in defence mechanisms and/or protect the lineage from diverse environmental conditions. When extrinsic differences in ATP availability between cells are additionally imposed, the population variability in efflux response will be amplified. Cells with lower energy levels will experience more variable expression of efflux components (amplifying intrinsic variability), slower responses to stressors, and hence slower removal of stressors from the environment. In future, experimental tests of this hypothesis could involve, for example, jointly tracking bacterial ATP levels and efflux gene expression in a culture prior to, during and post antibiotic stress using fluorescent reporters [11,12,16,59].

In our model, stress-free bacteria naturally express heterogeneous pulses of efflux genes (figure 3; electronic supplementary material, figure S7). This behaviour may be biologically valuable for several reasons. The products of acrAB have an extremely long half-life in E. coli [60], suggesting that the genes may only need to be transcribed in short pulses (which is predicted to be a less energetically demanding process compared to translation [61]) for the protein to be in an abundant supply. Dealing with external stressors is one role of the AcrAB-TolC efflux pump, but it also, for example, removes endogenous waste products [62], suggesting the pulse behaviour could aid physiological functions by expressing acrAB mRNA to be used in these processes. The pulse behaviour may also confer a natural resistance phenotype to a fraction of cells in stress-free environmental conditions at all times, as reported experimentally in S. enterica [13]. Therefore, the predicted pulsing behaviour may contribute to several essential intracellular functions.

Recent mathematical models have studied the mar network at the protein level, suggesting the motif is monostable, and supports a single deterministic steady-state in the absence of stress [63,64]. A theoretical study has also shown, at the protein level, the development of stochastic pulses in marRA products [64]. The cell-to-cell heterogeneity in mar expression that we predict agrees with ODE models of the mar operon—which account for the inherent stochasticity of biological systems—and experimental demonstrations of heterogeneous cell-to-cell steady-state mar expression within a population, in the presence and absence of stress [64-66] (figure 2; electronic supplementary material, figure S6). With the development of advanced experimental assays, future studies may resolve the behaviour at the transcript level to test our theoretical predictions.

Even for these relatively simple networks, our Boolean modelling approach has demonstrated general principles of energy influence without requiring characterization of a large parameter space. Developed with a view to scalability, this approach can be applied to many naturally occurring decision-making phenomena. Our previous work [14] has demonstrated the profound influence that energy variability can have on simple decision-making circuits; this GRN approach allows the expansion of this philosophy to a much wider range of architectures, including the broader regulatory context of these core networks (see the electronic supplementary material, S1–S2). As cell-to-cell variability in ATP is increasingly elucidated [18-20] we believe that our approach will help reveal the theoretical underpinnings of this important influence on regulatory dynamics.

It is important to note some limitations of our modelling approach. First, the coarse-grained nature of a Boolean model limits the precise quantitative agreement of individual observations with either more detailed models or biological reality. Instead, we focus on the relative magnitudes and directions of behaviour—for example, does X increase or decrease Y, and does Z do so more. Second, our models are deliberately simple. Model refinements, including allowing heterogeneity in the timescales [67] and ATP dependencies of different genes, could be invoked to allow a more detailed connection with biology, and doing so as more data becomes available is an important target for future work. Third, our models are not guaranteed to be unique. The edge weights in our models were manually chosen for compatibility with our training observations. As they encode logical interactions rather than physical quantities, there are ranges of values around our choices that will give networks that behave identically. The existence of non-overlapping sets of parameter values that would also describe and predict observations has not been addressed here. Parameter inference approaches including approximate Bayesian computation [68,69] could be used in future work to explore parameter space more fully.

Our hypothesized relationship between energy availability and efflux-pump activity, variability, and capacity to respond to antimicrobials has, as far as we are aware, not been previously considered. However, our hypothesis is not restricted to this AMR mechanism or motif. A connected phenomenon that may be mechanistically related is the formation of bacterial persisters in the presence of antibiotics [55,70,71]. Here, phenotypic variants display a transient tolerance to external stress [72,73], causing a leading factor in the recalcitrance of chronic infections [74-77]. An approach grounded in optimal control theory has found that the persister strategy is dependent upon stochastic fluctuations in both the proportion of persisters in the colony and the environment [78]. The results also predict that producing persisters in the absence of environmental volatility leads to a lower per capita growth rate, giving the population a possible evolutionary disadvantage. Increasing evidence has reported a mechanism linking persister cells and intracellular energy, or more specifically ATP [17,79-82]. One study has gone further, suggesting there may be a general low energy mechanism of persister formation in bacteria [83]. This relationship bears substantial similarities to our model where we predict high variability and noise in gene expression pulses in low energy cells, compared to intermediate and high energy cells (figures 2 and 5). High expression variability contributes to cell-to-cell variation within a clonal population, implying a different set of phenotypes may be expressed in the low energy sub-population. As energy, or more precisely ATP, increases, cell-to-cell variability decreases, implying a uniform phenotype displayed across the population. We would predict, in agreement with growing evidence, intracellular ATP is a factor in the mechanisms of persister generation, and future work will explore this using our approach when the governing GRNs can be identified.

Methods

Wiring diagrams

In the electronic supplementary material, figure S1, we present the wiring diagram form of each GRN (figure 1b,c). The wiring diagrams contain the regulatory interactions in each GRN (activatory or inhibitory) and edge weights, outputting a directed graph. A stress, such as antimicrobial exposure, is included through the ‘stress node’, S, within the wiring diagram. A stress input is modelled as repressing gene-gene regulatory edges, representing a ‘repressive’ modulation of transcription factors (MarR/RamR and AcrR). The capacity to expel substances by the efflux pump is integrated through a repressive regulation from acrAB to S.

Estimates of edge weights were gleaned from the relevant literature, integrating the natural behaviour of each network, and giving each regulatory actor the ability to operate. In E. coli, the marRA operon is separated into its individual components to allow each regulation of marRA to be captured, but are coupled in all iterative updates; the inability to model a node that both positively and negatively regulates itself is a limitation of the Boolean modelling framework. Experimental data have shown transcription of marA and acrA increase in marR mutant cells [84], and minimum inhibitory concentrations (MICs) of various antibiotics increase in an acrR mutant strain [44] (consistent with the limited repressor function of AcrR [41]). Together, these suggest, under normal conditions, edge weights of interactions from marR to marRA and acrR to acrAB should be dominant. Inactivation of the ramR gene upstream of ramA resulted in increased expression of ramA and the AcrAB efflux pump in Salmonella [85]. This requires the edge weight from ramR to ramA to outweigh ramA basal expression under normal conditions. Constitutive gene expression is captured using ‘ghost nodes’. These nodes are permanently ‘on’ in the network, supplying a constant positive influx to a node.

Interaction matrix J is constructed from the data contained within a wiring diagram. A weighted edge from node j to node i, is represented by element J within J. Each matrix element can take a positive, negative or zero value; edge weight sign corresponds to a negative, positive or absence of regulation respectively. Because edges are directed, J ≠ J, except for i = j or J = J = 0. Matrix J is implemented in numerical simulations of the Boolean model and has dimension N-by-N (N is the total number of nodes in a wiring diagram, including ghost and stress nodes).

Initial conditions

Two configurations of initial node states were considered in simulations of the E. coli and Salmonella wiring diagrams. These were named the ‘unstressed’ and ‘stressed’ initial condition, with the initial node states set as follows (tables 2 and 3).

Table 2

Escherichia coli node initial values. (Table defines the stressed and unstressed initial conditions for E. coli gene regulatory network components (figure 1b).)

	marR	marA	acrR	acrAB
unstressed	0	0	1	0
stressed	1	1	0	1

Table 3

Salmonella node initial values. (Table defines the stressed and unstressed initial conditions for Salmonella gene regulatory network components (figure 1c).)

	ramR	ramA	acrR	acrAB
unstressed	1	0	1	0
stressed	0	1	0	1

Modulation of transition rates allows energy variability to be captured in gene regulatory network models

Wiring diagrams, consisting of nodes and edges, schematically represent the interactions between elements within regulatory networks (as shown in the electronic supplementary material, figure S1). The edges often illustrate combined processes such as transcription, translation and phosphorylation. However, the information contained within edges in regulatory networks is coarse-grained, omitting a substantial amount of important biological detail. Several of the processes represented by these edges require ATP as an energy source [29,53], so there exists a core energy dependence in the dynamics modelled within each wiring diagram.

To address this and capture the fact that all the processes involved in gene expression depend on ATP, we constructed a modelling approach, modulating the strength of regulation from node X to node Y according to cellular energy levels. This reflects the fact that, for example, cells with low ATP will have a lower rate of transcription elongation per timestep [53]. We interpret this modulation as influencing the rate with which regulatory interactions can be applied from node X to node Y. The rationale for this picture is that if the expression of regulatory components is slowed by reduced energy levels, the corresponding biochemical interactions will occur at a lower rate.

In the absence of energy variability, our asynchronous Boolean modelling framework interprets regulatory interactions as stochastic events occurring with a characteristic mean rate of one per timestep per interaction. To include the influence of energy availability on the system, we allow this rate to be a function of energy level. To capture an order-of-magnitude range across cells, we use rates for a simulation timestep of 1 for ‘high’ energy (hence, all interactions are always realised), 0.5 for ‘intermediate’ energy, and 0.1 for ‘low’ energy. We thus modulate the probability of any given interaction being realized in a given simulation timestep by the level of available cellular energy.

Node update function

Nodes within the Boolean network are updated iteratively, using an update function. An update function defines the value of each node according to its regulatory nodes and remains fixed for all timesteps [86]. We consider a threshold update function to update the state of node i at timestep t = τ + 1, denoted σ(τ + 1). The argument of the update function consists of the regulatory nodes of node i, the states of these nodes at t = τ and the N-by-N interaction matrix J, where N is the total number of nodes in the wiring diagram. The summation value determines the state of node i at t = τ + 1.

During an active stress time frame, the state of stress node S (electronic supplementary material, figure S1) is imposed to be ‘on’ at the start of each timestep, regardless of the behaviour from the previous update. Once the final timestep within the stress period is reached, S is allowed to be switched ‘off’ from regulation within the network.

We consider the update function displayed in equation (4.1). Here, a node requires a non-zero sum of its inputs to change state from t = τ to t = τ + 1. In the case where the sum is zero, the node state remains unchanged:

Alternate threshold update functions are possible. We could, for example, consider the conditions that set a node as ‘on’ if there is an overall positive input, and ‘off’ otherwise. However, this would automatically remove a stress in the system by setting the stress node state to ‘off’ once reaching the endpoint of the stress time frame. This removes the capability of the efflux node to ‘expel’ the stress through feedback to the stress node S (electronic supplementary material, figure S1), behaviour which is not the focus of this study. By choosing update rule (4.1), it allows the efflux pump node to have the capacity to remove stress in the network, through regulatory interactions, depicting the behaviour observed in bacteria.

At each timestep, the binary vector containing all node states is the global, or network, state [87]. The set of all possible global states forms the state space of a network [86]. The state space for each network in this study is of size 2 for Salmonella, as we are interested in the dynamics of the M-component GRN network being modelled, and 2 for E. coli, owing to the coupling of marR and marA.

Node update methods

There are two methods for updating global states, synchronous and asynchronous. In the synchronous mode, all the nodes (including the stress node) are updated at every timestep, with the state of each node at t = τ + 1 depending on the states of its inputs at t = τ. The deterministic nature of synchronous Boolean models, combined with each model containing a finite state space (2 or 2), where M is the number of elements in the original GRN, guarantees that a synchronously updated time series simulation will eventually end at a steady-state attractor, or a cycle of repeating global states [86]. It also means the network will always reach the same state for a given initial condition and number of timesteps. We do not consider synchronous updating in this study owing to the artefacts that can arise from its imposition, and its inability to model stochastic processes. A more complex, but biologically realistic, protocol is the asynchronous mode, which gives a reasonable overview of the random dynamics in the cell and better captures transient behaviour. This method updates nodes of the system according to the last update of their regulator nodes, either from the previous or current iteration [88]. Previous studies have employed numerous different versions of asynchronous updating [27,67,89,90], including methods where all nodes are updated according to a random sequence, or one randomly selected node is updated at a timestep [91]. In this study, the asynchronous protocol updates η nodes per timestep, with η drawn from a Poisson distribution with a mean equal to the number of original elements in the regulatory network (M), plus the stress node. Additionally, each node can be updated more than once per timestep and the state of each node at t = τ + 1 depends on the states of its inputs at the most recent update. We used an asynchronous protocol because we believe it better captures the stochasticity in the cell.

Boolean analysis

Taken together, the wiring diagram, initial conditions, update functions and update method form the Boolean network model. Time-dependent numerical solutions of the presented Boolean model were generated in Python 3.8.2 using modules numpy, pandas and random. Numerical solutions were calculated from prescribed initial conditions covering the stressed and unstressed state of each network. Time-series figures, displaying the expression profiles of each node, were constructed from the mean of node expression and standard deviation per timestep, from 104 simulations. Visualizations of ‘single-cell’ expression profiles, and pulse statistics, were constructed from 20 and 5000 randomly selected simulations, respectively. Transition matrices were constructed from the simulation data and visualized through a heatmap structure, from 107 simulations. Each entry of the transition matrix is a non-negative real number representing the probability of transitioning between two network states from timestep t = τ to t = τ + 1. Clustering of global states was performed using the same data, with an Euclidean distance matrix from modules scipy.cluster and scipy.cluster.hierarchy. All scripts used in this study are openly accessible through https://github.com/StochasticBiology/boolean-efflux.git.