Sensitivity analysis of dynamic biological systems with timedelays
 Wu Hsiung Wu^{1},
 Feng Sheng Wang^{2}Email author and
 Maw Shang Chang^{1}
https://doi.org/10.1186/1471210511S7S12
© Wu et al; licensee BioMed Central Ltd. 2010
Published: 15 October 2010
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 timedelays. 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 directdecoupled algorithm is extended to solve the solution and dynamic sensitivities of timedelay 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 timedelays: 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 directcoupled 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 timedelays.
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 timedelays. 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]. Another approach which can be used to evaluate sensitivity equations is automatic differentiation. Automatic differentiation is a nonapproximate 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 roundoff 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 opensource code, automatic differentiation by overloading in C++ (ADOLC), for the automatic differentiation of C and C++ programs. The implementation of ADOLC utilizes the operator overloading capability of C/C++ compilers that accept userdefined data types, operators and functions. The implementation of either the source code transformation or the operator overloading is a compiletime solution. It allows one to generate derivatives from complicated existing code or userprovided model equations that expressed by userdefined 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 timedelay systems described by DDEs in this paper and named as extended adaptive modified collocation method (EAMCM). The EAMCM is implemented as a userfriendly program that facilitates the dynamic sensitivity analysis of DDE models through the implementation of adaptive directdecoupled method and automatic differentiation. EAMCM requires the user to supply only the model equations at runtime 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 directdecoupled 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 timedelays.
To evaluate the applicability of the extended algorithm, it is applied to two realistic models with timedelays: 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κ Bmediated survival pathway and the caspasemediated 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κ Bmediated 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κ Bmediated survival pathway and the caspasemediated 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 slowacting sympathetic nerves. This negativefeedback 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.
The value and definition of system parameters
Parameter  Definition  Value 

h _{0}  Uncontrolled heart rate  100 bpm 
p _{0}  Mean arterial blood pressure  100 mm Hg 
α  Sympathetic effect on peripheral resistance  15 
β  Sympathetic control of heart rate  10 
ν  Strength of vagal tone  9.63 
δ  Relaxation time  0.8 s^{1} 
γ  Damping effect of vagal activity on the sympathetic tone  0.2 
μ  3/(2 + α)  0.18 
A _{1}  Amplitude of the influence of respiration on blood pressure  0 
A _{2}  Amplitude of the influence of respiration on heart rate  0.003 
f _{r}  Breathing rate  0.17 Hz 
τ  Sympathetic time delay  3 s 
ϕ  Phase lag  3.14 s 
n  Hill exponent  8 
ε _{ h }  Relative coefficient for heart rate  1 
ε _{ p }  Relative coefficient for blood pressure  3 
The ranking of relative sensitivities of heart rate and blood pressure
Rank  1  2  3  4  5 

Heart rate  p _{0}  f _{ r }  τ  β  v 
10.667  4.791  4.791  1.284  1.237  
39.97%  18.61%  18.61%  4.81%  4.63%  
Blood pressure  p _{0}  v  β  f _{ r }  τ 
1.000  0.116  0.109  0.061  0.061  
70.56%  8.18%  7.68%  4.30%  4.30% 
Apoptosis signal 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.
The dynamic sensitivity profiles for all species with respect to k_{9}, the rate constant of the formation of survival complex, are nearly identical to that with respect to x_{10}(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 k_{15} /x_{17}(0), k_{18}/x_{20}(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.
The ranking of semirelative sensitivities for fragmented DNA
Rank  Parameter  Timeaveraged semirelative sensitivity  Percentage (%) 

1  k _{26}  51.349  25.49 
2  k _{21}  40.336  20.02 
3  k _{20}  33.558  16.66 
4  k _{15}  24.010  11.92 
5  k _{9}  23.755  11.79 
6  k _{17}  17.756  8.81 
7  k _{23}  7.492  3.72 
The ranking of semirelative sensitivities of NFκ B
Rank  Parameter  Timeaveraged semirelative sensitivity  Percentage (%) 

1  k _{29}  1.300E02  23.77 
2  k _{14}  1.082E02  19.79 
3  k _{11}  1.081E02  19.76 
4  k _{9}  6.709E03  12.27 
5  k _{15}  6.672E03  12.20 
6  k _{12}  2.255E03  4.12 
7  k _{28}  2.123E03  3.88 
The ranking of semirelative sensitivities of DISC
Rank  Parameter  Timeaveraged semirelative sensitivity  Percentage (%) 

1  k _{20}  0.1270  33.43 
2  k _{15}  0.0964  25.38 
3  k _{9}  0.0963  25.35 
4  k _{17}  0.0530  13.96 
Efficiency and accuracy
The number of evaluations of model equations
System  Finite difference method with dde23 solver  EAMCM method  

Cardiovascular  441,112(160)  60,000(8)  
Apoptosis  311,464(246)  177,329(175) 
Conclusions
We extend the applicability of the adaptive directdecoupled 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 runtime 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 directcoupled methods in theory, the adaptive directdecoupled 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 timedelays.
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
where x(t) ∈ ℝ^{ n }is a vector of state variables, θ ∈ ℝ^{ p }is a vector of parameters, τ_{ i }are positive timedelays, 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
where t_{ a }is the time point that x(t_{ a }) has been computed, t_{a+1}is the next time point, and the step size t_{a+1} t_{ a }is determined by the AMCM automatically.
Dynamic sensitivity analysis
where η_{ l }is the step size in time t_{ l }. Same as the AMCM algorithm, the extended algorithm (EAMCM) determines the step size automatically when computing the absolute dynamic sensitivity of x_{ i }.
where t_{ f }is the final time.
Automatic differentiation
where h indicates the function for each elementary arithmetic operator and o_{ i }is the i^{ th }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
 1.
A set of n delay differential equations $\dot{x}$ = f (x, y) with n dependent variables x_{ i }, i = 1, ..., n and m timedelay variables y_{ i }≡ x_{ j }(t  τ), i = 1, ..., m, j ∈ [1, n].
 2.
Two order sets x_{0} = {x_{ i }(0)i = 1, ..., n} and Φ_{0} = {ϕ_{ ij }(0)i = 1, ..., n, j = 1, ..., m}.
 3.
An order set T = {t_{1}, ..., t_{ k }} of sampling points, t_{ i }, 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 t_{ i }in T.
 1.
η_{ j }← t_{ i } t_{i1}, d_{ t }← 0, x^{ c }← x(t_{i1}), Φ^{ c }← Φ(t_{i1}).
 2.
Repeat the following steps until d_{ t }= t_{ i } t_{i1}.
 (a)
x^{ p }← x^{c} , y^{p} ← x(t_{i1}+ d_{ i } τ), Φ^{ p }← Φ^{ c }.
 (b)
If t_{i1}>τ , then Φ^{ d }← Φ(t_{i1}+ d_{ i } τ ); otherwise, Φ^{ d }← Φ_{0}.
 (c)
Evaluate the Jacobin matrix by automatic differentiation $A\leftarrow \left[\frac{\partial f\text{(}{x}^{p}\text{,}{y}^{p}\text{)}}{\partial x},\frac{\partial f\text{(}{x}^{p}\text{,}{y}^{p}\text{)}}{\partial y}\right]$.
 (d)
Compute the upper bound μ of the value of A_{2} by $\sqrt{n(m+n)}\parallel A{\parallel}_{\Delta},\parallel A{\parallel}_{\Delta}\equiv \underset{i,j}{\mathrm{max}}\left{a}_{ij}\right$.
 (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 x^{ c }stepped forward η_{ j }from x^{ p }.
 (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 x^{ c }and Φ^{ c }respectively, then d_{ t }← d_{ t }+ η_{ j }and η_{ j }← t_{ i } t_{i1}d_{ t }; otherwise exist this algorithm.
 3.
Save the value of the sensitivity matrix Φ^{ c }as Φ(t_{ i }).

return Φ(t_{ i }), i = 1, ..., k.
End of Algorithm EAMCM
Algorithm Iteration
 1.
A set of n delay differential equations $\dot{x}$ = f (x, y), y(t) ≡ x(t  τ ).
 2.
x(t), y(t), η_{ j }, and the iteration limitation.
 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
 1.
A set of n delay differential equations $\dot{x}$ = 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)).
 1.
Evaluate $A=\frac{\partial f(x(t),y(t))}{\partial x}$, $B=\frac{\partial f(x(t),y(t))}{\partial y}$, and $C=\frac{\partial f(x(t),y(t))}{\partial \theta}$ by automatic differentiation.
 2.
 3.
 4.
Evaluate $A=\frac{\partial f(x(t+\eta j),y(t+\eta j))}{\partial x}$, $B=\frac{\partial f(x(t+\eta j),y(t+\eta j))}{\partial y}$, and $C=\frac{\partial f(x(t+\eta j),y(t+\eta j))}{\partial \theta}$ 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 $\Phi $(x(t + η_{ j })) = AΦ(x(t + η_{ j })) + BΦ(y(t + η_{ j })) + C.
 (b)Evaluate the new value of Φ(x(t + η_{ j })) by$\Phi (x(t+{\eta}_{j}))=\Phi (x(t))+0.\text{5}{\eta}_{j}\ast (\dot{\Phi}(x(t+{\eta}_{j}))+\dot{\Phi}(x(t))).$
 6.
If the iteration limitation is reached, then exit this algorithm; otherwise, return Φ(x(t + η_{ j })).
End of Algorithm IterationOfSensitivity
Declarations
Acknowledgements
The financial support from the National Science Council, Taiwan, ROC (Grant NSC972221E194010MY3 and NSC992627B194001), 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/14712105/11?issue=S7.
Authors’ Affiliations
References
 Fowler AC, Mackey MC: Relaxation oscillations in a class of delay differential equations. SIAM Journal on Applied Mathematics 2002, 63: 299–323. 10.1137/S0036139901393512View ArticleGoogle Scholar
 Monk NAM: Oscillatory expression of Hes1, p53, and NF κ B driven by transcriptional time delays. Current Biology 2003, 13(16):1409–1413. 10.1016/S09609822(03)004949View ArticlePubMedGoogle Scholar
 Nelson PW, Perelson AS: Mathematical analysis of delay differential equation models of HIV1 infection. Mathematical Biosciences 2002, 179: 73–94. 10.1016/S00255564(02)000998View ArticlePubMedGoogle Scholar
 Bocharov GA, Rihan FA: Numerical modelling in biosciences using delay differential equations. Journal of Computational and Applied Mathematics 2000, 125(1–2):183–199. 10.1016/S03770427(00)004684View ArticleGoogle Scholar
 Enright WH, Hayashi H: A delay differential equation solver based on a continuous RungeKutta method with defect control. Numerical Algorithms 1997, 16(3–4):349–364. 10.1023/A:1019107718128View ArticleGoogle Scholar
 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. 10.1137/S1064827599363381View ArticleGoogle Scholar
 Paul CAH: Developing a delay differential equation solver. Applied Numerical Mathematics 1992, 9(3–5):403–414. 10.1016/01689274(92)90030HView ArticleGoogle Scholar
 Shampine L, Thompson S: Solving DDEs in Matlab. Appl Numer Math 2001, 37: 441–458. 10.1016/S01689274(00)000556View ArticleGoogle Scholar
 Shampine LF: Solving ODEs and DDEs with residual control. Applied Numerical Mathematics 2005, 52: 113–127. 10.1016/j.apnum.2004.07.003View ArticleGoogle Scholar
 Thompson S, Shampine LF: A friendly Fortran DDE solver. Applied Numerical Mathematics 2006, 56(3–4):503–516. 10.1016/j.apnum.2005.04.027View ArticleGoogle Scholar
 Willé DR, Baker CT: DELSOL  a numerical code for the solution of systems of delaydifferential equations. Applied Numerical Mathematics 1992, 9(3–5):223–234. 10.1016/01689274(92)900178View ArticleGoogle Scholar
 Zivaripiran H: DDVERK90: A userfriendly implementation of an effective DDE solver. Master's thesis, University of Toronto, Canada 2005.Google Scholar
 Daescu DN, Sandu A, Carmichael GR: Direct and adjoint sensitivity analysis of chemical kinetic systems with KPP: IInumerical validation and applications. Atmospheric Environment 2003, 37(36):5097–5114. 10.1016/j.atmosenv.2003.08.020View ArticleGoogle Scholar
 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. 10.1063/1.438530View ArticleGoogle Scholar
 Dunker AM: The decoupled direct method for calculating sensitivities coefficients in chemical kinetics. J Chem Phys 1984, 81: 2385–2393. 10.1063/1.447938View ArticleGoogle Scholar
 Rihan FA: Sensitivity analysis for dynamic systems with timelags. Journal of Computational and Applied Mathematics 2003, 151(2):445–462. 10.1016/S03770427(02)006593View ArticleGoogle Scholar
 Sandu A, Daescu DN, Carmichael GR: Direct and adjoint sensitivity analysis of chemical kinetic systems with KPP: Part Itheory and software tools. Atmospheric Environment 2003, 37(36):5083–5096. 10.1016/j.atmosenv.2003.08.019View ArticleGoogle Scholar
 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. 10.1109/99.537089Google Scholar
 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. 10.1016/S13522310(96)002403View ArticleGoogle Scholar
 Griewank A, Juedes D, Mitev H, Utke J, Vogel O, Walther A: ADOLC: A package for the automatic differentiation of algorithms written in C/C++. 1999.Google Scholar
 Wu WH, Wang FS, Chang MS: Dynamic sensitivity analysis of biological systems. BMC Bioinformatics 2008, 9: S17. 10.1186/147121059S12S17PubMed CentralView ArticlePubMedGoogle Scholar
 Fowler AC, McGuinness MJ: A delay recruitment model of the cardiovascular control system. Journal of Mathematical Biology 2005, 51(5):508–526. 10.1007/s0028500503391View ArticlePubMedGoogle Scholar
 Rangamani P, Sirovich L: Survival and apoptotic pathways initiated by TNF α : modeling and predictions. Biotechnol Bioeng 2007, 97(5):1216–1229. 10.1002/bit.21307View ArticlePubMedGoogle Scholar
 Ashe PC, Berry MD: Apoptotic signaling cascades. Prog Neuropsychopharmacol Biol Psychiatry 2003, 27: 199–214. 10.1016/S02785846(03)000162View ArticlePubMedGoogle Scholar
 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. 10.1038/376167a0View ArticlePubMedGoogle Scholar
 Beg AA, Baltimore D: An essential role for NF κ B in preventing TNF α induced cell death. Science 1996, 274: 782–784. 10.1126/science.274.5288.782View ArticlePubMedGoogle Scholar
 Bentele M, Lavrik I, Ulrich M, Stosser S, Heermann DW, Kalthoff H, Krammer PH, Eils R: Mathematical modeling reveals threshold mechanism in CD95induced apoptosis. J Cell Biol 2004, 166: 839–851. 10.1083/jcb.200404158PubMed CentralView ArticlePubMedGoogle Scholar
 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. 10.1177/0037549703040943View ArticleGoogle Scholar
 Fussenegger M, Bailey JE, Varner J: A mathematical model of caspase function in apoptosis. Nature Biotechnology 2000, 18: 768–774. 10.1038/81208View ArticlePubMedGoogle Scholar
 Hoffmann A, Levchenko A, Scott ML, Baltimore D: The I κ BNF κ B signaling module: temporal control and selective gene activation. Science 2002, 298: 1241–1245. 10.1126/science.1071914View ArticlePubMedGoogle Scholar
 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. 10.1016/j.jtbi.2004.01.001View ArticlePubMedGoogle Scholar
 Barrett CJ, Guild SJ, Ramchandra R, Malpas SC: Baroreceptor denervation prevents sympathoinhibition during angiotensin IIinduced hypertension. Hypertension 2005, 46: 168–172. 10.1161/01.HYP.0000168047.09637.d4View ArticlePubMedGoogle Scholar
 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. 10.1073/pnas.0811929106PubMed CentralView ArticlePubMedGoogle Scholar
 Ottesen JT: Modelling of the baroreflexfeedback mechanism with timedelay. Journal of Mathematical Biology 1997, (36):41–63. 10.1007/s002850050089Google Scholar
 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. full_textGoogle Scholar
 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. 10.1109/TBME.2003.808805View ArticlePubMedGoogle Scholar
 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. 10.1172/JCI26322PubMed CentralView ArticlePubMedGoogle Scholar
 Vuillard L, Nicholson J, Hay RT: A complex containing β TrCP recruits Ccd34 to catalyse ubiquitination of I κ B α . FEBS Letters 1999, 455: 311–314. 10.1016/S00145793(99)008959View ArticlePubMedGoogle Scholar
 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 NFB signaling control the dynamics of gene expression. Science 2004, 306: 704–708. 10.1126/science.1099962View ArticlePubMedGoogle Scholar
 Harper N, Hughes M, MacFarlane M, Cohen GM: Fasassociated death domain protein and caspase8 are not recruited to the tumor necrosis factor receptor 1 signaling complex during tumor necrosis factorinduced apoptosis. J Biol Chem 2003, 278: 25534–25541. 10.1074/jbc.M303399200View ArticlePubMedGoogle Scholar
 Sung MH, Simon R: In silico simulation of inhibitor drug effects on nuclear factor κ B pathway dynamics. Molecular Pharmacology 2004, 66: 70–75. 10.1124/mol.66.1.70View ArticlePubMedGoogle Scholar
 Haefner B: arresting a major culprit in cancer. Drug Discovery Today 2002, 7(12):653–663. 10.1016/S13596446(02)023097View ArticlePubMedGoogle Scholar
 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. 10.1074/jbc.274.38.27307View ArticlePubMedGoogle Scholar
 Yin MJ, Yamamoto Y, Gaynor RB: The antiinflammatory agents aspirin and salicylate inhibit the activity of I κ B kinase β . Nature 1998, 396(6706):77–80. 10.1038/23948View ArticlePubMedGoogle Scholar
 Baker CTH, Paul CAH, Willé DR: Issues in the numerical solution of evolutionary delay differential equations. Advances in Computational Mathematics 1995, 3: 171–196.View ArticleGoogle Scholar
 Varma A, Morbidelli M, Wu H: Parameter sensitivity in chemical systems. Cambridge: Cambridge University Press; 1999.View ArticleGoogle Scholar
Copyright
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.