Mathematical modelling reveals properties of TcdC required for it to be a negative regulator of toxin production in Clostridium difficile

The role of the protein TcdC in pathogenicity of the bacterium Clostridium difficile is currently unclear: conflicting reports suggest it is either a negative regulator of toxin production or, on the other hand, has no effect on virulence at all. We exploit a theoretical approach by taking what is known about the network of proteins surrounding toxin production by C. difficile and translating this into a mathematical model. From there it is possible to investigate a range of possible interactions (using numerical and asymptotic analyses), identifying properties of TcdC which would make it a realistic candidate as a toxin inhibitor. Our findings imply that if TcdC is really an inhibitor of toxin production then TcdC production should be at least as fast as that of the protein TcdR and TcdC should remain in the cells throughout growth. These are experimentally-testable hypotheses and are equally applicable to alternative candidates for toxin production inhibition.


Clostridium difficile and its pathogenicity
Clostridium difficile is a leading cause of hospital-associated infections in Europe, the United States and beyond. Infection generally occurs as a result of antibiotic use by a patient for a preceding infection: treatment causes disruption of the intestinal microflora, allowing colonisation by C. difficile in the gut (Rupnik et al. 2009). This Gram-positive organism has a range of effects on the host, from mild diarrhoea to pseudomembranous colitis (severe inflammation in areas of the colon) and thousands of cases every year in the United Kingdom alone are ultimately fatal (Statistical Bulletin 2011). Additionally, in surviving patients, infections often subsequently recur for months or years, even when treated with antibiotics. Worryingly, community outbreaks have begun to arise: C. difficile infection is no longer restricted to health care settings and a rise in what some believe are hypervirulent strains has been noted (Rupnik et al. 2009). Furthermore, antibiotic resistant strains are increasingly prevalent (French 2010). Evidently, as is the case with a growing number of bacterial species, novel drugs must be sought urgently.
One crucial obstacle to the development of targeted anti-C. difficile treatments is a general lack of knowledge about the virulence mechanisms used by this bacterium. Our understanding of the gene regulation networks governing pathogenesis in Clostridium species falls somewhat behind that of similar Bacillus species due to historical difficulties in genetically manipulating this bacterium. However, as more clostridial manipulation techniques are developed (Heap et al. 2007O'Connor et al. 2006), one may expect to see an increase in the level of detail uncovered.

The PaLoc
In this study, we focus our attention on the subset of the overall network controlling C. difficile virulence that is most directly responsible for determining toxin production: the Pathogenicity Locus (PaLoc) comprising the five genes tcdA, tcdB, tcdC, tcdE and tcdR (see Fig. 1 and Hundsberger et al. 1997, for example).
The gene products TcdA and TcdB are the precursors for Toxins A and B respectively, the two principal toxins involved in C. difficile infection (Lyras et al. 2009;Kuehne et al. 2010;Carter et al. 2010). 1 Transcription of tcdA and tcdB is positively controlled by TcdR, which also auto-regulates (Mani and Dupuy 2001;Mani et al. 2002). Conflicting reports exist surrounding the role of TcdE: it may (Tan et al. 2001;Govind and Dupuy 2012) or may not (Olling et al. 2012) be required for secretion of the toxins into the extracellular environment. Similarly, as we discuss below, recent publications have poured controversy on the role of TcdC. Before precise gene manipulation techniques were available for C. difficile, studies examining the role of TcdC had to be performed in a surrogate host, namely C. perfringens. Such work generated the hypothesis that TcdC negatively regulates toxin production via interference with the positive action of TcdR upon transcription of tcdA and tcdB (Dupuy et al. 2008;Matamouros et al. 2007). This was further supported with the discovery that ('hypervirulent') strains of C. difficile, believed to produce significantly higher toxin levels than previously seen strains, had a naturally occurring mutation in their tcdC gene, suggesting that these strains contain a non-functional TcdC protein that cannot temper toxin levels (Carter et al. 2011). However, two more recent studies in C. difficile itself (Cartman et al. 2012;Bakker et al. 2012) found that TcdC did not affect toxin levels in either the standard or 'hypervirulent' strains, drawing into question the role TcdC may have in regulating the amount of Toxin A and Toxin B produced. In addition, a mechanism of interaction between TcdC and TcdR has proved tricky to identify (van Leeuwen et al. 2013). Thus the role of TcdC, which could be crucial in C. difficile infection and in finding novel drug targets, remains unclear.
Mathematical modelling is increasingly being used as a means by which to glean information from gene regulation networks. Translating what is known about the system into a set of mathematical equations, and investigating various alternative assumptions about the more ambiguous elements, allows predictions to be made about the role of particular genes or proteins in the network and the effects of altering the expression levels of these genes and post-translational modifications to the system. We exploit this approach here to investigate the role of TcdC in the PaLoc and hence on toxin production by C. difficile.

Model development
A standard procedure to build deterministic ordinary-differential-equation models of gene regulation networks in either prokaryotic or eukaryotic systems is to employ conventional mass-action kinetic theory (Karlebach and Shamir 2008, for example): the rate of a reaction is proportional to the product of the concentrations involved in the reaction. Variable and parameter definitions are provided in Tables 1 and 2.
The following information and assumptions are used here to build (and in some cases simplify) a model of the PaLoc.
-The cells follow logistic growth (at rate r with carrying capacity K ) meaning that after an initial burst of exponential growth, the bacteria will settle into a stationary phase. This is appropriate for cells grown in batch culture over our time period of interest (roughly 24 h say). -We assume TcdE is responsible for toxin secretion causing cell death in the process as per Govind and Dupuy (2012). Toxins are released following Michaelis-Menten kinetics at maximum rate μ with Michaelis constant η.
-Protein levels are proportional to their respective mRNA levels in terms of their dynamics (evidence suggests that this is the case at least for TcdA, TcdB and TcdC- Dupuy and Sonenshein (1998), Govind et al. (2006)). This allows us to include mRNA in the model implicitly, thus reducing the number of equations. -Since we focus on qualitative implications, we consider only TcdA for simplicity (rather than both TcdA and TcdB), since they are believed to have the same dynamics but possibly at different quantitative levels (Dupuy and Sonenshein 1998). -TcdA and TcdR are produced initially at a basal level during early growth (at rates c l A and c l R respectively) and at a higher rate induced by TcdR (c h A and c h R ) when TcdR is present in sufficient levels later in growth (Dupuy and Sonenshein 1998;Mani et al. 2002). TcdE, on the other hand, is produced basally at rate c E throughout, while tcdC transcription appears during exponential growth and is Toxin secretion saturation constant Ratio of the separation to binding of TcdR from promoter site of tcdx suppressed upon entry to stationary phase (Hundsberger et al. 1997;Govind et al. 2006;Carter et al. 2011). Indeed the last of these cited studies indicates that TcdC is present in the cells for the first 12 h of growth; this corresponds to the onset of stationary phase. We mark the transition between exponential and stationary phase (i.e. as the number of cells begins to settle off to the carrying capacity K ), at time t * , with a step function F: but investigate also the possibility that TcdC is produced throughout growth. -All proteins are subject to natural degradation (at rate λ x for variable x) and dilution as a result of cell division during growth.
-If TcdC does indeed interfere with TcdR dynamics, the mechanism by which this would occur is unclear. In Carter et al. (2011), evidence suggests that TcdC operates at the transcriptional level, whereas van Leeuwen et al. (2013) indicate that TcdC is actually unlikely to bind the tcdR gene. Unless otherwise stated, we assume a functional TcdC protein inhibits toxin production by forming a heterodimer with TcdR (at rate β) thus blocking the positive action of TcdR on its own gene and on the toxin genes. This process whereby a protein is inactivated by heterodimerisation is known as molecular titration and is a common mechanism in regulatory networks to switch on or off a relevant response. See, for example, Buchler and Louis (2008) for a mathematical model demonstrating this. To model a strain with a nonfunctional TcdC protein, we simply use β = 0. We assume binding is irreversible to examine functional TcdC at maximum efficiency. We include some simulations to illustrate that our conclusions apply also to the case where TcdC instead blocks tcdR transcription but limit this to numerical solutions for brevity.
In the cases where TcdC acts at the transcriptional level, (6) and (7) are replaced respectively by Given the lack of relevant quantitative data for the PaLoc, rather than attempting to obtain precise estimates for the parameters, we exploit asymptotic methods which require only the relative sizes of groups of parameters to be estimated. For this, we must nondimensionalise the system so that groupings are compared like for like (i.e. their dimensions are the same). Nondimensional scalings are listed in Tables 1 and 2 and the resulting nondimensional system becomes (dropping primes) with replacing (13) and (14) if TcdC inhibits TcdR production at the transcriptional level. We include, in the above, scalings of the nondimensional parameters (denoted by hats) according to a small parameter which arises naturally in the ratio of constitutive to TcdR-induced transcription levels, the first of which is naturally much smaller than the second. The scalings for c = c h A /c l A and c R = c h R /c l R follow. To increase the possibility that TcdC has an effect upon the system dynamics, we choose also to use c R C = O(1/ ) and λ C = O( ). Finally, since we assume cells die through lysis when toxins are secreted, we take the rate of toxin secretion, μ, to be relatively slow to avoid the population dying out unrealistically quickly.

Intracellular dynamics
We solve the nonlinear system (10)-(15) numerically (using ode23 in Matlab) to yield a first level of insight into the dynamics of the components of the PaLoc. Since the time taken for the cells to enter stationary phase in our simulations is roughly τ = 10, we choose this as our cut-off point for TcdC transcription, i.e.  (10)-(15) using = 0.1 (i.e. small, as required) and a range of values of β, representing the rate of TcdC binding to TcdR. Large values of β lead to lowered TcdA and toxin levels while TcdC is present in the cells, but they ultimately cannot prevent toxin levels attaining those seen for β = 0, albeit at a later time point. Here, and henceforth, we use the arbitrary initial conditions N (0) = 0.01,

representing the cells not producing substantial levels of toxins initially
(though using a step function generates a rather sharp transition in TcdC levels, the qualitative conclusions are not affected by adopting a smoother function to represent this switch). In order to investigate the role of TcdC in toxin production we consider various values of β (the rate of TcdR inactivation by TcdC), see Fig. 2. We find -for sufficiently large values of β, toxin production can be lowered in the presence of functional TcdC, but as soon as TcdC disappears from the cells through natural degradation (around τ ≈ 20), toxin levels are restored to the level achieved if TcdC cannot bind TcdR (compare to the dotted line of Fig. 2). The strength of TcdR to activate toxin production in this phase is further increased by the fact it is positively auto-regulated, meaning that by the later stages of growth TcdR levels increase significantly. -When toxin levels are lowered (in τ < 20), they reach extremely low levels (consistent with the ultra-sensitive nature of switches that can be induced by mole-  -As TcdC approaches its non-zero steady state, TcdR proteins are impeded, preventing a long-term increase in TcdA and subsequently toxin levels. -Setting β = 0 and removing this TcdC-induced obstruction (dotted line) induces significantly higher toxin levels, both during the transient dynamics and in the steady state of TcdA (notice that toxin levels even for β = 1, given by the dot-dash line, are markedly lower than those for β = 0).
The fact that low toxin production cannot be maintained once TcdC vanishes from the cells suggests that the system is monostable in our parameter range and Fig. 4 confirms this for varying β (this steady-state plot was created using the bifurcation and continuation software XPPAUT).
We observe that if TcdC is produced constitutively throughout growth and β is greater than roughly 3 (for = 0.1) then TcdC has a significant negative effect on TcdA (and consequently on toxins). However, as alluded to earlier, if β is sufficiently large, it is possible that it actually reduces TcdA to unrealistically low levels, leaving  Fig. 4 A steady-state plot of the system in response to changes in β ( = 0.1) if TcdC is produced throughout growth (F = 1 for all τ ). We omit extracellular toxin levels from these solutions since they never achieve steady state under our assumptions (as we assume toxins are secreted as long as cells and TcdA are present), but any effect on toxin levels can be inferred from the solution for TcdA (since extracellular toxins are downstream of the rest of the system, this has no consequence on any other variable). Notice that the dip in TcdC levels occurs because TcdC is free (unbound) at low β, while for large β TcdR is sequestered sufficiently early for TcdC to be produced and remain free in the cell with no additional TcdR to bind. Intermediate values of β result in a decrease in TcdC as a result of binding to TcdR. Given that the system is monostable, our simple choice of initial conditions given in Fig. 2 does not affect the steady state of the system only a small interval of β that would yield results that might be observed in practice. We shall see in Sect. 3.3, where we shift our focus to the parameterĉ R C , that altering production rather than binding rates could produce more realistic levels.
These numerical results suggest that for non-functional TcdC to be the cause of any hypervirulence in C. difficile, it should be produced throughout growth. If TcdC is produced only in early growth then β must be sufficiently large for it to have any effect in the short to medium term, but eventually maximum toxin levels will be achieved once TcdC disappears from the cells, preventing functional TcdC from being the cause of lower toxin levels in the long term.
In all these simulations, we have assumed the mode of inhibition by TcdC is through binding to TcdR. In Figs. 5 and 6 we illustrate that equivalent conclusions hold if TcdC instead blocks tcdR transcription, though in this case toxin levels never reach amounts as low as in the previous case.
In the following sections we exploit asymptotic analyses to uncover the timescale on which any effect on toxin production by TcdC could first manifest itself (Sect. 3.2) and how fast TcdC must be produced (in relation to the other proteins in the PaLoc) for it to have a long-term effect on toxin levels (Sect. 3.3). In the interests of brevity, we limit  (12) and (15)-(17) lines illustrate an increasing inability for TcdC to bind and block tcdR transcription (i.e. this represents the hypervirulent case). Toxin production is only mildly suppressed while TcdC is present in the cells this to the model representing TcdC-TcdR heterodimerisation but the methods applied therein are equally applicable to the model whereby TcdC acts via transcription.

Time-dependent asymptotic analysis
In this section we split the full transient solution into distinct timescales where different reactions dominate in each case (see Jabbari et al. 2010 for an additional example of such time-dependent asymptotic analysis in the context of gene regulation networks). This enables mathematical visualisation of when TcdC first has an effect on toxin levels, providing an estimate for the minimum time TcdC must exist in the cells for it to exert its hypothesised effect. In addition, the simple analytical expressions will facilitate parameter inference from time course measurements of protein levels or gene expression in future work. Three timescales emerge from the asymptotic analysis; the final two must be split into three cases depending on the value of what we will show to be a key combination of parameters: Fĉ R C −ĉ R (in fact when this is negative the system resolves to steady state already on the second timescale). We summarise the analysis in Table 3 and include full details in the Appendix.
Given that secreted toxin levels never achieve steady state under our assumptions, for simplicity we drop this variable from the remainder of the study and infer all results from TcdA, the precursor to secreted toxins. Time-dependent numerical solution to the full system, assuming TcdC inhibits TcdR production at the transcriptional level and is produced throughout growth (Eqs. (10)- (12) and (15)-(17) lines illustrate an increasing inability for TcdC to bind and block tcdR transcription (i.e. this represents the hypervirulent case). As with the case where TcdC binds TcdR, in this scenario toxin production can be lowered in the long term Table 3 A summary of the scalings required for the time-dependent asymptotic analysis of Sect. 3.2 and the large time behaviour of each variable on each timescale (since this facilitates derivation of the scalings for the subsequent timescale). The final two timescales differ depending on the value of Fĉ R C −ĉ R . When Fĉ R C −ĉ R < 0, the system settles to steady state already on the second timescale. We have that ρ =λ C +μ/(1 + ηλ E ) − FU R β/ĉ R and θ = Fλ R /β. We replace the steady state approximations for A, R and C by the relevant equation number for the final timescale of Case III in the interests of space

Timescale 1: τ = τ
On this first timescale, all variables except τ are O(1). The system (10)- (14) becomes meaning that the O(1) (i.e. the largest) terms yield the following leading-order behaviour: Only TcdA and TcdR dynamics evolve on this very early timescale, both seeing the positive influence of TcdR upon their production rates. Notably, TcdC remains at its initial value. These approximations to the behaviour on an early timescale are depicted in Figs. 7, 8, and 9 (division into parameter-led subcases on subsequent timescales requires multiple plots).

Timescale
The influence of TcdR upon TcdA and its own production results in A and R increasing on this timescale. The resulting rescaled problem reads Notice that the equations for N , E and A decouple from the remaining system at leading order and can be solved to give where for largeτ . At leading order therefore, cell number is governed solely by the logistic equation, toxin production not being significant enough to negatively impact on growth. Adopting (35) and combining the leading-order terms of (30) and (31) enables us to obtain the largeτ behaviour of C andR: The dynamics of these two variables are determined by whether Fĉ R C −ĉ R is negative, zero or positive. We examine each case separately (note that (35) holds for all three scenarios outlined below).
Case I: Fĉ R C <ĉ R If Fĉ R C <ĉ R , the leading-order system evolves to steady state on this second timescale, namelyR see Fig. 7. TcdA is not affected by any parameters relating to TcdC at leading order since TcdR is produced at a faster rate than TcdC, thus rendering TcdC ineffective. Instead, throughout the whole time course of the simulation TcdA is governed solely by its own production rate, balanced by its loss rate.
Case II: Fĉ R C =ĉ R Following the analysis detailed in the Appendix, we find for largeτ . In this instance,ĉ R C is influencing TcdR levels on this intermediary timescale but the effect has not yet filtered through to TcdA.
Case III: Fĉ R C >ĉ R Under this parameter regime, hold for largeτ . As above,ĉ R C is not yet affecting TcdA on this timescale, but it does have a significant effect on both TcdC and TcdR, with the difference betweenĉ R C andĉ R being key to determining the rates at which these proteins are produced. In addition, in contrast to Case II where an increase in the rate of binding of TcdR and TcdC negatively impacts on the levels of both of these proteins, here (where TcdC production occurs faster than that of TcdR) β affects only TcdR at leading order because TcdC is produced in sufficient quantities for it to be replenished following loss through complex formation. This will prove to be a key property of TcdC for it to be a negative regulator of toxin production in this study.

Timescale 3
Case I: Fĉ R C <ĉ R When Fĉ R C <ĉ R , the system has already evolved to steady state on Timescale 2 and therefore no further discussion is needed here.

Case II: Fĉ
On the third timescale, if the rates of TcdC and TcdR production are identical, we have The first three of these equations give (the O( ) term being needed to obtain the leading-order balance in (46)), meaning that, again, TcdA is not affected by TcdC at leading order at any point in the time course. We find (see Appendix) where The resulting approximations are depicted in Fig. 8. Remember that in this scenario Fĉ R C =ĉ R , which in terms of the dimensional parameters requires that the rate of TcdC production is exactly the same as the fastest rate of TcdR production. Though it is nigh impossible that this occur in reality, we include the analysis here for mathematical completeness. Given that the rate of protein production is unlikely to be absolutely constant, this could in theory cater for the scenario where the production rates of TcdC and TcdR waver around similar values.
Case III: Fĉ R C >ĉ R ,τ = −1τ , C = −1Č ,R = Ř In this case, where TcdC production occurs at a faster rate than that of TcdR, we have on the final timescale: so that (47) holds also in this case. The remaining approximations are given by The approximations forŘ and A follow. At steady state therefore we havē A =ĉĉ c R C has a direct influence on A at leading order on this timescale (τ = O(1/ )) meaning that TcdC must remain in the cells for this length of time for it to have any effect upon toxin levels (even when, given that Fĉ R C >ĉ R , TcdC production is faster than that of TcdR). We note from Fig. 9 that this is beyond the onset of stationary phase.
The above workings make it explicit that the value of Fĉ R C −ĉ R is crucial for determining the overall behaviour of the PaLoc. TcdC can only be effective in lowering toxin production if it is produced at a faster rate than TcdR (Fĉ R C >ĉ R ). Otherwise the system will settle to a steady state with TcdA present in high levels. If TcdC production is switched off at some stage of cell growth (F = 0 for some τ ) then only Case I can ever arise, and TcdC does not affect toxin levels.
Though providing useful suggestions for experiments already, it would be interesting to extend this study to consider the possibility that the value of Fĉ R C −ĉ R is not constant throughout growth, but rather is subject to variations. If high-frequency experimental time course data were available to suggest that this were the case, this would provide a fascinating extension to the asymptotic analysis.

Steady-state approximations
Nonlinear systems of the kind presented in this study rarely have explicit steady-state solutions from which information can be gleaned about the influence of specific parameters on the system as a whole. We have already seen in the previous section how the leading-order terms of each equation (these vary depending upon the parameter regime and timescale of interest) can be extracted to obtain tractable asymptotic approximations to the full steady state. These expressions enable clear insight into which parameters, variables and reactions are dominating the behaviour of the full system in the long-term. We extend these approximations here to include lower order terms. We thus summarise our findings of the previous section whilst also gleaning further information about the system, in addition to tracking the steady state of the system in response to variations to the key parameter identified earlier,ĉ R C : the ratio of TcdC to TcdR production-see Fig. 10.
Following the results of Sect. 3.2 which illustrate that TcdC must be present throughout growth for it to affect toxin levels, we here further exploit the aforementioned asymptotic techniques to infer properties of TcdC which would enable it to affect toxin production at steady state if this were indeed the case (i.e., F = 1 for all τ ). For this we require expressions for C, R and A (toxin levels follow directly from intracellular TcdA levels and we have seen from Sect. 3.2 that E and N are given at steady state by the appropriate expressions in (47) and (48)). We are able to unravel the magnitude of TcdC's contribution to the steady-state of the system as a whole in the regimes outlined in the previous section.
Taking c R C = O(1/ ) as before (and c R = O(1/ ) to be fixed), we must consider two regimes:ĉ R C < 1 andĉ R C > 1 (we omitĉ R C = 1 for brevity: the leading-order steady states can be obtained from Sect. 3.2.3). In both regimes, we derive the asymptotic expansions of the variables to second order to gain a full picture of the effect of varyinĝ c R C on the system and to glean an estimate of the size of TcdC's contribution on toxin levels.

Regime 1:ĉ R C < 1
When c R C = 1 ĉ R C andĉ R C < 1 both TcdC and TcdR production rates are fast but that of TcdC is lower (similar to Case I in Sect. 3.2). The following scalings and approximations hold (bars denote variables at steady state): where (notice that the leading-order terms, r −1 , c 0 and a −1 , concur with (35) and (37) from the dynamic solutions of Sect. 3.2). These approximations are marked by crosses in Fig. 10. For each variable, the scaling of the leading-order terms with respect to reflects its magnitude. Hence in this parameter range, levels of TcdR and TcdA are significantly higher than those of TcdC. Significantly larger amounts of TcdR over TcdC mean that the former has more influence over TcdA levels, preventing TcdC from exerting a negative influence and resulting in high levels of TcdA. Crucially,ĉ R C only appears in the approximation for A at second order, thus rendering any effect TcdC has to be small relative to that of R via the parameterĉ (the rate of TcdA production).

Regime 2:ĉ R C > 1
In this scenario TcdC is produced in higher amounts than TcdR (corresponding with Case III of Sect. 3.2) and hence their orders of magnitude switch around and TcdC is able to exert an influence over TcdA. The scalings and approximations are given bȳ where (69)-(75) are depicted in Fig. 10 with circles; it is evident that these track very closely the true steady state of the full system. In this instance, TcdC levels are notably higher than those of TcdR, and although it is still O(1/ ), TcdA is now affected byĉ R C at leading order (notice the dependence of r 0 onĉ R C , and consequently of a −1 ). With an increase toĉ R C rather than β (as investigated in Sect. 3.1 and Figs. 2, 3 where TcdA levels quickly become negligible), TcdA levels can be lowered but do remain at significant levels in the cell, thus a negative regulator acting via binding to TcdR and produced in much higher amounts than TcdR could feasibly be responsible for the difference between toxin levels in so-called hypervirulent and non-hypervirulent strains (we stress again, however, that the negative regulator must be present throughout growth).
Since TcdA remains at O(1/ ), somehow inducing tcdC transcription would not be successful in eliminating toxin production, even if TcdC is indeed a negative regulator. An implication being that this approach is unlikely to be a viable option for a novel therapeutic strategy.
In summary, forĉ R C to influence TcdA levels at leading order, we see thatĉ R C must be greater than one. Remembering that we have already scaled c R C to be O(1/ ), this amounts to the unscaled nondimensional c R C being at least O(1/ ). Translating this into dimensional parameters, we have that From our original parameter choice, we also have Taking (76) and (77) together, it is evident that the rate of TcdC production (c C ) must be at least as high as the highest rate of TcdR production (c h R ) for it to have a noticeable negative effect upon toxin levels at steady state. This could be investigated experimentally in order to identify if TcdC is produced in sufficient quantities to exert its proposed negative effect. Likewise, this property could be applied to any other protein that is identified as a possible candidate for inhibiting toxin production in this manner: it must be produced at a faster rate than TcdR.

Discussion
Different experimental approaches have yielded conflicting results regarding the role of TcdC in the PaLoc. In this complementary theoretical study, we have shown how a mathematical model can assist in unravelling the purpose of particular proteins in such networks. While in silico conclusions cannot be assumed automatically to translate into in vivo knowledge, they can be used to recommend experiments for corroboration or simply to provide guidance on the optimal way to obtain more information experimentally. We have indicated that TcdC production should be faster than that of TcdR for it to be able to inhibit TcdR action effectively (if it acts via binding), in addition to the duration of time over which TcdC must actually be produced for it to have a long-term effect (longer than the onset of stationary phase). Both of these things can be investigated experimentally to determine whether TcdC really is a potential candidate for toxin inhibition. These results are equally applicable to alternative proteins, and thus could be used as a generic indicator of potential candidates for toxin inhibition via TcdR-binding.
Current experimental results suggest that TcdC is produced only in early growth (Hundsberger et al. 1997;Govind et al. 2006;Carter et al. 2011) in a laboratory setting (i.e. in batch culture) but it should be considered a possibility that an in vivo infection could consist of cells entirely in exponential growth and thus TcdC could be produced throughout infection. However, this cannot explain the discrepancy in experimental results which suggest that when TcdC is produced only in early growth, it can have a long-term effect on toxin levels. Our findings imply that if TcdC does block TcdR action, then toxin production should only be delayed in this setting, agreeing therefore with Bakker et al. (2012) and Cartman et al. (2012) that TcdC is unlikely to be responsible for lowered toxin levels.
The range of β (the rate of TcdC binding to TcdR) which allows for a potentially realistic response to TcdC if it were a negative regulator of toxin production (i.e. some inhibition, but not a complete block) is fairly narrow under our parameter regime, suggesting that this would be a very sensitive mechanism (and perhaps therefore somewhat risky if a total switch is not what is required) for C. difficile to employ. However, any complete block could be caused by our assumption that, once bound, TcdC and TcdR do not separate (in doing so we were assuming TcdC operates at maximum efficiency, thus if it cannot lower toxin production in the long-term under this assumption we can be more certain it cannot do this under more slack conditions). Importantly, this result could mean that a novel drug which acts by binding TcdR and consequently blocking its action could be effective in quenching C. difficile infection via molecular titration. Overexpressing production of a negative regulator, whether it acts by binding the protein or blocking tcdR transcription, appears to be a less promising strategy since in our simulations toxin production was not abolished in these cases.
Of course, a number of additional assumptions were made in the formulation of the model, but we stress that these do not affect our results. For instance, TcdE may not actually be responsible for toxin secretion and therefore the formation of extracellular toxins (as argued by Olling et al. 2012), but all of our conclusions regarding toxin levels can be drawn from the variables representing intracellular TcdA (extracellu-lar quantities were really only considered in Sect. 3.1 for completeness). Similarly, including TcdB in the model would have given identical qualitative results to those seen for TcdA.
This study highlights the benefits of using a combined numerical and asymptotic approach to investigate gene regulation networks. While both enable a visualisation of the dynamics of the system, including investigating a range of parameter values, the asymptotic analysis permits significant insight into which reactions, variables and parameters are actually determining the behaviour of the system. The generation of smaller and less complicated subsystems which capture the behaviour of the full system on timescales or parameter regimes of interest are useful not only to improve our understanding of the system as a whole and to simplify the underlying mathematics, but also to pinpoint elements worthy of further experimental investigation.

Summary
If elements of the pathogenesis-related regulatory network of C. difficile are to be used as targets of novel drugs (in order to minimise the likelihood of the bacteria developing resistance to the drugs), it is crucial that efforts across a range of disciplines be adopted to ensure we fully understand the implications of interfering with the life-cycle of this important bacterium.
TcdC has long been assumed to play a pivotal role in toxin production with only more recent studies suggesting that its contribution is in fact negligible. Our findings here support the hypothesis that TcdC negatively regulates toxin production in conventional (i.e. non-hypervirulent) strains and that a lack of its functionality could be responsible for any increase in toxins in hypervirulent strains only if it is produced throughout growth and in higher quantities than TcdR. We thus present model-derived properties which can be investigated in a laboratory environment to lend support to either of the biological hypotheses. This will yield experimental data which can further inform the model and strengthen the reliability of future predictions.
A.1 Timescale 1: τ = τ On this first timescale, all variables except τ are O(1). The scaling τ = τ enables the appearance of the fastest reactions to appear here at leading order. We assume all initial conditions are O(1) implying that the bacteria do not produce toxins in significant quantities initially and do not require scaling. The system (78)-(82) becomes meaning that the O(1) terms yield the following leading-order behaviour: In order to link to the following timescale we must derive the long-term behaviour of each variable and match these to the subsequent short-term dynamics: we have that R ∼ĉ Rτ and A ∼ĉτ for largeτ (91) on this timescale, while the other variables remain constant. This informs our choice of variable scalings on the subsequent timescale: since both R and A behave likeτ for largeτ , they must be scaled with the same power of as isτ to move to the second timescale. The relevant power of is chosen to bring the mathematically appropriate terms into the leading-order behaviour on the second timescale.
A.2 Timescale 2:τ = −1τ , A = −1Ã , R = −1R When identifying the scalings required to move to a next timescale, it is helpful to examine the comparisons between the approximations on the previous timescale and the full system (Figs. 7,8,9). It is clear that the approximations for A and R grow too fast on the first timescale. Similarly, those for N , E and C do not grow at all. These facilitate the derivation of the above scalings for the second timescale. The rescaled problem reads Notice that the equations for N , E and A decouple from the remaining system at leading order and can be solved to give where for largeτ . Adopting (100) and combining the leading-order terms of (95) and (96) enables us to obtain the largeτ behaviour of C andR: The dynamics of these two variables are determined by whether Fĉ R C −ĉ R is negative, zero or positive. We examine each case separately (note that (100) holds for all three scenarios outlined below).

A.2.1 Case I: Fĉ R C <ĉ R
If Fĉ R C <ĉ R , the leading-order system evolves to steady state on this second timescale, namelyR (see Fig. 7) and therefore no further timescales are required for this parameter regime.
Substituting (103) back into the leading order terms of (96) (again, using N ∼ 1) gives for largeτ However, because C 1, From (103) which gives and hence from (95) and (100)  When Fĉ R C <ĉ R , the system has already evolved to steady state on Timescale 2 and therefore no further discussion is needed here.

A.3.2 Case II: Fĉ
On the previous timescale C ∼τ 1 2 and R ∼τ − 1 2 (see (105) and (106)). Thus, ifτ is scaled with α then C must be scaled with α 2 and R with − α 2 . The above scalings (that set α = −1 to achieve a new balance by bringing to leading-order TcdR production and TcdC loss-see Fig. 8) abide to this rule. On this third timescale, therefore, if the rates of TcdC and TcdR production are identical, we have The first three of these equations give (the O( ) term being needed to obtain the leading-order balance in (114)), Combining (113) and (114) and taking leading-order terms, we obtain the following expressionĉ but from (114) and the condition that Hence, withČ → 0 asτ → 0. This can be solved to givě where The resulting approximations are depicted in Fig. 8.
A.3.3 Case III: Fĉ R C >ĉ R ,τ = −1τ , C = −1Č ,R = Ř On the previous timescale C ∼τ andR ∼τ −1 for this Case (see (108) and (109)). Therefore, (following the reasoning outlined for Case II) the above scalings facilitate the transition between the two timescales. In this case, where TcdC production occurs at a faster rate than that of TcdR, we have on this final timescale: so that (115) holds also in this case. The remaining approximations are given by (129) can be solved to givě (matching toČ ∼ ((Fĉ R C −ĉ R )/ĉ R C )τ asτ → 0). The solutions forŘ and A follow. At steady-state therefore we havē A =ĉĉ