Abstract
Background
In this paper the dynamics of the transcriptiontranslation system for XlnR regulon in Aspergillus niger is modeled. The model is based on Hill regulation functions and uses ordinary differential equations. The network response to a trigger of Dxylose is considered and stability analysis is performed. The activating, repressive feedback, and the combined effect of the two feedbacks on the network behavior are analyzed.
Results
Simulation and systems analysis showed significant influence of activating and repressing feedback on metabolite expression profiles. The dynamics of the Dxylose input function has an important effect on the profiles of the individual metabolite concentrations. Variation of the time delay in the feedback loop has no significant effect on the pattern of the response. The stability and existence of oscillatory behavior depends on which proteins are involved in the feedback loop.
Conclusions
The dynamics in the regulation properties of the network are dictated mainly by the transcription and translation degradation rate parameters, and by the Dxylose consumption profile. This holds true with and without feedback in the network. Feedback was found to significantly influence the expression dynamics of genes and proteins. Feedback increases the metabolite abundance, changes the steady state values, alters the time trajectories and affects the response oscillatory behavior and stability conditions. The modeling approach provides insight into network behavioral dynamics particularly for smallsized networks. The analysis of the network dynamics has provided useful information for experimental design for future in vitro experimental work.
Background
The filamentous fungus Aspergillus niger is an important organism in the production of enzymes and precursors for the food and chemical industries. Industrial citric acid production by A. niger represents one of the most efficient, highest yield bioprocesses in use by industry. The xylanolytic activator gene xlnR is a main controlling gene in the XlnR regulon of A. niger.
The XlnR regulon is activated by Dxylose in the culturing media [1]. The current description of this system is based on static interpretation of the system. Experiments [2] showed, however, that the expression of genes in the XlnR regulon is dynamic. Therefore, to advance the application of A. niger by better understanding of the XlnR regulon, the dynamics of the regulon needs to be quantified. For this purpose time course experiments are scheduled. However, planning of the experiments is improved by quantifying the behavior by a simulation and analysis study prior to the experiments.
For the XlnR regulon, literature information on the network structure was used as a basis for the simulations. To our knowledge, currently very little has been done on modeling the dynamics of the XlnR regulon and also on time course profiling of the genes that constitute the XlnR regulon in A. niger. The challenge with genetic network modeling lies with determining a specific equation formalism to represent the network structure. One of the suggested strategies of modeling using differential equations is to fix the form of the equation [3,4]. Prior knowledge on the network structure is essential to develop a quantitative model [5]. The descriptive information on the XlnR regulon [1] enables us to hypothesize models for the interaction between the different network components.
Generally, in the study of biological networks, positive feedback (PFB), negative feedback (NFB) [69], feedforward loops and time delay [10] have been shown to be influential. NFB loops cause oscillatory behavior if the signal propagation around the feedback loop is low and PFBs can lead to a bistability behavior [11]. Overall, feedback plays an important part in biological networks by allowing the cell to adjust to the repertoire of functional proteins to current needs. Other examples of biological systems in which the effect of feedback and time delay were extensively studied can be found in the developmental regulator Hes1, which inhibits its own transcription [12,13] or in the NFB loops in the p53 response [14,15]. Bliss et al. [16] investigated conditions on parameters that ensure stability of the unique steady states in an operon using differential equations modeling. These authors then chose parameters that allowed the model to describe the tryptophan operon of E. coli. Their investigations focus on network stability in steady state. They showed that with parameters corresponding to a mutant with reduced repression, stability conditions were violated leading to oscillations. The analysis of the condition on parameters that ensure steady state stability also lead to insight into direct repression of a gene by its own product [16].
In computational systems biology, numerous studies have been done on genetic network reconstruction using time course data but little attention has been given to understanding the network dynamics. It is crucial to understand or at worst have a fuzzy idea of a biological network dynamics if one is to gain deeper insight into the biological network dynamics and functionality.
This paper concerns the analysis of the network dynamic behavior, the effect of feedback loops and the conditions under which oscillatory responses in metabolite expression may be exhibited are investigated. Modeling of the XlnR regulon is explored by using nonlinear differential equations and Hill functions for the transcription and linear reaction kinetics for the translation process. To ensure that detailed aspects of the transcriptiontranslation model formalism are captured, some assumptions are incorporated in the modeling. In silico perturbation experiments were performed by triggering the genetic network at steady state. The advantages of using ordinary differential equations (ODEs) are vast since they are capable of modeling degradation effects and causal effects in a network [17].
Applications of dynamical systems in modeling transcription regulatory networks can be found in [1821]. Many of these studies used the continuoustime domain to model gene expression as biochemical processes using in ODEs. The modeling and analysis aims to identify which factors determine the dynamics to aid and guide future time course experimental studies. In our work, we highlight the need to understand the dynamics of biological networks with the hypothesis that modeling and experimentation should go hand in hand.
Methods
Regulation mechanism for the XlnR regulon
In the model organism Aspergillus niger, transcription of genes encoding xylanolytic and cellulolytic enzymes take place [1]. Activation enables the degradation of the cellulose and hemicellulose from plant cell walls. XlnR is a zinc binuclear cluster protein consisting of about 875 amino acids. It is suspected that XlnR binds as a monomer [1]. The xlnR gene is induced in the presence of Dxylose in the culturing media and repressed in the presence of the carbon catabolic repressor, CreA [22].
Gene regulation can take place at different stages of the central dogma of molecular biology (DNA→ RNA → Protein). These stages include among others transcription, translation and posttranslational modifications (PTMs) of the associated protein. In Figure 1 a scheme of the activities in the XlnR regulon is given. The xlnR gene is induced by Dxylose. At induction the xlnR gene produces messenger RNA (mRNA) which is translated in proteins. These proteins then activate the target genes (TG). For the XlnR regulon, the number of target genes are estimated to be in the order of 20 to 40. In Figure 1 all target genes are represented by TG. After transcription and translation of the target genes, the target proteins (TP) are obtained. Protein from PTMs can be involved in the regulation of the xlnR gene trough a feedback loop. At each step in transcription and translation mRNA and proteins can be degraded and/or used for other processes (D1–D4).
Figure 1. The XlnR regulon scheme The XlnR regulon induced by Dxylose in the presence or absence of CreA. The representations P1 and TP are the proteins from the xlnR gene and target genes, respectively. The terms mRNA1 and TGmRNA represent the transcription products from the xlnR gene and target genes, respectively. FB represents the feedback protein in which posttranslational modifications might take place.
Transcription model
Commonly, hyperbolic functions and the sigmoid class of functions are used to represent the kinetics of gene regulation [23]. Such functions mimic the nonlinearity in gene regulation, by assuming that a critical amount of protein has to be accumulated before a gene can be considered regulated or repressed. The most common form of function used for modeling gene transcription is the Hill function [24,25].
Let the vector z = [z_{1},…,z_{n}]^{Τ} represent the concentrations of the translated proteins corresponding to the genes 1,…,n; where n is the number of genes involved. Throughout this work, the notation, e.g. z_{i} is used to represent the time dependent variable z_{i}(t) (where t ∈ [0, ∞)). Then the activating and repressing functions are given by
where ψ^{–}(z_{i}, θ_{i}) = 1 – ψ^{+}(z_{i}, θ_{i}), θ_{i} is the gene specific halfsaturation parameter and the positive number h represents the Hill coefficient. The regulation mechanism for each target gene i is captured by the function Ψ(z_{i}, θ_{i}) in (1).
According to Hasper et al. [26] there is evidence that although most zinc binuclear cluster proteins bind as a dimer, it seems that XlnR binds as a monomer  therefore, a Hill coefficient with h = 1 is used. Given the availability of structural prior knowledge and that the master regulator activates the target genes, the nonlinear system that models the transcription process is given by
where k_{i}_{1} = 1/θ_{i} for h = 1, x_{i}  mRNA concentration from gene i, ρ_{i} is the basal (or leaky) transcription rate for gene i and is associated with very low levels of mRNA, k_{i}_{1}  effective affinity constant for gene 1 activating gene i (i = 2,…,n), k_{is}  maximum synthesis parameter for gene i, k_{id}  first order degradation rate (or consumption rate) for gene i, x_{0}  vector of initial mRNA concentration, z_{i}  concentration of translated protein from gene i, b = [b_{1},…,b_{n}]^{Τ} is the input matrix and u = [u_{1},…,u_{n}]^{Τ} is the input vector (gene triggering compounds). The model formulation with no feedback later on aids the assessment of the marginal contribution of a feedback loop in the network dynamics.
Translation model
A system of linear differential equations to model the protein abundance (translation process) is then considered. The linear model representations (3) are used to capture the dynamics of the translation process with both the production and degradation terms being linear.
where r_{i}  specific translation rate for gene i, η_{i}  degradation rate for protein i. The z_{i}’s represent the target proteins in the scheme in Figure 1. At steady state the response rate and degradation rate balance, i.e. ẋ_{1} ≈ ẋ_{2} ≈ … ≈ ẋ_{n} ≈ 0. By setting the transcription rates ẋ_{i} = 0 for all i in (2), it follows that
The steady state values in (4) are based on the assumption that, for a small time window the change in concentration of the input stimulus and metabolite concentrations remain nearly unchanged. The model specifications for the transcription and translation process describe the rates of change of concentration of the genes and proteins. Overall, the system of 2n coupled differential equations in (2) and (3) describe the network dynamics.
System stability
The interesting case to analyze is the systems behavior in the absence of the inhibitor, CreA. Let us denote the equilibrium concentrations of mRNA and protein quantities by the vectors and respectively. Using (3) the steady states lead to the relationships for all i. The stability of each steady state (from (2) and (3)) can be analyzed using Hopf Bifurcation, an analytic approach that has been widely used in investigating stability conditions in gene expression networks [2730].
Let F : ℝ^{2n} → ℝ^{2n} be a set of smooth functions (with F = (F_{1},…,F_{2n})) that capture the XlnR regulon system dynamics. In this case we have F_{1} = ẋ_{1},…,F_{n} = ẋ_{n}, and F_{n+1} = ż_{1},…,F_{2n} = ż_{n} in (2) and (3) respectively. By definition, the Jacobian matrix is given by
This Jacobian matrix is used to assess the regulon stability and to identify which parameters dictate the transcript abundance. First, consider a case of three genes and three proteins, n = 3. The Jacobian is given by
where x = [x_{1}, x_{2}, x_{3}]^{Τ} and z = [z_{1}, z_{2}, z_{3}]^{Τ} . The corresponding steady states of the vectors and can be computed, accordingly. Using expressions (2) and (3) in (5) we obtain
where
for i = 2, 3. A similar generalized expression for can be obtained given a regulon with a known number of transcripts n. The characteristic polynomial obtained from (7) is given by
In the case of this regulon, the derived characteristic polynomial turns out to be the same as the determinant, i.e. . Using the expression (9), conditions that ensure global stability can be established on the parameters.
The formulations of the Jacobian matrix and the eigenvalue spectra can be extended to an ndimensional network system. The generalization for the eigenvalue spectra using a similar analysis as above leads to the expression
for n ∈ ℤ^{+}, a positive integer. In the network without feedback, it turns out that the trace of the Jacobian matrix is equal to the sum of all the eigenvalues (i.e. holds true). The above analysis shows that for the open loop system there is no possibility for oscillatory behavior to occur, and that the time constant only depends on the degradation coefficients.
According to Aro et al. [31], van Peij et al. [1] and Hasper et al. [32], the A. niger genes eglA, eglB, eglC, cbhA, cbhB, xlnB, xlnC and xlnD contain binding sequences (GGCTAAA) to the XlnR protein as well as binding sequences to CreA, a repressor protein acting in the presence of monomeric sugars (i.e., glucose) as an autoregulating mechanism. This property ensures that most target genes have similar expression dynamics in time.
Feedback in the network
Numerous transcription systems are known to include genes that regulate their own expression values [33]. In our analysis, to model the effect of feedback we hypothesize that the TPs and PTMs in the feedback loop in the scheme in Figure 1 only act on the xlnR gene. Therefore, only the equation that captures the dynamics for the first gene (x_{1}) has to be modified accordingly. The adapted equation is given by
where
is the repressor Hill function and C_{A}  quantitative activity state for CreA, k_{A}  inverse of the Hill constant of CreA. The term τ > 0 represents a time delay in the feedback loop. The sets S_{1} = {j  j = 1,…,m} and S_{2} = {l  l = m + 1,…,n – 1} where S_{1} ⋃ S_{2} = {1,…,n – 1} i.e. collection of all the target proteins in the regulon. All the supposed repressing and activating proteins are lumped in the sets S_{1} and S_{2}, respectively. The effect of the Dxylose and the feedback loop is modeled as additive. Equation (11) also specifies the build up of proteins and repression or activation of the xlnR gene through the feedback loop. Through the sequence of PTMs the protein availability in the feedback loop is delayed. All the other components representative of the target genes in the network models (2) and (3) remain unchanged.
xlnR gene promoter activity under feedback
Let us define the promoter activities by the expressions (13) and (14). The promoter activity corresponding to the case when an activating protein is involved in the feedback loop is represented by the term Γ_{A} and that for the case of a repressing feedback effect denoted by Γ_{R}.
The extracts from the denominator functions are given by (15) and (16), respectively.
These terms are used in the calculations for the activating and repressing promoter activity for the XlnR regulon. For the sake of illustrations, two target genes were considered (i.e. values of j = 1 and l = 2) in the simulation with one as an activator and the other as a repressor (the values k_{RL} = k_{AL} = 1 were used). We consider the sets S_{1} and S_{2} of unit elements which index the proteins that are responsible for regulating the xlnR gene through a series of mechanisms in the PTM channel.
Existence of oscillatory behavior
The eigenvalue spectra from the derived Jacobian matrix can be used for this analysis. The presence of at least a pair of eigenvalues with complex parts implies the existence of oscillatory behavior. The PTMs in the feedback loop may produce oscillatory behavior depending on the individual attributes of the target genes and the consequent proteins in the feedback loop.
We observed that in the absence of a feedback loop, the system dynamics is dictated by the degradation parameters. Active degradation of proteins or mRNA is a major part of many metabolic and stress response systems [11]. This may not necessarily hold true for the system with a feedback loop because of the structural change in the Jacobian matrix. The metabolites involved in this oscillatory dynamics are presumably determined by the individual biochemical and mechanistic attributes of the individual molecules. Therefore, no single hard rule for classifying which metabolites are responsible for the overall oscillatory behavior of the expression profiles of the genes and proteins can be stipulated. This might however be possible for some specific network pathways for which extensive information is available. Consider the Jacobian matrix corresponding to three genes and three proteins (17). We now analyze the effect of having a protein in the feedback loop. These proteins are considered to regulate the xlnR gene and their possible positions in the Jacobian matrix are indicated by ”×” as shown in (17).
Using the adapted model (11), the computed entry in the (1, n + i)th cell (i = j, l for all values of j and l) of the Jacobian matrix (5) is given by (18) and/or (19) depending on which proteins in the feedback loop are involved in the regulation of the xlnR gene. By taking partial derivatives of the function F_{1} with respect to the variable of interest within each of the sets S_{1} and S_{2}, we then have the more compact expression
This term corresponds to the repressing proteins, and the terms given by
for the activating proteins. The parameters k_{AL} and k_{RL} represent the lumped affinity constants for the activating and repressing proteins, respectively. Suppose that the XlnR protein (i = 1) is the only metabolite responsible for autoregulation in the feedback loop. By representing the corresponding entry at the (1, 4)th position in the matrix (5) by a nonzero parameter ω_{1} ∈ ℝ\ {0}, a parameter that intrinsically represents the autoregulation effect of the xlnR gene. The computed eigenvalue spectrum from (7) using the expression is given by
where
are the conjugate roots. The eigenvalues λ_{5}(·) and λ_{6}(·) may take on values from the real space, ℝ or the complex space, ℂ. From (21) and (22) we observe that oscillation can only be obtained if the condition ω_{1} < –(η_{1} – k_{1d})^{2}/4r_{1} for r_{1} > 0 is satisfied. This condition on ω_{1} signifies a contribution from a NFB loop in the XlnR regulon network. Notice that the expression (η_{1} – k_{1d})^{2} > 0 for all values of η_{1} and k_{1d}. This finding adds to consolidate the findings by Tiana et al. [15] from a theoretical analysis of three eukaryotic genetic regulatory network in which they attributed the existence of oscillation to a common design of a NFB with underlying time delay. Considering the expressions for the eigenvalues in (21) and (22), we observe that for a PFB effect of the XlnR protein there is no possibility for oscillatory behavior. This result does not necessarily hold true for the proteins in the feedback loop corresponding to the cells at positions (1, n + i) for i = 2,…,n.
The presence or absence of oscillatory behavior is insufficient for drawing conclusions about stability in system responses. Stability, using (21) and (22) exists if the conditions Re(λ_{5}(·)) < 0 and Re(λ_{6}(·)) < 0 are simultaneously fulfilled. Hence, the inequality
Solving the inequality (23) for ω_{1} leads to the condition ω_{1} <η_{1}k_{1d}/r_{1}. A similar analysis for the existence of oscillatory behavior and stability dynamics can be done for the other proteins in the feedback loop, for example at position (1, 5) and (1, 6) or combinations in the matrix (17). However, although such analysis is conceptually simple, the analytic expressions are very complex to work with. Information about the stability and oscillatory behavior is obtained by numerical solutions. An example is considered to investigate the time evolution of gene activity and protein abundance in the XlnR regulon.
Bifurcation analysis
Bifurcation analysis relates to stability on the system parameters. Stability properties for the system without feedback are given by (10), where it was shown that the roots of the characteristic polynomial correspond to the degradation rate constants for the mRNA expression and protein abundance. As these constants are positive, the system is always stable. In the case that one of them equals to zero, then the system is critically stable.
For the network with feedback loop, consider the conjugate roots in (21) and (22) denoted by λ_{i}(·) for i = 5, 6. In the event that the conditions Re(λ_{i}(·)) = 0 and  Im(λ_{i}(·)) ≠ 0 are simultaneously fulfilled  then there exists a Hopf bifurcation for the corresponding genes and proteins. Such a bifurcation occurs when the root of the positive discriminant function in (23) equates to the sum of the degradation parameters for the xlnR gene and XlnR protein, i.e.
or after working out becomes ω_{1} = η_{1}k_{1d}/r_{1}. This example illustrates the case of a feedback in the cell at position(1, 4) of the matrix in (17). The analysis for the other entries of ω results in highly complex expressions, therefore a numerical analysis is preferred.
Results
System specification
The analysis is illustrated by an example case. Consider a regulon network of three genes given a perturbation of Dxylose. The pulse perturbation takes place at time t = 0. During fermentation, the Dxylose is consumed and the Dxylose concentration follows the expression u(t) = u(0)(1/(β + e^{Kt})), where u(t) ≡ [Dxylose] and β > 0, with K = 0.3 and u(0) = 50 mM as the initial Dxylose concentration. The parameters used for the simulation are: b_{1} = 1, ρ_{1} = 2e – 3, ρ_{2} = 2.5e – 3, ρ_{3} = 1e – 3, k_{1d} = 0.5, k_{2d} = 0.4, k_{3d} = 0.3, k_{2s} = 5, k_{3s} = 6, k_{21} = 0.1, k_{31} = 0.1, r_{1} = r_{2} = r_{3} = 0.5, η_{1} = 1, η_{2} = 1 and η_{3} = 1.
Stability and response analysis  without feedback
The expression for the characteristic polynomial, P(·) in (10) is independent of the translation rate parameters r_{i}, the gene synthesis coefficient k_{is} and the terms in the expression (8). From (10) it can be observed that without feedback, the system is globally stable (i.e. the conditions and are satisfied for all and ). The system stability behavior is dictated by how fast the translation and transcription rates are (i.e. magnitudes of k_{id} and η_{i}).
In Figure 2 both the gene expressions in plot (B) and protein abundance plot (C) show similar behavioral dynamics. Moreover, with the chosen input pattern of Dxylose the target genes show phase plots similar in patterns but with variations that are dictated by individual gene or protein kinetic parameters. A relaxation time of τ_{R1} = 1/k_{1d} ≈ 2 hours is noticed for the master regulator and for the target genes, τ_{R1} <τ_{R2}, τ_{R3}. The relaxation time is an approximation for the time required for the system to relax into steady state. This represents the time it takes a system to react to a persistent external input (Dxylose).
Figure 2. Dxylose consumption, gene expression, protein abundance and phase plots: without feedback (A): The simulated trajectory for Dxylose consumption. (B): Gene expressions profiles. (C): Proteins abundance plots. (D): Phase plot for gene expression showing variation of mRNA concentrations of the xlnR gene and the other target genes, x_{i}. (E): Corresponding protein abundance phase plot.
Feedback in the network
Since the presence of CreA is a strong repressor that inhibits the xlnR gene activity by blocking the promoter binding site, we chose to model this influence by considering a switchlike function with H ∈ {0, 1}. Here the assignment of H = 0 and H = 1 means CreA is present and absent respectively. In the absence of CreA the protein products from the target genes are involved in regulating the activity of the master regulator. These protein products may either inhibit or activate the xlnR gene.
A comparison of the metabolite expression dynamics for the network with and without feedback loops in the absence of CreA is shown in Figure 3. The same parameter values in the section System specification were used for the simulation with the extra parameters from (11) being k_{RL} = 1 and k_{AL} = 1 and the lumped synthesis parameter from (11) chosen as k_{ls} = 1. Figure 3 indicates the enhanced metabolite expression as a result of incorporating a feedback loop with delay (with τ = 1) in the model  a result that is similar to what was observed by Maithreye et al. [34] during a theoretical kinetics analysis of the concentration of green florescent protein (GFP) in time.
Figure 3. Dxylose consumption, gene expression, protein abundance and phase plots: with feedback (A): The simulated trajectory for Dxylose. (B): Gene expression profiles with solid lines (–) showing the expression profiles for the genes in the absence of CreA. The corresponding dotted lines (⋯) show the simulated effect of competitive feedback (with τ = 1). (C): Protein abundance profiles (solid lines).
Activating and repressing feedback
The expressions (18) and (19) have the potential to yield oscillatory behavior in the metabolite response profiles. The oscillatory behavior (if and when it exists) is purely governed by the values of the system mechanistic parameters. Such oscillatory behavioral patterns of gene expression may vary from organism to organism, and can be detected from time series data if enough samples are taken.
To assess the effect of time delays in the transcription and translation processes, some cases were simulated. The results of the expression timedynamics for both the genes and proteins are shown in Figure 4. The simulations were performed for specific cases of τ = 1 hour and τ = 5 hours and the subsequent outputs compared. The metabolite expression patterns from the two cases are nearly similar with the main difference occurring at the maximum level. Overall, longer time delays lead to a small and non significant reduction in expression values.
Figure 4. Effect of time delay on expression (A)(B): Plots showing the effect of variation in time delay in the feedback loops corresponding to the transcription and translation processes, respectively. The observed effect on the responses is small except for the slight deviation at the peak of the expression profiles.
xlnR gene promoter site activity
The competitive effect of the activators and repressors for the promoter binding sites was also simulated. The effect of which transcription factor (TF) (either an activator or a repressor) wins occupancy of a promoter binding site depends partly on the strength of the synthesis parameter k_{ls} (Figure 5).
Figure 5. xlnR promoter region activity Plot of the xlnR promoter region activities defined by Γ_{A}, Γ_{R} and Γ_{A}Γ_{R} depending on the regulator. The term Γ_{A}Γ_{R}  is the combined affect of competitive binding to promoter region by activators and repressors. Plots (A) and (B) show the influence of weak (k_{ls} = 1) and strong (k_{ls} = 5) synthesis parameters respectively.
The promoter is most active (activity around 50 – 80%) when the regulon is fully active. This corresponds with the time window at which the network is fully responsive to the external perturbation. We observe that the xlnR gene activator has a tendency of occupying most of the promoter sites at any given time (Figure 5).
Bifurcation and oscillatory behavior analysis
A range of values of activating and repressing parameters ω_{1}, ω_{2}, ω_{3}, respectively on entries (1, 4), (1, 5) and (1, 6) in expression (17) was considered for analyzing the stability behavior of the network. It was observed that NFB on ω_{1} gives a stable system and values of ω_{2} and ω_{3} below –977 results in an unstable systems, Figures 6 (A)(C). The PFB effect of the XlnR protein on the xlnR gene leads to unstable system dynamics for ω_{1} > 1. This can be seen from Figure 6 (A). This result also leads to the conclusion that the XlnR regulon is unstable if the xlnR gene has a PFB from its own protein.
Figure 6. Stability index curves (A)(C): Plots of the stability indices for various values of ω_{j} ∈ (–1000,1000) for j = 1, 2, 3. A stability index value of –1 and +1 indicates stable and unstable responses of gene and protein expression in time, respectively. The value ω_{j} < 0 represents repressing feedback and ω_{j} > 0 is the activating feedback effect. Simulations performed using parameters from System specification and the pseudo steady state expression values at which the Jacobian is estimated.
An analysis of how the various feedback parameters affect the oscillatory behavior of the gene and protein expression was also considered. The results show that there exist threshold values (or a range of parameter values) for the feedback parameters ω_{j}’s for which their oscillatory behavior may or may not occur (see Figure 7 (A)(C)). The value ω_{j} = 0 corresponds to no feedback in the system and according to the previous analysis (under the subsection: stability and response analysis  without feedback), no oscillation occurs in this case. A transient oscillatory behavior is observed for values of the parameter ω_{j} ≈ 0 for all j, Figure 7 (B)(C). The observed stability curve corresponding to the XlnR protein in the feedback loop (ω_{1}) is a near reflection of the corresponding resultant oscillation curve (Figure 6 (A) and Figure 7 (A)).
Figure 7. Oscillation index curves (A)(C): The plots indicate possible oscillatory behavior for various values of ω_{j} ∈ (–1000,1000) for j = 1, 2, 3. Each ω_{j} is considered independently while the others are set to zero. The oscillation indices –1 and +1 represent nonexistence and existence of oscillatory behavior in the response behavior of the metabolites, respectively.
Discussion
The model gives a better understanding of the rate limiting steps in the process of activating the XlnR regulon and therefore, helps to define the biological control points. Similarly this knowledge can be used to obtain strains that have enhanced xylanolytic enzyme production. These enzymes are industrially of importance as food and feed additives, but are also part of a system that is used to bleach paper pulp. Given that the transcription rate and degradation rates have been shown to be the key parameters that dictate the systems dynamics for the XlnR regulon; this information is important for designing and sampling of time course experiments. Once the transcripts are unstable, the proteins get quickly degraded; otherwise they remain stable. This observation is linked to the Dxylose uptake in fermentation experiments. The consumption of Dxylose also indirectly controls the regulation of the target genes and therewith the breakdown of sugars.
Simulations showed that the dynamics of the Dxylose input function considered in the examples has an important effect on the profiles of the individual metabolite concentrations. This is particularly dictated by the value of the parameters in the external input function u(t). The larger the value of the K, the faster the consumption of Dxylose. This depends on the chemical reactions taking place in any given cell, or the saturation levels of the individual compounds in a cell.
Feedback significantly affects the response of the output profiles for the metabolites and changes the final steady state values (Figure 3). Further simulations showed that variations of the time delay in the feedback loop (τ = 1, 2,…, 5 hours) have a small effect on the pattern of the response (Figure 4). The stability analysis subsection shows that the metabolite response dynamics exhibits no oscillatory behavior for a network without feedback loop. For a network with feedback loops, the results from numerical analysis showed that feedback conditions for which the system is stable or for which the system exhibits oscillatory behavior can be obtained (Figure 7 (A)(C)). In modeling the feedback loop, time delay was accounted for and included in the model.
According to Bliss et al. [16] and Thomas and d’Ari. [35], including time delay in modeling biological networks is considered important because many biological systems exhibit some delays in their feedback loops. However, according to our finding (Figure 4), incorporating the time delay had no strong effect on the overall dynamics of the metabolite expression profiles.
The analysis shows that the existence or absence of oscillatory behavior is dictated by the numerical values of the individual mechanistic parameters. The conditions for oscillatory behavior follow from the eigenvalue spectra. The eigenvalue spectra analysis like that in (20) and the corresponding conditions for which all the eigenvalues are less than zero, gives also indication for the stability properties for the XlnR network with feedback loop. If all the eigenvalues satisfy the condition for all entries, then the system is stable, otherwise it is unstable.
Two scenarios can be considered: one in which the proteins involved in the feedback loop are activating and the other in which the proteins are repressing. The details of the expected behavioral dynamics from such a system requires a case by case analysis (like in Figure 6 (A)(C)) of the effects of the proteins in the feedback loop. A similar analysis can be extended to study the stability in case of a combined effect of any two or more proteins of interest. When the number of network components become large, obtaining explicit analytic solutions and expressions for the eigenvalues and other quantities of interest increasingly become complex  in which case the alternative of numerical methods can be used (see Figure 6 (A)(C) ). Thomas and d’Adri [36], and Thomas et al. [35] investigated the properties of mathematical Boolean net Modeling Genetic Networks works  investigations that provided significant insight into genetic network dynamics. In their work they showed the importance of NFB loops for maintaining homeostasis in levels of gene products. Our analysis leads to the observation that having a NFB loop stabilizes the response of the metabolite expressions (Figure 6 (A)(C)). However, there exists certain ranges of values of the strength of feedback effects that make the system unstable. This sets constraints to the feasible parameter for the system if instability is not observed. In some cases having a PFB loop yields a stable network response, Figures 6 (B) and (C). This result is in agreement with that found in a study by Maithreye et al. [34]. In their investigations they found that NFB loops provide stability and withstand considerable variations and random perturbations of biochemical parameters.
The effect of time delay on stability can be analyzed from a transfer function of the model in the ”s” domain, or by a transformation to the ”z” domain. In these cases the delay time is considered as a finite dimensional system. Stability analysis can be done by searching for stability properties in the ”s” domain or ”z” domain. Examples of other methods that deal with the delay times are given for state estimation in the work by Liu et al. [37] and Yu et al. [38].
The adaptive filtering approach developed in [38] is based on the adaptive synchronization setting, for estimating unknown delayed genetic regulatory networks with noise disturbance. Using this approach, no exclusive knowledge of system parameters is required, e.g. those lacking in the XlnR regulon and many other biological networks. Liu et al. [37] proposed an adaptive feedback control approach for simultaneously identifying unknown (or uncertain) network topological structure, unknown parameters of uncertain general complex networks with time delay from available mRNA data and estimation of protein concentration. The effectiveness and applicability of their approach was shown using in silico numerical simulations. In contrast to [37] and [38], we study the XlnR regulon dynamics and do not focus on the system structure identification and parameter estimation.
According to BalsaCanto et al. [39], powerful mathematical analytic tools highlight the value for successful study of many biological systems. However, such success can mainly be attributed to the unrelenting endeavors for an indepth understanding of both computational methods and the biological problems of interest. For the case of the XlnR regulon, our analysis provides a basis for understanding the behavioral dynamics of genes and proteins after network perturbation. This will form a basis for future wetlab experiments, particularly with the genes from the XlnR regulon. Given that the metabolite expression dynamics are known, this study provides a basis for strategic thinking in line with experimental design. The modeling approach used in this paper provides good information for understanding network behavioral dynamics particularly for smallsized networks. This is illustrated by the XlnR regulon in which even the simplest of structures can yield interestingly complex dynamics. Therefore, a reasons for having limited our focus to the regulon dynamics. Having detailed information regarding the basal parameters and the other mechanistic parameters might further improve the analysis and investigations into the network dynamics. Nevertheless, with informed parameter guesses, simulation studies provide good information into the systems behavior.
Conclusions
The investigations in this paper considers the XlnR regulon as a dynamic system instead of a static system. Our study provides insight into the dynamic properties of the XlnR regulon. By studying this system, it has become more clear that the transcription and translation degradation rate parameters and the Dxylose consumption profile dictate most of the dynamics in the regulation properties of the network. The existence of oscillatory behavior depends on the conditions of the mechanistic parameters in the feedback loop  conditions that cannot always be generalized analytically and therefore, must be treated by numerical analysis. The role played by feedback in the network dynamics was found to be significant on the expression dynamics of genes and proteins. This means that the effect of the feedback should be considered in the model if there is sufficient supportive biological need or evidence from data. Just like for most biological systems, this is no doubt an important piece of information for the accurate modeling of biological network.
The analysis of the network dynamics has provided useful information for future in vitro experimental work. Particularly the potential for hypothesis testing basing on this work and the design of related perturbation experiments to generate time series data. Once there are available techniques for the network structural identification and parameter estimation for the XlnR regulon can be investigated.
Authors’ contributions
LHdG provided the biological knowledge that was used for the model formulations. JO performed the modeling, data analysis and wrote the manuscript. GvS and AJBvB also contributed in calculations and critical review of the methods used in the analysis. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Acknowledgements
This work is supported by the graduate school VLAG and the IPOP program of Wageningen University.
This article has been published as part of BMC Systems Biology Volume 5 Supplement 1, 2011: Selected articles from the 4th International Conference on Computational Systems Biology (ISB 2010). The full contents of the supplement are available online at http://www.biomedcentral.com/17520509/5?issue=S1.
References

van Peij NNME, Gielkens MMC, de Vries RP, Visser J, de Graaff LH: The Transcriptional Activator XlnR Regulates Both Xylanolytic and Endoglucanase Gene Expression in Aspergillus niger.
Appl Environ. Microbiol 1998, 64:36153619. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

van der Veen D, Oliveira JM, van den Berg WAM, de Graaff LH: Varience components analysis reveals contribution of sample processing to transcript variation.
Appl Environ. Microbiol 2009, 75:24142422. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Sakamoto E, Iba H: Inferring a System of Differential Equations for a Gene Regulatory Network by using Genetic Programming. In Proc. Congress on Evolutionary Computation. IEEE Press; 2001:720726.

Tyson J, Othmer HG: The dynamics of feedback control circuits in biochemical pathways.

Lee TI, Rinaldi NJ, Robert F, Odom DT, BarJoseph Z, Gerber GK, Hannett NM, Harbison CT, Thompson CM, Simon I, et al.: Transcriptional regulatory networks in Saccharomyces cerevisiae.
Science 2002, 298:799804. PubMed Abstract  Publisher Full Text

Becskei A, Serrano L: Engineering stability in gene networks by autoregulation.
Nature 2000, 405:590593. PubMed Abstract  Publisher Full Text

Kepler TB, Elston TC: Stochasticity in transcriptional regulation: Origins, consequences, and mathematical representations.
Biophys J 2001, 81:31163136. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Simpson ML, Cox CD, Sayler GS: Frequency domain analysis of noise in autoregulated gene circuits.
Proc Natl. Acad. Sci 2003, 100:45514556. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Tao Y: Intrinsic and external noise in an autoregulatory genetic network.
J Theor. Biol 2004, 229:147156. PubMed Abstract  Publisher Full Text

Goutsias J, Kim S: Stochastic Transcriptional Regulatory Systems with Time Delays: A MeanField Approximation.

Sneppen K, Krishna S, Semsey S: Simplified Models of Biological Networks.
Annu Rev. Biophys 2010, 39:4359. PubMed Abstract  Publisher Full Text

Hirata H, Yoshiura S, Ohtsuka T, Bessho Y, Harada T, Yoshikawa K, Kageyama R: Oscillatory expression of the bHLH factor Hes1 regulated by a negative feedback loop.
Science 2002, 298:840843. PubMed Abstract  Publisher Full Text

Jensen M, Tiana G, Sneppen K: Sustained oscillations and time delays in gene expression of protein Hes1.
FEBS Lett 2003, 541:176177. PubMed Abstract  Publisher Full Text

Tiana G, Jensen M, Sneppen K: Time delay as a key to apoptosis induction in the p53 network.
Eur Phys 2002, 29:135140. Publisher Full Text

Tiana G, Krishna S, Pigolotti S, Jensen MH, Sneppen K: Oscillations and temporal signalling in cells.
Phys Biol 2007, 4:R1R17. PubMed Abstract  Publisher Full Text

Bliss RD, Painter PR, Marr AG: Role of feedback inhibition in stabilizing the classical operon.
J Theor. Biol 1982, 97:177193. PubMed Abstract  Publisher Full Text

Mendes P, Sha W, Ye K: Artificial gene networks for objective comparison of analysis algorithms.
Bioinformatics 2003, 19:122129. Publisher Full Text

de Jong H: Modeling and simulation of genetic regulatory systems: a literature review.
J Comput. Biol 2002, 9:67103. PubMed Abstract  Publisher Full Text

Isaacs FJ, Hasty J, Cantor CR, Collins JJ: Prediction and measurement of an autoregulatory genetic module.
Proc Natl. Acad. Sci 2003, 100:77147719. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

SayyedAhmad A, Tuncay K, Ortoleva PJ: Transcriptional regulatory network refinement and quantification through kinetic modeling, gene expression microarray data and information theory.
BMC Bioinformatics 2007, 8:20. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Agrawal S, Archer C, Schaffer DV: Computational Models of the Notch Network Elucidate Mechanisms of Contextdependent Signaling.
PLoS Comput Biol 2009, 5:e1000390. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

de Vries RP, V J, de Graaff LH: CreA modulates the XlnR induced expression on xylose of Aspergillus nigerAspergillus niger genes involved in xylan degradation.
Res Microbiol 1999, 150:281285. PubMed Abstract  Publisher Full Text

Yagil G, Yagil E: On the relation between effector concentration and the rate of induced enzyme synthesis.
Biophysical Journal 1971, 11:1127. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hill AV: The possible effect of the aggregation of the molecules of haemoglobin on its dissociation curves.

Polynikis A, Hogan SJ, di Bernardo M: Comparing different ODE modeling approaches of gene regulatory networks.
Journal of Theoretical Biology 2009, 261:511530. PubMed Abstract  Publisher Full Text

Hasper L: Function and mode of regulation of the transcriptional activator XlnR from Aspergillus. PhD thesis. Wageningen University, Microbiology Department; 2003.

Mahaffy JM: Genetic control models with diffusion and delays.
Math Biosci 1988, 90:519533. Publisher Full Text

Mahaffy JM, Jorgensen DA, Vanderheyden RL: Oscillations in a model of repression with external control.
J Math. Biol 1992, 30:669691. PubMed Abstract  Publisher Full Text

Mocek WT, Rudbicki R, Voit EO: Approximation of delays in biochemical systems.
Math Biosci 2005, 198:190216. PubMed Abstract  Publisher Full Text

Verdugo A, Rand R: Hopf bifurcation in a DDE model of gene expression.

Aro N, Pakula T, Penttila M: Transcriptional regulation of plant cell wall degradation by filamentous fungi.
FEMS Microbiol Rev 2005, 29:719739. PubMed Abstract  Publisher Full Text

Hasper AA, Dekkers E, van Mil M, van de Vondervoort PJI, de Graaff LH: EglC, A New Endoglucanase from Aspergillus niger with major activity towards Xyloglucan.
Appl. Environ. Microbiol 2002, 68:15561560. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Keseler I, ColladoVides J, GamaCastro S, Ingraham J, Paley S, Paulsen I, PeraltaGil M, Karp P: EcoCyc: a comprehensive database resource for Eschercichia coli.
Nucleic Acids Res 2005, 33:D334337. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Maithreye R, Sarkar RR, Parnaik VK, Sinha S: DelayInduced Transient Increase and Heterogeneity in Gene Expression in Negatively AutoRegulated Gene Circuits.
PLoS ONE 2008, 3:e2972. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Thomas R, Thieffry D, Kauffman M: Dynamical behaviour of biological regulatory networksI. Biological role of feedback loops and practical use of the concept of the loopcharacteristic state.
Bull. Math. Biol 1995, 57:247276. PubMed Abstract

Liu H, Lu JA, Lu J, Hill D: Structure identification of uncertain general complex dynamical networks with time delay.
Automatica 2009, 45:17991807. Publisher Full Text

Yu W, Lu J, Chen G, Duan Z, Zhou Q: Estimating Uncertain Delayed Genetic Regulatory Networks: An Adaptive Filtering Approach.

BalsaCanto E, Alonso A, Banga JR: Computational procedures for optimal experimental design in biological systems.
IET Syst. Biol 2008, 2:163172. PubMed Abstract  Publisher Full Text