Email updates

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

This article is part of the supplement: Ninth International Conference on Bioinformatics (InCoB2010): Bioinformatics

Open Access Proceedings

Sensitivity analysis of dynamic biological systems with time-delays

Wu Hsiung Wu1, Feng Sheng Wang2* and Maw Shang Chang1

Author Affiliations

1 Department of Computer Science and Information Engineering, National Chung Cheng University, Chiayi 62102, Taiwan

2 Department of Chemical Engineering, National Chung Cheng University, Chiayi 62102, Taiwan

For all author emails, please log on.

BMC Bioinformatics 2010, 11(Suppl 7):S12  doi:10.1186/1471-2105-11-S7-S12

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2105/11/S7/S12


Published:15 October 2010

© 2010 Wu 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

Mathematical modeling has been applied to the study and analysis of complex biological systems for a long time. Some processes in biological systems, such as the gene expression and feedback control in signal transduction networks, involve a time delay. These systems are represented as delay differential equation (DDE) models. Numerical sensitivity analysis of a DDE model by the direct method requires the solutions of model and sensitivity equations with time-delays. The major effort is the computation of Jacobian matrix when computing the solution of sensitivity equations. The computation of partial derivatives of complex equations either by the analytic method or by symbolic manipulation is time consuming, inconvenient, and prone to introduce human errors. To address this problem, an automatic approach to obtain the derivatives of complex functions efficiently and accurately is necessary.

Results

We have proposed an efficient algorithm with an adaptive step size control to compute the solution and dynamic sensitivities of biological systems described by ordinal differential equations (ODEs). The adaptive direct-decoupled algorithm is extended to solve the solution and dynamic sensitivities of time-delay systems describing by DDEs. To save the human effort and avoid the human errors in the computation of partial derivatives, an automatic differentiation technique is embedded in the extended algorithm to evaluate the Jacobian matrix. The extended algorithm is implemented and applied to two realistic models with time-delays: the cardiovascular control system and the TNF-α signal transduction network. The results show that the extended algorithm is a good tool for dynamic sensitivity analysis on DDE models with less user intervention.

Conclusions

By comparing with direct-coupled methods in theory, the extended algorithm is efficient, accurate, and easy to use for end users without programming background to do dynamic sensitivity analysis on complex biological systems with time-delays.

Background

Mathematical modeling has been applied to the study and analysis of complex biological systems for a long time. Many mathematical models for dynamic biological systems are formulated as nonlinear ordinary differential equations (ODEs). Some processes in biological systems, such as the gene expression and feedback control in signal transduction networks, involve a time delay. These systems are represented as delay differential equation (DDE) models. Many DDE models have been proposed in the last decade [1-3]. Bocharov et al. [4] reviewed various applications of DDE models in population dynamics, epidemiology, immunology, neural networks, and cell kinetics. Sensitivity analysis can shed light on the dynamic behavior of biological systems and assist the modeling process by identifying the critical parameters that determine the essential behavior of the system. Numerical sensitivity analysis of a DDE model by the direct method requires to obtain the solutions of model and sensitivity equations with time-delays. To do dynamic sensitivity analysis on DDE models, an efficient and accurate approach to compute the solution of DDEs is the basic requirement. Many sophisticated DDE solvers are available recently [5-12]. The major effort is the computation of Jacobian matrix when computing the solution of sensitivity equations. The computation of partial derivatives of complicated equations either by the analytic method or by symbolic manipulation is time consuming, inconvenient, and prone to introduce human errors. To surmount this problem, an automatic approach to obtain the derivatives of complex functions efficiently and accurately is necessary.

Dynamic sensitivity analysis is an important tool for assessing dynamic behavior of biological systems. The common used approach for sensitivity analysis is the numerical differentiation by finite difference approximation. The main drawback of this approach is that the accuracy is hard to analyze. Due to the efficiency and accuracy, a variety of direct methods are proposed [13-15]. Rihan [16] derives a general theory for sensitivity analysis of DDE models by using adjoint equations and direct methods to estimate the sensitivity equations with variable and constant parameters, respectively. The kinetic preprocessor (KPP) numerical library is a comprehensive set of software tools for direct and adjoint sensitivity analysis [17]. An-other approach which can be used to evaluate sensitivity equations is automatic differentiation. Automatic differentiation is a non-approximate differentiation technique that permits the fast and accurate evaluation of partial derivatives in Jacobian matrix. The values for the derivatives obtained by automatic differentiation are exact if we do not take account of the unavoidable round-off error due to the finite precision arithmetic of a computer. The theoretical exactness of the automatic differentiation comes from the fact that it uses the same rules of differentiation as in differential calculus, but these rules are applied to an algorithmic evaluation of the function rather than to a formula. The implementation of automatic differentiation can be divided into two different classes: source code transformation and operator overloading. The most widely used source code transformation program is ADIFOR 2.0 [18]. This program, like as the preprocessor of a compiler, accepts and analyzes Fortran 77 source code and produces code to evaluate the derivatives of the function in Fortran 77 standard. The output code is optimized by eliminating unnecessary arithmetic operations and temporary variables and then compiled with a standard compiler into an object code that can simultaneously evaluate derivatives and function values. Hwang et al. [19] demonstrated that ADIFOR is a powerful tool for ODE models from three aspects of performance: accuracy, efficiency, and scaled capability. Griewank et al. [20] developed an open-source code, automatic differentiation by overloading in C++ (ADOL-C), for the automatic differentiation of C and C++ programs. The implementation of ADOL-C utilizes the operator overloading capability of C/C++ compilers that accept user-defined data types, operators and functions. The implementation of either the source code transformation or the operator overloading is a compile-time solution. It allows one to generate derivatives from complicated existing code or user-provided model equations that expressed by user-defined data types, operators and functions. These available codes are implemented for ODEs and is suitable for users with programming background.

We have proposed an efficient algorithm with an adaptive step size control, called adaptive modified collocation method (AMCM), to compute the solution and dynamic sensitivities of biological systems described by ODEs [21]. The algorithm is extended to solve the solution and dynamic sensitivities of time-delay systems described by DDEs in this paper and named as extended adaptive modified collocation method (EAMCM). The EAMCM is implemented as a user-friendly program that facilitates the dynamic sensitivity analysis of DDE models through the implementation of adaptive direct-decoupled method and automatic differentiation. EAMCM requires the user to supply only the model equations at run-time in a form of mathematical expression to compute dynamic sensitivities of DDE models. The evaluation of sensitivity equations is done automatically by automatic differentiation technique along with the inevitable evaluation of model equations. By combining the adaptive direct-decoupled AMCM algorithm with automatic differentiation technique and extending its usage to DDE models, the extended algorithm, EAMCM, is efficient, accurate, and easy to use for end users without programming background to do dynamic sensitivity analysis on complex biological systems with time-delays.

To evaluate the applicability of the extended algorithm, it is applied to two realistic models with time-delays: the cardiovascular control system and the TNF-α signal transduction network. The first DDE model for human cardiovascular control system was developed by Fowler et al. [22] to explore the interactions between the heart rate and blood pressure under the baroreflex control. The time delay arises from the slow activity of sympathetic nervous system. Sensitivity analysis is applied to this DDE model through the EAMCM program to identify the key parameters that could provide useful diagnostic information about the state of the cardiovascular system. The second DDE model for TNF-α signal transduction network built by Rangamani and Sirovich [23] considers both the NF-κB-mediated survival pathway and the caspase-mediated apoptosis pathway simultaneously. Due to the delay of translocation of NF-κB to the nucleus, the transcription processes of cIAP and IκB in the NF-κB-mediated survival pathway were represented by DDEs. The EAMCM is applied to this model to analyze its dynamic sensitivities and decipher the relationship between the NF-κB-mediated survival pathway and the caspase-mediated apoptosis pathway.

Results and discussion

Cardiovascular disease is the major cause of human death. A detailed understanding of cardiovascular systems is important for the diagnosis of cardiovascular disease, ultimately leading to improved treatment. The EAMCM program can be used to do dynamic sensitivity analysis on the cardiovascular control system to investigate the hemodynamics and regulation control of cardiovascular systems.

Apoptosis is central to the development of cancer. In the recent years, the prevalent model to explain why cancer therapies fail has been that cell resistant to or inhibition of apoptosis [24]. So now, the new treatment goal is how to control apoptosis that brings on cell death of the cancer cells. NF-κB is a transcription factor family that activating numerous genes that are related to cell survival pathways. Most commonly, NF-κB activation inhibits apoptosis pathways, as evidenced by several knockout mouse models [25,26]. Based on these findings, the goal to design more effective cancer therapies can be initiated by apoptosis induction and inhibition of NF-κB. Many mathematical models describing the dynamics of apoptosis and NF-κB pathways have been published in last decade [27-31]. Most of the models to date have concentrated on only one of signaling pathways and do not consider the delayed transcription processes. The EAMCM program is applied to a TNF-α-induced signaling network considering both signaling pathways simultaneously to investigate how these two pathways work together to regulate cell fate in response to TNF-α.

Cardiovascular control system

Mathematical models are useful to investigate the hemodynamics and regulation control of cardiovascular systems. In general, cardiovascular models consist of two major systems: the hemodynamic system and the autonomic nerve control system. The hemodynamic system is a systemic circulatory blood distribution network to deliver oxygen, nutrients, and hormones to cells and remove carbon dioxide and metabolic wastes. The left ventricle pumps blood to arteries, capillaries, veins, and back to the heart. The blood hemodynamics of this circulation can be represented by the relationship between blood pressure and heart rate in the cardiovascular system. The control of the blood pressure is crucial to human health due to that the blood pressure is the energy, generated by the heart, to push blood around the body. Although the endogenous regulation of arterial pressure is not completely understood, there are evidences that baroreceptors are important for minimizing changes in blood pressure. Animal studies have shown that blood pressure is much more variable if the influence of baroreceptors is removed [32,33]. Baroreceptors detect changes in arterial pressure and send signals ultimately to the medulla of the brain stem. The medulla, by way of the autonomic nerve control system, adjusts the mean arterial pressure by altering the heart rate and the total peripheral resistance. The autonomic nerve control system includes the sympathetic and parasympathetic nervous systems. When blood pressure starting to fall, the baroreceptor stimulation decreasing and the reflex response causes the peripheral resistance increasing and the heart to beat faster and harder by slow-acting sympathetic nerves. This negative-feedback mechanism largely restores the blood pressure. Conversely, if blood pressure increases, stimulation of baroreceptors gives rise to nerve impulses which run to the brain and slow down the heartbeat through the fast activity in the parasympathetic pathway.

Fowler et al. [22] developed a lumped DDE model of the integrated cardiovascular system combined with a baroreflex feedback control of blood pressure to describe the interactions between heart rate, blood pressure, and respiration. This DDE model is a simple model without considering the pulmonary part of the cardiovascular system and is derived from the model introduced by Ottesen [34] by adding an intrinsically controlled heart rate and baroreflex control of peripheral resistance. This simple model consists of only 2 delay differential equations which include 16 parameters and is expressed as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M1">View MathML</a>

(1)

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M2">View MathML</a>

(2)

where x1 is the heart rate, x2 is the arterial pressure, and t* = t/τ is the dimensionless time. The functions g1 and g2 are defined by

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M3">View MathML</a>

where x2(t* - 1) and x2 (t*) denote the blood pressure with and without a time delay, respectively. The values and definitions of system parameters are given in Table 1. The state variables are the heart rate and blood pressure. A sinusoidal function is added to the model equations to mimic respiration. The Hill function gi with an exponent n is used to model the baroreflex feedback control of heart rate. This simple model has shown to be able to reproduce many of the empirical observations such as respiratory sinus arrhythmia (RSA), Mayer waves, and synchronization [35].

Table 1. The value and definition of system parameters

The EAMCM program is used to do sensitivity analysis on the lumped cardiovascular model. The non-constant exponent of Hill function and sinusoidal functions in differential model equations complicate the evaluation of Jacobian matrix for computing the solution of sensitivity equations. By the help of automatic differentiation embedded in the EAMCM program, user can provide the model equations only at run-time for solving the dynamic sensitivities of the cardiovascular system. The dynamic sensitivities of heart rate and blood pressure with respect to all system parameters and initial conditions are computed. All relative parameter sensitivities are presented by 100% stacked column chart and shown in Figure 1. It is easy to find which parameter makes more effects on heart rate and blood pressure than the others from Figure 1. The values of top five sensitivities for the heart rate and blood pressure are shown in Table 2. The uncontrolled average arterial blood pressure (p0), breathing rate (fr), sympathetic delay (τ), sympathetic control of heart rate (β), and strength of vagal tone (ν) are identified as the key parameters for the control of heart rate and blood pressure. The relative sensitivities of heart rate and blood pressure with respect to the uncontrolled average arterial blood pressure are shown in Figure 2. The dynamic sensitivities of heart rate with respect to p0 oscillate symmetrically between positive and negative values. This result indicates that the uncontrolled average arterial blood pressure amplifies the variation of heart rate. In contrast, the dynamic sensitivities of blood pressure with respect to p0 oscillate but are all positive. This means that an increase of the uncontrolled average arterial blood pressure shifts the blood pressure to a higher value but does not change the variation of blood pressure. As shown in Table 2, the effect of uncontrolled average arterial blood pressure on the variation of average heart rate is tenfold larger than the variation of blood pressure. There is evidence that the slow-acting sympathetic nerves and the fast-acting vagal nerves compete with each other to increase and decrease the heart rate, respectively [36]. The relative sensitivities of heart rate with respect to parameters for slow sympathetic control (β) and fast vagal control (ν) are investigated and shown in Figure 3. Figure 3 shows the sympathovagal balance in physiology and both sympathetic control and vagal control amplify the variation of heart rate. The relative sensitivities of blood pressure with respect to parameters for slow sympathetic control (β) and fast vagal control (ν) are shown in Figure 4. The slow-acting sympathetic control upregulates the blood pressure, but does not change its variation. The relative sensitivity of blood pressure with respect to the sympathetic control is positive over the time. In contrast, the fast-acting vagal control downregulates the blood pressure and has a negative relative sensitivity over the time.

thumbnailFigure 1. Stacked 100% column chart for individual state variables. Each column in the stack column chart shows all relative parameter sensitivities for a state variable. The proportion of a parameter sensitivity to the total sensitivity for a state variable is displayed as a color area in each column. The values of time-averaged relative parameter sensitivities are used as the data.

Table 2. The ranking of relative sensitivities of heart rate and blood pressure

thumbnailFigure 2. The relative sensitivities of heart rate and blood pressure with respect to p0. The relative sensitivities of heart rate and blood pressure with respect to the uncontrolled average arterial blood pressure. The time is in dimensionless scale.

thumbnailFigure 3. The relative sensitivities of heart rate with respect to β and ν. The relative sensitivities of heart rate with respect to parameters for slow sympathetic control and fast vagal control. The time is in dimensionless scale.

thumbnailFigure 4. The relative sensitivities of blood pressure with respect to β and ν. The relative sensitivities of blood pressure with respect to parameters for slow sympathetic control and fast vagal control. The time is in dimensionless scale.

Apoptosis signal network

The TNF-α signal transduction network was developed by Rangamani and Sirovich [23], which considers both the NF-κB-mediated survival pathway and the caspase-mediated apoptosis pathway simultaneously. These two pathways involve 31 species in 19 reactions and the schematic diagram is shown in Figure 5. The formation of this network involves binding reactions between ligand and death receptor, protein-protein reactions, enzymatic reactions, translocations, and transcription processes. The network is induced by ligation of TNF-α to the cell surface receptor TNFR1. The ligation of TNFR1 by TNF results in the recruitment of the adaptor proteins such as TNFR-associated death domain (TRADD), TNFR-associated factor 2 (TRAF2), receptor-interacting protein 1 (RIP1), and possibly other yet unidentified proteins to form the early complex. In the NF-κB-mediated transcription pathway, the inactive inhibitor kappa B kinase (IKK) binds to the early complex leads to the activation of IKK, IκB phosphorylation, and release of NF-κB. The free NF-κB translocates to the nucleus, binds to DNA, and leads to the transcription of IκB and cellular inhibitor of apoptosis protein (cIAP) that protects cells from TNF-induced apoptosis by binding to activated caspase-3 [37]. The newly synthesized free IκB enters the nucleus and binds to nuclear NF-κB and this complex is exported to the cytoplasm [38]. This NF-κB-IκB complex is the target for IκB phosphorylation by active IKK and the liberating NF-κB will translocate to the nucleus again. Nelson et al. [39] proposed that this oscillatory feedback behavior of NF-κB regulates the expression of cIAP. In the caspase-mediated apoptosis pathway, TRADD, RIP1, and TRAF2 are dissociated from TNFR1 and recruit Fas Associated Death Domain (FADD) and caspase-8 to form a protein complex called the death-inducing signaling complex (DISC) [40]. As a result of DISC formation, caspase-8 is cleaved at the DISC resulting in the activation of caspase-8. The activated caspase-8 in turn activates effector caspases, such as caspase-3, causing the cell to undergo apoptosis.

thumbnailFigure 5. Schematic diagram of TNF-α signal transduction network. The solid lines indicate reversible reactions and the dash-dot lines denote irreversible reactions. The dash lines indicate the delayed transcription processes. The reactions and components of the survival pathway are shown in green. The reactions and components of the apoptotic pathway are shown in blue. The boxes with red border denote the components with nonzero initial value in the network.

The transcription processes of cIAP and IκB due to the translocation of NF-κB to the nucleus are represented as delayed reactions. The delay time used for transcription is 20 minutes as suggested by Sung and Simon [41]. Based on material balance, this model consists of 31 delay differential equations which include 29 parameters. The state variables are the concentration of the molecules in the survival and apoptosis pathways and the input variable is the concentration of TNF-α that stimulates the signal transduction pathways. The output variable is the concentration of fragmented DNA, which can be used as a marker for apoptosis. The fragmented DNA is defined as the fraction of DNA sites that have been attacked by the activity of effector caspase. The set of delay differential equations, all of the relevant definitions of variables, and parameters appearing in the DDE model, together with the nominal values can be found in Additional file 1. The reason for representing the model equations here and not just referring to the article by Rangamani et al. [23] is that some of the parameters and state variables have different names as that in the original model.

Additional file 1. TNF-α signal transduction model. This file includes the set of delay differential equations of the TNF-α signal transduction model, all of the relevant definitions of state variables, and the nominal values of parameters appearing in the delay differential equations.

Format: PDF Size: 49KB Download file

This file can be viewed with: Adobe Acrobat ReaderOpen Data

The EAMCM program is applied to the TNF-α signal transduction model using the initial conditions as described in Table S1 of Additional file 1. All dynamic sensitivities with respect to 29 parameters and 31 initial conditions are computed simultaneously without any difficulty. All time-averaged semi-relative parameter sensitivities for each state variable are shown in Figure 6. Most of the time-averaged semi-relative parameter sensitivities for each state variable are too small compared with the largest and can be ignored. It is easy to find from Figure 6 that only some few parameter sensitivities get significant percentage of the total sensitivity for each state variable.

thumbnailFigure 6. Stacked 100% column chart for individual state variables. Each column in the stack column chart shows all semi-relative parameter sensitivities for a state variable. The proportion of a parameter sensitivity to the total sensitivity for a state variable is displayed as a color area in each column. The values of time-averaged semi-relative parameter sensitivities are used as the data.

The dynamic sensitivity profiles for all species with respect to k9, the rate constant of the formation of survival complex, are nearly identical to that with respect to x10(0), the initial value of IKK (data not shown here). This is not surprising because the kinetic order is set to one for each flux in the model. So each relative effect on the output with respect to the rate constant is the same as that with respect to the initial concentration of the corresponding species. The same situation can be found for each pair of k15 /x17(0), k18/x20(0), etc. In the following, we analyze the dynamic parameter sensitivities only, because the same results for the corresponding dynamic initial sensitivities can be found from the dynamic parameter sensitivities.

All of the dynamic sensitivities with respect to k9, the rate constant of the formation of survival complex, and k15 , the rate constant of the formation of death complex, are symmetric with respect to the time axis. This means that if we have plotted the sensitivity profile of a species with respect to k9, the corresponding sensitivity profile with respect to k15 can be obtained simply by reflecting about the time axis. To elucidate the effects of IKK (x10) and FADD (x17) on the oscillatory behavior of NF-κB (x16) and IκB (x31) in the survival pathway, the semi-relative sensitivities of NF-κB and IκB with respect to the rate constants k9 and k15 , of the formation of survival complex and death complex are shown in Figure 7. We observe, from Figure 7, the negative regulation of oscillatory behavior of NF-κB and IκB when the rate constant (k15) or the initial concentration of FADD increases. The reverse effect is seen by increasing the rate constant (k9) or the initial concentration of IKK. We also investigate the responses of the apoptosis pathway to the variances of FADD and IKK. The activation of effector caspase-3 (x23) is the finial reaction of the apoptosis pathway, so the concentration of the active caspase-3 (x25) can be used as the response of the apoptosis pathway. The semi-relative sensitivities of active caspase-3 with respect to k9 and k15 are shown in Figure 7. The negative values of the semi-relative sensitivities of active caspase-3 with respect to k9 show that the active caspase-3 decreases when the rate constant (k9) or the initial concentration of IKK increases. By contrast, increasing the rate constant (k15) or the initial concentration of FADD induces the increase of the active caspase-3. These results are in agreement with observations by mutant studies [23] and show the interplay of the apoptotic and survival pathways in response to the variations of IKK and FADD. As Figure 7 shows, a combination of increasing the initial concentration of FADD and decreasing the initial concentration of IKK gets more effects on the DNA fragmentation.

thumbnailFigure 7. The symmetry of semi-relative sensitivities with respect to k9 and k15. The solid lines are the semi-relative sensitivities with respect to the rate constant k9 of the formation of survival complex and the short dash lines are the semi-relative sensitivities with respect to the rate constant k15 of the formation of death complex. The semi-relative sensitivities of NF-κB (x16) are shown in red, IκB (x31) in green, activated caspase-3 (x25) in blue, and fragmented DNA (x26) in pink.

The fragmented DNA (x26) is the output of the signal transduction model in response to TNF-α stimulus. Our goal is to identify the essential reactions that have significant effect on the system output. The semi-relative sensitivities of fragmented DNA are used to achieve this goal. Since dynamic sensitivities vary with time, it is hard to determine the most important reaction that has the largest effect on the system output. We consider the usual used significance measure, time-averaged semi-relative sensitivity defined similar to equation (12), as the strength of effects on fragmented DNA for individual reactions. Table 3 shows the ranking of dynamic sensitivities of fragmented DNA with respect to all parameters based on time-averaged semi-relative sensitivities. The values of time-averaged semi-relative sensitivities and the corresponding ratios to the total semi-relative sensitivity are shown if the ratios are greater than 2%. Six out of the top seven key parameters belong to the apoptosis pathway. This means that the apoptosis pathway dominates the cell fate in response to TNF-α. The rate constants for the formation of phosphorylated caspase-8 (k20) and phosphorylated caspase-3 (k23) and the regulation of activating caspase-3 by phosphorylated caspase-8 (k21 ) are identified as important parameters. This result indicates that the activation cascades of caspase-8 and caspase-3 are important reactions in the apoptosis signal transduction pathway. Figure 8 shows the top seven semi-relative sensitivities of fragmented DNA. All key parameters positively regulate the activity of DNA fragmentation except parameter k9.

Table 3. The ranking of semi-relative sensitivities for fragmented DNA

thumbnailFigure 8. The semi-relative sensitivities of fragmented DNA. The semi-relative sensitivities of fragmented DNA (x26) with respect to the rate constants of the formation of survival complex (k9), the formation of death complex (k15), the formation of DISC without TNFR1 (k17), the caspase-8 activation (k20), the cleavage of procaspase-3 (k21), the caspase-3 activation (k23), and the fragmentation of DNA (k26).

The transcription factor NF-κB (x16) plays an important role in the survival pathway. To further understand the control of NF-κB, we analyze the semi-relative sensitivities of NF-κB with respect to all parameters. The ranking of parameter sensitivities of NF-κB based on the time-averaged semi-relative sensitivities is shown in Table 4. A parameter is referred to as an important parameter if the ratio of its time-averaged semi-relative parameter sensitivity to the total parameter sensitivity of NF-κB is greater than 2%. Ignoring the rate constants for the reverse reactions, all parameters -k9, k11, k12, k14, k28, and k29 - in the survival pathway are identified as important parameters. The parameter k15 is also identified as an important parameter, although NF-κB is negatively sensitive to it. The top seven sensitivities of NF-κB are shown in Figure 9. The dynamic sensitivity of NF-κB with respect to the IKK activation (k11) is similar to that with respect to the NF-κB activation (k14) as Figure 9 shows. This means that the inhibitor of IKK acts the same function as the inhibitor of NF-κB to inhibit the activity of NF-κB. The activation of NF-κB in most types of cells leads to the inhibition of apoptosis, accelerates cell proliferation, and promotes tumorigenesis. To inhibit cell growth, some small-molecule inhibitors targeting IKK have already been developed to treat certain type tumors [42]. Several synthetic drugs that are be able to inhibit the activities of IKK and NF-κB have been shown to have the same effects on tumor development [43,44]. Larger the rate constant for the transcription of IκB (k29), more IκB is generated to deactivate NF-κB and causes lower active NF-κB. A large rate constant for the formation of death complex (k15) leads to less signals to activate IKK that then phosphorylates IκB and triggers the activation of NF-κB. Both rate constants regulate negatively the activation of NF-κB as shown in Figure 9. The sensitivities of NF-κB with respect to the rates of IKK activation (k11), the formation of NF-κB/IκB/IKK* (k12), the NF-κB activation (k14), and the deactivation of NF-κB (k28) alternate between positive and negative values. These oscillations of sensitivities may be caused by the negative feedback regulation of IκB.

Table 4. The ranking of semi-relative sensitivities of NF-κB

thumbnailFigure 9. The semi-relative sensitivities of NF-κB. The semi-relative sensitivities of NF-κB (x16) with respect to the rate constants of the formation of survival complex (k9), the IKK activation (k11), the formation of NF-κB/IκB/IKK* (k12), the NF-κB activation (k14), the formation of death complex (k15), the deactivation of NF-κB (k28), and the transcription of cIAP and IκB (k29).

Following similar procedures mentioned above, we investigate the regulation of the apoptosis pathway. The DISC complex is essential for TNF-induced apoptosis and it is required for casepase-8 activation. To investigate the regulation of apoptosis, we identify the important reactions that regulate the formation of DISC by sensitivity analysis. The ranking of dynamic sensitivities of DISC based on the time-averaged semi-relative sensitivities is shown in Table 5. The key parameters -k9, k15, k17, and k20 - are identified and the dynamic sensitivities of DISC with respect to these four parameters are shown in Figure 10. The reaction of dissociation of DISC from the death receptor TNFR1 is essential for the following caspase-8 activation and its corresponding rate constant k17 is identify as an important parameter. This result is in agreement with the observation in an in vitro binding assay by Harper et al. [40].

thumbnailFigure 10. The semi-relative sensitivities of DISC. The semi-relative sensitivities of DISC (x21) with respect to the rate constants of the formation of survival complex (k9), the formation of death complex (k15), the formation of DISC without TNFR1 (k17 ), the caspase-8 activation (k20).

Table 5. The ranking of semi-relative sensitivities of DISC

Efficiency and accuracy

To verify the result obtained by the EAMCM algorithm, it is compared with the finite difference method using the dde23 as the DDE solver. The dde23 solver is available in MATLAB 6.5 and later. Forward difference is considered in the finite difference method. The dynamic sensitivities of these two systems mentioned above are solved by the finite difference method with spacing ratio 0.1 and 0.01, respectively. The relative parameter sensitivities of heart rate and blood pressure with respect to β obtained by the finite difference method with spacing ratio 0.1 and 0.01, respectively, and the EAMCM method are shown in Figure 11 as an illustration (another data is similar and not shown here). According to the definition of relative sensitivity, the theoretical value of relative sensitivity is obtained when the spacing ratio is approaching to zero. From Figure 11, the relative sensitivities obtained by the EAMCM are close to the theoretical values. We analyze the performance of the finite difference method and the EAMCM method for computing the dynamic sensitivities by measuring the number of evaluations of model equations. The results are shown in Table 6. The CPU time in second running by a 1.86 GHz Intel Xeon CPU with 4 GMB RAM is shown in the parenthesis for reference. Based on the comparison, the EAMCM program surely outperforms the finite difference method using dde23 solver. The EAMCM program can be accessed from http://www.che.ccu.edu.tw/~bioproc/index_english.files/page00064.htm webcite and a brief manual can be found in the Additional file 2.

thumbnailFigure 11. The relative sensitivities obtained by the finite difference method and the EAMCM method. a) The relative sensitivities of heart rate with respect to the uncontrolled average arterial blood pressure (β); b) The relative sensitivities of blood pressure with respect to β. The green and red lines are obtained by the finite difference method with spacing ratio 0.1 and 0.01, respectively. The blue line is obtained by the EAMCM method. The time is in dimensionless scale.

Table 6. The number of evaluations of model equations

Additional file 2. Manual of EAMCM program. This document provides a brief introduction to the extended adaptive modified collocation method (EAMCM) that we implemented to compute the solution and dynamic sensitivities of ordinal differential equation (ODE) and delay differential equation (DDE) systems.

Format: PDF Size: 226KB Download file

This file can be viewed with: Adobe Acrobat ReaderOpen Data

Conclusions

We extend the applicability of the adaptive direct-decoupled algorithm for ODE models to DDE models and include the implementation of automatic differentiation technique among it. The most attractive feature of the EAMCM program is minimal user intervention that can reduce the human effort required for solving the dynamic sensitivities of complex biological systems and reduce the number of human errors introduced. EAMCM requires the user to supply only the model equations at run-time to compute dynamic sensitivities of DDE models. The evaluation of sensitivity equations is done automatically by automatic differentiation technique along with the inevitable evaluation of model equations. The computations of partial derivatives and values of model equations simultaneously induce less overhead cost of computer time. The exact accuracy of the computed derivatives is achieved by the property of automatic differentiation. By compared with direct-coupled methods in theory, the adaptive direct-decoupled EAMCM algorithm is efficient, accurate, and easy to use for end users without programming background to do dynamic sensitivity analysis on complex biological systems with time-delays.

We illustrate the use of the EAMCM program in the sensitivity analysis of two DDE models: the cardiovascular control system and the TNF-α signal transduction network. The parameters for sympathetical and vagal control of heart rates are identified as key parameters in the cardiovascular control system. From the symmetry of dynamic effects of sympathetical and vagal control on heart rate obtained by sensitivity analysis, it reflects the sympathovagal balance in physiology. The TNF-α signal transduction network is a more complicated system than the first model and symbolic differentiation is unaffordable in this case to obtain the sensitivity equations. By using the EAMCM program, users can provide the model equations only for solving the dynamic sensitivities of the model. The formation of survival and death complexes are identified as the key reactions for the fragmentation of DNA via sensitivity analysis. This result reveals that the interplay between the components of the survival and apoptosis pathways plays an important role in the TNF-α signal transduction network.

Methods

Delay differential equations (DDEs) arise in either natural or technological control problems for a large and important class of dynamical systems. This type of dynamical system now occupies a place of central importance in all areas of science and particularly in the biological sciences [45]. There are various kinds of delay differential equations. Here, we focus on equations with fixed, discrete delays, namely those of the form

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M4">View MathML</a>

(3)

where x(t) ∈ ℝn is a vector of state variables, θ ∈ ℝp is a vector of parameters, τi are positive time-delays, and r is the number of multiple delays. DDE models are similar to ODE models, but their evolution involves past values of the state variables. When giving initial conditions for ODE systems, we only need to specify the initial values of the state and input variables. In order to solve DDE systems, we have to look back to earlier values of x at every time step. Therefore, it is necessary to provide an initial function to specify the value of the solution before time t = 0. This function has to cover a period at least as long as the longest delay since we look back in time that far.

DDE solver

Most DDEs do not have analytic solutions, so it is generally necessary to resort to numerical methods. We have presented an adaptive modified collocation method (AMCM) for computing the solution of autonomous ODE systems [21]. This method is easy to extend to compute the solution of simple scale DDE systems. To simplify discussion, we assume the DDE model consists of a set of DDEs with a single delay and expressed as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M5">View MathML</a>

(4)

where the delay τ is a positive constant. The initial function x(t) defined on the interval [-τ, 0] is set to x(0). Equation (4) can be reduced to ODEs by introducing a new input variable y for each delayed variablen x(t -τ) as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M6">View MathML</a>

and can be solved by our ODE solver in the AMCM. The AMCM algorithm with piecewise linear polynomials as the shape functions solves the ODEs transformed from equation (4) for each subinterval [tj-1, tj], ti-1 tj-1 <tj ti by

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M7">View MathML</a>

(5)

where ηj is the step size in time tj . The adaptive ODE solver in the AMCM controls the step size automatically. The earlier values of x defined or computed on the interval [t - τ, t] have to be stored for y at every computation time t. Due to the automatic step size control, the solution at time t - τ may be not computed by the AMCM. To solve this problem, interpolation is used to generate the solution for the input variable y as follows:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M8">View MathML</a>

(6)

where ta is the time point that x(ta) has been computed, ta+1 is the next time point, and the step size ta+1 - ta is determined by the AMCM automatically.

Dynamic sensitivity analysis

For a model described by equation (4), the absolute parameter sensitivity s(xi , θj ) of dependent variable xi x with respect to a change in parameter θj θ is defined as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M9">View MathML</a>

(7)

where xi(t; θj + Δθj) is the ith component of the solution of equation (4) with a change Δθj on the jth parameter and the others fixed. The absolute parameter sensitivity s(xi , θj ) is also defined as the first-order local sensitivity of xi with respect to θj [46]. The term local refers to the fact that the value of s(xi , θj ) depends on the given set of values for the parameters θ. Under the assumption that the system responds linearly for small perturbations, s(xi , θj ) measures the ratio between the effect on xi and the variation of θj . It is useful to consider the ratio between the relative effect on the output and the relative variation of a parameter when comparing different parameter sensitivities with respect to different parameters. The relative parameter sensitivity S(xi , θj ) of xi with respect to θj , a dimensionless quantity, is defined as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M10">View MathML</a>

(8)

Sometimes the use of relative parameter sensitivities has the difficulty of numerical instability caused by the division by zero when xi approaching zero. To address this problem, the semi-relative parameter sensitivity <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M11">View MathML</a> of xi with respect to θj is used and defined as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M12">View MathML</a>

(9)

Once the local sensitivity is known, the calculation of the relative sensitivity is straightforward. For brevity's sake, we limit our explanation on the absolute sensitivity only below. For an autonomous system describing by equation (4), the sensitivity equations are given as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M13">View MathML</a>

(10)

where fi is the ith element of f . The AMCM algorithm is extended to compute the solution of equations (10), i.e., dynamic sensitivities of DDE systems. When the solution of a DDE system is obtained, the absolute dynamic sensitivity of xi with respect to θj is computed by

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M14">View MathML</a>

(11)

where ηl is the step size in time tl . Same as the AMCM algorithm, the extended algorithm (EAMCM) determines the step size automatically when computing the absolute dynamic sensitivity of xi .

The dynamic sensitivity in equation (8) reflects a relative relationship between the magnitudes of a parameter and a state variable at any time. We can define the integral value for the relative dynamic sensitivity over the whole time as a measure for ranking all sensitivities to identify the bottleneck reactions of the system. The time-averaged relative sensitivity is therefore defined as:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M15">View MathML</a>

(12)

where tf is the final time.

Automatic differentiation

The main challenge in computing the solution of equation (10) is the evaluation of partial derivatives of f with respect to all state variables x, delay variables x(t - τ), and system parameters θ. One way to obtain the partial derivatives is to use symbolic differentiation tools, such as Maple and Mathematica, to perform the algebra of the differentiation. The explicit expression of the partial derivatives of f can be generated automatically. This is very useful because it saves the human effort and avoids the human errors in the analytical differentiation process. In principle, this approach gives exact values of the partial derivatives of f at the expense of high computation cost. Symbolic differentiation tools always generate lengthy formulas containing many common subexpressions that require considerable computation to evaluate. If only the values of the partial derivatives rather than the explicit expressions of the partial derivatives of f are needed, the simplest and common used approach is the numerical differentiation by finite difference approximation

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M16">View MathML</a>

(13)

The main drawback of this approach is that the accuracy is hard to analyze. Another approach which can be used to evaluate partial derivatives is automatic differentiation. Automatic differentiation is a numerical computation of exact values of the partial derivatives without generating a formula for the partial derivatives and is much more effcient than symbolic differentiation. This approach is based on the fact that every function, no matter how complicated, can be represented by a well-formed expression that is a finite combination of elementary arithmetic operators, such as addition (+), subtraction (-), multiplication (*), division (/), and power (^), primitives, such as a constant or a state variable, and intrinsic functions, such as sin, cos, etc. Each elementary arithmetic operation involves at most two operands which either have been computed in a previous step or are primitives. The chain rule can be applied to each of elementary arithmetic operators as follows:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M17">View MathML</a>

(14)

where h indicates the function for each elementary arithmetic operator and oi is the ith operand. The derivative of each primitive and each intrinsic function would have to know. By applying the chain rule recursively to each elementary arithmetic operator and each intrinsic function, the partial derivatives of a function can be computed exactly and in a completely mechanical fashion.

EAMCM algorithm

The proposed algorithm EAMCM is shown as follows:

Algorithm EAMCM

Input:

1. A set of n delay differential equations <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M18">View MathML</a> = f (x, y) with n dependent variables xi , i = 1, ..., n and m time-delay variables yi xj (t - τ), i = 1, ..., m, j ∈ [1, n].

2. Two order sets x0 = {xi (0)|i = 1, ..., n} and Φ0 = {ϕij (0)|i = 1, ..., n, j = 1, ..., m}.

3. An order set T = {t1, ..., tk} of sampling points, ti, 1 ≤ i k is the sampling time of the solution of each DDE, k is the number of sampling points.

4. A tolerance ε.

Output: The dynamic sensitivities of dependent variables at each sampling time.

• For each sampling time ti in T.

1. ηj ti - ti-1, dt ← 0, xc x(ti-1 ), Φc Φ(ti-1).

2. Repeat the following steps until dt = ti - ti-1.

(a) xp xc , yp x(ti-1 + di - τ), Φp Φc .

(b) If ti-1 >τ , then Φd Φ(ti-1+ di - τ ); otherwise, Φd Φ0.

(c) Evaluate the Jacobin matrix by automatic differentiation <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M19">View MathML</a>.

(d) Compute the upper bound μ of the value of ||A||2 by <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M20">View MathML</a>.

(e) If μ * ε ≥ 1, it means the DDEs are stiff, then exit this algorithm.

(f) If μ * ηj > 1, then ηj ← 0.9/μ.

(g) Call the Iteration algorithm to compute the value of xc stepped forward ηj from xp .

(h) Call the IterationOfSensitivity algorithm to compute the value of the sensitivity matrix Φc stepped forward ηj from Φp.

(i) If the Iteration and IterationOfSensitivity algorithms succeed in computing xc and Φc respectively, then dt dt + ηj and ηj ti - ti-1-dt ; otherwise exist this algorithm.

3. Save the value of the sensitivity matrix Φc as Φ(ti).

• return Φ(ti), i = 1, ..., k.

End of Algorithm EAMCM

Algorithm Iteration

Input:

1. A set of n delay differential equations <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M18">View MathML</a> = f (x, y), y(t) ≡ x(t - τ ).

2. x(t), y(t), ηj , and the iteration limitation.

Output: x(t + ηj ).

1. Evaluate the value of f (x(t), y(t)).

2. x(t + ηj) ← x(t) + f (x(t), y(t)) * ηj.

3. y(t + ηj) ← x(t + ηj - τ).

4. Repeat the following steps until the iteration limitation is reached or the value of x(t+ηt) is convergent.

(a) Evaluate the value of f (x(t + ηj), y(t + ηj )).

(b) x(t + ηj) ← x(t) + 0.5 * ηj * (f (x(t), y(t)) + f (x(t + ηj), y(t + ηj))).

5. If the iteration limitation is reached, then exit this algorithm; otherwise, return x(t + ηj).

End of Algorithm Iteration

Algorithm IterationOfSensitivity

Input:

1. A set of n delay differential equations <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M18">View MathML</a> = f (x, y), y(t) ≡ x(t - τ).

2. ηj and the iteration limitation.

3. Two vectors of dependent variables x(t) at time t and x(t + ηj) at time t + ηj .

4. Two sensitivity matrices Φ(x(t)) and Φ(y(t)).

Output: The sensitivity matrix Φ(x(t + ηj )).

1. Evaluate <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M21">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M22">View MathML</a>, and <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M23">View MathML</a> by automatic differentiation.

2. Evaluate the derivative <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M24">View MathML</a>(x(t)) = (x(t)) + (y(t)) + C.

3. Φ(x(t + ηj)) = Φ(x(t)) + <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M24">View MathML</a>(x(t)) * ηj .

4. Evaluate <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M25">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M26">View MathML</a>, and <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M27">View MathML</a> by automatic differentiation.

5. Repeat the following steps until the iteration limitation is reached or the value of Φ(x(t + ηj)) is convergent.

(a) Evaluate the derivative <a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M24">View MathML</a>(x(t + ηj)) = (x(t + ηj)) + (y(t + ηj)) + C.

(b) Evaluate the new value of Φ(x(t + ηj)) by

<a onClick="popup('http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/11/S7/S12/mathml/M28">View MathML</a>

6. If the iteration limitation is reached, then exit this algorithm; otherwise, return Φ(x(t + ηj)).

End of Algorithm IterationOfSensitivity

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

WHW developed and implemented the algorithm and drafted the manuscript. FSW conceived of the study, participated in its design and coordination, and helped to draft the manuscript. MSC assisted in developing the algorithm and finalizing the manuscript. All authors read and approved the final manuscript.

Acknowledgements

The financial support from the National Science Council, Taiwan, ROC (Grant NSC97-2221-E-194-010-MY3 and NSC99-2627-B-194-001), is highly appreciated.

This article has been published as part of BMC Bioinformatics Volume 11 Supplement 7, 2010: Ninth International Conference on Bioinformatics (InCoB2010): Bioinformatics. The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2105/11?issue=S7.

References

  1. Fowler AC, Mackey MC: Relaxation oscillations in a class of delay differential equations.

    SIAM Journal on Applied Mathematics 2002, 63:299-323. Publisher Full Text OpenURL

  2. Monk NAM: Oscillatory expression of Hes1, p53, and NF-κB driven by transcriptional time delays.

    Current Biology 2003, 13(16):1409-1413. PubMed Abstract | Publisher Full Text OpenURL

  3. Nelson PW, Perelson AS: Mathematical analysis of delay differential equation models of HIV-1 infection.

    Mathematical Biosciences 2002, 179:73-94. PubMed Abstract | Publisher Full Text OpenURL

  4. Bocharov GA, Rihan FA: Numerical modelling in biosciences using delay differential equations.

    Journal of Computational and Applied Mathematics 2000, 125(1-2):183-199. Publisher Full Text OpenURL

  5. Enright WH, Hayashi H: A delay differential equation solver based on a continuous Runge-Kutta method with defect control.

    Numerical Algorithms 1997, 16(3-4):349-364. Publisher Full Text OpenURL

  6. Engelborghs K, Luzyanina T, Hout KJ, Roose D: Collocation methods for the computation of periodic solutions of delay differential equations.

    SIAM Journal on Scientific Computing 2001, 22(5):1593-1609. Publisher Full Text OpenURL

  7. Paul CAH: Developing a delay differential equation solver.

    Applied Numerical Mathematics 1992, 9(3-5):403-414. Publisher Full Text OpenURL

  8. Shampine L, Thompson S: Solving DDEs in Matlab.

    Appl Numer Math 2001, 37:441-458. Publisher Full Text OpenURL

  9. Shampine LF: Solving ODEs and DDEs with residual control.

    Applied Numerical Mathematics 2005, 52:113-127. Publisher Full Text OpenURL

  10. Thompson S, Shampine LF: A friendly Fortran DDE solver.

    Applied Numerical Mathematics 2006, 56(3-4):503-516. Publisher Full Text OpenURL

  11. Willé DR, Baker CT: DELSOL - a numerical code for the solution of systems of delay-differential equations.

    Applied Numerical Mathematics 1992, 9(3-5):223-234. Publisher Full Text OpenURL

  12. Zivaripiran H: DDVERK90: A user-friendly implementation of an effective DDE solver.

    Master's thesis, University of Toronto, Canada 2005. OpenURL

  13. Daescu DN, Sandu A, Carmichael GR: Direct and adjoint sensitivity analysis of chemical kinetic systems with KPP: II-numerical validation and applications.

    Atmospheric Environment 2003, 37(36):5097-5114. Publisher Full Text OpenURL

  14. Dougherty EP, Hwang JT, Rabitz H: Further developments and applications of the Green's function method of sensitivity analysis in chemical kinetics.

    J Chem Phys 1979, 71:1794-1808. Publisher Full Text OpenURL

  15. Dunker AM: The decoupled direct method for calculating sensitivities coefficients in chemical kinetics.

    J Chem Phys 1984, 81:2385-2393. Publisher Full Text OpenURL

  16. Rihan FA: Sensitivity analysis for dynamic systems with time-lags.

    Journal of Computational and Applied Mathematics 2003, 151(2):445-462. Publisher Full Text OpenURL

  17. Sandu A, Daescu DN, Carmichael GR: Direct and adjoint sensitivity analysis of chemical kinetic systems with KPP: Part I-theory and software tools.

    Atmospheric Environment 2003, 37(36):5083-5096. Publisher Full Text OpenURL

  18. Bischof C, Khademi P, Mauer A, Carle A: Adifor 2.0: Automatic differentiation of Fortran 77 programs.

    Computing in Science and Engineering 1996, 3(3):18-32. Publisher Full Text OpenURL

  19. Hwang D, Byun DW, Odman MT: An automatic differentiation technique for sensitivity analysis of numerical advection schemes in air quality models.

    Atmospheric Environment 1997, 31(6):879-888. Publisher Full Text OpenURL

  20. Griewank A, Juedes D, Mitev H, Utke J, Vogel O, Walther A: ADOL-C: A package for the automatic differentiation of algorithms written in C/C++.

    1999.

  21. Wu WH, Wang FS, Chang MS: Dynamic sensitivity analysis of biological systems.

    BMC Bioinformatics 2008, 9:S17. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  22. Fowler AC, McGuinness MJ: A delay recruitment model of the cardiovascular control system.

    Journal of Mathematical Biology 2005, 51(5):508-526. PubMed Abstract | Publisher Full Text OpenURL

  23. Rangamani P, Sirovich L: Survival and apoptotic pathways initiated by TNF-α: modeling and predictions.

    Biotechnol Bioeng 2007, 97(5):1216-1229. PubMed Abstract | Publisher Full Text OpenURL

  24. Ashe PC, Berry MD: Apoptotic signaling cascades.

    Prog Neuropsychopharmacol Biol Psychiatry 2003, 27:199-214. PubMed Abstract | Publisher Full Text OpenURL

  25. Beg AA, Sha WC, Bronson RT, Ghosh S, Baltimore D: Embryonic lethality and liver degeneration in mice lacking the RelA component of NF-κB.

    Nature 1995, 376:167-170. PubMed Abstract | Publisher Full Text OpenURL

  26. Beg AA, Baltimore D: An essential role for NF-κB in preventing TNF-α-induced cell death.

    Science 1996, 274:782-784. PubMed Abstract | Publisher Full Text OpenURL

  27. Bentele M, Lavrik I, Ulrich M, Stosser S, Heermann DW, Kalthoff H, Krammer PH, Eils R: Mathematical modeling reveals threshold mechanism in CD95-induced apoptosis.

    J Cell Biol 2004, 166:839-851. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  28. Cho KH, Shin SY, Kolch W, Wolkenhauer O: Experimental design in systems biology, based on parameter sensitivity analysis using a Monte Carlo method: a case study for the TNFα-mediated NF-κB signal transduction pathway.

    Simulation 2003, 79:726-739. Publisher Full Text OpenURL

  29. Fussenegger M, Bailey JE, Varner J: A mathematical model of caspase function in apoptosis.

    Nature Biotechnology 2000, 18:768-774. PubMed Abstract | Publisher Full Text OpenURL

  30. Hoffmann A, Levchenko A, Scott ML, Baltimore D: The IκB-NF-κB signaling module: temporal control and selective gene activation.

    Science 2002, 298:1241-1245. PubMed Abstract | Publisher Full Text OpenURL

  31. Lipniacki T, Paszek P, Brasier AR, Luxon B, Kimmel M: Mathematical model of NF-κB regulatory module.

    Journal of Theoretical Biology 2004, 228(2):195-215. PubMed Abstract | Publisher Full Text OpenURL

  32. Barrett CJ, Guild SJ, Ramchandra R, Malpas SC: Baroreceptor denervation prevents sympathoinhibition during angiotensin II-induced hypertension.

    Hypertension 2005, 46:168-172. PubMed Abstract | Publisher Full Text OpenURL

  33. Ramchandra R, Hood SG, Denton DA, Woods RL, McKinley MJ, McAllen RM, May CN: Basis for the preferential activation of cardiac sympathetic nerve activity in heart failure.

    Proc Natl Acad Sci U S A 2009, 106(3):924-928. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  34. Ottesen JT: Modelling of the baroreflex-feedback mechanism with time-delay.

    Journal of Mathematical Biology 1997, (36):41-63. PubMed Abstract | Publisher Full Text OpenURL

  35. McSharry PE, McGuinness MJ, Fowler AC: Confronting a cardiovascular system model with heart rate and blood pressure data.

    Computers in Cardiology 2005, 32:587-590. Publisher Full Text OpenURL

  36. McSharry PE, Clifford GD, Tarassenko L, Smith LA: A dynamical model for generating synthetic electrocardiogram signals.

    IEEE Transactions on Biomedical Engineering 2003, 50(3):289-294. PubMed Abstract | Publisher Full Text OpenURL

  37. Luo JL, Kamata H, Karin M: IKK/NF-κB signaling: balancing life and death - a new approach to cancer therapy.

    J Clin Invest 2005, 115:2625-2632. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  38. Vuillard L, Nicholson J, Hay RT: A complex containing βTrCP recruits Ccd34 to catalyse ubiquitination of IκBα.

    FEBS Letters 1999, 455:311-314. PubMed Abstract | Publisher Full Text OpenURL

  39. Nelson DE, Ihekwaba AEC, Elliott M, Johnson JR, Gibney CA, Foreman BE, Nelson G, See V, Horton CA, Spiller DG, Edwards SW, McDowell HP, Unitt JF, Sullivan E, Grimley R, Benson N, Broomhead D, Kell DB, White MRH: Oscillations in NF-B signaling control the dynamics of gene expression.

    Science 2004, 306:704-708. PubMed Abstract | Publisher Full Text OpenURL

  40. Harper N, Hughes M, MacFarlane M, Cohen GM: Fas-associated death domain protein and caspase-8 are not recruited to the tumor necrosis factor receptor 1 signaling complex during tumor necrosis factor-induced apoptosis.

    J Biol Chem 2003, 278:25534-25541. PubMed Abstract | Publisher Full Text OpenURL

  41. Sung MH, Simon R: In silico simulation of inhibitor drug effects on nuclear factor-κB pathway dynamics.

    Molecular Pharmacology 2004, 66:70-75. PubMed Abstract | Publisher Full Text OpenURL

  42. Haefner B: arresting a major culprit in cancer.

    Drug Discovery Today 2002, 7(12):653-663. PubMed Abstract | Publisher Full Text OpenURL

  43. Yamamoto Y, Yin MJ, Lin KM, Gaynor RB: Sulindac inhibits activation of the NF-κB pathway.

    Journal of Biological Chemistry 1999, 274(38):27307-27314. PubMed Abstract | Publisher Full Text OpenURL

  44. Yin MJ, Yamamoto Y, Gaynor RB: The anti-inflammatory agents aspirin and salicylate inhibit the activity of IκB kinase-β.

    Nature 1998, 396(6706):77-80. PubMed Abstract | Publisher Full Text OpenURL

  45. Baker CTH, Paul CAH, Willé DR: Issues in the numerical solution of evolutionary delay differential equations.

    Advances in Computational Mathematics 1995, 3:171-196. OpenURL

  46. Varma A, Morbidelli M, Wu H: Parameter sensitivity in chemical systems. Cambridge: Cambridge University Press; 1999. OpenURL