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. and describe how other genes affect gene i. They are the functions of and , 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 , where t+_{i} and t−_{i} are indices of the threshold values, , and . Ψ_{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
with defined similarly. and may be set to 0 or 1, or different forms when appropriate threshold values are chosen. For example, and . and describe how the drug u affect gene i. and are the synthesis and degradation factors of the drug on gene i. and are used when the drug is activating or repressing certain genes, respectively. Since most drugs are used to repress genes, only is considered in the examples of this article. Note that γ^{u} is defined as a drugeffect 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 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 . The resulting 2gene network is given by
where β_{1}, β_{2} are synthesis factors, γ_{2} is a degradation factor, and are threshold values. is a drugeffect factor. Using dynamical systems theory, the statetrajectory schematic diagrams of this 2gene 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, , the system converges to a new steady state, .
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., , 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 statespace 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, 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 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 , 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.
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 (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 is the ratio between the drugeffect factor and the effective drug concentration 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 and its effect saturates when its concentration level exceeds an upper threshold . 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 . 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 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 is the drugeffect 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
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.
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 . 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 and the lower threshold . The third scenario is the case that and the drug is not effective.
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 . 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 x_{1}((k+1)τ) < x_{1}(kτ), 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 . The same methodology can be applied to a simpler scenario where . 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 stagebystage. 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.
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 . 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.
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 , 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 statespace trajectories follow the state trajectory schematic in Figure 10, as predicted by the analytical results.
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 . The resulting 2gene network is shown in Figure 13 and is given by
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 is a drugeffect factor. η_{1} > 0 and η_{2} > 0 are the parameters of the positive feedback loop between the two genes. The 2gene network under drug perturbation model, Equation (27), can be rewritten as a secondorder 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 , and d = β_{1}γ_{2} + β_{2}η_{1}, the two eigenvalues are given by . It is easy to verify that . 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 . 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 , then for the case λ_{1} ≠ λ_{2}, we have
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 () when the drug is taken periodically.
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 and , the above secondorder ODE has no closedform solution. In this case, the solution can be obtained numerically.
• 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 . Since does not depend on explicitly, x_{1} has a closedform solution given by Equation (29).
• 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 in this stage.
• 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 (), 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 xaxis of Figure 14e is 20 times bigger than that of Figure 14b.
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 are threshold values. is a drugeffect factor. We assume periodic drug intake and drug concentration level follows exponential decay during each period, i.e., , where kτ_{i} ≤ t ≤ (k+1)τ_{i}. A graphical model of the 3gene network is given in Figure 15. There exist two positive feedback loops between x_{1} and x_{3}.
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 , , and , 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 semistable limit cycle is shown in Figure 16b.
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 , Equation (36) becomes , and thus x_{3} is negative regulated by x_{1} and converge to . However, when initial condition , Equation (36) becomes , 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 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 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 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 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
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