 Proceedings
 Open access
 Published:
Sensitivity analysis of dynamic biological systems with timedelays
BMC Bioinformatics volume 11, Article number: S12 (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.
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
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
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 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 n is used to model the baroreflex feedback control of heart rate. This simple model has shown to be able to reproduce many of the empirical observations such as respiratory sinus arrhythmia (RSA), Mayer waves, and synchronization [35].
The EAMCM program is used to do sensitivity analysis on the lumped cardiovascular model. The nonconstant 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 runtime for solving the dynamic sensitivities of the cardiovascular system. The dynamic sensitivities of heart rate and blood pressure with respect to all system parameters and initial conditions are computed. All relative parameter sensitivities are presented by 100% stacked column chart and shown in Figure 1. It is easy to find which parameter makes more effects on heart rate and blood pressure than the others from Figure 1. The values of top five sensitivities for the heart rate and blood pressure are shown in Table 2. The uncontrolled average arterial blood pressure (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 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 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 slowacting sympathetic nerves and the fastacting vagal nerves compete with each other to increase and decrease the heart rate, respectively [36]. The relative sensitivities of heart rate with respect to parameters for slow sympathetic control (β) and fast vagal control (ν) are investigated and shown in Figure 3. Figure 3 shows the sympathovagal balance in physiology and both sympathetic control and vagal control amplify the variation of heart rate. The relative sensitivities of blood pressure with respect to parameters for slow sympathetic control (β) and fast vagal control (ν) are shown in Figure 4. The slowacting 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 fastacting vagal control downregulates the blood pressure and has a negative relative sensitivity over the time.
Apoptosis signal network
The TNFα signal transduction network was developed by Rangamani and Sirovich [23], which considers both the NFκ Bmediated survival pathway and the caspasemediated apoptosis pathway simultaneously. These two pathways involve 31 species in 19 reactions and the schematic diagram is shown in Figure 5. The formation of this network involves binding reactions between ligand and death receptor, proteinprotein 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 TNFRassociated death domain (TRADD), TNFRassociated factor 2 (TRAF2), receptorinteracting protein 1 (RIP1), and possibly other yet unidentified proteins to form the early complex. In the NFκ Bmediated 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 TNFinduced apoptosis by binding to activated caspase3 [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κ BIκ 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 caspasemediated apoptosis pathway, TRADD, RIP1, and TRAF2 are dissociated from TNFR1 and recruit Fas Associated Death Domain (FADD) and caspase8 to form a protein complex called the deathinducing signaling complex (DISC) [40]. As a result of DISC formation, caspase8 is cleaved at the DISC resulting in the activation of caspase8. The activated caspase8 in turn activates effector caspases, such as caspase3, 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 timeaveraged semirelative parameter sensitivities for each state variable are shown in Figure 6. Most of the timeaveraged semirelative parameter sensitivities for each state variable are too small compared with the largest and can be ignored. It is easy to find from Figure 6 that only some few parameter sensitivities get significant percentage of the total sensitivity for each state variable.
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 semirelative 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 7. We observe, from Figure 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 caspase3 (x_{23}) is the finial reaction of the apoptosis pathway, so the concentration of the active caspase3 (x_{25}) can be used as the response of the apoptosis pathway. The semirelative sensitivities of active caspase3 with respect to k_{9} and k_{15} are shown in Figure 7. The negative values of the semirelative sensitivities of active caspase3 with respect to k_{9} show that the active caspase3 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 caspase3. These results are in agreement with observations by mutant studies [23] and show the interplay of the apoptotic and survival pathways in response to the variations of IKK and FADD. As Figure 7 shows, a combination of increasing the initial concentration of FADD and decreasing the initial concentration of IKK gets more effects on the DNA fragmentation.
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 semirelative 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, timeaveraged semirelative sensitivity defined similar to equation (12), as the strength of effects on fragmented DNA for individual reactions. Table 3 shows the ranking of dynamic sensitivities of fragmented DNA with respect to all parameters based on timeaveraged semirelative sensitivities. The values of timeaveraged semirelative sensitivities and the corresponding ratios to the total semirelative 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 caspase8 (k_{20}) and phosphorylated caspase3 (k_{23}) and the regulation of activating caspase3 by phosphorylated caspase8 (k_{21} ) are identified as important parameters. This result indicates that the activation cascades of caspase8 and caspase3 are important reactions in the apoptosis signal transduction pathway. Figure 8 shows the top seven semirelative 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 semirelative sensitivities of NFκ B with respect to all parameters. The ranking of parameter sensitivities of NFκ B based on the timeaveraged semirelative sensitivities is shown in Table 4. A parameter is referred to as an important parameter if the ratio of its timeaveraged semirelative 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 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 9 shows. This means that the inhibitor of IKK acts the same function as the inhibitor of NFκ B to inhibit the activity of NFκ B. The activation of NFκ B in most types of cells leads to the inhibition of apoptosis, accelerates cell proliferation, and promotes tumorigenesis. To inhibit cell growth, some smallmolecule 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 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 TNFinduced apoptosis and it is required for casepase8 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 timeaveraged semirelative sensitivities is shown in Table 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 10. The reaction of dissociation of DISC from the death receptor TNFR1 is essential for the following caspase8 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].
Efficiency and accuracy
To verify the result obtained by the EAMCM algorithm, it is compared with the finite difference method using the dde23 as the DDE solver. The dde23 solver is available in MATLAB 6.5 and later. Forward difference is considered in the finite difference method. The dynamic sensitivities of these two systems mentioned above are solved by the finite difference method with spacing ratio 0.1 and 0.01, respectively. The relative parameter sensitivities of heart rate and blood pressure with respect to β obtained by the finite difference method with spacing ratio 0.1 and 0.01, respectively, and the EAMCM method are shown in Figure 11 as an illustration (another data is similar and not shown here). According to the definition of relative sensitivity, the theoretical value of relative sensitivity is obtained when the spacing ratio is approaching to zero. From Figure 11, the relative sensitivities obtained by the EAMCM are close to the theoretical values. We analyze the performance of the finite difference method and the EAMCM method for computing the dynamic sensitivities by measuring the number of evaluations of model equations. The results are shown in Table 6. The CPU time in second running by a 1.86 GHz Intel Xeon CPU with 4 GMB RAM is shown in the parenthesis for reference. Based on the comparison, the EAMCM program surely outperforms the finite difference method using dde23 solver. The EAMCM program can be accessed from http://www.che.ccu.edu.tw/~bioproc/index_english.files/page00064.htm and a brief manual can be found in the Additional file 2.
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
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
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
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
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
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_{j1}, t_{ j }], t_{i1}≤ t_{j1}<t_{ j }≤ t_{ i }by
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:
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
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
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 firstorder 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
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 semirelative parameter sensitivity \tilde{S}({x}_{i},{\theta}_{j}) of x_{ i }with respect to θ_{ j }is used and defined as
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
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
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 timeaveraged relative sensitivity is therefore defined as:
where t_{ f }is the final time.
Automatic differentiation
The main challenge in computing the solution of equation (10) is the evaluation of partial derivatives of f with respect to all state variables x, delay variables x(t  τ), and system parameters θ. One way to obtain the partial derivatives is to use symbolic differentiation tools, such as Maple and Mathematica, to perform the algebra of the differentiation. The explicit expression of the partial derivatives of f can be generated automatically. This is very useful because it saves the human effort and avoids the human errors in the analytical differentiation process. In principle, this approach gives exact values of the partial derivatives of f at the expense of high computation cost. Symbolic differentiation tools always generate lengthy formulas containing many common subexpressions that require considerable computation to evaluate. If only the values of the partial derivatives rather than the explicit expressions of the partial derivatives of f are needed, the simplest and common used approach is the numerical differentiation by finite difference approximation
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 wellformed 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:
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
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 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
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 }) ← x(t + η_{ j }  τ).

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

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

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

5.
If the iteration limitation is reached, then exit this algorithm; otherwise, return x(t + η_{ j }).
End of Algorithm Iteration
Algorithm IterationOfSensitivity
Input:

1.
A set of n delay differential equations \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)).
Output: The sensitivity matrix Φ(x(t + η_{ j })).

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.
Evaluate the derivative (x(t)) = AΦ(x(t)) + BΦ(y(t)) + C.

3.
Φ(x(t + η_{ j })) = Φ(x(t)) + (x(t)) * η_{ j } .

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
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/S0036139901393512
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)004949
Nelson PW, Perelson AS: Mathematical analysis of delay differential equation models of HIV1 infection. Mathematical Biosciences 2002, 179: 73–94. 10.1016/S00255564(02)000998
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)004684
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:1019107718128
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/S1064827599363381
Paul CAH: Developing a delay differential equation solver. Applied Numerical Mathematics 1992, 9(3–5):403–414. 10.1016/01689274(92)90030H
Shampine L, Thompson S: Solving DDEs in Matlab. Appl Numer Math 2001, 37: 441–458. 10.1016/S01689274(00)000556
Shampine LF: Solving ODEs and DDEs with residual control. Applied Numerical Mathematics 2005, 52: 113–127. 10.1016/j.apnum.2004.07.003
Thompson S, Shampine LF: A friendly Fortran DDE solver. Applied Numerical Mathematics 2006, 56(3–4):503–516. 10.1016/j.apnum.2005.04.027
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)900178
Zivaripiran H: DDVERK90: A userfriendly 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: IInumerical validation and applications. Atmospheric Environment 2003, 37(36):5097–5114. 10.1016/j.atmosenv.2003.08.020
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.438530
Dunker AM: The decoupled direct method for calculating sensitivities coefficients in chemical kinetics. J Chem Phys 1984, 81: 2385–2393. 10.1063/1.447938
Rihan FA: Sensitivity analysis for dynamic systems with timelags. Journal of Computational and Applied Mathematics 2003, 151(2):445–462. 10.1016/S03770427(02)006593
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.019
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.537089
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)002403
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.
Wu WH, Wang FS, Chang MS: Dynamic sensitivity analysis of biological systems. BMC Bioinformatics 2008, 9: S17. 10.1186/147121059S12S17
Fowler AC, McGuinness MJ: A delay recruitment model of the cardiovascular control system. Journal of Mathematical Biology 2005, 51(5):508–526. 10.1007/s0028500503391
Rangamani P, Sirovich L: Survival and apoptotic pathways initiated by TNF α : modeling and predictions. Biotechnol Bioeng 2007, 97(5):1216–1229. 10.1002/bit.21307
Ashe PC, Berry MD: Apoptotic signaling cascades. Prog Neuropsychopharmacol Biol Psychiatry 2003, 27: 199–214. 10.1016/S02785846(03)000162
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/376167a0
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.782
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.200404158
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/0037549703040943
Fussenegger M, Bailey JE, Varner J: A mathematical model of caspase function in apoptosis. Nature Biotechnology 2000, 18: 768–774. 10.1038/81208
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.1071914
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.001
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.d4
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.0811929106
Ottesen JT: Modelling of the baroreflexfeedback mechanism with timedelay. Journal of Mathematical Biology 1997, (36):41–63. 10.1007/s002850050089
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_text
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.808805
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/JCI26322
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)008959
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.1099962
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.M303399200
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.70
Haefner B: arresting a major culprit in cancer. Drug Discovery Today 2002, 7(12):653–663. 10.1016/S13596446(02)023097
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.27307
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/23948
Baker CTH, Paul CAH, Willé DR: Issues in the numerical solution of evolutionary delay differential equations. Advances in Computational Mathematics 1995, 3: 171–196.
Varma A, Morbidelli M, Wu H: Parameter sensitivity in chemical systems. Cambridge: Cambridge University Press; 1999.
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.
Author information
Authors and Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
WHW developed and implemented the algorithm and drafted the manuscript. FSW conceived of the study, participated in its design and coordination, and helped to draft the manuscript. MSC assisted in developing the algorithm and finalizing the manuscript. All authors read and approved the final manuscript.
Electronic supplementary material
12859_2010_4250_MOESM1_ESM.pdf
Additional file 1: TNF α signal transduction model. This file includes the set of delay differential equations of the TNFα signal transduction model, all of the relevant definitions of state variables, and the nominal values of parameters appearing in the delay differential equations. (PDF 49 KB)
12859_2010_4250_MOESM2_ESM.pdf
Additional file 2: Manual of EAMCM program. This document provides a brief introduction to the extended adaptive modified collocation method (EAMCM) that we implemented to compute the solution and dynamic sensitivities of ordinal differential equation (ODE) and delay differential equation (DDE) systems. (PDF 226 KB)
Rights and permissions
Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Wu, W.H., Wang, F.S. & Chang, M.S. Sensitivity analysis of dynamic biological systems with timedelays. BMC Bioinformatics 11 (Suppl 7), S12 (2010). https://doi.org/10.1186/1471210511S7S12
Published:
DOI: https://doi.org/10.1186/1471210511S7S12