Drug discovery today is a complex, expensive, and time-consuming process with high attrition rate. A more systematic approach is needed to combine innovative approaches in order to lead to more effective and efficient drug development. This article provides systematic mathematical analysis and dynamical modeling of drug effect under gene regulatory network contexts. A hybrid systems model, which merges together discrete and continuous dynamics into a single dynamical model, is proposed to study dynamics of the underlying regulatory network under drug perturbations. The major goal is to understand how the system changes when perturbed by drugs and give suggestions for better therapeutic interventions. A realistic periodic drug intake scenario is considered, drug pharmacokinetics and pharmacodynamics information being taken into account in the proposed hybrid systems model. Simulations are performed using MATLAB/SIMULINK to corroborate the analytical results.
Keywords:Drug effect; Hybrid systems; PK/PD; Gene regulatory network (GRN); Dosing regimens
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 . This approach is possible because the phenotypic behavior of the cell reflects the dynamics of the gene and protein-based regulatory network. When it comes to drug therapeutics and disease modeling, the major goal is to understand how the system changes when perturbed and how to modify the system to achieve a desired outcome. To understand and exploit the complicated mapping between genome and phenome, especially in the context of drug discovery, it is critical to evaluate the regulatory interactions between the genes and proteins that form the gene regulatory network (GRN). To date, the hope of the rapid translation of “genes to drugs” has foundered on the reality that disease biology is complex and drug development must be driven by insights into biological responses . 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 . 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 . 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 . 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.  demonstrate that BCR–ABL inhibitor dasatinib, which has greater potency and a short half-life, can achieve deep clinical remission in CML patients by achieving transient potent BCR–ABL inhibition, while traditionally approved tyrosine kinase inhibitors usually have prolonged half-lives that result in continuous target inhibition. A similar study of whether short pulses of higher dose or persistent dosing with lower doses has the most favorable outcomes has been carried out by Amin et al.  in the setup of inactivation of HER2–HER3 signaling. Finding an optimal dosing regimen based on the dynamics of biological systems and relevant PK/PD information is critically important.
System modeling is emerging as a valuable tool in therapeutics to address these challenges [3,10-12]. The process begins with building a quantitative model of a biological system. Consequences of particular perturbations, such as optimal dosing regimens, optimal drug targets, or combinational therapy, can be simulated in time courses using such models. In this study, we propose a hybrid systems model for GRNs and incorporate a drug’s PK and PD information by using a state-space approach. We first study drug effect assuming the drug target to be a gene or protein in the proposed drug perturbation model using dynamical system theory, considering the case of periodic drug intake and analytically deriving the conditions for the drug to be effective. We extend the analysis to the 2-gene case and then to the case of a network with multiple coupled genes and positive feedback loops. Simulations are performed using MATLAB/SIMULINK to supplement our analytical results.
While discrete modeling leaves out many details, continuous modeling includes so many details that computational demands preclude their applications to many larger systems. Hybrid systems, which aim to merge ideas from both continuous and discrete modeling into one paradigm, are appealing for GRN modeling under drug perturbations because biological systems are naturally nonlinear, have highly varied regulatory requirements, and possess a wide range of control strategies for meeting their needs. While some simple, local, feedback control methods can provide sufficient regulation of many more-or-less continuous cellular processes, the regulation of discontinuous processes possessing the character of computational decision making requires more elaborate regulatory methods . In particular, some genes display regulation in a thresholded switch-like manner .
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 . Other biological hybrid systems analyzed in similar ways include the Delta-Notch decision process [18,19], GRNs of carbon starvation , and nutritional stress response  in Escherichia coli. As far as we know, the only hybrid systems modeling concerning treatment or drug effects is contained in our earlier work .
where xi≥0 corresponds to the concentrations of proteins encoded by genes in the network and can be interpreted as the gene expression level. fi(.) is a general nonlinear function and represents the rate of synthesis. It can be approximated by a sigmoidal function or a unit step function, and unit step function is used in this article. γixi is the rate of degradation. To use hybrid systems and incorporate drug effect, we propose the following model for a GRN of N genes under drug perturbation:
where the last two terms on the right-hand side of Equation (2) correspond to drug perturbation u. β > 0 and γ > 0 are synthesis and degradation rates, respectively. Ki ≥ 1 is an integer representing the number of activation/synthesis terms. 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 drug-effect factor, which is closely related to the drug pharmacology model discussed in the following section.
It should be kept in mind that the focus of this article is studying the effect of dosing, in particular, dosing regimens, on the expression of genes involved in a pathology by using hybrid systems theory. Whereas the simpler Equation (1) is widely accepted, it does not contain drug-effect terms. Equation (2) extends Equation (1) by including such terms. While the structure is intuitively reasonable and somewhat general, the actual details of the drug-effect terms are unknown. Finding the specific form of Equation (2) for a specific disease is a system identification problem, which is quite distinct from the analysis problem addressed in this article. We are addressing optimization of treatment intervention, given the system. The details of our analysis might change when the details of Equation (2) are clarified, but we expect that the hybrid systems approach taken in the article will go through with appropriate modifications in the mathematical details.
We consider a 2-gene example to illustrate the feasibility of using hybrid systems for modeling drug effect. Specifically, we assume that there are two interactive genes x1,x2 that repress each other, and x1 is a disease gene which loses its self-regulation. We also assume that a drug targets x1 by reducing its expression level and providing a negative feedback term . The resulting 2-gene network is given by
where β1, β2 are synthesis factors, γ2 is a degradation factor, and are threshold values. is a drug-effect factor. Using dynamical systems theory, the state-trajectory schematic diagrams of this 2-gene network without and with drug input are obtained and plotted in Figures 1 and 2, respectively. It is observed that without drug input, the gene expression level of x1 increases unbounded, while with proper drug input, , the system converges to a new steady state, .
Figure 1. State trajectory schematic of 2-gene example without drug intake.
Figure 2. State trajectory schematic of 2-gene example with drug intake.
We assume periodic drug intake and the drug concentration level in the effect-site follows exponential decay during each period τi, i.e., , where kτi ≤ t ≤ (k + 1)τi and λd is the degradation factor. The response of gene expression levels of the two genes under periodic drug intake is shown in Figure 3. The state-space trajectory of gene expression level of x1 vs. the drug concentration level u is given in Figure 4. A comparison of trajectory of the gene expression level x1 vs. x2 with and without drug are provided in Figures 5 and 6, respectively. It is observed that the drug is quite effective in bringing down the expression level of x1. The simulation study matches the theoretical analysis, as in Figures 1 and 2, that with proper drug intervention the system will converge to a new steady state, and x2 = β2 / γ2 = 1, while x2 → 0 and x1 → ∞ without drug input.
Figure 3. The state response under periodic drug intake.
The basis of clinical pharmacology is the fact that the intensities of many pharmacological effects are functions of the amount of drug in the body and, more specifically, the concentration of drug at the effect-site . Historically, PK and PD were considered as separate disciplines; however, the information provided by these disciplines is limited if regarded in isolation . A drug-effect factor γu is included in our proposed model (Equation 2), which is related to drug’s PD characteristic (concentration–response) and its PK information (dose–concentration). In order to describe the time course of drug effect in response to different dosing regimens, the integrated PK/PD model is indispensable because it builds the bridge between these two classical disciplines of pharmacology . Following each dosing regimen, instead of a two-dimensional PK and PD relationship, the proposed approach enables a description of a three-dimensional dose–concentration–effect relationship. Specifically, PK and PD are linked through γu by a state-space approach to facilitate the description and prediction of the time course of drug effects resulting from different drug administration regimens.
Drug concentration–response curve: PD model
In general, the magnitude of a pharmacological effect increases monotonically with increased dose, eventually reaching a plateau level where further increase in dose has little additional effect . The classic and most commonly used concentration–response model is the Hill equation , also known as the sigmoidal Emax model  or logistic model . The relationship between the concentration of the drug and its effect is most often nonlinear. In this study, we use hybrid systems to approximate PD curves. A common method is to replace certain slowly changing variables by their piecewise linear approximations (see Figure 7). For example, the PD model used in our study can approximate the popular sigmoidal Emax model (see Figure 8). The Emax model has the general form , where Emax is the maximum effect, C is the concentration, EC50 is the concentration necessary to produce 50% of Emax, and m represents a sigmoidity factor or steepness of the curve.
Figure 7. The PD model: concentration-response curve used in this study.
Figure 8. Sigmoidal Emax model (m=4), and approximation by our PD model.
We assume a threshold of concentration below which the drug candidate is ineffective, the minimum effective dose (MinED), and another threshold value, called maximum effective dose (MaxED), above which there is no clinically significant increase in pharmacological effect in this study. As an example, we use a linear curve to approximate the concentration–response curve between MinED and MaxED. It is assumed that the drug-effect coefficient (the drug target is x1) is related to the concentration u through a sigmoid function and can be approximated by the curve shown in Figure 7. The corresponding relationship can be expressed as
where is the ratio between the drug-effect 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 Emax model can be well approximated by the proposed PD model. By taking the derivative of E with respect to C and evaluating it at EC50, we obtain the slope as . The upper and lower bounds should satisfy . An example of the sigmoidal Emax model when m=4 and our proposed PD model are plotted together in Figure 8, where the proposed model closely resembles the sigmoidal Emax model.
Periodic drug intake: PK model
Drug concentration at the effect-site is critical for its pharmacological effect. Currently, plasma drug concentrations are markers that serve as surrogates for drug concentration at the effect-site for beneficial and adverse effects; however, markers not grounded on a sound theoretical foundation and therapeutic mechanism-based intervention can limit the usefulness of PK/PD modeling to drug development. For example, it has been demonstrated that the intracellular PK of a drug is quite different from plasma drug concentration [29,30]. As observed in the study by Kuh et al. , the intracellular concentration of a drug will exponentially increase as the drug is absorbed after each drug intake. The drug concentration may change very slowly (in our model, we approximate that as a flat curve) when the intracellular and extracellular drug concentration approach equilibrium. In time, drug concentration will exponentially decrease as the rate at which it is eliminated is more than the rate at which it enters the effect-site and, as a result, effects diminish.
Based on the study by Kuh et al. , a general model for drug concentration-time profile is given in Figure 9. Drug concentration is plotted on a logarithmic scale against time following each periodic drug intake. λa denotes the exponential increase quotient; λd is the exponential decrease quotient; τ is the interval between each drug intake; and p1, p2, and p3 denote the time stayed in the increase, equilibrium, and decrease stage, respectively. Different drugs work in different ways and the proposed model is general enough to cover various cases. Drug concentration may increase very quickly and, as a result, the increase stage may be neglected, or the equilibrium stage may be very short and can be ignored for simplicity. By adjusting the parameters in the proposed model, specific drug characteristics can be represented. In the case when the proposed model cannot approximate a drug’s PK profile, extensive simulations can be performed based on the drug’s actual PK profile. In this article, we consider a periodic drug intake scenario. Specifically, we are interested in investigating and comparing the following two potential scenarios: large dose with a longer interval versus small dose with a shorter interval.
Figure 9. A general drug concentration-time profile.
Mathematical analysis of drug effect
In this section, we study the time course of drug effect for different dosage and schedule arrangements where the drug is designed to repress a “target gene”. The case with a special PK profile (drug concentration only has exponential decay) was analytically studied in our previous work . In this study, we extend the analysis considering a general PK profile given in Figure 9 and PD model given in Figure 8. Closed-form analytical solution is provided and simulations are performed to validate theoretical analysis. In later sections, we show that the same methodology can be applied to interactive genes, where not only will the drug affect the gene expression level, but the target gene is also coupled with other genes.
It is assumed that the disease gene has lost part of its self-regulation capacity and the dynamical equation of the expression level x1 is given by
There is a steady state x1 = β1 / γ1 in such a system, however, if the synthesis rate is much bigger than the self-degradation rate, β1 ≫ γ1, then the gene expression level will be too high. A drug is used as a control input to repress the target gene expression level. The corresponding dynamical equation after drug intake is changed to
where is the drug-effect factor defined in the previous section. After incorporating a drug’s PK/PD (Figures 7 and 9) into our proposed hybrid system model Equation (8), considering the scenario that the patient is taking the drug periodically, the resulting model is given by
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.
The state-space and a sample trajectory schematic of the state (target gene expression and drug concentration level) under periodic drug intake are shown in Figure 10. As is common in hybrid systems, there are both continuous quantitative changes and discrete transitions in our proposed model. The entire state space may be divided into different domains according to the value of discrete state. As shown in Figure 10, there are five domains in the state space, with D1,D3,D5 not being transient domains. The figure shows the case when the drug is effective and the drug dosage is large enough that ζi is higher than the upper threshold . 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.
When the state transits in each period under periodic drug intake, it may pass through different domains (depending on the changes of drug concentration along time). During the transit time through domains D5 and D3, the gene expression level is pushed lower (to the left), while the driving strength will depend on the drug’s PD characteristic. During the transit time through D1, the expression level will rise (to the right), since the drug concentration is lower than . For the drug to be effective, the reduction of the expression level in D5 and D3 has to be larger than the increase of the expression level in D1. In summary, we should have x1((k+1)τ) < x1(kτ), so that after each treatment the expression level x1 will decrease.
State trajectory analysis
We analyze the drug effect considering the scenario shown in Figure 10, where . 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 stage-by-stage. The time notations used in the derivation are given by
• t1: the traveling time from the initial state to the boundary between D3 and D1.
• t2: from initial state to boundary between D5 and D3.
• t3: from initial state to the end of stage c, t3 = kτi+p1.
• t4: time at which the drug concentration starts to decrease, t4 = kτi+p1+p2.
• t5: from the initial state to the end of stage e.
• t6: from the initial state to the end of stage f.
For kτi ≤ t ≤ (k+1)τi, where i is the index for different dosing regimens, the corresponding equations and solutions for each stage are given by:
• Stage (a) - D1 (kτi ≤ t ≤ t1):
• Stage (b) - D3 (t1≤ t ≤ t2):
• Stage (c) - D5 (t2 ≤ t ≤ t3 = kτi+p1):
• Stage (d) - D5 (t3 ≤ t ≤ t4 = kτi + p1 + p2):
• Stage (e) - D5 (t4 ≤ t ≤ t5):
• Stage (f) - D3 (t5 ≤ t ≤ t6):
• Stage (g) - D1 (t6 ≤ t ≤ (k+1)τi):
We can deduce the necessary and sufficient condition for the effectiveness of the drug by expressing the inequality x1((k+1)τ) < x1(kτ) in terms of dosing period τ and unit dose, assuming the dosage is proportional to the maximum drug concentration ζi reached after taking the drug. When the initial conditions are x1 = x1(kτi), the equations governing the state trajectory from time kτi to time (k+1)τi are given by
For the drug to be effective, we need the disease gene expression level to decrease following each period of drug intake. Hence, we can express x1((k+1)τ) < x1(kτ) in terms of dosage and frequency schedule and derive the region where the drug is effective using the above listed equations.
Results and analysis
Based on the theoretical analysis in previous two sections, we demonstrate that the drug efficacy depends on total drug intake, different dosages, and frequencies. The density of drug intake is defined as . 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 . 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 D3 and D1. Figure 12d shows that the state-space trajectory settles in a small limit cycle and disease gene expression level settles at 11.5 at the end of each period of treatment. On the other hand, the big infrequent drug intake case results in a big limit cycle as in Figure 12f. Although the dosage is high, the long period between dosages means that the period stayed in D1 is getting longer (where drug concentration is below , hence not effective), and disease gene expression level settles at 12.1 at the end of each period of treatment. As a comparison, the drug effect for the case with intermediate dosage and frequency shown in Figure 12b,e is superior to the other two cases. Disease gene expression level settles at 8.8 at the end of each treatment period. If we check the curve in Figure 11 with α = 0.4, the intermediate dosage case with τ = 12 is located near the bottom of the “U” shape. Lastly, we observe that all state-space trajectories follow the state trajectory schematic in Figure 10, as predicted by the analytical results.
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 2-gene networks
We extend the theoretical analysis to 2-gene networks and show that the same framework applies to the modeling and analysis of drug effect in more complex gene networks. We assume that x1 is the target gene, there exists a positive feedback loop between x1 and another gene x2, and that a drug targets x1 by reducing its expression level and providing a negative feedback term . The resulting 2-gene network is shown in Figure 13 and is given by
Figure 13. A 2-gene network with positive feedback loop and drug input.
where β1, β2 are synthesis factors, γ1 and γ2 are degradation factors, and is a drug-effect factor. η1 > 0 and η2 > 0 are the parameters of the positive feedback loop between the two genes. The 2-gene network under drug perturbation model, Equation (27), can be rewritten as a second-order ODE:
The solution of this equation is given by
where λ1 and λ2 are the two eigenvalues of Equation 28 and k1 and k2 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 a2−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 x1 due to the strength of the drug is faster than the increase of x1 due to the positive feedback loop, both eigenvalues are negative, the system is stable and x1 will experience exponential decay; on the other hand, if the effect of the positive feedback loop is dominant, then one of the eigenvalues will be positive and x1 will increase exponentially.
Given initial condition x1(t0) and , then for the case λ1 ≠ λ2, we have
Now with the baseline analysis of the second-order system, we provide detailed state trajectory analysis by taking into account the practical form of PK/PD ( ) when the drug is taken periodically.
State trajectory analysis
We analyze the drug-effect following the same framework given in the subsection “State trajectory analysis” under the main section “Mathematical analysis of drug effect”. For kτi ≤ t ≤ (k+1)τi, i = 1,2,…, the corresponding equations and solutions for each stages are given as follows:
• Stage (a) - D1 (kτi ≤ t ≤ t1):
• The solution of x1(t) is given by Equation (29), with k1 and k2 given by Equations (32) and (33), and t0 = kτi.
• Stage (b) - D3(t1 ≤ t ≤ t2):
• where a, b, d are defined as before. When incorporating the practical form of and , the above second-order ODE has no closed-form solution. In this case, the solution can be obtained numerically.
• Stage (c) - D5 (t2 ≤ t ≤ t3 = kτi + p1): The set of equations are the same as in Stage (b) except that . Since does not depend on explicitly, x1 has a closed-form solution given by Equation (29).
• Stage (d) - D5 (t3 ≤ t ≤ t4 = kτi + p1 + p2): The solution of x1 is the same as that in Stage (c) except the start and end times, and the equation of u, which is in this stage.
• Stage (e) - D5 (t4 ≤ t ≤ t5): The solution of x1 is the same as in Stage (c) except the start and end times, and the equation of u, which now is .
• Stage (f) - D3 (t5 ≤ t ≤ t6): The solution of x1 is the same as in Stage (b) except the start and end times, and the equation of u, which now is .
• Stage (g) - D1 (t6 ≤ t ≤ (k+1)τi): The solution of x1 is the same as in Stage (a) except the start and end times, and the equation of u, which now is .
We can deduce the necessary and sufficient condition for the effectiveness of the drug by expressing the inequality x1((k+1)τ) < x1(kτ) in terms of dosing period τ and unit dose. In the 2-gene case, no explicit closed-form expression can be deduced for the solutions in stages (b) and (f) and numerical methods have to be applied. However, through such analysis, it is observed that the same methodology for analyzing drug effect can be extended to GRNs with multiple interactive genes, although the mathematics involved will become more complicated and sometimes numerical methods must be applied when there is no closed-form solution.
Simulation results and analysis
When drug input is not present, the disease gene expression will grow unbounded owing to the positive feedback loop between the two genes. Here, we study response of the disease gene expression to drug input and compare two different schedules for τ = 20 and τ = 30, keeping α = 0.8 fixed. The response and state trajectories in 2D and 3D are given in Figure 14a–f, with the first and second rows corresponding to τ = 20 and τ = 30, respectively. We observe that both cases have periodic responses, but the disease gene expression is much better controlled when τ = 20. This is because the drug concentration is high enough in both cases compared to the threshold ( ), while the decay of the drug concentration is shorter in the case when τ = 20. In Figure 14c,f, the 3D state (disease gene expression) trajectories show that the trajectory settles to an inner circle when τ = 20, whereas the trajectory settles to an outer circle when τ=30. Similar observations apply to Figure 14b,e. Note the scale of x-axis of Figure 14e is 20 times bigger than that of Figure 14b.
Figure 14. Drug response follows two different schedules but same drug intake (α = 0.8): (a–c) the state response, state space trajectory, and 3D state-space trajectory at τ = 20, respectively; (d–f) the state response, state space trajectory, and 3D state space trajectory at τ = 30, respectively. Other parameters are .
Extension and discussion
In previous sections, we have considered the drug-effect on one-gene and a 2-gene case. In this section, we will consider the drug-effect on a target gene in a more sophisticated GRN context.
3-gene network with multiple feedback loops
Suppose a 3-gene network is given by
where η1 is a perturbation parameter (from x3 to x1), β1, β2, β3 are activation factors, γ2, γ3 are degradation factors, and are threshold values. is a drug-effect 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 3-gene network is given in Figure 15. There exist two positive feedback loops between x1 and x3.
Figure 15. The 3-gene GRN model. The solid line is the real interaction between genes. The dashed line is derived to show the positive feedback loop for certain conditions.
When the target gene is in GRN context, not only its expression level is related to drug perturbation, but also depends on network contexts. Several interesting phenomena are observed through our simulations study:
1. Drug response is related to disease stage. Simulations are performed with different initial target gene expression level (x1(0)). Figure 16a–c shows the system responses with x1(0) = 20, which is not too high (corresponding to early disease state). As shown in Figure 16a, x1 expression level reduces to the range [7.7, 8.4] under periodic drug intake, while x2 and x3, the two other interactive genes settle at 1.0 and 4.0, respectively. The system reaches a new steady state (a semi-stable limit cycle, to be exact), with , , and , where x1 is well controlled. The trajectories of x1 vs. u and x1 versus x3 are given in Figure 16b,c, respectively. The semi-stable limit cycle is shown in Figure 16b.
Figure 16. Simulation results for the 3-gene network under drug perturbation with different initial conditions: (a–c) with initial condition (x1(0) = 20), (d–f) with initial condition (x1(0) = 40). Other parameter settings are: .
System responses with x1(0) = 40 (corresponding to late disease state) are shown in Figure 16d–f for comparison. Although the other parameter settings are exactly the same, the drug will not repress the disease gene x1 (Figure 16d) owing to the interaction between the disease gene x1 and gene x3. When , Equation (36) becomes , and thus x3 is negative regulated by x1 and converge to . However, when initial condition , Equation (36) becomes , and thus x3 is positively regulated by x2 and its expression level will keep increasing. As a result, x1 will keep increasing as well, and a positive feedback loop is formed between x1 and x3. This is confirmed by the trajectories of x1 versus u and x1 versus x3 given in Figure 16e,f, respectively.
2. Under certain conditions, single drug perturbation may not be enough. A drug is usually designed to a specific target. In this example, the drug tries to provide negative feedback to the regulation of x1 (tries to repress x1); however, since the target gene is interactive (or, in a more general setting, pathways have crosstalk), only repressing the target gene (or blocking the signal of one pathway) may not prevent the target gene from expressing itself through interactions with other genes (or through inter-connected pathways). In our case, x1 is interactive with x3. To continue with previous simulation (results shown in Figure 16d–f, we try to increase the drug dosage tenfold from u(kτ) = 24 to u(kτ) = 240 with the same dosing period τ = 8 trying to bring down the expression level of x1. However, from system responses shown in Figure 17a–c, it is observed that the drug is not effective although the dosage is increased tenfold.
One step further, not only we increase dosage to u(kτ) = 240, but also to increase the dosing frequency (dosing period is decreased from τ = 8 to τ = 2), systems responses are shown in Figure 18a–d, where Figure 18c shows the left part of the trajectory shown in Figure 18b. It can be observed that although the drug perturbation is very strong, and the drug concentration is always staying in domain D5, drug is still not effective.
Figure 18. Simulation results for the 3-gene network under drug perturbation (high dosage and short interval): (a-d) state response, trajectory of x1 versus u, detailed initial trajectory of x1 versus u, and trajectory of x1 versusx3, respectively. Parameter settings are the same with Figure 16 except u(kτ) = 240 and τ = 2.
From the nonlinear dynamical system perspective, the equation represents a semi-stable limit cycle. If the initial condition is from the inside of the limit cycle, then the system will converge to the limit cycle; however, if the initial condition is from the outside of the limit cycle, then the system will diverge from the limit cycle. Such simulation results demonstrate the heterogeneity of the drug’s responses due to the nonlinearities in complex systems, where multiple inputs affect each output and the underpinning structure may include parallel, redundant, and feedback loop processes, it is likely that some cases will not respond to a single drug perturbation no matter how strong it is. As a result, innovative perturbation methods, such as finding a better target or combinatorial therapy, are necessary.
Simulation of effects of different drugs and a drug combination on NF−κB pathway
In this article, the models and examples are selected such that they are mathematically tractable and important insights can be obtained, and we can verify the theoretical results with simulation results. For large-scale networks and multiple drugs/drug targets, the proposed model is still applicable; however, analytical results may not be attainable even for this simplistic model. In that case, simulations can be carried out case-by-case. To illustrate this point of view, we carried out a simulation study of the NF−κB pathway under two different drugs and each drug with different drug targets.
NF−κB signaling regulates inflammation, cell proliferation, and apoptosis by increasing the expression of specific cellular genes in response to a variety of extracellular ligands. How to explore therapeutic strategies to prevent the prolonged activation of the NF−κB pathway attracts lots of attention [32,33]. An ODE model of the NF−κB pathway  is adopted and the two drugs under consideration are drug X (drug A in ) and FDA approved drug proteasome inhibitor Bortezomib (Velcade) . 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 . 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 . 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 . We first validate the results in [34,35]. In Figure 20a, we show that oscillatory behaviors occur for NF−κB pathway with constant stimulus. Under this constant stimulus, it is observed in Figure 20b,c that only very high dose of drug X can effectively block NF−κB nuclear translocation. Similar observation is obtained for Bortezomib from Figure 20d,e, where low drug dosage (65% inhibition) is not effective, while the side effects are unacceptable when the drug is effective (95% inhibition). All the above simulation results are consistent with those in [34,35]. In this article, we go a step further and consider the combination of these two drugs. It is interesting to see in Figure 20f that some combinations of the drugs with non-overlapping toxicities, e.g., combined Bortezomib (65% inhibition) and drug X (0.2 μM), might provide enormous benefit by keeping the level of nuclear NF−κB low while having tolerable toxicities.
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 , where drug effect was only studied for a specific PK profile (drug concentration only has exponential decay stage) when the drug is targeted to a single gene, three major extensions are provided in this article: (i) we provide analytical results of drug effect under a very general PK profile, where three stages of drug concentration change (increase, equilibrium, and decrease) are considered; (ii) the proposed methodology is applied to interactive genes in a GRN context, with detailed analytical derivations for both one-gene and two-gene cases; and (iii) we perform extensive simulations for a more complicated GRN setting and explain several interesting observations due to multiple feedback loops and the existence of limit cycles.
It is expected that the theoretical framework proposed in this article, when correlated to real biological networks, can help improve drug development productivity and make drug discovery more systematic. During such process, cross disciplinary effort is indispensable. For example, application of such a framework will require experiments designed to elucidate model parameters, such as protein concentration levels and synthesis and degradation speeds. While some parameters may be relatively easy to obtain, others may be difficult to get based on current techniques and model simplification may be necessary; nonetheless, the basic hybrid systems model and the conclusions drawn from it, such as the nature of DERs and the role of limit cycles, will remain valid, only their particular forms being changed to represent experimental instantiation of the model.
The authors declare that they have no competing interests.
This study was supported in part by the National Cancer Institute (2 R25CA090301-06) and the National Science Foundation (NSF-1238918).
S Huang, Rational drug discovery: what can we learn from regulatory networks. Drug Discov. Today 7, 163–169 (2002). Publisher Full Text
N Shah, C Kasap, C Weier, M Balbas, J Nicoll, E Bleickardt, C Nicaise, C Sawyers, Transient potent BCR-ABL inhibition is sufficient to commit chronic myeloid leukemia cells irreversibly to apoptosis. Cancer cell 14(6), 485–493 (2008). PubMed Abstract | Publisher Full Text
X Li, L Qian, J Hua, M Bittner, E Dougherty, Assessing the efficacy of molecularly targeted agents on cell line-based platforms by using system identification. BMC Genomics 13(6), S11 (2012). PubMed Abstract | Publisher Full Text | PubMed Central Full Text
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
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, C Tomlin, Symbolic reachable set computation of piecewise affine hybrid automata and its application to biological modelling: delta-notch protein signalling. IET Syst. Biol 1, 170–183 (2004). Publisher Full Text
S Drulhe, G Ferrari-Trecate, H de Jone, A Viari, Reconstruction of switching thresholds in piecewise-affine models of genetic regulatory networks. Hybrid Syst.: Comput. Control LNCS 3927, 184–199 (2006)
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
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
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