A mathematical model of flagellar gene regulation and construction in Salmonella enterica

Millions of people worldwide develop foodborne illnesses caused by Salmonella enterica (S. enterica) every year. The pathogenesis of S. enterica depends on flagella, which are appendages that the bacteria use to move through the environment. Interestingly, populations of genetically identical bacteria exhibit heterogeneity in the number of flagella. To understand this heterogeneity and the regulation of flagella quantity, we propose a mathematical model that connects the flagellar gene regulatory network to flagellar construction. A regulatory network involving more than 60 genes controls flagellar assembly. The most important member of the network is the master operon, flhDC, which encodes the FlhD4C2 protein. FlhD4C2 controls the construction of flagella by initiating the production of hook basal bodies (HBBs), protein structures that anchor the flagella to the bacterium. By connecting a model of FlhD4C2 regulation to a model of HBB construction, we investigate the roles of various feedback mechanisms. Analysis of our model suggests that a combination of regulatory mechanisms at the protein and transcriptional levels induce bistable FlhD4C2 levels and heterogeneous numbers of flagella. Also, the balance of regulatory mechanisms that become active following HBB construction is sufficient to provide a counting mechanism for controlling the total number of flagella produced.


Introduction
Millions of people worldwide develop foodborne illnesses caused by Salmonella enterica (S. enterica) every year [1]. The pathogenesis of S. enterica depends on flagella, which are appendages that the bacteria use to move through the environment. Flagella are located peritrichously on the bacterium's surface and generally number between 6 and 10 per cell [2]. However, not all bacteria produce flagella. Partridge and Harshey describes the distributions of the number of flagella per bacterium under different growth conditions [3]. It was observed under each growth condition that some bacteria did not produce any flagella, while the number of flagella on bacteria that produced them roughly follows a normal distribution, where the mean depended on the growth condition. This phenotypic heterogeneity enables the evasion of the host immune system during acute infection [4]. It is not yet clear what mechanisms determine flagella heterogeneity and how the bacteria regulate the number of flagella produced.
Structurally, a flagellum is divided into three parts: basal body, hook, and filament. The basal body anchors the flagellum to the cell membrane, the hook is a flexible joint that links the rigid basal body with the filament, and the filament is a rigid helical structure which rotates to drive the cell forward. Assembly of functional flagella involves the temporally coordinated expression of more than 60 genes [5]. In particular, these genes are divided into three classes: Class I, Class II, and Class III depending on the timing of their activation. The flhDC operon is the only member of Class I. This operon encodes the master regulator proteins FlhD and FlhC, which form a transcriptional activation complex FlhD 4 C 2 [6]. The FlhD 4 C 2 complex is responsible for activation of Class II genes. Class II promoters control the expression of all the proteins required for the assembly of a functional hook-basal body (HBB) structure [7]. The HBB includes the flagellar type III secretion system, which exports flagellar proteins from the cytoplasm through the growing structure during assembly [8]. Once the HBB assembly is complete, FlhD 4 C 2 activates expression of flgM, fliA, and other Class III genes. fliA encodes the flagella-specific sigma factor, σ 28 , which is essential for Class III promoter activation. flgM encodes an anti-sigma factor FlgM. Prior to HBB completion FlgM and σ 28 form a complex that inhibits the activity of σ 28 . After the HBB is completed, the specificity of the secretion apparatus switches, so filament subunits and FlgM are secreted from the cell. σ 28 facilitates the secretion of FlgM through the HBB, thereby acting as the FlgM type III secretion chaperone [9]. The secretion of FlgM frees σ 28 inside the cell, which leads to the activation of Class III genes. Class III genes control the assembly of the filament, the chemotaxis system, and motor proteins.
While the components of the gene network have been identified, understanding the additional details and dynamics of the system requires further study. In particular, genetically identical bacteria grown in the same environment exhibit bistable expression of multiple flagellar genes [3]. Expression of Class II and Class III genes, flhB and fliC, respectively, is bimodal, and this bimodality is likely due to bistability upstream in the network [10][11][12][13]. Other factors also impact the flagella gene network, including nutrient availability, which enhances the expression of flagellar and motility genes [14]. Nutrients repress RflP, formerly YdiV, expression through the protein CsrA [15]. Low nutrient levels correspond to high RflP expression, and high nutrient levels correspond to low RflP expression. The tuning of RflP expression influences the flagellar gene network because RflP promotes the degradation of FlhD 4 C 2 through the protease ClpXP [16].
A series of mathematical models have been proposed to explore the dynamics of the flagellar gene network [17][18][19]. While these models investigated the roles of some of the feedback mechanisms in the network, they fail to describe the heterogeneity observed in bacterial populations, since all bacteria are predicted to grow flagella. Further, the number of HBBs was modeled as a continuous variable, although HBBs are discrete entities.
Meanwhile, Koirala et al. presented a model that attempted to explain how nutrients tune flagellar gene expression dynamics in Salmonella enterica [13]. Using experimental observations and a deterministic model, the authors demonstrated that environmental nutrient levels can regulate RflP expression and that this expression can lead to bistable expression of Class II and Class III genes in the bacteria [13]. The model does not include HBB construction, or, consequently, the additional feedback mechanisms that are activated following HBB completion.
To investigate the heterogeneous flagella quantities within populations of clonal bacteria, we propose two mathematical models of the flagellar gene regulatory network. We first consider a sub-system of the network to focus on FlhD 4 C 2 production and degradation. We use bifurcation analysis of an ordinary differential equation model to explore the role of flhDC and rflP expression rates on overall FlhD 4 C 2 concentration. We find bistability in FlhD 4 C 2 concentration as a function of both expression rates. We then extend the model to consider HBB construction and additional regulatory mechanisms. We retain the bistability observed in the subsystem model and investigate factors that impact the bistable regime. We also propose a mechanism for how S. enterica controls the number of flagella produced via FliT:FliD dynamics and Class III regulation of FliZ.

Sub-system model
We first consider a model of the key components that regulate flhDC expression and FlhD 4 C 2 protein concentrations. The model includes the FliZ-RflP-FlhD 4 C 2 feedback loop and RflM activity. Expression of fliZ is up-regulated by FlhD 4 C 2 and FliZ represses expression of rflP [20]. RflP inhibits the activity of FlhD 4 C 2 by binding to the FlhD subunit and targeting FlhD 4 C 2 for proteolysis via ClpXP [16]. Thus, FliZ is involved in a negative-negative feedback loop through RflP on the FlhD 4 C 2 protein. FlhD 4 C 2 up-regulates expression of rflM, and RflM represses the transcription of FlhD 4 C 2 [21][22][23]. Fig 1 summarizes the mechanisms included in this model.
We assume that on the time scale of flagellar gene regulation and construction, the bacterial cytoplasm is well-mixed. The diffusion coefficients for proteins in the bacterial cytoplasm range from 1-10 μm 2 s −1 , and since S. enterica are 2-5 microns long, proteins can diffuse across the bacterium on the order of seconds [24,25]. Since HBBs are completed 30 minutes after induction of flhDC expression and nascent flagellar filaments become visible at 45 minutes, spatial non-uniformity of the proteins will not be significant [2].
Using the well-mixed assumption, we set We model transcription and translation processes as Michaelis-Menten functions and protein binding, unbinding, and degradation processes linearly. FlhD 4 C 2 is produced at a basal rate of ρ 1 , and the production rate is divided by the concentration of RflM, scaled by η 1 . FlhD 4 C 2 binds to RflP at rate κ 1 . FliZ is produced as a function of FlhD 4 C 2 concentration. RflP is produced at a basal rate of ρ 2 , and the production is repressed by FliZ. External nutrient availability also represses rflP expression, so ρ 2 is a proxy for nutrient availability [15]. RflP forms a complex with FlhD 4 C 2 and is recycled at rate δ 1 where FlhD 4 C 2 is degraded by ClpXP. We assume that the complex does not dissociate spontaneously. RflM is produced as a function of FlhD 4 C 2 concentration. All components, including the RflP:FlhD 4 C 2 complex, are assumed to degrade at a small, basal rate, γ i , where i = 1, 2, . . ., 5. Table 1 lists the model parameters with descriptions and values.
Bifurcation analysis of sub-system model. Since some bacteria produce flagella and others do not, we suspect that FlhD 4 C 2 concentration is bistable. We use resultant analysis of the model to determine what parameter regimes support bistability [26]. To apply resultant analysis, we take the model to steady state, reduce the system to a single polynomial in one variable, and perform a rescaling to reduce the number of parameters. We then solve the equation for a rescaled parameter of interest. To find double roots of the first derivative, we then take the resultant of the first and second derivatives of the right side of the equation. The curve in parameter space on which the first and second derivatives have a common root corresponds to the separation between monotonic and triphasic behavior. We first take the model to steady state and find : We reduce the number of parameters by taking the rescaling x = k 1 w and letting K ¼ k 2 k 1 ,

PLOS COMPUTATIONAL BIOLOGY
We solve for P 2 to find : We next use resultant analysis of f 0 (w) and f @(w) to determine when f(w) is monotone and when it is triphasic. The triphasic regions correspond to bistability, and the monotone regions correspond to monostability. We first consider the case in which η 1 = 0 (A 2 = 0), which means that RlfM does not impact flhDC expression. f(w) simplifies to We next calculate the resultant of the numerators of the first and second derivatives of f A 2 =0(w) using Maple (2019; Maplesoft, Waterloo, ON), and we find The function R f,A 2 =0 = 0 separates parameter space into regions for which f A 2 =0(w) is monotone increasing (R f,A 2 =0 < 0) and triphasic (R f,A2 =0 > 0). In the triphasic region, there is bistability, which will be verified later. The solid curve in Fig 2 shows the zero level contour for R f, is negative to the left of the curve and is positive to the right of the curve. f(w) is monotone for values of P 1 and A 1 that fall to the left of the curve, whereas for values that fall to the right of the curve, f(w) is triphasic, and there is bistability. Thus, bistability is possible under constitutive expression of flhDC (i.e. when A 2 = 0).
We note that in the case that D = 0, the resultant reduces to For positive values of A 1 and P 1 , R f,A 2 =0,D=0 is always negative, and either A 1 = 0 or P 1 = 0 satisfy R f = 0. This indicates that we must have D > 0 to have bistability. In terms of the original parameters, this means that γ 4 > 0 is required to have bistability. We next consider the case in which η 1 6 ¼ 0 (A 2 6 ¼ 0), meaning RflM represses the expression of flhDC. We calculate the resultant of the numerators of the first and second derivatives of f (w), but it is too long to show here. The dashed curve in Fig 2 is the zero level contour for R f when RflM activity is included. We see that the R f = 0 curve moves up and to the right as A 2 increases, which decreases the size of the parameter space in which there is bistability. In particular, the repressive activity of RflM requires an increase in the baseline production rates of FlhD 4 C 2 (P 1 ) and FliZ (A 1 ) to have bistability.
We next use the parameters in Table 1 and XPPAUT (8.0) to visualize the bistable regime identified by the resultant analysis. Fig 3 shows the steady state concentration of FlhD 4 C 2 plotted as a function of the parameters ρ 1 and ρ 2 , demonstrating bistability for certain regions of  Depending on the values of ρ 1 and ρ 2 , the steady state value of FlhD 4 C 2 concentration has one solution or three solutions. As shown in Fig 3A, the FlhD 4 C 2 concentration is low for small ρ 1 values. This means that the bacteria do not produce enough FlhD 4 C 2 to initiate flagella production. With ρ 1 increased sufficiently, a saddle node bifurcation occurs, and bistability arises. The lower solution represents the case in which the bacterium produces no flagella, and the upper solution represents the case in which the bacterium produces flagella. Increasing ρ 1 further results in another saddle node bifurcation, and a return to a single steady state value for FlhD 4 C 2 concentration. The expectation is that for high flhDC production rates, the bacterium will have a high FlhD 4 C 2 concentration, and, consequently, produce flagella. Saddle node bifurcations are also observed for ρ 2 , as shown in Fig 3B. In this case, small ρ 2 values result in a high steady state FlhD 4 C 2 concentration, intermediate ρ 2 values result in bistability, and large ρ 2 values result in a low FlhD 4 C 2 concentration.
Many factors control flhDC expression, so it is reasonable to expect that clonal bacteria may have different values for ρ 1 . In this case, the variation in ρ 1 values could contribute to heterogeneity in flagella number, assuming the bacteria on the lower branch have no flagella and bacteria on the upper branch have enough FlhD 4 C 2 to produce flagella. Also, since nutrient availability represses rflP expression, small variations in nutrient availability could result in differing values for ρ 2 . Variability in ρ 2 would result in different concentrations of FlhD 4 C 2 , and, consequently, the presence or absence of flagella. We note that, as with the case A 1 = 0, γ 4 > 0 is required for bistability.
We also consider the η 1 = 0 case, which corresponds to a knockout of RflM or constitutive expression of flhDC via an inducible promoter. Under this condition, the FlhD 4 C 2 production term reduces to ρ 1 . Fig 4 shows the steady state concentration of FlhD 4 C 2 is bistable for certain regions of ρ 1 and ρ 2 , but the FlhD 4 C 2 production term is not multivalued.
In this case, constitutive expression of flhDC can lead to bistable FlhD 4 C 2 levels. In particular, the η 1 = 0 case shows bistability in FlhD 4 C 2 even when the production of flhDC is not multivalued. Therefore, we suggest that RflM knockout cells should still exhibit bistability in the protein level of FlhD 4 C 2 . In anticipation of a need to understand the role of FlhD 4 C 2 degradation, we next investigate the impact of the basal degradation rate of FlhD 4 C 2 on FlhD 4 C 2 concentration. Fig 5A shows bistability in FlhD 4 C 2 concentration as a function of the parameter γ 1 .
A saddle node bifurcation is observed for FlhD 4 C 2 concentration as a function of γ 1 . Small γ 1 values result in bistable steady state FlhD 4 C 2 concentration, and increasing γ 1 sufficiently results in a low FlhD 4 C 2 concentration. We next investigate the impact of γ 1 on the bistability observed in FlhD 4 C 2 as a function of ρ 2 . Fig 5B shows a cusp bifurcation in the γ 1 -ρ 2 plane. This demonstrates that the bistability previously observed in FlhD 4 C 2 as a function of ρ 2 is sensitive to γ 1 .

Extended model
To explore how the bacteria regulate the total number of flagella constructed, we extend the model to consider the construction of HBBs and additional regulatory mechanisms that become active following HBB completion. In particular, FlhD 4 C 2 up-regulates the fliDST operon, which encodes the FliT and FliD proteins. Before HBB completion, FliT and FliD are bound together, which inhibits the activity of FliT. Following HBB completion, FliD is secreted through the HBB and is assembled as the cap of the flagellum. The secretion of FliD frees FliT in the cell to bind to FlhD 4 C 2 and target FlhD 4 C 2 for proteolysis via ClpXP [27][28][29]. FliT is recycled during FlhD 4 C 2 degradation [27]. FlhD 4 C 2 also up-regulates production of FlgM and σ 28 , which form a complex, and complete HBBs secrete FlgM, resulting in free σ 28 in the bacterium. σ 28 up-regulates the expression of Class III genes, including fliC. FliC is also secreted through the HBB to become the flagellar filament. σ 28 also up-regulates fliZ, which means fliZ is under both Class II and Class III regulation. The gene network for the full model is shown in Fig 6. For simplicity, we assume that FliT and FliD are produced as a complex (FliT:FliD), rather than modeling the production of each protein and formation of the FliT:FliD complex. Similarly, σ 28 and FlgM are assumed to be produced as the complex σ 28 :FlgM.
As before, we assume that the system is well-mixed, and we model transcription and translation processes as Michaelis-Menten functions and protein binding, unbinding, and degradation processes linearly. J is the number of complete HBBs that are capable of secreting FlgM, FliD, and FliC, and we assume that the secretion of the three components, FlgM, FliD, and FliC, through the HBBs is competitive and saturating. We assume that these three components have the same secretion affinity, σ. As in the sub-system model, we assume that every compound, including the complexes, degrade at some basal rate.
Model reduction. The model involves both protein level dynamics, including binding and degradation, and transcription and translation. We assume that the protein level dynamics occur on a faster time scale than transcription and translation, so we take a quasi-steady-state approximation for the protein level dynamics. We assume that c 1 = [RflP:FlhD 4 C 2 ], y = [RflP], and c 2 = [FliT:FlhD 4 C 2 ] are fast equilibrating, so that they are at equilibrium. We therefore assume that We also assume that We now let u = v + c 2 , the sum of [FliT] and [FliT:FlhD 4 C 2 ], so we have We then assume that Using the quasi-steady-state approximation, the model no longer includes equations for c 1 = [RflP:FlhD 4 C 2 ], y = [RflP], and c 2 = [FliT:FlhD 4 C 2 ]. The reduced model with competitive, saturating secretion, is HBB model. Since little is understood about HBB initiation and construction, we propose a simple model that assumes the environment within the bacteria is well-mixed. We assume that a general HBB protein is produced at a rate that depends on the concentration of FlhD 4 C 2 in the bacterium. Mouslim and Hughes showed that flagella construction requires a threshold of flhDC expression, and, consequently, Class II expression [21]. This finding suggests that the concentration of HBB proteins must reach some threshold before initiation of HBB construction. We assume that once the concentration of the HBB protein reaches a threshold, nucleation of a HBB occurs, and the free HBB protein concentration drops to zero. The newly nucleated HBB grows by recruiting free HBB proteins until it reaches a threshold level, at which point the HBB is considered complete [30,31]. Complete HBBs do not recruit free HBB proteins but secrete FlgM, FliD, and FliC, which results in free σ 28 and FliT within the bacterium. Depending on the HBB protein production rate and recruitment rate of incomplete HBBs, the bacteria are capable of growing multiple HBBs simultaneously. Let h be the concentration of free HBB proteins in the bacterium, and b i be the number of HBB proteins in the i th incomplete HBB. We then have where α 5 is the production rate of free HBB proteins, k 7 is the Michaelis constant for HBB protein expression, r is the recruitment rate of free HBB proteins to incomplete HBBs, k h is the Michaelis constant for HBB protein recruitment, γ 12 is the degradation rate of free HBB proteins, and g i is zero or one, depending on whether the i th HBB is growing. K is the number of incomplete HBBs, K = ∑ i g i . The nucleation and completion thresholds are h nuc and h max , respectively. Simulation methods. We use the parameters in Tables 1 and 2 to simulate the model in MATLAB (2019a; MathWorks, Natick, MA).
We use all zero initial conditions, unless specified otherwise, and set K = 0 and J = 0. We then numerically integrate Eqs 17-26 until h reaches the nucleation threshold, h nuc , which indicates the first nucleation event. Using the state of the system just prior to nucleation and setting h ! 0, b 1 ! h nuc , g 1 ! 1, K ! 1, we numerically integrate the system again until h = h nuc or b i = h max . If h = h nuc , another HBB is nucleated, meaning h ! 0, b j ! h nuc , g j ! 1, K ! K + 1, where the jth HBB is newly nucleated. If b i = h max , we complete construction of the ith HBB, meaning b i ! 0, g i ! 0, K ! K − 1, J ! J + 1. The process continues until the system reaches steady-state.

Results
The bistability observed in the sub-system model is maintained in the extended model. Depending on the initial conditions, HBBs are either produced or not. Since we do not explicitly model the construction of the filament, the number of HBBs is a proxy for the number of flagella. We first assume an initial condition of zero for all components of the gene network. Due to basal production of FlhD 4 C 2 , there is some expression of Class II genes (Fig 7A), but an insufficient concentration of FlhD 4 C 2 is accumulated to initiate HBB production (Fig 7B).
In contrast, when we assume a small initial concentration of FliZ and FlhD 4 C 2 , FlhD 4 C 2 accumulates sufficiently to initiate HBB production (Fig 8A), resulting in the production of seven HBBs (Fig 8B). This value falls within the experimentally observed range of 6-10 per bacterium [2].
The differing outcomes confirm that the FlhD 4 C 2 -FliZ-RflP feedback loop, as predicted by the sub-system model, determines if HBBs are produced. We then compute the separatrix in [FlhD 4 C 2 ]-[FliZ] initial condition space to determine what initial conditions are necessary to stimulate HBB production (solid curve in Fig 9). We next consider the impact of nutrient availability on HBB production by changing the RflP production rate, ρ 2 . We use ρ 2 as a proxy for environmental nutrient conditions in which increased nutrient availability decreases ρ 2 and decreased nutrient availability increases ρ 2 . A 20% decrease in ρ 2 results in the production of eight HBBs (Fig 10A), which is one more in comparison to the baseline case, whereas increasing ρ 2 by 20% results in no HBB production ( Fig 10B).
Since changing ρ 2 results in a different number of HBBs produced, we know that ρ 2 impacts the separatrix between no HBB production and HBB production. Increasing ρ 2 by 20% shifts the curve to require larger initial concentrations of FlhD 4 C 2 and FliZ to produce HBBs (dashed curve in Fig 9). A 20% decrease in ρ 2 lowers the curve out of the positive parameter quadrant, so that HBBs are produced under any biologically reasonable initial conditions. We next investigate the impact of the FliT:FliD production rate, α 3 , on HBB production. Decreasing α 3 results in more HBBs (Fig 11A), while increasing α 3 results in fewer HBBs ( Fig  11B), in comparison to the case with baseline α 3 .

PLOS COMPUTATIONAL BIOLOGY
The change in HBB production in response to changes in α 3 suggests that the FliT:FliD concentration, and, consequently, the free FliT concentration, is a key factor in controlling the termination of HBB production. This is exemplified by setting α 3 to zero (Fig 12). In this case, HBB construction continues indefinitely, which means the repressive activity of FliT on FlhD 4 C 2 is required to stop HBB production.
To further consider the impact of FliT on FlhD 4 C 2 concentration, we recall the sub-system model (Eqs 1-5). Even though FliT is not explicitly included in the model, increasing the basal rate of degradation of FlhD 4 C 2 , γ 1 , can be interpreted as a proxy for the repressive activity of FliT on FlhD 4 C 2 . As described earlier, Fig 5A shows that FlhD 4 C 2 concentration is bistable as a function of γ 1 . This means that as HBBs are completed and free FliT levels increase, γ 1 increases, which results in a saddle-node bifurcation and collapse of the system to a single, low steady state value. Similarly, Fig 5B shows that the bistability of FlhD 4 C 2 concentration as a function of ρ 2 is sensitive to γ 1 . As γ 1 increases, the system moves outside the region in parameter space that supports bistability. These results suggest that there is a threshold level of free FliT required to shut-down HBB production. Finally, we investigate the role of Class III regulation of FliZ. In particular, we set the Class III production rate of FliZ, β 1 , to zero. Without Class III regulation of FliZ, the FliZ concentration drops dramatically following the first HBB completion (Fig 13A), and, as a result, only five HBBs are produced (Fig 13B). The reduction in the number of HBBs produced, in comparison to the baseline case, is due to the lack of HBB nucleation following the first HBB completion event. Following HBB completion, FliT activity decreases the concentration of FlhD 4 C 2 . This slows the production of FliZ, and, without Class III regulation, the FliZ concentration drops dramatically. The drop in FliZ concentration further contributes to the decrease in FlhD 4 C 2 concentration. FlhD 4 C 2 then drops to a level that does not support further HBB nucleation. Together, these results suggest that the FlhD 4 C 2 -FliZ-RflP feedback loop determines whether flagella production will be initiated, and the balance of FliT:FliD dynamics and Class III regulation of FliZ control the total number of flagella produced.

Discussion
S. enterica is a common foodborne pathogen that impacts millions of people each year, and the flagellum is an essential component of S. enterica pathogenesis. Motility is critical in enabling bacteria to reach the site of infection and establish disease. Therefore, understanding how S. enterica regulate and construct flagella is of interest for disease prevention and treatment. An interesting feature of flagella regulation in S. enterica is the presence of a heterogeneous quantity of flagella per bacterium in a clonal population [3].
Here, we developed mathematical models of the gene network that regulates flagella construction to improve our understanding of flagella heterogeneity. Analysis of a sub-system model (Eqs 1-5) of the network identified the possibility of bistable FlhD 4 C 2 concentrations. The bistability was determined by the production rates of FlhD 4 C 2 and RflP, and could explain why some bacteria within a clonal population grow flagella while others do not. This finding could be tested by measuring the expression of flhDC in a population of rflM knockout bacteria.
We then extended the sub-system model to include the construction of HBBs and additional regulatory mechanisms. Simulations of this model (Eqs 17-26) suggest that FliT activity and Class III regulation of FliZ provide a counting mechanism by which the bacteria can control how many flagella are produced. Class III production of FliZ prevents underproduction of flagella, and FliT prevents overproduction. Hence, FlhD 4 C 2 is regulated by multiple negative feedback loops to ensure that protein levels in the bacteria can be tightly controlled, and both under-and overproduction of flagella is avoided.
The flhDC operon is primarily transcribed from two promoters, P1 flhDC and P5 flhDC , and many additional factors are shown to positively and negatively regulate flhDC expression, including HilD, RcsDBC, RtsB, and LrhA [21,[32][33][34]. Our model includes the main components involved in the expression of flhDC from P1 flhDC . Expression from the P5 flhDC , which is regulated via a FliZ-HilD positive feedback loop [21], plays a role in virulence rather than motility. Therefore, we focus on expression from P1 flhDC . Future models could include expression from P5 flhDC and additional regulatory mechanisms.
The limited availability of data meant that parameter values were selected manually, rather than measured or estimated. This study investigated the possible behaviors of the gene network, and, thus, does not rely on accurate parameter estimates. As additional data becomes available, the parameter values could be updated. Further, due to a limited understanding of HBB construction, we have proposed a simple model of the process. As we learn more about HBB construction, the model could also be updated.
Our results explain the bimodality of flagella count data and how the bacteria regulate the total number of flagella produced. However, since the model is deterministic, it cannot capture the spread of 6-10 flagella in the count data. Future models could incorporate sources of stochasticity, such as cell division, to understand the variability in the number of flagella. Since the system is assumed to be well-mixed but HBBs are discrete entities, cell division could be modeled by equally dividing the proteins in the parent bacterium and randomly splitting the HBBs between the two daughter cells. The randomness of HBB inheritance could be sufficient to generate the observed flagella distributions. This paper has proposed two mathematical models of the flagellar gene network in S. enterica. We predict that the FlhD 4 C 2 -FliZ-RflP feedback loop determines whether flagella are produced, and the combination of FliT activity and Class III regulation of FliZ determines the total number of flagella produced.