Home | About | Journals | Submit | Contact Us | Français |

**|**BMC Bioinformatics**|**v.11(Suppl 7); 2010**|**PMC2957680

Formats

Article sections

- Abstract
- Background
- Results and discussion
- Conclusions
- Methods
- EAMCM algorithm
- Competing interests
- Authors' contributions
- Supplementary Material
- References

Authors

Related links

BMC Bioinformatics. 2010; 11(Suppl 7): S12.

Published online 2010 October 15. doi: 10.1186/1471-2105-11-S7-S12

PMCID: PMC2957680

Wu Hsiung Wu: wt.ude.ucc.sc@hww; Feng Sheng Wang: wt.ude.ucc@wsfmhc; Maw Shang Chang: wt.ude.ucc.sc@gnahcsm

Ninth International Conference on Bioinformatics (InCoB2010): Bioinformatics

Christian Schönbach, Kenta Nakai, Tin Wee Tan and Shoba
Ranganathan

http://www.biomedcentral.com/content/pdf/1471-2105-11-S7-info.pdfAsia Pacific Bioinformatics Network (APBioNet) Ninth International Conference on
Bioinformatics (InCoB2010)

26-28 September 2010

Tokyo, Japan

Copyright ©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.

This article has been cited by other articles in PMC.

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.

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.

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.

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.

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

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

$$\frac{\text{d}{x}_{1}}{\text{d}{t}^{*}}=\frac{{h}_{0}}{{\epsilon}_{h}}\left[\frac{\beta {g}_{1}}{1+\gamma {g}_{2}}-\nu {g}_{2}\right]+\frac{\delta}{{\epsilon}_{h}}({h}_{0}-{x}_{1}),$$

(1)

$$\frac{\text{d}{x}_{2}}{\text{d}{t}^{*}}=\frac{\mu {p}_{0}}{{\epsilon}_{p}{h}_{0}}{x}_{1}-\frac{{x}_{2}}{{\epsilon}_{p}(1+\alpha {g}_{1})},$$

(2)

where *x*_{1 }is the heart rate, *x*_{2 }is the
arterial pressure, and *t** = *t*/*τ *is the dimensionless
time. The functions *g*_{1 }and *g*_{2 }are defined
by

$$\begin{array}{lll}{g}_{1}\hfill & =\hfill & \frac{1}{1+{[{x}_{2}({t}^{*}-1)/{p}_{0}+{r}_{1}]}^{n}},\hfill \\ {g}_{2}\hfill & =\hfill & 1-\frac{1}{1+{[{x}_{2}({t}^{*})/{p}_{0}+{r}_{2}]}^{n}},\hfill \\ {r}_{1}\hfill & =\hfill & {A}_{1}\mathrm{sin}\{2\pi {f}_{r}\tau ({t}^{*}-1)-\varphi \},\hfill \\ {r}_{2}\hfill & =\hfill & {A}_{2}\mathrm{sin}\{2\pi {f}_{r}\tau {t}^{*}-\varphi \},\hfill \end{array}$$

where *x*_{2}(*t** - 1) and *x*_{2
}(*t**) denote the blood pressure with and without a time delay,
respectively. The values and definitions of system parameters are given in Table
Table1.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 *g _{i }*with an exponent

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 Figure1.1. It is easy to find which parameter makes more effects on heart rate and
blood pressure than the others from Figure Figure1.1. The values of top
five sensitivities for the heart rate and blood pressure are shown in Table Table2.2. The uncontrolled average arterial blood pressure
(*p*_{0}), breathing rate (*f*_{r}),
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 Figure2.2. The
dynamic sensitivities of heart rate with respect to *p*_{0 }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
*p*_{0 }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
Table2,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 Figure3.3. Figure Figure33 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 Figure4.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.

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

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 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 Figure6.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 Figure66 that only some few
parameter sensitivities get significant percentage of the total sensitivity for each
state variable.

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.

All of the dynamic sensitivities with respect to *k*_{9}, the rate
constant of the formation of survival complex, and *k*_{15 }, 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 *k*_{9}, the corresponding sensitivity profile with respect
to *k*_{15 }can be obtained simply by reflecting about the time axis.
To elucidate the effects of IKK (*x*_{10}) and FADD
(*x*_{17}) on the oscillatory behavior of NF-*κ*B
(*x*_{16}) and I*κ*B (*x*_{31}) in the
survival pathway, the semi-relative sensitivities of NF-*κ*B and
I*κ*B with respect to the rate constants *k*_{9 }and
*k*_{15 }, of the formation of survival complex and death complex
are shown in Figure Figure7.7. We observe, from Figure Figure7,7, the negative regulation of oscillatory behavior of NF-*κ*B
and I*κ*B when the rate constant (*k*_{15}) or the initial
concentration of FADD increases. The reverse effect is seen by increasing the rate
constant (*k*_{9}) 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 (*x*_{23}) is the finial reaction
of the apoptosis pathway, so the concentration of the active caspase-3
(*x*_{25}) can be used as the response of the apoptosis pathway.
The semi-relative sensitivities of active caspase-3 with respect to *k*_{9
}and *k*_{15 }are shown in Figure Figure7.7. The
negative values of the semi-relative sensitivities of active caspase-3 with respect
to *k*_{9 }show that the active caspase-3 decreases when the rate
constant (*k*_{9}) or the initial concentration of IKK increases. By
contrast, increasing the rate constant (*k*_{15}) 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 Figure7
shows,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.

The fragmented DNA (*x*_{26}) 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 Table33 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 (*k*_{20}) and phosphorylated caspase-3
(*k*_{23}) and the regulation of activating caspase-3 by
phosphorylated caspase-8 (*k*_{21 }) 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 Figure88 shows the top seven semi-relative sensitivities of
fragmented DNA. All key parameters positively regulate the activity of DNA
fragmentation except parameter *k*_{9}.

The transcription factor NF-*κ*B (*x*_{16}) 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 Table4.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
-*k*_{9}, *k*_{11}, *k*_{12},
*k*_{14}, *k*_{28}, and *k*_{29 }- in
the survival pathway are identified as important parameters. The parameter
*k*_{15 }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 Figure9.9. The dynamic
sensitivity of NF-*κ*B with respect to the IKK activation
(*k*_{11}) is similar to that with respect to the
NF-*κ*B activation (*k*_{14}) as Figure Figure99 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
(*k*_{29}), 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 (*k*_{15}) 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 Figure9.9. The
sensitivities of NF-*κ*B with respect to the rates of IKK activation
(*k*_{11}), the formation of
NF-*κ*B/I*κ*B/IKK* (*k*_{12}), the
NF-*κ*B activation (*k*_{14}), and the deactivation of
NF-*κ*B (*k*_{28}) alternate between positive and
negative values. These oscillations of sensitivities may be caused by the negative
feedback regulation of I*κ*B.

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 Table5.5. The key
parameters -*k*_{9}, *k*_{15}, *k*_{17},
and *k*_{20 }- are identified and the dynamic sensitivities of DISC
with respect to these four parameters are shown in Figure Figure10.
The10.
The reaction of dissociation of DISC from the death receptor TNFR1 is essential for
the following caspase-8 activation and its corresponding rate constant
*k*_{17 }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].

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 Figure1111 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 Figure11,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 Table6.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
and a brief manual can be found in the Additional file 2.

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.

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

$$\frac{\text{d}x}{\text{d}t}=f(x(t),x(t-{\tau}_{1}),x(t-{\tau}_{2}),\mathrm{...},x(t-{\tau}_{r});\theta ),$$

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

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

$$\frac{\text{d}x}{\text{d}t}=f(x(t),x(t-\tau );\theta ),$$

(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

$$y(t)=\{\begin{array}{ll}x(0)\hfill & \text{if}\text{\hspace{0.17em}}t\le \tau ,\hfill \\ x(t-\tau )\hfill & \text{if}\text{\hspace{0.17em}}t>\tau ,\hfill \end{array}$$

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 [*t*_{j-1},
*t*_{j}], *t*_{i-1 }≤
*t*_{j-1 }<*t*_{j }≤
*t*_{i }by

$$x({t}_{j})=x({t}_{j-1})+\frac{1}{2}{\eta}_{j}\{f(x({t}_{j}),y({t}_{j});\theta )+f(x({t}_{j-1}),y({t}_{j-1});\theta )\},$$

(5)

where *η*_{j }is the step size in time
*t*_{j }. 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:

$$y(t)=x(t-\tau )=x({t}_{a})+\frac{t-\tau -{t}_{a}}{{t}_{a+1}-{t}_{a}}[x({t}_{a+1})-x({t}_{a})],\text{\hspace{0.17em}}{t}_{a}<t-\tau <{t}_{a+1},$$

(6)

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.

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

$$s({x}_{i},{\theta}_{j})=\underset{\Delta {\theta}_{j}\to 0}{\mathrm{lim}}\frac{{x}_{i}(t;{\theta}_{j}+\Delta {\theta}_{j})-{x}_{i}(t;{\theta}_{j})}{\Delta {\theta}_{j}}=\frac{\partial {x}_{i}(t;\theta )}{\partial {\theta}_{j}},$$

(7)

where *x*_{i}(*t*; *θ*_{j
}+ Δ*θ*_{j}) is the
*i*^{th }component of the solution of equation (4) with
a change Δ*θ*_{j }on the *j*^{th
}parameter and the others fixed. The absolute parameter sensitivity
*s*(*x*_{i }, *θ*_{j
}) is also defined as the first-order *local sensitivity *of
*x*_{i }with respect to *θ*_{j
}[46]. The term *local
*refers to the fact that the value of *s*(*x*_{i
}, *θ*_{j }) depends on the given set of
values for the parameters *θ*. Under the assumption that the system
responds linearly for small perturbations, *s*(*x*_{i
}, *θ*_{j }) measures the ratio between the
effect on *x*_{i }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*(*x*_{i
}, *θ*_{j }) of *x*_{i
}with respect to *θ*_{j }, a dimensionless
quantity, is defined as

$$S({x}_{i},{\theta}_{j})=\frac{\partial \mathrm{ln}{x}_{i}(t;\theta )}{\partial \mathrm{ln}{\theta}_{j}}=\frac{{\theta}_{j}}{{x}_{i}}s({x}_{i},{\theta}_{j}).$$

(8)

Sometimes the use of relative parameter sensitivities has the difficulty of numerical
instability caused by the division by zero when *x*_{i
}approaching zero. To address this problem, the *semi*-*relative
parameter sensitivity *$\tilde{S}({x}_{i},{\theta}_{j})$ of
*x*_{i }with respect to *θ*_{j
}is used and defined as

$$\tilde{S}({x}_{i},{\theta}_{j})={\theta}_{j}s({x}_{i},{\theta}_{j}).$$

(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

$$\frac{\text{d}s({x}_{i}(t),{\theta}_{j})}{\text{d}t}={\displaystyle \sum _{k=1}^{n}\frac{\partial {f}_{i}}{\partial {x}_{k}(t)}}s({x}_{k}(t),{\theta}_{j})+{\displaystyle \sum _{k=1}^{n}\frac{\partial {f}_{i}}{\partial {x}_{k}(t-\tau )}}s({x}_{k}(t-\tau ),{\theta}_{j})+\frac{\partial {f}_{i}}{\partial {\theta}_{j}},$$

(10)

where *f*_{i }is the *i*^{th }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 *x*_{i }with
respect to *θ*_{j }is computed by

$$\begin{array}{c}s({x}_{i}({t}_{l}),{\theta}_{j})=s({x}_{i}({t}_{l-1}),{\theta}_{j})\\ +\frac{1}{2}{\eta}_{l}\left\{{\displaystyle \sum _{k=1}^{n}\frac{\partial {f}_{i}({t}_{l})}{\partial {x}_{k}({t}_{l})}}s({x}_{k}({t}_{l}),{\theta}_{j})+{\displaystyle \sum _{k=1}^{n}\frac{\partial {f}_{i}({t}_{l})}{\partial {x}_{k}({t}_{l}-\tau )}}s({x}_{k}({t}_{l}-\tau ),{\theta}_{j})+\frac{\partial {f}_{i}({t}_{l})}{\partial {\theta}_{j}}\right\}\\ +\frac{1}{2}{\eta}_{l}\left\{{\displaystyle \sum _{k=1}^{n}\frac{\partial {f}_{i}({t}_{l-1})}{\partial {x}_{k}({t}_{l-1})}}s({x}_{k}({t}_{l-1}),{\theta}_{j})+{\displaystyle \sum _{k=1}^{n}\frac{\partial {f}_{i}({t}_{l-1})}{\partial {x}_{k}({t}_{l-1}-\tau )}}s({x}_{k}({t}_{l-1}-\tau ),{\theta}_{j})+\frac{\partial {f}_{i}({t}_{l-1})}{\partial {\theta}_{j}}\right\},\end{array}$$

(11)

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

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:

$${\overline{S}}_{ij}=\frac{1}{{t}_{f}}{\displaystyle {\int}_{0}^{{t}_{f}}|S({x}_{i},{\theta}_{j})|\text{d}t,}$$

(12)

where *t*_{f }is the final time.

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

$$\frac{\partial f}{\partial x}=\underset{\Delta \text{x}\to 0}{\mathrm{lim}}\frac{f(x+\Delta x)-f(x)}{\Delta x}.$$

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

$$\frac{\partial h({o}_{1},{o}_{2})}{\partial x}=\frac{\partial h}{\partial {o}_{1}}\frac{\partial {o}_{1}}{\partial x}+\frac{\partial h}{\partial {o}_{2}}\frac{\partial {o}_{2}}{\partial x},$$

(14)

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.

The proposed algorithm EAMCM is shown as follows:

Algorithm EAMCM

**Input**:

1. A set of *n *delay differential equations $\dot{x}$ = **f
**(**x, y**) with *n *dependent variables *x*_{i
}, *i *= 1, ..., *n *and *m *time-delay 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*_{i-1}, *d*_{t
}← 0, **x**^{c }←
**x**(*t*_{i-1 }), **Φ**^{c
}← **Φ**(*t*_{i-1}).

2. Repeat the following steps until *d*_{t }=
*t*_{i }- *t*_{i-1}.

(a) **x**^{p }← **x**^{c },
**y**^{p }← **x**(t_{i-1 }+ *d*_{i
}- *τ*), **Φ**^{p }←
**Φ**^{c }.

(b) If *t*_{i-1 }>*τ *, then
**Φ**^{d }←
**Φ**(*t*_{i-1}+ *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*_{i-1}*-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

Input:

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.

**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}*) ←

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}*),

(b) **x**(*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 $\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

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

**Output**: The sensitivity matrix **Φ**(**x**(*t *+
*η*_{j })).

1. Evaluate $A=\frac{\partial f(x(t),y(t))}{\partial x}$, $B\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\partial f(x(t),y(t))}{\partial y}$, and $C=\frac{\partial f(x(t),y(t))}{\partial \theta}$ by automatic differentiation.

2. Evaluate the derivative $\dot{\Phi}$(**x**(*t*))
= **AΦ**(**x**(*t*)) + **BΦ**(**y**(*t*)) +
**C**.

3. **Φ**(**x**(*t *+ *η*_{j}))
= **Φ**(**x**(*t*)) + $\dot{\Phi}$(**x**(*t*))
* *η _{j }*.

4. Evaluate $A=\frac{\partial f(x(t+\eta j),y(t+\eta j))}{\partial x}$, $B\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\partial f(x(t+\eta j),y(t+\eta j))}{\partial y}$, and $C\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\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 $\dot{\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))\text{}+\text{}0.\text{5}\ast {\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

The authors declare that they have no competing interests.

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.

**TNF- α signal transduction model**. This file includes the
set of delay differential equations of the TNF-

Click here for file^{(49K, pdf)}

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

Click here for file^{(226K, pdf)}

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.

- Fowler AC, Mackey MC. Relaxation oscillations in a class of delay differential equations. SIAM Journal on Applied Mathematics. 2002;11:299–323. doi: 10.1137/S0036139901393512. [Cross Ref]
- Monk NAM. Oscillatory expression of Hes1, p53, and NF-
*κ*B driven by transcriptional time delays. Current Biology. 2003;11(16):1409–1413. doi: 10.1016/S0960-9822(03)00494-9. [PubMed] [Cross Ref] - Nelson PW, Perelson AS. Mathematical analysis of delay differential equation models of HIV-1 infection. Mathematical Biosciences. 2002;11:73–94. doi: 10.1016/S0025-5564(02)00099-8. [PubMed] [Cross Ref]
- Bocharov GA, Rihan FA. Numerical modelling in biosciences using delay differential equations. Journal of Computational and Applied Mathematics. 2000;11(1-2):183–199. doi: 10.1016/S0377-0427(00)00468-4. [Cross Ref]
- Enright WH, Hayashi H. A delay differential equation solver based on a continuous Runge-Kutta method with defect control. Numerical Algorithms. 1997;11(3-4):349–364. doi: 10.1023/A:1019107718128. [Cross Ref]
- 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;11(5):1593–1609. doi: 10.1137/S1064827599363381. [Cross Ref]
- Paul CAH. Developing a delay differential equation solver. Applied Numerical Mathematics. 1992;11(3-5):403–414. doi: 10.1016/0168-9274(92)90030-H. [Cross Ref]
- Shampine L, Thompson S. Solving DDEs in Matlab. Appl Numer Math. 2001;11:441–458. doi: 10.1016/S0168-9274(00)00055-6. [Cross Ref]
- Shampine LF. Solving ODEs and DDEs with residual control. Applied Numerical Mathematics. 2005;11:113–127. doi: 10.1016/j.apnum.2004.07.003. [Cross Ref]
- Thompson S, Shampine LF. A friendly Fortran DDE solver. Applied Numerical Mathematics. 2006;11(3-4):503–516. doi: 10.1016/j.apnum.2005.04.027. [Cross Ref]
- Willé DR, Baker CT. DELSOL - a numerical code for the solution of systems of delay-differential equations. Applied Numerical Mathematics. 1992;11(3-5):223–234. doi: 10.1016/0168-9274(92)90017-8. [Cross Ref]
- Zivaripiran H. DDVERK90: A user-friendly implementation of an effective DDE solver. Master's thesis, University of Toronto, Canada. 2005.
- 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;11(36):5097–5114. doi: 10.1016/j.atmosenv.2003.08.020. [Cross Ref]
- 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;11:1794–1808. doi: 10.1063/1.438530. [Cross Ref]
- Dunker AM. The decoupled direct method for calculating sensitivities coefficients in chemical kinetics. J Chem Phys. 1984;11:2385–2393. doi: 10.1063/1.447938. [Cross Ref]
- Rihan FA. Sensitivity analysis for dynamic systems with time-lags. Journal of Computational and Applied Mathematics. 2003;11(2):445–462. doi: 10.1016/S0377-0427(02)00659-3. [Cross Ref]
- 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;11(36):5083–5096. doi: 10.1016/j.atmosenv.2003.08.019. [Cross Ref]
- Bischof C, Khademi P, Mauer A, Carle A. Adifor 2.0: Automatic differentiation of Fortran 77 programs. Computing in Science and Engineering. 1996;11(3):18–32. doi: 10.1109/99.537089. [Cross Ref]
- Hwang D, Byun DW, Odman MT. An automatic differentiation technique for sensitivity analysis of numerical advection schemes in air quality models. Atmospheric Environment. 1997;11(6):879–888. doi: 10.1016/S1352-2310(96)00240-3. [Cross Ref]
- 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.
- Wu WH, Wang FS, Chang MS. Dynamic sensitivity analysis of biological systems. BMC Bioinformatics. 2008;11:S17. doi: 10.1186/1471-2105-9-S12-S17. [PMC free article] [PubMed] [Cross Ref]
- Fowler AC, McGuinness MJ. A delay recruitment model of the cardiovascular control system. Journal of Mathematical Biology. 2005;11(5):508–526. doi: 10.1007/s00285-005-0339-1. [PubMed] [Cross Ref]
- Rangamani P, Sirovich L. Survival and apoptotic pathways initiated by TNF-
*α*: modeling and predictions. Biotechnol Bioeng. 2007;11(5):1216–1229. doi: 10.1002/bit.21307. [PubMed] [Cross Ref] - Ashe PC, Berry MD. Apoptotic signaling cascades. Prog Neuropsychopharmacol Biol Psychiatry. 2003;11:199–214. doi: 10.1016/S0278-5846(03)00016-2. [PubMed] [Cross Ref]
- 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;11:167–170. doi: 10.1038/376167a0. [PubMed] [Cross Ref] - Beg AA, Baltimore D. An essential role for NF-
*κ*B in preventing TNF-*α*-induced cell death. Science. 1996;11:782–784. doi: 10.1126/science.274.5288.782. [PubMed] [Cross Ref] - 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;11:839–851. doi: 10.1083/jcb.200404158. [PMC free article] [PubMed] [Cross Ref]
- 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;11:726–739. doi: 10.1177/0037549703040943. [Cross Ref]
- Fussenegger M, Bailey JE, Varner J. A mathematical model of caspase function in apoptosis. Nature Biotechnology. 2000;11:768–774. doi: 10.1038/81208. [PubMed] [Cross Ref]
- Hoffmann A, Levchenko A, Scott ML, Baltimore D. The I
*κ*B-NF-*κ*B signaling module: temporal control and selective gene activation. Science. 2002;11:1241–1245. doi: 10.1126/science.1071914. [PubMed] [Cross Ref] - Lipniacki T, Paszek P, Brasier AR, Luxon B, Kimmel M. Mathematical model of NF-
*κ*B regulatory module. Journal of Theoretical Biology. 2004;11(2):195–215. doi: 10.1016/j.jtbi.2004.01.001. [PubMed] [Cross Ref] - Barrett CJ, Guild SJ, Ramchandra R, Malpas SC. Baroreceptor denervation prevents sympathoinhibition during angiotensin II-induced hypertension. Hypertension. 2005;11:168–172. doi: 10.1161/01.HYP.0000168047.09637.d4. [PubMed] [Cross Ref]
- 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;11(3):924–928. doi: 10.1073/pnas.0811929106. [PubMed] [Cross Ref]
- Ottesen JT. Modelling of the baroreflex-feedback mechanism with time-delay. Journal of Mathematical Biology. 1997. pp. 41–63. [PubMed] [Cross Ref]
- McSharry PE, McGuinness MJ, Fowler AC. Confronting a cardiovascular system model with heart rate and blood pressure data. Computers in Cardiology. 2005;11:587–590. full_text.
- McSharry PE, Clifford GD, Tarassenko L, Smith LA. A dynamical model for generating synthetic electrocardiogram signals. IEEE Transactions on Biomedical Engineering. 2003;11(3):289–294. doi: 10.1109/TBME.2003.808805. [PubMed] [Cross Ref]
- Luo JL, Kamata H, Karin M. IKK/NF-
*κ*B signaling: balancing life and death - a new approach to cancer therapy. J Clin Invest. 2005;11:2625–2632. doi: 10.1172/JCI26322. [PMC free article] [PubMed] [Cross Ref] - Vuillard L, Nicholson J, Hay RT. A complex containing
*β*TrCP recruits Ccd34 to catalyse ubiquitination of I*κ*B*α*FEBS Letters. 1999;11:311–314. doi: 10.1016/S0014-5793(99)00895-9. [PubMed] [Cross Ref] - 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;11:704–708. doi: 10.1126/science.1099962. [PubMed] [Cross Ref]
- 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;11:25534–25541. doi: 10.1074/jbc.M303399200. [PubMed] [Cross Ref]
- Sung MH, Simon R. In silico simulation of inhibitor drug effects on nuclear factor-
*κ*B pathway dynamics. Molecular Pharmacology. 2004;11:70–75. doi: 10.1124/mol.66.1.70. [PubMed] [Cross Ref] - Haefner B. arresting a major culprit in cancer. Drug Discovery Today. 2002;11(12):653–663. doi: 10.1016/S1359-6446(02)02309-7. [PubMed] [Cross Ref]
- Yamamoto Y, Yin MJ, Lin KM, Gaynor RB. Sulindac inhibits activation of the NF-
*κ*B pathway. Journal of Biological Chemistry. 1999;11(38):27307–27314. doi: 10.1074/jbc.274.38.27307. [PubMed] [Cross Ref] - Yin MJ, Yamamoto Y, Gaynor RB. The anti-inflammatory agents aspirin and salicylate inhibit the activity of
I
*κ*B kinase-*β*Nature. 1998;11(6706):77–80. doi: 10.1038/23948. [PubMed] [Cross Ref] - Baker CTH, Paul CAH, Willé DR. Issues in the numerical solution of evolutionary delay differential equations. Advances in Computational Mathematics. 1995;11:171–196.
- Varma A, Morbidelli M, Wu H. Parameter sensitivity in chemical systems. Cambridge: Cambridge University Press; 1999.

Articles from BMC Bioinformatics are provided here courtesy of **BioMed Central**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |