Email updates

Keep up to date with the latest news and content from BMC Systems Biology and BioMed Central.

This article is part of the supplement: Selected articles from the 4th International Conference on Computational Systems Biology (ISB 2010)

Open Access Report

Modeling and analysis of the dynamic behavior of the XlnR regulon in Aspergillus niger

Jimmy Omony12*, Leo H de Graaff2, Gerrit van Straten1 and Anton J B van Boxtel1

Author Affiliations

1 Systems and Control group, Wageningen University, Wageningen, P.O. Box 17, 6700 AA, The Netherlands

2 Laboratory of Systems and Synthetic Biology, Wageningen University, Wageningen, Dreijenplein 10, 6703 HB, The Netherlands

For all author emails, please log on.

BMC Systems Biology 2011, 5(Suppl 1):S14  doi:10.1186/1752-0509-5-S1-S14

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1752-0509/5/S1/S14


Published:20 June 2011

© 2011 Omony et al; licensee BioMed Central Ltd.

This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

In this paper the dynamics of the transcription-translation 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 D-xylose 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 D-xylose 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 D-xylose 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 small-sized 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 bio-processes 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 D-xylose 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) [6-9], 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 bi-stability 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 transcription-translation 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 [18-21]. Many of these studies used the continuous-time 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 D-xylose 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 post-translational 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 D-xylose. 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).

thumbnailFigure 1. The XlnR regulon scheme The XlnR regulon induced by D-xylose 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 post-translational 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 = [z1,…,zn]Τ 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. zi is used to represent the time dependent variable zi(t) (where t ∈ [0, ∞)). Then the activating and repressing functions are given by

(1)

where ψ(zi, θi) = 1 – ψ+(zi, θi), θi is the gene specific half-saturation parameter and the positive number h represents the Hill coefficient. The regulation mechanism for each target gene i is captured by the function Ψ(zi, θ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

(2)

where ki1 = 1/θi for h = 1, xi - 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, ki1 - effective affinity constant for gene 1 activating gene i (i = 2,…,n), kis - maximum synthesis parameter for gene i, kid - first order degradation rate (or consumption rate) for gene i, x0 - vector of initial mRNA concentration, zi - concentration of translated protein from gene i, b = [b1,…,bn]Τ is the input matrix and u = [u1,…,un]Τ 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.

(3)

where ri - specific translation rate for gene i, ηi - degradation rate for protein i. The zi’s represent the target proteins in the scheme in Figure 1. At steady state the response rate and degradation rate balance, i.e. 12 ≈ … ≈ n ≈ 0. By setting the transcription rates i = 0 for all i in (2), it follows that

(4)

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 [27-30].

Let F : ℝ2n → ℝ2n be a set of smooth functions (with F = (F1,…,F2n)) that capture the XlnR regulon system dynamics. In this case we have F1 = 1,…,Fn = n, and Fn+1 = ż1,…,F2n = żn in (2) and (3) respectively. By definition, the Jacobian matrix is given by

(5)

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

(6)

where x = [x1, x2, x3]Τ and z = [z1, z2, z3]Τ . The corresponding steady states of the vectors and can be computed, accordingly. Using expressions (2) and (3) in (5) we obtain

(7)

where

(8)

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

(9)

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 n-dimensional network system. The generalization for the eigenvalue spectra using a similar analysis as above leads to the expression

(10)

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 auto-regulating 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 (x1) has to be modified accordingly. The adapted equation is given by

(11)

where

(12)

is the repressor Hill function and CA - quantitative activity state for CreA, kA - inverse of the Hill constant of CreA. The term τ > 0 represents a time delay in the feedback loop. The sets S1 = {j | j = 1,…,m} and S2 = {l | l = m + 1,…,n – 1} where S1S2 = {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 S1 and S2, respectively. The effect of the D-xylose 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.

(13)

(14)

The extracts from the denominator functions are given by (15) and (16), respectively.

(15)

(16)

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 kRL = kAL = 1 were used). We consider the sets S1 and S2 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).

(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 F1 with respect to the variable of interest within each of the sets S1 and S2, we then have the more compact expression

(18)

This term corresponds to the repressing proteins, and the terms given by

(19)

for the activating proteins. The parameters kAL and kRL 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 auto-regulation 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 auto-regulation effect of the xlnR gene. The computed eigenvalue spectrum from (7) using the expression is given by

(20)

where

(21)

(22)

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 < –(η1k1d)2/4r1 for r1 > 0 is satisfied. This condition on ω1 signifies a contribution from a NFB loop in the XlnR regulon network. Notice that the expression (η1k1d)2 > 0 for all values of η1 and k1d. 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

(23)

Solving the inequality (23) for ω1 leads to the condition ω1 <η1k1d/r1. 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.

(24)

or after working out becomes ω1 = η1k1d/r1. 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 D-xylose. The pulse perturbation takes place at time t = 0. During fermentation, the D-xylose is consumed and the D-xylose concentration follows the expression u(t) = u(0)(1/(β + eKt)), where u(t) ≡ [D-xylose] and β > 0, with K = 0.3 and u(0) = 50 mM as the initial D-xylose concentration. The parameters used for the simulation are: b1 = 1, ρ1 = 2e – 3, ρ2 = 2.5e – 3, ρ3 = 1e – 3, k1d = 0.5, k2d = 0.4, k3d = 0.3, k2s = 5, k3s = 6, k21 = 0.1, k31 = 0.1, r1 = r2 = r3 = 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 ri, the gene synthesis coefficient kis 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 kid 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 D-xylose 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/k1d ≈ 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 (D-xylose).

thumbnailFigure 2. D-xylose consumption, gene expression, protein abundance and phase plots: without feedback (A): The simulated trajectory for D-xylose 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, xi. (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 switch-like 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 kRL = 1 and kAL = 1 and the lumped synthesis parameter from (11) chosen as kls = 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.

thumbnailFigure 3. D-xylose consumption, gene expression, protein abundance and phase plots: with feedback (A): The simulated trajectory for D-xylose. (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 time-dynamics 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.

thumbnailFigure 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 kls (Figure 5).

thumbnailFigure 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 (kls = 1) and strong (kls = 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.

thumbnailFigure 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)).

thumbnailFigure 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 non-existence 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 D-xylose uptake in fermentation experiments. The consumption of D-xylose also indirectly controls the regulation of the target genes and therewith the breakdown of sugars.

Simulations showed that the dynamics of the D-xylose 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 D-xylose. 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 Balsa-Canto 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 in-depth 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 wet-lab 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 small-sized 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 D-xylose 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/1752-0509/5?issue=S1.

References

  1. 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:3615-3619. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  2. 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:2414-2422. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  3. 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:720-726. OpenURL

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

    Prog Theor. Biol 1978, 5:1-62. OpenURL

  5. Lee TI, Rinaldi NJ, Robert F, Odom DT, Bar-Joseph Z, Gerber GK, Hannett NM, Harbison CT, Thompson CM, Simon I, et al.: Transcriptional regulatory networks in Saccharomyces cerevisiae.

    Science 2002, 298:799-804. PubMed Abstract | Publisher Full Text OpenURL

  6. Becskei A, Serrano L: Engineering stability in gene networks by autoregulation.

    Nature 2000, 405:590-593. PubMed Abstract | Publisher Full Text OpenURL

  7. Kepler TB, Elston TC: Stochasticity in transcriptional regulation: Origins, consequences, and mathematical representations.

    Biophys J 2001, 81:3116-3136. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  8. Simpson ML, Cox CD, Sayler GS: Frequency domain analysis of noise in autoregulated gene circuits.

    Proc Natl. Acad. Sci 2003, 100:4551-4556. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  9. Tao Y: Intrinsic and external noise in an auto-regulatory genetic network.

    J Theor. Biol 2004, 229:147-156. PubMed Abstract | Publisher Full Text OpenURL

  10. Goutsias J, Kim S: Stochastic Transcriptional Regulatory Systems with Time Delays: A Mean-Field Approximation.

    Phys Biol 2006, 13:1049-1076. OpenURL

  11. Sneppen K, Krishna S, Semsey S: Simplified Models of Biological Networks.

    Annu Rev. Biophys 2010, 39:43-59. PubMed Abstract | Publisher Full Text OpenURL

  12. 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:840-843. PubMed Abstract | Publisher Full Text OpenURL

  13. Jensen M, Tiana G, Sneppen K: Sustained oscillations and time delays in gene expression of protein Hes1.

    FEBS Lett 2003, 541:176-177. PubMed Abstract | Publisher Full Text OpenURL

  14. Tiana G, Jensen M, Sneppen K: Time delay as a key to apoptosis induction in the p53 network.

    Eur Phys 2002, 29:135-140. Publisher Full Text OpenURL

  15. Tiana G, Krishna S, Pigolotti S, Jensen MH, Sneppen K: Oscillations and temporal signalling in cells.

    Phys Biol 2007, 4:R1-R17. PubMed Abstract | Publisher Full Text OpenURL

  16. Bliss RD, Painter PR, Marr AG: Role of feedback inhibition in stabilizing the classical operon.

    J Theor. Biol 1982, 97:177-193. PubMed Abstract | Publisher Full Text OpenURL

  17. Mendes P, Sha W, Ye K: Artificial gene networks for objective comparison of analysis algorithms.

    Bioinformatics 2003, 19:122-129. Publisher Full Text OpenURL

  18. de Jong H: Modeling and simulation of genetic regulatory systems: a literature review.

    J Comput. Biol 2002, 9:67-103. PubMed Abstract | Publisher Full Text OpenURL

  19. Isaacs FJ, Hasty J, Cantor CR, Collins JJ: Prediction and measurement of an autoregulatory genetic module.

    Proc Natl. Acad. Sci 2003, 100:7714-7719. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  20. Sayyed-Ahmad 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 OpenURL

  21. Agrawal S, Archer C, Schaffer DV: Computational Models of the Notch Network Elucidate Mechanisms of Context-dependent Signaling.

    PLoS Comput Biol 2009, 5:e1000390. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  22. 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:281-285. PubMed Abstract | Publisher Full Text OpenURL

  23. Yagil G, Yagil E: On the relation between effector concentration and the rate of induced enzyme synthesis.

    Biophysical Journal 1971, 11:11-27. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

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

    J Physiol 1910, 40:4-7. OpenURL

  25. Polynikis A, Hogan SJ, di Bernardo M: Comparing different ODE modeling approaches of gene regulatory networks.

    Journal of Theoretical Biology 2009, 261:511-530. PubMed Abstract | Publisher Full Text OpenURL

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

  27. Mahaffy JM: Genetic control models with diffusion and delays.

    Math Biosci 1988, 90:519-533. Publisher Full Text OpenURL

  28. Mahaffy JM, Jorgensen DA, Vanderheyden RL: Oscillations in a model of repression with external control.

    J Math. Biol 1992, 30:669-691. PubMed Abstract | Publisher Full Text OpenURL

  29. Mocek WT, Rudbicki R, Voit EO: Approximation of delays in biochemical systems.

    Math Biosci 2005, 198:190-216. PubMed Abstract | Publisher Full Text OpenURL

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

    Science Direct 2008, 13:235-242. OpenURL

  31. Aro N, Pakula T, Penttila M: Transcriptional regulation of plant cell wall degradation by filamentous fungi.

    FEMS Microbiol Rev 2005, 29:719-739. PubMed Abstract | Publisher Full Text OpenURL

  32. 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:1556-1560. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  33. Keseler I, Collado-Vides J, Gama-Castro S, Ingraham J, Paley S, Paulsen I, Peralta-Gil M, Karp P: EcoCyc: a comprehensive database resource for Eschercichia coli.

    Nucleic Acids Res 2005, 33:D334-337. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  34. Maithreye R, Sarkar RR, Parnaik VK, Sinha S: Delay-Induced Transient Increase and Heterogeneity in Gene Expression in Negatively Auto-Regulated Gene Circuits.

    PLoS ONE 2008, 3:e2972. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  35. Thomas R, Thieffry D, Kauffman M: Dynamical behaviour of biological regulatory networks-I. Biological role of feedback loops and practical use of the concept of the loop-characteristic state.

    Bull. Math. Biol 1995, 57:247-276. PubMed Abstract OpenURL

  36. Thomas R, d’Ari R: Biological Feedback.

    CRC, in press. OpenURL

  37. Liu H, Lu JA, Lu J, Hill D: Structure identification of uncertain general complex dynamical networks with time delay.

    Automatica 2009, 45:1799-1807. Publisher Full Text OpenURL

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

    IEEE Trans. Automat. Contr 2009, 54:792-897. OpenURL

  39. Balsa-Canto E, Alonso A, Banga JR: Computational procedures for optimal experimental design in biological systems.

    IET Syst. Biol 2008, 2:163-172. PubMed Abstract | Publisher Full Text OpenURL