- Report
- Open access
- Published:

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

*BMC Systems Biology*
**volume 5**, Article number: S14 (2011)

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

### 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 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 Ψ(*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 2*n* 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* = (*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 *n*-dimensional 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 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 (*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 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}.

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 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

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}/4*r*_{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 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/(*β* + *e*^{Kt})), 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: *b*_{1} = 1, *ρ*_{1} = 2*e* – 3, *ρ*_{2} = 2.5*e* – 3, *ρ*_{3} = 1*e* – 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 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 *τ*_{R 1} = 1/*k*_{1d} ≈ 2 hours is noticed for the master regulator and for the target genes, *τ*_{R 1} <*τ*_{R 2}, *τ*_{R 3}. 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).

### 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 *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.

#### 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.

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

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.

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

## 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.

## 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: 3615-3619.

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. 10.1128/AEM.02270-08.

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

Tyson J, Othmer HG: The dynamics of feedback control circuits in biochemical pathways. Prog Theor. Biol. 1978, 5: 1-62.

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. 10.1126/science.1075090.

Becskei A, Serrano L: Engineering stability in gene networks by autoregulation. Nature. 2000, 405: 590-593. 10.1038/35014651.

Kepler TB, Elston TC: Stochasticity in transcriptional regulation: Origins, consequences, and mathematical representations. Biophys J. 2001, 81: 3116-3136. 10.1016/S0006-3495(01)75949-8.

Simpson ML, Cox CD, Sayler GS: Frequency domain analysis of noise in autoregulated gene circuits. Proc Natl. Acad. Sci. 2003, 100: 4551-4556. 10.1073/pnas.0736140100.

Tao Y: Intrinsic and external noise in an auto-regulatory genetic network. J Theor. Biol. 2004, 229: 147-156. 10.1016/j.jtbi.2004.03.011.

Goutsias J, Kim S: Stochastic Transcriptional Regulatory Systems with Time Delays: A Mean-Field Approximation. Phys Biol. 2006, 13: 1049-1076.

Sneppen K, Krishna S, Semsey S: Simplified Models of Biological Networks. Annu Rev. Biophys. 2010, 39: 43-59. 10.1146/annurev.biophys.093008.131241.

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. 10.1126/science.1074560.

Jensen M, Tiana G, Sneppen K: Sustained oscillations and time delays in gene expression of protein Hes1. FEBS Lett. 2003, 541: 176-177. 10.1016/S0014-5793(03)00279-5.

Tiana G, Jensen M, Sneppen K: Time delay as a key to apoptosis induction in the p53 network. Eur Phys. 2002, 29: 135-140. 10.1140/epjb/e2002-00271-1.

Tiana G, Krishna S, Pigolotti S, Jensen MH, Sneppen K: Oscillations and temporal signalling in cells. Phys Biol. 2007, 4: R1-R17. 10.1088/1478-3975/4/2/R01.

Bliss RD, Painter PR, Marr AG: Role of feedback inhibition in stabilizing the classical operon. J Theor. Biol. 1982, 97: 177-193. 10.1016/0022-5193(82)90098-4.

Mendes P, Sha W, Ye K: Artificial gene networks for objective comparison of analysis algorithms. Bioinformatics. 2003, 19: 122-129. 10.1093/bioinformatics/btg1069.

de Jong H: Modeling and simulation of genetic regulatory systems: a literature review. J Comput. Biol. 2002, 9: 67-103. 10.1089/10665270252833208.

Isaacs FJ, Hasty J, Cantor CR, Collins JJ: Prediction and measurement of an autoregulatory genetic module. Proc Natl. Acad. Sci. 2003, 100: 7714-7719. 10.1073/pnas.1332628100.

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-10.1186/1471-2105-8-20.

Agrawal S, Archer C, Schaffer DV: Computational Models of the Notch Network Elucidate Mechanisms of Context-dependent Signaling. PLoS Comput Biol. 2009, 5: e1000390-10.1371/journal.pcbi.1000390.

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. 10.1016/S0923-2508(99)80053-9.

Yagil G, Yagil E: On the relation between effector concentration and the rate of induced enzyme synthesis. Biophysical Journal. 1971, 11: 11-27. 10.1016/S0006-3495(71)86192-1.

Hill AV: The possible effect of the aggregation of the molecules of haemoglobin on its dissociation curves. J Physiol. 1910, 40: 4-7.

Polynikis A, Hogan SJ, di Bernardo M: Comparing different ODE modeling approaches of gene regulatory networks. Journal of Theoretical Biology. 2009, 261: 511-530. 10.1016/j.jtbi.2009.07.040.

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

Mahaffy JM: Genetic control models with diffusion and delays. Math Biosci. 1988, 90: 519-533. 10.1016/0025-5564(88)90081-8.

Mahaffy JM, Jorgensen DA, Vanderheyden RL: Oscillations in a model of repression with external control. J Math. Biol. 1992, 30: 669-691. 10.1007/BF00173263.

Mocek WT, Rudbicki R, Voit EO: Approximation of delays in biochemical systems. Math Biosci. 2005, 198: 190-216. 10.1016/j.mbs.2005.08.001.

Verdugo A, Rand R: Hopf bifurcation in a DDE model of gene expression. Science Direct. 2008, 13: 235-242.

Aro N, Pakula T, Penttila M: Transcriptional regulation of plant cell wall degradation by filamentous fungi. FEMS Microbiol Rev. 2005, 29: 719-739. 10.1016/j.femsre.2004.11.006.

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. 10.1128/AEM.68.4.1556-1560.2002.

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. 10.1093/nar/gki108.

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-10.1371/journal.pone.0002972.

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.

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

Liu H, Lu JA, Lu J, Hill D: Structure identification of uncertain general complex dynamical networks with time delay. Automatica. 2009, 45: 1799-1807. 10.1016/j.automatica.2009.03.022.

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.

Balsa-Canto E, Alonso A, Banga JR: Computational procedures for optimal experimental design in biological systems. IET Syst. Biol. 2008, 2: 163-172. 10.1049/iet-syb:20070069.

## 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.

## Author information

### Authors and Affiliations

### Corresponding author

## Additional information

### 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.

## Rights and permissions

This article is published under license to 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.

## About this article

### Cite this article

Omony, J., de Graaff, L.H., van Straten, G. *et al.* Modeling and analysis of the dynamic behavior of the XlnR regulon in Aspergillus niger.
*BMC Syst Biol* **5**
(Suppl 1), S14 (2011). https://doi.org/10.1186/1752-0509-5-S1-S14

Published:

DOI: https://doi.org/10.1186/1752-0509-5-S1-S14