Dynamical modeling of drug effect using hybrid systems

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.


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 *Correspondence: xiangfangli@ieee.org 1 Department of Electrical and Computer Engineering, Texas A&M University, College Station, TX 77843, USA Full list of author information is available at the end of the article 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 http://bsb.eurasipjournals.com/content/2012/1 /19 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][11][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 moreor-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 where x i ≥ 0 corresponds to the concentrations of proteins encoded by genes in the network and can be interpreted as the gene expression level. f i (.) 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. γ i x i 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: 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. K i ≥ 1 is an integer representing the number of activation/synthesis terms.
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 T i j , where t+ i and t− i are indices of the threshold values, 0 ≤ θ t+ i j ∈ T i j , and 0 ≤ θ t− i j ∈ T i j . 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  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 x 1 , x 2 that repress each other, and x 1 is a disease gene which loses its self-regulation. We also assume that a drug targets x 1 by reducing its expression level and providing a negative feedback term −γ u 1 x 1 . The resulting 2-gene network is given byẋ where β 1 , β 2 are synthesis factors, γ 2 is a degradation factor, and θ 2 1 , θ 1 2 are threshold values. γ u 1 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 x 1 increases unbounded, while  with proper drug input, β 1 /γ u 1 < θ 2 1 , the system converges to a new steady state, (β 1 /γ u 1 , β 2 /γ 2 ). We assume periodic drug intake and the drug concentration level in the effect-site follows exponential decay during each period τ i , i.e., u i (t) = ζ i e −λ d (t−kτ i ) , 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 x 1 vs. the drug concentration level u is given in Figure 4. A comparison of trajectory of the gene expression level x 1 vs. x 2 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 x 1 . 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, x 1 = β 1 /γ u 1 and x 2 = β 2 /γ 2 = 1, while x 2 → 0 and x 1 → ∞ without drug input.

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 (doseconcentration). 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 E max 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 E max model (see Figure 8). The E max model has the general form E = E max C m EC m 50 +C m , where E max is the maximum effect, C is the concentration, EC 50 is the concentration necessary to produce 50% of E max , and m represents a sigmoidity factor or steepness of the curve.
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 http://bsb.eurasipjournals.com/content/2012/1/19

Figure 7
The PD model: concentration-response curve used in this study.
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 γ u 1 (the drug target is x 1 ) 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 where q u 1 is the ratio between the drug-effect factor γ u 1 and the effective drug concentration u − θ u 1 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 θ u 1 and its effect saturates when its concentration level exceeds an upper threshold θ u 1 . Note that the sigmoidal E max model can be well approximated by the proposed PD model. By taking the derivative of E with respect to C and evaluating it at EC 50 , we obtain the slope as q u 1 = mE max 4EC 50 . The upper and lower bounds should satisfy An example of the sigmoidal E max model when m = 4 and our proposed PD model are plotted together in Figure 8, where the proposed model closely resembles the sigmoidal E max 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 p 1 , p 2 , and p 3 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. http://bsb.eurasipjournals.com/content/2012/1/19

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 x 1 is given bẏ There is a steady state x 1 = β 1 /γ 1 in such a system, however, if the synthesis rate is much bigger than the selfdegradation 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 tȯ where γ u 1 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 bẏ where for kτ i ≤ t ≤ (k + 1)τ i and i = 1, 2, . . . denoting the index of different dosing regimens, 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 D 1 , D 3 , D 5 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 θ u 1 . 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 θ u 1 and the lower threshold θ u 1 . The third scenario is the case that ζ i < θ u 1 and the drug is not effective. 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 D 5 and D 3 , 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 D 1 , the expression level will rise (to the right), since the drug concentration is lower than θ u 1 . For the drug to be effective, the reduction of the expression level in D 5 and D 3 has to be larger than the increase of the expression level in D 1 . In summary, we should have so that after each treatment the expression level x 1 will decrease.

State trajectory analysis
We analyze the drug effect considering the scenario shown in Figure 10, where ζ i > θ u 1 . The same methodology can be applied to a simpler scenario where θ u 1 < ζ i < θ u 1 . 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 • t 1 : the traveling time from the initial state to the boundary between D 3 and D 1 . • t 2 : from initial state to boundary between D 5 and D 3 . • t 3 : from initial state to the end of stage c, t 3 = kτ i +p 1 .
• t 4 : time at which the drug concentration starts to decrease, t 4 = kτ i + p 1 + p 2 . • t 5 : from the initial state to the end of stage e. • t 6 : from the initial state to the end of stage f.
where i is the index for different dosing regimens, the corresponding equations and solutions for each stage are given by: • Stage (d) -D 5 (t 3 ≤ t ≤ t 4 = kτ i + p 1 + p 2 ): • Stage (e) -D 5 (t 4 ≤ t ≤ t 5 ): • Stage (f ) -D 3 (t 5 ≤ t ≤ t 6 ): • Stage (g) -D 1 (t 6 ≤ t ≤ (k + 1)τ i ): We can deduce the necessary and sufficient condition for the effectiveness of the drug by expressing the inequality x 1 ((k + 1)τ ) < x 1 (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 x 1 = x 1 (kτ i ), the http://bsb.eurasipjournals.com/content/2012/1/19 equations governing the state trajectory from time kτ i to time (k + 1)τ i are given by 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 x 1 ((k + 1)τ ) < x 1 (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 α = Dose1 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.
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 12af, 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 D 3 and D 1 . 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 D 1 is getting longer (where drug concentration is below θ

Figure 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: q u 1 = 0.1, 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.

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 x 1 is the target gene, there exists a positive feedback loop between x 1 and another gene x 2 , and that a drug targets x 1 by reducing its expression level and providing a negative feedback term −γ u 1 x 1 . The resulting 2-gene network is shown in Figure 13 and is given byẋ where β 1 , β 2 are synthesis factors, γ 1 and γ 2 are degradation factors, and γ u 1 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: The solution of this equation is given by where λ 1 and λ 2 are the two eigenvalues of Equation 28 and k 1 and k 2 are parameters depending on the initial conditions. Letting a = γ 1 + γ u 1 + γ 2 , b = (γ 1 + γ u 1 )γ 2 − η 1 η 2 and d = β 1 γ 2 + β 2 η 1 , the two eigenvalues are given by It is easy to verify that a 2 − 4b = (γ 1 + γ u 1 − γ 2 ) 2 + 4η 1 η 2 > 0. Since a 2 − 4b > 0, we conclude that λ 1 = λ 2 and both eigenvalues are real. Furthermore, one of the eigenvalues, λ 2 , is always negative since λ 2 = −a− √ a 2 −4b 2 < 0. The sign of λ 1 will be determined by the sign of b: In other words, The above equation has an important biological interpretation: when the degradation of x 1 due to the strength of the drug is faster than the increase of x 1 due to the positive feedback loop, both eigenvalues are negative, the system is stable and x 1 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 x 1 will increase exponentially.
Given initial condition x 1 (t 0 ) andẋ 1 (t 0 ), then for the case λ 1 = λ 2 , we have 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 (γ u 1 ) 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) -D 1 (kτ i ≤ t ≤ t 1 ): The solution of x 1 (t) is given by Equation (29), with k 1 and k 2 given by Equations (32) and (33), and t 0 = kτ i . : where a, b, d are defined as before. When incorporating the practical form of γ u 1 = q u 1 (u − θ u 1 ) and u = e λ a (t−kτ i ) − 1, the above second-order ODE has no closed-form solution. In this case, the solution can be obtained numerically.
The set of equations are the same as in Stage (b ) except that Since γ u 1 does not depend on u = e λ a (t−kτ i ) − 1 explicitly, x 1 has a closed-form solution given by Equation (29).
solution of x 1 is the same as that in Stage (c) except the start and end times, and the equation of u, which is The solution of x 1 is the same as in Stage (c) except the start and end times, and the equation of u, which now is The solution of x 1 is the same as in Stage (b ) except the start and end times, and the equation of u, which now is is the same as in Stage (a) except the start and end times, and the equation of u, which now is We can deduce the necessary and sufficient condition for the effectiveness of the drug by expressing the inequality x 1 ((k + 1)τ ) < x 1 (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 http://bsb.eurasipjournals.com/content/2012/1/19  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 (θ u 1 ), 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.

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.
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 (x 1 (0)). Figure 16a-c shows the system responses with x 1 (0) = 20, which is not too high (corresponding to early disease state). As shown in Figure 16a, x 1 expression level reduces to the range [ 7.7, 8.4] under periodic drug intake, while x 2 and x 3 , 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 x s 1 = x s 3 = β 3 γ 3 , where x 1 is well controlled. The trajectories of x 1 vs. u and x 1 versus x 3 are given in Figure 16b,c, respectively. The semi-stable limit cycle is shown in Figure 16b. System responses with x 1 (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 x 1 (Figure 16d) owing to the interaction between the disease gene x 1 and gene x 3 . When 2 ) − γ 3 x 3 , and thus x 3 is negative regulated by x 1 and converge to x s 3 = β 3 γ 3 . However, when initial condition x 1 (0) = 40 > θ 3 1 = 21, Equation (36) becomesẋ 3 = β 3 s − (x 2 , θ 3 2 ), and thus x 3 is positively regulated by x 2 and its expression level will keep increasing. As a result, x 1 will keep increasing as well, and a positive feedback loop is formed between x 1 and x 3 . This is confirmed by the trajectories of x 1 versus u and x 1 versus x 3 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 x 1 (tries to repress x 1 ); 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, x 1 is interactive with x 3 . To continue with previous simulation (results shown in Figure 16d trying to bring down the expression level of x 1 . However, from system responses shown in Figure 17a-c, it is observed that the drug is not effective although the dosage is increased tenfold. 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 D 5 , drug is still not effective.
From the nonlinear dynamical system perspective, the equation represents a semistable 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 http://bsb.eurasipjournals.com/content 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. 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.

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. 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 http://bsb.eurasipjournals.com/content/2012/1/19 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.

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