Abstract
Drug discovery today is a complex, expensive, and timeconsuming 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 regimensIntroduction
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 proteinbased 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 halflife, can achieve deep clinical remission in CML patients by achieving transient potent BCR–ABL inhibition, while traditionally approved tyrosine kinase inhibitors usually have prolonged halflives 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,1012]. 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 statespace 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 2gene 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 moreorless 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 switchlike 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 DeltaNotch 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 righthand 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.
with
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 drugeffect terms. Equation (2) extends Equation (1) by including such terms. While the structure is intuitively reasonable and somewhat general, the actual details of the drugeffect 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 2gene 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 selfregulation. We also assume that a drug targets
x_{1} by reducing its expression level and providing a negative feedback term
where β_{1}, β_{2} are synthesis factors, γ_{2} is a degradation factor, and
Figure 1. State trajectory schematic of 2gene example without drug intake.
Figure 2. State trajectory schematic of 2gene example with drug intake.
We assume periodic drug intake and the drug concentration level in the effectsite
follows exponential decay during each period τ_{i}, i.e.,
Figure 3. The state response under periodic drug intake.
Figure 4. The statespace trajectory under periodic drug intake. Parameter setting of 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 effectsite [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 drugeffect 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 twodimensional PK and PD relationship, the proposed approach enables a description of a threedimensional dose–concentration–effect relationship. Specifically, PK and PD are linked through γ^{u} by a statespace 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
Figure 7. The PD model: concentrationresponse curve used in this study.
Figure 8. Sigmoidal E_{max }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 drugeffect coefficient
where
Periodic drug intake: PK model
Drug concentration at the effectsite is critical for its pharmacological effect. Currently, plasma drug concentrations are markers that serve as surrogates for drug concentration at the effectsite for beneficial and adverse effects; however, markers not grounded on a sound theoretical foundation and therapeutic mechanismbased 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 effectsite and, as a result, effects diminish.
Based on the study by Kuh et al. [29], a general model for drug concentrationtime 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.
Figure 9. A general drug concentrationtime 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. Closedform 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 selfregulation capacity and the dynamical equation of the expression level x_{1} is given by
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 to
where
where for kτ_{i} ≤ t ≤ (k+1)τ_{i} and i = 1,2,… denoting the index of different dosing regimens,
Statespace analysis
The statespace 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
Figure 10. The state trajectory schematic (target gene expression versus drug concentration)
under PK profile (Figure 9) assuming
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
State trajectory analysis
We analyze the drug effect considering the scenario shown in Figure 10, where
• 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.
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)  D_{1} (kτ_{i} ≤ t ≤ t_{1}):
• Stage (b)  D_{3} (t_{1}≤ t ≤ t_{2}):
• Stage (c)  D_{5} (t_{2 } ≤ t ≤ t_{3} = kτ_{i}+p_{1}):
• 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 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
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:
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 D_{3} and D_{1}. Figure 12d shows that the statespace 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 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
Analysis of 2gene networks
We extend the theoretical analysis to 2gene 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
Figure 13. A 2gene network with positive feedback loop and drug input.
where β_{1}, β_{2} are synthesis factors, γ_{1} and γ_{2} are degradation factors, and
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
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
Now with the baseline analysis of the secondorder system, we provide detailed state
trajectory analysis by taking into account the practical form of PK/PD (
State trajectory analysis
We analyze the drugeffect 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}.
• Stage (b)  D_{3}(t_{1} ≤ t ≤ t_{2}):
• where a, b, d are defined as before. When incorporating the practical form of
• Stage (c)  D_{5} (t_{2} ≤ t ≤ t_{3} = kτ_{i} + p_{1}): The set of equations are the same as in Stage (b) except that
• Stage (d)  D_{5} (t_{3} ≤ t ≤ t_{4} = kτ_{i} + p_{1} + p_{2}): The 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
• Stage (e)  D_{5} (t_{4} ≤ t ≤ t_{5}): 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
• Stage (f)  D_{3} (t_{5} ≤ t ≤ t_{6}): 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
• Stage (g)  D_{1} (t_{6} ≤ t ≤ (k+1)τ_{i}): The solution of x_{1} 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 2gene case, no explicit closedform 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 closedform 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 (
Figure 14. Drug response follows two different schedules but same drug intake (α = 0.8): (a–c)
the state response, state space trajectory, and 3D statespace trajectory at τ = 20,
respectively; (d–f) the state response, state space trajectory, and 3D state space
trajectory at τ = 30, respectively. Other parameters are
Extension and discussion
In previous sections, we have considered the drugeffect on onegene and a 2gene case. In this section, we will consider the drugeffect on a target gene in a more sophisticated GRN context.
3gene network with multiple feedback loops
Suppose a 3gene network is given by
where η_{1} is a perturbation parameter (from x_{3} to x_{1}), β_{1}, β_{2}, β_{3} are activation factors, γ_{2}, γ_{3} are degradation factors, and
Figure 15. The 3gene 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 (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 semistable limit cycle, to be exact), with
Figure 16. Simulation results for the 3gene network under drug perturbation with different initial
conditions: (a–c) with initial condition (x_{1}(0) = 20), (d–f) with initial condition (x_{1}(0) = 40). Other parameter settings are:
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. 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 interconnected pathways). In our case, x_{1} is interactive with x_{3}. 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 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.
Figure 17. Simulation results for the 3gene network under drug perturbation (high dosage), parameter settings are the same with Figure 16d–f except dosage u(kτ) = 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 D_{5}, drug is still not effective.
Figure 18. Simulation results for the 3gene network under drug perturbation (high dosage and short interval): (ad) state response, trajectory of x_{1 }versus u, detailed initial trajectory of x_{1 }versus u, and trajectory of x_{1 }versusx_{3}, respectively. Parameter settings are the same with Figure 16 except u(kτ) = 240 and τ = 2.
From the nonlinear dynamical system perspective, the equation
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 largescale 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 casebycase. 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.
Figure 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 nonoverlapping 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.
Figure 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 onegene and twogene 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
ODE model of the NF − κB pathway
The ODE model of the NF − κB pathway is adopted from [34,35].
Competing interests
The authors declare that they have no competing interests.
Acknowledgements
This study was supported in part by the National Cancer Institute (2 R25CA09030106) and the National Science Foundation (NSF1238918).
References

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

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

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

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

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

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

S Undevia, G GomezAbuin, M Ratain, Pharmacokinetic variability of anticancer agents. Nat. Rev. Cancer 5, 447–458 (2005). PubMed Abstract  Publisher Full Text

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

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 HER2HER3 tumorigenic driver. Sci. Transl. Med 2(16), 1–9 (2010)

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

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

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

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

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

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

B MullerHill, The Lac Operon: A Short History of a Genetic Paradigm (New York: Walter de Gruyter, 1996)

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

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

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

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

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)

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

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

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

J PrezUrizar, V GranadosSoto, FJ FloresMurrieta, G CastadaHernndez, Pharmacokinetic–pharmacodynamic modeling: why? Arch. Med. Res 31(6), 539–545 (2000). PubMed Abstract  Publisher Full Text

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

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

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

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

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

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

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

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

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

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

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