SpringerOpen Newsletter

Receive periodic news and updates relating to SpringerOpen.

Open Access Research

Dynamical modeling of drug effect using hybrid systems

Xiangfang Li1*, Lijun Qian2 and Edward R Dougherty134

Author Affiliations

1 Department of Electrical and Computer Engineering, Texas A&M University, College Station, TX 77843, USA

2 Department of Electrical and Computer Engineering, Prairie View A&M University, Prairie View, TX 77446, USA

3 Computational Biology Division, Translational Genomics Research Institution, Phoenix, AZ 85004, USA

4 Department of Bioinformatics and Computational Biology, University of Texas M.D. Anderson Cancer Center, Houston, TX 77030, USA

For all author emails, please log on.

EURASIP Journal on Bioinformatics and Systems Biology 2012, 2012:19  doi:10.1186/1687-4153-2012-19

The electronic version of this article is the complete one and can be found online at: http://bsb.eurasipjournals.com/content/2012/1/19


Received:30 April 2012
Accepted:12 November 2012
Published:26 December 2012

© 2012 Li et al.; licensee Springer.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Drug discovery today is a complex, expensive, and time-consuming process with high attrition rate. A more systematic approach is needed to combine innovative approaches in order to lead to more effective and efficient drug development. This article provides systematic mathematical analysis and dynamical modeling of drug effect under gene regulatory network contexts. A hybrid systems model, which merges together discrete and continuous dynamics into a single dynamical model, is proposed to study dynamics of the underlying regulatory network under drug perturbations. The major goal is to understand how the system changes when perturbed by drugs and give suggestions for better therapeutic interventions. A realistic periodic drug intake scenario is considered, drug pharmacokinetics and pharmacodynamics information being taken into account in the proposed hybrid systems model. Simulations are performed using MATLAB/SIMULINK to corroborate the analytical results.

Keywords:
Drug effect; Hybrid systems; PK/PD; Gene regulatory network (GRN); Dosing regimens

Introduction

The ultimate goal of drug therapy is to modulate the phenotypic behavior of cells by altering the behavior of the gene and protein components of the cell [1]. This approach is possible because the phenotypic behavior of the cell reflects the dynamics of the gene and protein-based regulatory network. When it comes to drug therapeutics and disease modeling, the major goal is to understand how the system changes when perturbed and how to modify the system to achieve a desired outcome. To understand and exploit the complicated mapping between genome and phenome, especially in the context of drug discovery, it is critical to evaluate the regulatory interactions between the genes and proteins that form the gene regulatory network (GRN). To date, the hope of the rapid translation of “genes to drugs” has foundered on the reality that disease biology is complex and drug development must be driven by insights into biological responses [2]. A systems approach is crucial for moving biology from a descriptive to a predictive science [3,4]. This calls for appropriate modeling to establish a functional understanding of disease–drug interaction, in order to better predict drug effects and make drug discovery a faster and more systematic process.

Pharmacokinetics (PK) is the study of what the body does to the drug, i.e., the absorption, distribution, metabolism, and excretion of the drug, and pharmacodynamics (PD) seeks to study what the drug does to the body. A salient challenge is to link a drug’s PK information with PD characteristics to provide a better understanding of the time course of drug effect (PK/PD) after drug administration [5]. Modeling and simulation tools are required to integrate PK and PD data and optimize drug regimens.

A salient problem is finding a dosing regimen of a drug candidate that is both efficacious and safe [6]. Traditionally, drugs have been administered on an experimental basis, but it is virtually impossible to optimize dosing regimens using strictly empirical methods, especially since different patients may respond differently to the same drug dosage [7]. Moreover, traditionally designing the dosing regimen to achieve some desired target goal such as relatively constant serum concentration may not be optimal because of underlying dynamic biological networks. For example, Shah et al. [8] demonstrate that BCR–ABL inhibitor dasatinib, which has greater potency and a short half-life, can achieve deep clinical remission in CML patients by achieving transient potent BCR–ABL inhibition, while traditionally approved tyrosine kinase inhibitors usually have prolonged half-lives that result in continuous target inhibition. A similar study of whether short pulses of higher dose or persistent dosing with lower doses has the most favorable outcomes has been carried out by Amin et al. [9] in the setup of inactivation of HER2–HER3 signaling. Finding an optimal dosing regimen based on the dynamics of biological systems and relevant PK/PD information is critically important.

System modeling is emerging as a valuable tool in therapeutics to address these challenges [3,10-12]. The process begins with building a quantitative model of a biological system. Consequences of particular perturbations, such as optimal dosing regimens, optimal drug targets, or combinational therapy, can be simulated in time courses using such models. In this study, we propose a hybrid systems model for GRNs and incorporate a drug’s PK and PD information by using a state-space approach. We first study drug effect assuming the drug target to be a gene or protein in the proposed drug perturbation model using dynamical system theory, considering the case of periodic drug intake and analytically deriving the conditions for the drug to be effective. We extend the analysis to the 2-gene case and then to the case of a network with multiple coupled genes and positive feedback loops. Simulations are performed using MATLAB/SIMULINK to supplement our analytical results.

Model formulation

While discrete modeling leaves out many details, continuous modeling includes so many details that computational demands preclude their applications to many larger systems. Hybrid systems, which aim to merge ideas from both continuous and discrete modeling into one paradigm, are appealing for GRN modeling under drug perturbations because biological systems are naturally nonlinear, have highly varied regulatory requirements, and possess a wide range of control strategies for meeting their needs. While some simple, local, feedback control methods can provide sufficient regulation of many more-or-less continuous cellular processes, the regulation of discontinuous processes possessing the character of computational decision making requires more elaborate regulatory methods [13]. In particular, some genes display regulation in a thresholded switch-like manner [14].

Hybrid systems include a broad space of models and systems. Several hybrid systems models have been developed for biological networks. Some of these have been used to perform reachability analysis to elucidate biologically meaningful properties. For example, the Lac operon system has been well studied both experimentally and using continuous models [15,16]. A hybrid model and use of a reachability algorithm were validated by comparison with experimental data and continuous models [17]. Other biological hybrid systems analyzed in similar ways include the Delta-Notch decision process [18,19], GRNs of carbon starvation [20], and nutritional stress response [21] in Escherichia coli. As far as we know, the only hybrid systems modeling concerning treatment or drug effects is contained in our earlier work [22].

Gene regulation can be modeled by rate equations expressing the difference between rate of production and rate of degradation [23,24]. We adopt the general model

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M1">View MathML</a>

(1)

where xi≥0 corresponds to the concentrations of proteins encoded by genes in the network and can be interpreted as the gene expression level. fi(.) is a general nonlinear function and represents the rate of synthesis. It can be approximated by a sigmoidal function or a unit step function, and unit step function is used in this article. γixi is the rate of degradation. To use hybrid systems and incorporate drug effect, we propose the following model for a GRN of N genes under drug perturbation:

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M2">View MathML</a>

(2)

where the last two terms on the right-hand side of Equation (2) correspond to drug perturbation u. β > 0 and γ > 0 are synthesis and degradation rates, respectively. Ki ≥ 1 is an integer representing the number of activation/synthesis terms. <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M3">View MathML</a> and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M4">View MathML</a> describe how other genes affect gene i. They are the functions of <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M5">View MathML</a> and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M6">View MathML</a>, where s+(.) is the unit step function, s(.) = 1−s+(.), and the θ terms are the corresponding threshold values. For each gene j, the set of threshold values related to gene i is denoted by <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M7">View MathML</a>, where t+i and ti are indices of the threshold values, <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M8">View MathML</a>, and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M9">View MathML</a>. Ψi and Φi represent the two sets of genes that affect the expression of gene i in different manners. Specifically, in this article, we consider

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M10">View MathML</a>

(3)

with <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M11">View MathML</a> defined similarly. <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M12">View MathML</a> and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M13">View MathML</a> may be set to 0 or 1, or different forms when appropriate threshold values are chosen. For example, <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M14">View MathML</a> and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M15">View MathML</a>. <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M16">View MathML</a> and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M17">View MathML</a> describe how the drug u affect gene i. <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M18">View MathML</a> and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M19">View MathML</a> are the synthesis and degradation factors of the drug on gene i. <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M20">View MathML</a> and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M21">View MathML</a> are used when the drug is activating or repressing certain genes, respectively. Since most drugs are used to repress genes, only <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M22">View MathML</a> is considered in the examples of this article. Note that γu is defined as a drug-effect factor, which is closely related to the drug pharmacology model discussed in the following section.

It should be kept in mind that the focus of this article is studying the effect of dosing, in particular, dosing regimens, on the expression of genes involved in a pathology by using hybrid systems theory. Whereas the simpler Equation (1) is widely accepted, it does not contain drug-effect terms. Equation (2) extends Equation (1) by including such terms. While the structure is intuitively reasonable and somewhat general, the actual details of the drug-effect terms are unknown. Finding the specific form of Equation (2) for a specific disease is a system identification problem, which is quite distinct from the analysis problem addressed in this article. We are addressing optimization of treatment intervention, given the system. The details of our analysis might change when the details of Equation (2) are clarified, but we expect that the hybrid systems approach taken in the article will go through with appropriate modifications in the mathematical details.

We consider a 2-gene example to illustrate the feasibility of using hybrid systems for modeling drug effect. Specifically, we assume that there are two interactive genes x1,x2 that repress each other, and x1 is a disease gene which loses its self-regulation. We also assume that a drug targets x1 by reducing its expression level and providing a negative feedback term <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M23">View MathML</a>. The resulting 2-gene network is given by

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M24">View MathML</a>

(4)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M25">View MathML</a>

(5)

where β1, β2 are synthesis factors, γ2 is a degradation factor, and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M26">View MathML</a> are threshold values. <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M27">View MathML</a> is a drug-effect factor. Using dynamical systems theory, the state-trajectory schematic diagrams of this 2-gene network without and with drug input are obtained and plotted in Figures 1 and 2, respectively. It is observed that without drug input, the gene expression level of x1 increases unbounded, while with proper drug input, <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M28">View MathML</a>, the system converges to a new steady state, <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M29">View MathML</a>.

thumbnailFigure 1. State trajectory schematic of 2-gene example without drug intake.

thumbnailFigure 2. State trajectory schematic of 2-gene example with drug intake.

We assume periodic drug intake and the drug concentration level in the effect-site follows exponential decay during each period τi, i.e., <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M30">View MathML</a>, where kτi ≤ t ≤ (k + 1)τi and λd is the degradation factor. The response of gene expression levels of the two genes under periodic drug intake is shown in Figure 3. The state-space trajectory of gene expression level of x1 vs. the drug concentration level u is given in Figure 4. A comparison of trajectory of the gene expression level x1 vs. x2 with and without drug are provided in Figures 5 and 6, respectively. It is observed that the drug is quite effective in bringing down the expression level of x1. The simulation study matches the theoretical analysis, as in Figures 1 and 2, that with proper drug intervention the system will converge to a new steady state, <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M31">View MathML</a> and x2 = β2 / γ2 = 1, while x2 → 0 and x1 →  without drug input.

thumbnailFigure 3. The state response under periodic drug intake.

thumbnailFigure 4. The state-space trajectory under periodic drug intake. Parameter setting of Figures 3 and 4: <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M32">View MathML</a>.

thumbnailFigure 5. The state-space trajectory with drug:<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M34">View MathML</a>. The rest parameter settings are the same with Figures 3 and 4.

thumbnailFigure 6. The state-space trajectory without drug intake. The rest parameter settings are the same with Figures 3 and 4.

Pharmacology model

The basis of clinical pharmacology is the fact that the intensities of many pharmacological effects are functions of the amount of drug in the body and, more specifically, the concentration of drug at the effect-site [5]. Historically, PK and PD were considered as separate disciplines; however, the information provided by these disciplines is limited if regarded in isolation [25]. A drug-effect factor γu is included in our proposed model (Equation 2), which is related to drug’s PD characteristic (concentration–response) and its PK information (dose–concentration). In order to describe the time course of drug effect in response to different dosing regimens, the integrated PK/PD model is indispensable because it builds the bridge between these two classical disciplines of pharmacology [25]. Following each dosing regimen, instead of a two-dimensional PK and PD relationship, the proposed approach enables a description of a three-dimensional dose–concentration–effect relationship. Specifically, PK and PD are linked through γu by a state-space approach to facilitate the description and prediction of the time course of drug effects resulting from different drug administration regimens.

Drug concentration–response curve: PD model

In general, the magnitude of a pharmacological effect increases monotonically with increased dose, eventually reaching a plateau level where further increase in dose has little additional effect [6]. The classic and most commonly used concentration–response model is the Hill equation [26], also known as the sigmoidal Emax model [27] or logistic model [28]. The relationship between the concentration of the drug and its effect is most often nonlinear. In this study, we use hybrid systems to approximate PD curves. A common method is to replace certain slowly changing variables by their piecewise linear approximations (see Figure 7). For example, the PD model used in our study can approximate the popular sigmoidal Emax model (see Figure 8). The Emax model has the general form <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M35">View MathML</a>, where Emax is the maximum effect, C is the concentration, EC50 is the concentration necessary to produce 50% of Emax, and m represents a sigmoidity factor or steepness of the curve.

thumbnailFigure 7. The PD model: concentration-response curve used in this study.

thumbnailFigure 8. Sigmoidal Emax model (m=4), and approximation by our PD model.

We assume a threshold of concentration below which the drug candidate is ineffective, the minimum effective dose (MinED), and another threshold value, called maximum effective dose (MaxED), above which there is no clinically significant increase in pharmacological effect in this study. As an example, we use a linear curve to approximate the concentration–response curve between MinED and MaxED. It is assumed that the drug-effect coefficient <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M36">View MathML</a> (the drug target is x1) is related to the concentration u through a sigmoid function and can be approximated by the curve shown in Figure 7. The corresponding relationship can be expressed as

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M37">View MathML</a>

(6)

where <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M38">View MathML</a> is the ratio between the drug-effect factor <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M39">View MathML</a> and the effective drug concentration <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M40">View MathML</a> in the linear range. This reflects the fact that the drug only starts to take effect when its concentration level is above a lower threshold <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M41">View MathML</a> and its effect saturates when its concentration level exceeds an upper threshold <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M42">View MathML</a>. Note that the sigmoidal Emax model can be well approximated by the proposed PD model. By taking the derivative of E with respect to C and evaluating it at EC50, we obtain the slope as <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M43">View MathML</a>. The upper and lower bounds should satisfy <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M44">View MathML</a>. An example of the sigmoidal Emax model when m=4 and our proposed PD model are plotted together in Figure 8, where the proposed model closely resembles the sigmoidal Emax model.

Periodic drug intake: PK model

Drug concentration at the effect-site is critical for its pharmacological effect. Currently, plasma drug concentrations are markers that serve as surrogates for drug concentration at the effect-site for beneficial and adverse effects; however, markers not grounded on a sound theoretical foundation and therapeutic mechanism-based intervention can limit the usefulness of PK/PD modeling to drug development. For example, it has been demonstrated that the intracellular PK of a drug is quite different from plasma drug concentration [29,30]. As observed in the study by Kuh et al. [29], the intracellular concentration of a drug will exponentially increase as the drug is absorbed after each drug intake. The drug concentration may change very slowly (in our model, we approximate that as a flat curve) when the intracellular and extracellular drug concentration approach equilibrium. In time, drug concentration will exponentially decrease as the rate at which it is eliminated is more than the rate at which it enters the effect-site and, as a result, effects diminish.

Based on the study by Kuh et al. [29], a general model for drug concentration-time profile is given in Figure 9. Drug concentration is plotted on a logarithmic scale against time following each periodic drug intake. λa denotes the exponential increase quotient; λd is the exponential decrease quotient; τ is the interval between each drug intake; and p1, p2, and p3 denote the time stayed in the increase, equilibrium, and decrease stage, respectively. Different drugs work in different ways and the proposed model is general enough to cover various cases. Drug concentration may increase very quickly and, as a result, the increase stage may be neglected, or the equilibrium stage may be very short and can be ignored for simplicity. By adjusting the parameters in the proposed model, specific drug characteristics can be represented. In the case when the proposed model cannot approximate a drug’s PK profile, extensive simulations can be performed based on the drug’s actual PK profile. In this article, we consider a periodic drug intake scenario. Specifically, we are interested in investigating and comparing the following two potential scenarios: large dose with a longer interval versus small dose with a shorter interval.

thumbnailFigure 9. A general drug concentration-time profile.

Mathematical analysis of drug effect

In this section, we study the time course of drug effect for different dosage and schedule arrangements where the drug is designed to repress a “target gene”. The case with a special PK profile (drug concentration only has exponential decay) was analytically studied in our previous work [22]. In this study, we extend the analysis considering a general PK profile given in Figure 9 and PD model given in Figure 8. Closed-form analytical solution is provided and simulations are performed to validate theoretical analysis. In later sections, we show that the same methodology can be applied to interactive genes, where not only will the drug affect the gene expression level, but the target gene is also coupled with other genes.

It is assumed that the disease gene has lost part of its self-regulation capacity and the dynamical equation of the expression level x1 is given by

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M45">View MathML</a>

(7)

There is a steady state x1 = β1 / γ1 in such a system, however, if the synthesis rate is much bigger than the self-degradation rate, β1 ≫ γ1, then the gene expression level will be too high. A drug is used as a control input to repress the target gene expression level. The corresponding dynamical equation after drug intake is changed to

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M46">View MathML</a>

(8)

where <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M47">View MathML</a> is the drug-effect factor defined in the previous section. After incorporating a drug’s PK/PD (Figures 7 and 9) into our proposed hybrid system model Equation (8), considering the scenario that the patient is taking the drug periodically, the resulting model is given by

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M48">View MathML</a>

(9)

where for kτi ≤ t ≤ (k+1)τi and i = 1,2,… denoting the index of different dosing regimens, <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M49">View MathML</a>, <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M50">View MathML</a>, <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M51">View MathML</a> for any i, since we assume that the same drug is taken in different dosage and schedule settings. ζi is the highest concentration level reached after taking the drug.

State-space analysis

The state-space and a sample trajectory schematic of the state (target gene expression and drug concentration level) under periodic drug intake are shown in Figure 10. As is common in hybrid systems, there are both continuous quantitative changes and discrete transitions in our proposed model. The entire state space may be divided into different domains according to the value of discrete state. As shown in Figure 10, there are five domains in the state space, with D1,D3,D5 not being transient domains. The figure shows the case when the drug is effective and the drug dosage is large enough that ζi is higher than the upper threshold <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M52">View MathML</a>. The sample trajectory of the state corresponds to two periods of drug intake (numbers 1–6 corresponding to the junctions of the drug concentration and the boundaries of the domains, also marked in Figure 9). Another possible scenario is that the drug dosage is not large enough that ζi is between the upper threshold <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M53">View MathML</a> and the lower threshold <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M54">View MathML</a>. The third scenario is the case that <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M55">View MathML</a> and the drug is not effective.

thumbnailFigure 10. The state trajectory schematic (target gene expression versus drug concentration) under PK profile (Figure 9) assuming<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M57">View MathML</a>.

When the state transits in each period under periodic drug intake, it may pass through different domains (depending on the changes of drug concentration along time). During the transit time through domains D5 and D3, the gene expression level is pushed lower (to the left), while the driving strength will depend on the drug’s PD characteristic. During the transit time through D1, the expression level will rise (to the right), since the drug concentration is lower than <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M58">View MathML</a>. For the drug to be effective, the reduction of the expression level in D5 and D3 has to be larger than the increase of the expression level in D1. In summary, we should have x1((k+1)τ) < x1(kτ), so that after each treatment the expression level x1 will decrease.

State trajectory analysis

We analyze the drug effect considering the scenario shown in Figure 10, where <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M59">View MathML</a>. The same methodology can be applied to a simpler scenario where <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M60">View MathML</a>. We divide the state trajectory in a period kτi ≤ t ≤ (k+1)τi into stages a, b, c, d, e, f, and g as marked in Figure 10 and examine the drug effect stage-by-stage. The time notations used in the derivation are given by

t1: the traveling time from the initial state to the boundary between D3 and D1.

t2: from initial state to boundary between D5 and D3.

t3: from initial state to the end of stage c, t3 = kτi+p1.

t4: time at which the drug concentration starts to decrease, t4 = kτi+p1+p2.

t5: from the initial state to the end of stage e.

t6: from the initial state to the end of stage f.

For kτi ≤ t ≤ (k+1)τi, where i is the index for different dosing regimens, the corresponding equations and solutions for each stage are given by:

Stage (a) - D1 (kτi ≤ t ≤ t1):

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M61">View MathML</a>

(10)

Stage (b) - D3 (t1 t t2):

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M62">View MathML</a>

(11)

Stage (c) - D5 (t2  ≤ t ≤ t3 = kτi+p1):

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M63">View MathML</a>

(12)

Stage (d) - D5 (t3 ≤ t ≤ t4 = kτi + p1 + p2):

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M64">View MathML</a>

(13)

Stage (e) - D5 (t4 ≤ t ≤ t5):

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M65">View MathML</a>

(14)

Stage (f) - D3 (t5 ≤ t ≤ t6):

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M66">View MathML</a>

(15)

Stage (g) - D1 (t6 ≤ t ≤ (k+1)τi):

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M67">View MathML</a>

(16)

We can deduce the necessary and sufficient condition for the effectiveness of the drug by expressing the inequality x1((k+1)τ) < x1(kτ) in terms of dosing period τ and unit dose, assuming the dosage is proportional to the maximum drug concentration ζi reached after taking the drug. When the initial conditions are x1 = x1(kτi), the equations governing the state trajectory from time kτi to time (k+1)τi are given by

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M68">View MathML</a>

(17)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M69">View MathML</a>

(18)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M70">View MathML</a>

(19)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M71">View MathML</a>

(20)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M72">View MathML</a>

(21)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M73">View MathML</a>

(22)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M74">View MathML</a>

(23)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M75">View MathML</a>

(24)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M76">View MathML</a>

(25)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M77">View MathML</a>

(26)

For the drug to be effective, we need the disease gene expression level to decrease following each period of drug intake. Hence, we can express x1((k+1)τ) < x1(kτ) in terms of dosage and frequency schedule and derive the region where the drug is effective using the above listed equations.

Results and analysis

Based on the theoretical analysis in previous two sections, we demonstrate that the drug efficacy depends on total drug intake, different dosages, and frequencies. The density of drug intake is defined as <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M78">View MathML</a>. It is proportional to the total drug intake, and hence, related to the drug toxicity level in practice. First, we demonstrate the effect of total drug intake (equivalently, α) towards drug efficacy. For each fixed total drug intake, we plot the target gene expression reduction based on Equations 17 to 26 for different dosing period τ as a curve in Figure 11. It is observed that the curve is “U” shaped and there exist “sweet spot” for certain dosages and schedules given a fixed α. If we define drug efficacy region (DER) as the drug drives down the target gene expression by more than a desired percentage (say 60% in this case), it is demonstrated that the DER is related to the total drug intake, dosing period τ and dosage. DER is marked by the shaded area in Figure 11 for the case that α ≤ 0.5. It is also observed that when α gets bigger, which indicating more toxicity, DER is getting bigger accordingly. We would like to emphasize that toxicity is one of the primary causes for drug attrition and long development cycle times [31]. If a drug’s toxicity profile is available, for example, the maximum dosage and maximum exposure (α), we can find a good compromise between toxicity and drug efficacy based on such study, and determine the “sweet spot” (a good dosage and schedule balance) given the obtained α, and hence provide valuable suggestions of the dosing regimens to the clinical study.

thumbnailFigure 11. The percentage change of the disease gene expression versus the period of drug intake τ for α = 0.4,0.5,0.6,0.7,0.8, respectively. The change must be significant (say at least 60%) and α must be smaller than the acceptable toxic level. Parameters used: <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M79','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M79">View MathML</a>.

Second, we test the analytical results via numerical simulation using MATLAB/SIMULINK. Given a fixed total drug intake, or equivalently, a fixed density of drug intake (α = 0.4), three typical scenarios are studied by simulation: small frequent drug intake with τ = 7 and dose = 2.8; big infrequent drug intake with τ = 22 and dose = 8.8; and intermediate dosage and frequency with τ = 12 and dose = 4.8. The results are shown in Figure 12a–f, with the first row corresponding to the state responses and the second row corresponding to the state space trajectory. Although the three cases have the same total drug intake and initial condition (initial gene expression level x(0) = 20), the drug efficacy is different. In the small frequent intake case, the dosage is small and the drug concentration is mostly changing between domains D3 and D1. Figure 12d shows that the state-space trajectory settles in a small limit cycle and disease gene expression level settles at 11.5 at the end of each period of treatment. On the other hand, the big infrequent drug intake case results in a big limit cycle as in Figure 12f. Although the dosage is high, the long period between dosages means that the period stayed in D1 is getting longer (where drug concentration is below <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M80">View MathML</a>, hence not effective), and disease gene expression level settles at 12.1 at the end of each period of treatment. As a comparison, the drug effect for the case with intermediate dosage and frequency shown in Figure 12b,e is superior to the other two cases. Disease gene expression level settles at 8.8 at the end of each treatment period. If we check the curve in Figure 11 with α = 0.4, the intermediate dosage case with τ = 12 is located near the bottom of the “U” shape. Lastly, we observe that all state-space trajectories follow the state trajectory schematic in Figure 10, as predicted by the analytical results.

thumbnailFigure 12. Drug response at three different schedules but same drug intake: (a–c) the state response at τ = 7,12,22, respectively; (d–f) the state space trajectory at τ = 7,12,22, respectively. Other parameters are <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M81','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M81">View MathML</a>.

Analysis of 2-gene networks

We extend the theoretical analysis to 2-gene networks and show that the same framework applies to the modeling and analysis of drug effect in more complex gene networks. We assume that x1 is the target gene, there exists a positive feedback loop between x1 and another gene x2, and that a drug targets x1 by reducing its expression level and providing a negative feedback term <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M82','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M82">View MathML</a>. The resulting 2-gene network is shown in Figure 13 and is given by

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M83','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M83">View MathML</a>

(27)

thumbnailFigure 13. A 2-gene network with positive feedback loop and drug input.

where β1, β2 are synthesis factors, γ1 and γ2 are degradation factors, and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M84','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M84">View MathML</a> is a drug-effect factor. η1 > 0 and η2 > 0 are the parameters of the positive feedback loop between the two genes. The 2-gene network under drug perturbation model, Equation (27), can be rewritten as a second-order ODE:

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M85','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M85">View MathML</a>

(28)

The solution of this equation is given by

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M86">View MathML</a>

(29)

where λ1 and λ2 are the two eigenvalues of Equation 28 and k1 and k2 are parameters depending on the initial conditions. Letting <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M87">View MathML</a>, <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M88">View MathML</a> and d = β1γ2 + β2η1, the two eigenvalues are given by <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M89','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M89">View MathML</a>. It is easy to verify that <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M90','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M90">View MathML</a>. Since a2−4b > 0, we conclude that λ1 ≠ λ2 and both eigenvalues are real. Furthermore, one of the eigenvalues, λ2, is always negative since <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M91">View MathML</a>. The sign of λ1 will be determined by the sign of b:

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M92','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M92">View MathML</a>

(30)

In other words,

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M93">View MathML</a>

(31)

The above equation has an important biological interpretation: when the degradation of x1 due to the strength of the drug is faster than the increase of x1 due to the positive feedback loop, both eigenvalues are negative, the system is stable and x1 will experience exponential decay; on the other hand, if the effect of the positive feedback loop is dominant, then one of the eigenvalues will be positive and x1 will increase exponentially.

Given initial condition x1(t0) and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M94','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M94">View MathML</a>, then for the case λ1 ≠ λ2, we have

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M95','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M95">View MathML</a>

(32)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M96">View MathML</a>

(33)

Now with the baseline analysis of the second-order system, we provide detailed state trajectory analysis by taking into account the practical form of PK/PD (<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M97','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M97">View MathML</a>) when the drug is taken periodically.

State trajectory analysis

We analyze the drug-effect following the same framework given in the subsection “State trajectory analysis” under the main section “Mathematical analysis of drug effect”. For kτi ≤ t ≤ (k+1)τi, i = 1,2,…, the corresponding equations and solutions for each stages are given as follows:

Stage (a) - D1 (kτi ≤ t ≤ t1):

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M98','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M98">View MathML</a>

(34)

• The solution of x1(t) is given by Equation (29), with k1 and k2 given by Equations (32) and (33), and t0 = kτi.

Stage (b) - D3(t1 ≤ t ≤ t2):

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M99','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M99">View MathML</a>

(35)

• where a, b, d are defined as before. When incorporating the practical form of <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M100','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M100">View MathML</a> and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M101','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M101">View MathML</a>, the above second-order ODE has no closed-form solution. In this case, the solution can be obtained numerically.

Stage (c) - D5 (t2 ≤ t ≤ t3 = kτi + p1): The set of equations are the same as in Stage (b) except that <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M102','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M102">View MathML</a>. Since <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M103','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M103">View MathML</a> does not depend on <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M104','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M104">View MathML</a> explicitly, x1 has a closed-form solution given by Equation (29).

Stage (d) - D5 (t3 ≤ t ≤ t4 = kτi + p1 + p2): The solution of x1 is the same as that in Stage (c) except the start and end times, and the equation of u, which is <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M105','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M105">View MathML</a> in this stage.

Stage (e) - D5 (t4 ≤ t ≤ t5): The solution of x1 is the same as in Stage (c) except the start and end times, and the equation of u, which now is <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M106','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M106">View MathML</a>.

Stage (f) - D3 (t5 ≤ t ≤ t6): The solution of x1 is the same as in Stage (b) except the start and end times, and the equation of u, which now is <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M107','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M107">View MathML</a>.

Stage (g) - D1 (t6 ≤ t ≤ (k+1)τi): The solution of x1 is the same as in Stage (a) except the start and end times, and the equation of u, which now is <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M108','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M108">View MathML</a>.

We can deduce the necessary and sufficient condition for the effectiveness of the drug by expressing the inequality x1((k+1)τ) < x1(kτ) in terms of dosing period τ and unit dose. In the 2-gene case, no explicit closed-form expression can be deduced for the solutions in stages (b) and (f) and numerical methods have to be applied. However, through such analysis, it is observed that the same methodology for analyzing drug effect can be extended to GRNs with multiple interactive genes, although the mathematics involved will become more complicated and sometimes numerical methods must be applied when there is no closed-form solution.

Simulation results and analysis

When drug input is not present, the disease gene expression will grow unbounded owing to the positive feedback loop between the two genes. Here, we study response of the disease gene expression to drug input and compare two different schedules for τ = 20 and τ = 30, keeping α = 0.8 fixed. The response and state trajectories in 2D and 3D are given in Figure 14a–f, with the first and second rows corresponding to τ = 20 and τ = 30, respectively. We observe that both cases have periodic responses, but the disease gene expression is much better controlled when τ = 20. This is because the drug concentration is high enough in both cases compared to the threshold (<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M109','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M109">View MathML</a>), while the decay of the drug concentration is shorter in the case when τ = 20. In Figure 14c,f, the 3D state (disease gene expression) trajectories show that the trajectory settles to an inner circle when τ = 20, whereas the trajectory settles to an outer circle when τ=30. Similar observations apply to Figure 14b,e. Note the scale of x-axis of Figure 14e is 20 times bigger than that of Figure 14b.

thumbnailFigure 14. Drug response follows two different schedules but same drug intake (α = 0.8): (a–c) the state response, state space trajectory, and 3D state-space trajectory at τ = 20, respectively; (d–f) the state response, state space trajectory, and 3D state space trajectory at τ = 30, respectively. Other parameters are <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M110','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M110">View MathML</a>.

Extension and discussion

In previous sections, we have considered the drug-effect on one-gene and a 2-gene case. In this section, we will consider the drug-effect on a target gene in a more sophisticated GRN context.

3-gene network with multiple feedback loops

Suppose a 3-gene network is given by

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M111','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M111">View MathML</a>

(36)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M112','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M112">View MathML</a>

(37)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M113','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M113">View MathML</a>

(38)

where η1 is a perturbation parameter (from x3 to x1), β1, β2, β3 are activation factors, γ2, γ3 are degradation factors, and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M114','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M114">View MathML</a> are threshold values. <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M115','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M115">View MathML</a> is a drug-effect factor. We assume periodic drug intake and drug concentration level follows exponential decay during each period, i.e., <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M116','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M116">View MathML</a>, where kτi ≤ t ≤ (k+1)τi. A graphical model of the 3-gene network is given in Figure 15. There exist two positive feedback loops between x1 and x3.

thumbnailFigure 15. The 3-gene GRN model. The solid line is the real interaction between genes. The dashed line is derived to show the positive feedback loop for certain conditions.

When the target gene is in GRN context, not only its expression level is related to drug perturbation, but also depends on network contexts. Several interesting phenomena are observed through our simulations study:

1. Drug response is related to disease stage. Simulations are performed with different initial target gene expression level (x1(0)). Figure 16a–c shows the system responses with x1(0) = 20, which is not too high (corresponding to early disease state). As shown in Figure 16a, x1 expression level reduces to the range [7.7, 8.4] under periodic drug intake, while x2 and x3, the two other interactive genes settle at 1.0 and 4.0, respectively. The system reaches a new steady state (a semi-stable limit cycle, to be exact), with <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M117','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M117">View MathML</a>, <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M118','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M118">View MathML</a>, and <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M119','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M119">View MathML</a>, where x1 is well controlled. The trajectories of x1 vs. u and x1 versus x3 are given in Figure 16b,c, respectively. The semi-stable limit cycle is shown in Figure 16b.

thumbnailFigure 16. Simulation results for the 3-gene network under drug perturbation with different initial conditions: (a–c) with initial condition (x1(0) = 20), (d–f) with initial condition (x1(0) = 40). Other parameter settings are: <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M120','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M120">View MathML</a>.

System responses with x1(0) = 40 (corresponding to late disease state) are shown in Figure 16d–f for comparison. Although the other parameter settings are exactly the same, the drug will not repress the disease gene x1 (Figure 16d) owing to the interaction between the disease gene x1 and gene x3. When <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M121','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M121">View MathML</a>, Equation (36) becomes <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M122','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M122">View MathML</a>, and thus x3 is negative regulated by x1 and converge to <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M123','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M123">View MathML</a>. However, when initial condition <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M124','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M124">View MathML</a>, Equation (36) becomes <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M125','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M125">View MathML</a>, and thus x3 is positively regulated by x2 and its expression level will keep increasing. As a result, x1 will keep increasing as well, and a positive feedback loop is formed between x1 and x3. This is confirmed by the trajectories of x1 versus u and x1 versus x3 given in Figure 16e,f, respectively.

2. Under certain conditions, single drug perturbation may not be enough. A drug is usually designed to a specific target. In this example, the drug tries to provide negative feedback to the regulation of x1 (tries to repress x1); however, since the target gene is interactive (or, in a more general setting, pathways have crosstalk), only repressing the target gene (or blocking the signal of one pathway) may not prevent the target gene from expressing itself through interactions with other genes (or through inter-connected pathways). In our case, x1 is interactive with x3. To continue with previous simulation (results shown in Figure 16d–f, we try to increase the drug dosage tenfold from u(kτ) = 24 to u(kτ) = 240 with the same dosing period τ = 8 trying to bring down the expression level of x1. However, from system responses shown in Figure 17a–c, it is observed that the drug is not effective although the dosage is increased tenfold.

thumbnailFigure 17. Simulation results for the 3-gene network under drug perturbation (high dosage), parameter settings are the same with Figure 16d–f except dosage u() = 240.

One step further, not only we increase dosage to u(kτ) = 240, but also to increase the dosing frequency (dosing period is decreased from τ = 8 to τ = 2), systems responses are shown in Figure 18a–d, where Figure 18c shows the left part of the trajectory shown in Figure 18b. It can be observed that although the drug perturbation is very strong, and the drug concentration is always staying in domain D5, drug is still not effective.

thumbnailFigure 18. Simulation results for the 3-gene network under drug perturbation (high dosage and short interval): (a-d) state response, trajectory of x1 versus u, detailed initial trajectory of x1 versus u, and trajectory of x1 versusx3, respectively. Parameter settings are the same with Figure 16 except u(kτ) = 240 and τ = 2.

From the nonlinear dynamical system perspective, the equation <a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M126','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M126">View MathML</a> represents a semi-stable limit cycle. If the initial condition is from the inside of the limit cycle, then the system will converge to the limit cycle; however, if the initial condition is from the outside of the limit cycle, then the system will diverge from the limit cycle. Such simulation results demonstrate the heterogeneity of the drug’s responses due to the nonlinearities in complex systems, where multiple inputs affect each output and the underpinning structure may include parallel, redundant, and feedback loop processes, it is likely that some cases will not respond to a single drug perturbation no matter how strong it is. As a result, innovative perturbation methods, such as finding a better target or combinatorial therapy, are necessary.

Simulation of effects of different drugs and a drug combination on NFκB pathway

In this article, the models and examples are selected such that they are mathematically tractable and important insights can be obtained, and we can verify the theoretical results with simulation results. For large-scale networks and multiple drugs/drug targets, the proposed model is still applicable; however, analytical results may not be attainable even for this simplistic model. In that case, simulations can be carried out case-by-case. To illustrate this point of view, we carried out a simulation study of the NFκB pathway under two different drugs and each drug with different drug targets.

NFκB signaling regulates inflammation, cell proliferation, and apoptosis by increasing the expression of specific cellular genes in response to a variety of extracellular ligands. How to explore therapeutic strategies to prevent the prolonged activation of the NFκB pathway attracts lots of attention [32,33]. An ODE model of the NFκB pathway [34] is adopted and the two drugs under consideration are drug X (drug A in [35]) and FDA approved drug proteasome inhibitor Bortezomib (Velcade) [36]. The detailed simulation setup is available in the Appendix and the SIMULINK model is given in Figure 19. The specificity of some drugs to inhibit several of these components of the NFκB pathway is one of the concerns. For example, the proteasome which is responsible for the IκBα degradation has many other important functions. Thus, Bortezomib modulates a variety of cellular processes that may contribute to toxicity if the dosage is too high [36]. Hence, we design combination therapy to induce a better effect and at the same time to contain toxicity to a certain threshold.

thumbnailFigure 19. Simulation results for the NF−κB pathway under various drug perturbations with different drug administration. (a) Effect of continuous stimulus, no drug input. (b) Effect of drug X (0.2 μM). (c) Effect of drug X (1.0 μM). (d) Effect of Bortezomib (65% inhibition). (e) Effect of Bortezomib (95% inhibition). (f) Effect of combined Bortezomib (65% inhibition) and drug X (0.2 μM). The detailed parameters are available in the Appendix.

To achieve this, drug X is a protein kinase inhibitor, which competitively inhibit IKK with the binding kinetics the same as that of the natural reaction involving NFκB:IκBα and IKK [35]. While Bortezomib is a proteasome inhibitor that will inhibit IκBα degradation, its effect is adjusted through the parameter setting related to individual terms for IκBα and NFκB:IκBα molecules rescued from inhibition of IκBα degradation [35]. We first validate the results in [34,35]. In Figure 20a, we show that oscillatory behaviors occur for NFκB pathway with constant stimulus. Under this constant stimulus, it is observed in Figure 20b,c that only very high dose of drug X can effectively block NFκB nuclear translocation. Similar observation is obtained for Bortezomib from Figure 20d,e, where low drug dosage (65% inhibition) is not effective, while the side effects are unacceptable when the drug is effective (95% inhibition). All the above simulation results are consistent with those in [34,35]. In this article, we go a step further and consider the combination of these two drugs. It is interesting to see in Figure 20f that some combinations of the drugs with non-overlapping toxicities, e.g., combined Bortezomib (65% inhibition) and drug X (0.2 μM), might provide enormous benefit by keeping the level of nuclear NFκB low while having tolerable toxicities.

thumbnailFigure 20. SIMULINK model for the NF−κB pathway.

Conclusions and future work

This article provides systematic mathematical analysis and dynamical modeling of drug effect in the GRN context, where a drug functions as a control input to reduce the elevated target gene expression level. A hybrid systems model is proposed to study the dynamics of the underlying regulatory network under drug perturbation. Drug pharmacology information is incorporated into drug therapeutic response modeling to demonstrate the significant difference in drug effect for different dosing regimens. Considering the complicated nature of gene regulation, this study is a small step towards quantitative modeling of therapeutic effect. We have kept the examples mathematically tractable so that valuable insights and reasonable predictions can be obtained from theoretical analysis.

Compared to our previous work [22], where drug effect was only studied for a specific PK profile (drug concentration only has exponential decay stage) when the drug is targeted to a single gene, three major extensions are provided in this article: (i) we provide analytical results of drug effect under a very general PK profile, where three stages of drug concentration change (increase, equilibrium, and decrease) are considered; (ii) the proposed methodology is applied to interactive genes in a GRN context, with detailed analytical derivations for both one-gene and two-gene cases; and (iii) we perform extensive simulations for a more complicated GRN setting and explain several interesting observations due to multiple feedback loops and the existence of limit cycles.

It is expected that the theoretical framework proposed in this article, when correlated to real biological networks, can help improve drug development productivity and make drug discovery more systematic. During such process, cross disciplinary effort is indispensable. For example, application of such a framework will require experiments designed to elucidate model parameters, such as protein concentration levels and synthesis and degradation speeds. While some parameters may be relatively easy to obtain, others may be difficult to get based on current techniques and model simplification may be necessary; nonetheless, the basic hybrid systems model and the conclusions drawn from it, such as the nature of DERs and the role of limit cycles, will remain valid, only their particular forms being changed to represent experimental instantiation of the model.

Appendix

See Tables 1 and 2.

Table 1. Definition of variables

Table 2. Parameter values

ODE model of the NF − κB pathway

The ODE model of the NF − κB pathway is adopted from [34,35].

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M127','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M127">View MathML</a>

(39)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M128','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M128">View MathML</a>

(40)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M129','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M129">View MathML</a>

(41)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M130','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M130">View MathML</a>

(42)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M131','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M131">View MathML</a>

(43)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M132','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M132">View MathML</a>

(44)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M133','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M133">View MathML</a>

(45)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M134','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M134">View MathML</a>

(46)

<a onClick="popup('http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M135','MathML',630,470);return false;" target="_blank" href="http://bsb.eurasipjournals.com/content/2012/1/19/mathml/M135">View MathML</a>

(47)

Competing interests

The authors declare that they have no competing interests.

Acknowledgements

This study was supported in part by the National Cancer Institute (2 R25CA090301-06) and the National Science Foundation (NSF-1238918).

References

  1. S Huang, Rational drug discovery: what can we learn from regulatory networks. Drug Discov. Today 7, 163–169 (2002). Publisher Full Text OpenURL

  2. E Butcher, E Berg, E Kunkel, System biology in drug discovery. Nat. Biotechnol 22, 1253–1259 (2004). PubMed Abstract | Publisher Full Text OpenURL

  3. N Kumar, B Hendriks, K Janes, D Graaf, D Lauffenburger, Applying computational modeling to drug discovery and development. Drug Discov. Today 11, 806–811 (2006). PubMed Abstract | Publisher Full Text OpenURL

  4. J Chen, E Dougherty, S Demir, C Friedman, C Li, S Wong, Grand challenges for multimodal bio-medical systems. IEEE Circuits Syst. Mag 5(2), 46–52 (2005)

  5. H Derendorf, B Meibohm, Modeling of pharmacokinetic/pharmacodynamic (PK/PD) relationships: concepts and perspectives. Pharm. Res 16(2), 176–185 (1999). PubMed Abstract | Publisher Full Text OpenURL

  6. N Ting, Introduction and new drug development process. Dose Finding in Drug Development, ed. by Ting N (New York: Springer, 2006)

  7. S Undevia, G Gomez-Abuin, M Ratain, Pharmacokinetic variability of anticancer agents. Nat. Rev. Cancer 5, 447–458 (2005). PubMed Abstract | Publisher Full Text OpenURL

  8. N Shah, C Kasap, C Weier, M Balbas, J Nicoll, E Bleickardt, C Nicaise, C Sawyers, Transient potent BCR-ABL inhibition is sufficient to commit chronic myeloid leukemia cells irreversibly to apoptosis. Cancer cell 14(6), 485–493 (2008). PubMed Abstract | Publisher Full Text OpenURL

  9. D Amin, N Sergina, D Ahuja, M McMahon, J Blair, D Wang, B Hann, K Koch, K Shokat, M Moasser, Resiliency and vulnerability in the HER2-HER3 tumorigenic driver. Sci. Transl. Med 2(16), 1–9 (2010)

  10. P Rajasethupathy, S Vayttaden, U Bhalla, Systems modeling: a pathway to drug discovery. Curr. Opin. Chem. Biol 9, 400–406 (2005). PubMed Abstract | Publisher Full Text OpenURL

  11. E Butcher, Can cell systems biology rescue drug discovery? Nat. Rev. Drug Discov 4, 461–467 (2005). PubMed Abstract | Publisher Full Text OpenURL

  12. X Li, L Qian, J Hua, M Bittner, E Dougherty, Assessing the efficacy of molecularly targeted agents on cell line-based platforms by using system identification. BMC Genomics 13(6), S11 (2012). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  13. E Dougherty, M Brun, J Trent, M Bittner, Conditioning-based modeling of contextual genomic regulation. IEEE/ACM Trans. Comput. Biol. Bioinf 6, 310–320 (2009)

  14. A Kringstein, F Rossi, A Hofmann, H Blau, Graded transcriptional response to different concentrations of a single transactivator. Proc. Natl Acad. Sci. USA 95, 13670–13675 (1998). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  15. M Santillan, M Mackey, Quantitative approaches to the study of bistability in the lac operon of Escherichia coli. J. R. Soc. Interface 5(Suppl 1), S29–S39 (2008). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  16. B Muller-Hill, The Lac Operon: A Short History of a Genetic Paradigm (New York: Walter de Gruyter, 1996)

  17. A Halasz, V Kumar, M Imielinski, C Belta, O Sokolsky, S Pathak, H Rubin, Analysis of lactose metabolism in E. coli using reachability analysis of hybrid systems. IET Syst. Biol 1, 120–148 (2007). PubMed Abstract | Publisher Full Text OpenURL

  18. R Ghosh, A Tiwari, C Tomlin, Automated symtolic reachability analysis; with application to delta-notch signaling automata. Hybrid Syst.: Comput. Control LNCS 2623, 233–248 (2003)

  19. R Ghosh, C Tomlin, Symbolic reachable set computation of piecewise affine hybrid automata and its application to biological modelling: delta-notch protein signalling. IET Syst. Biol 1, 170–183 (2004). Publisher Full Text OpenURL

  20. S Drulhe, G Ferrari-Trecate, H de Jone, A Viari, Reconstruction of switching thresholds in piecewise-affine models of genetic regulatory networks. Hybrid Syst.: Comput. Control LNCS 3927, 184–199 (2006)

  21. G Batt, D Ropers, H de Jone, J Geiselmann, M Page, D Schneider, Qualitative analysis and verification of hybrid models of genetic regulatory networks: nutritional stress response in Escherichia coli. Hybrid Syst.: Comput. Control LNCS 3414, 134–150 (2005)

  22. X Li, L Qian, M Bittner, E Dougherty, Characterization of drug efficacy regions based on dosage and frequency schedules. IEEE Trans. Biomed. Eng 58(3), 488–498 (2011). PubMed Abstract | Publisher Full Text OpenURL

  23. H de Jong, Modeling and simulation of genetic regulatory systems: a literature review. J. Comput. Biol 9, 67–103 (2002). PubMed Abstract | Publisher Full Text OpenURL

  24. L Glass, S Kauffman, The logical analysis of continuous non-linear biochemical control networks. J. Theor. Biol 39, 103–129 (1973). PubMed Abstract | Publisher Full Text OpenURL

  25. J Prez-Urizar, V Granados-Soto, FJ Flores-Murrieta, G Castada-Hernndez, Pharmacokinetic–pharmacodynamic modeling: why? Arch. Med. Res 31(6), 539–545 (2000). PubMed Abstract | Publisher Full Text OpenURL

  26. AV Hill, The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves. J. Physiol 40, iv–vii (1910)

  27. NH Holford, LB Sheiner, Understanding the dose–effect relationship: clinical application of pharmacokinetic–pharmacodynamic models. Clin. Pharmacokinet 6(6), 429–453 (1981). PubMed Abstract | Publisher Full Text OpenURL

  28. DR Waud, RB Parker, Pharmacological estimation of drug–receptor dissociation constants. Statistical evaluation. II. Competitive antagonists. J. Pharmacol. Exp. Ther 177, 13–24 (1971). PubMed Abstract | Publisher Full Text OpenURL

  29. H Kuh, S Jang, M Wientjes, J Au, Computational model of intracellular pharmacokinetics of paclitaxel. J. Pharmacol. Exp. Ther 293(3), 761–770 (2000). PubMed Abstract | Publisher Full Text OpenURL

  30. AR Tzafriri, ER Edelman, Endosomal receptor kinetics determine the stability of intracellular growth factor signalling complexes. Biochem. J 402, 537–549 (2007). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  31. J Kramer, J Sagartz, D Morris, The application of discovery toxicology and pathology towards the design of safer pharmaceutical lead candidates. Nat. Rev. Drug Discov 6, 636–649 (2007). PubMed Abstract | Publisher Full Text OpenURL

  32. Y Yamamoto, R Gaynor, Therapeutic potential of inhibition of the NF-κB pathway in the treatment of inflammation and cancer. J. Clin. Invest 107(2), 135–142 (2001). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  33. V Baud, M Karin, Is NF-κB a good target for cancer therapy? Hopes and pitfalls. Nat. Rev. Drug Discov 8, 33–40 (2009). PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  34. A Hoffmann, A Levchenko, ML Scott, D Baltimore, The IkBNF-kB signaling module: temporal control and selective gene activation. Science 298(8), 1241–1245 (2002). PubMed Abstract | Publisher Full Text OpenURL

  35. M Sung, R Simon, In silico simulation of inhibitor drug effects on nuclear factor-kB pathway dynamics. Mol. Pharmacol 66, 70–75 (2004). PubMed Abstract | Publisher Full Text OpenURL

  36. R Orlowski, D Kuhn, Proteasome inhibitors in cancer therapy: lessons from the first decade. Clin. Cancer Res 14(6), 1649–1657 (2008). PubMed Abstract | Publisher Full Text OpenURL