Time-Based Switching Control of Genetic Regulatory Networks: Toward Sequential Drug Intake for Cancer Therapy.

As cancer growth and development typically involves multiple genes and pathways, combination therapy has been touted as the standard of care in the treatment of cancer. However, drug toxicity becomes a major concern whenever a patient takes 2 or more drugs simultaneously at the maximum tolerable dosage. A potential solution would be administering the drugs in a sequential or alternating manner rather than concurrently. This study therefore examines the feasibility of such an approach from a switched system control perspective. Particularly, we study how genetic regulatory systems respond to sequential (switched) drug inputs using the time-based switching mechanism. The design of the time-driven drug switching function guarantees the stability of the genetic regulatory system and the repression of the diseased genes. Simulation results using proof-of-concept models and the proliferation and survival pathways with sequential drug inputs show the effectiveness of the proposed approach.

Introduction

The perturbation of genetic regulatory networks (GRNs) is usually tied to oncogenesis because multiple genetic mutations and pathway breaches are commonly implicated in the growth and development of cancer. Therefore, an important step toward an effective cancer therapy is to understand the effects of drugs on GRNs. Many cancer treatments employ combination therapies in which 2 or more drugs are administered in targeting the cancer cells and those drugs function in a synergistic manner to disrupt specific phases of the cell reproduction cycles.1 Such combinatorial targeted therapies have been touted as the standard of care in preventing genetic mutation and drug-related resistance.

Toxicity is a major concern when drug combinations are taken simultaneously at the maximum tolerated dosage (MTD)2 even though combination therapies typically take into account the toxicity of drug combinations and how effective they are against the cancer cells.1 Toxicity is generally a function of the reduction in the immune system performance of the patient, bodyweight loss, pain, and other side effects experienced by the patients. Such adverse events may range from grade 0 to grade 5 (grades 0, 1, 2, 3, 4, and 5 imply the absence of adverse events, presence of mild, moderate, severe, life-threatening, and death-related adverse events, respectively). The vital question is, “Could the patient take the drugs in a sequential manner, rather than simultaneously, such that the undesired biological signals are blocked and yet the toxicity is low?” In this study, we try to answer this question from a mathematical modeling perspective using switched systems control theory.

As a result of the increased toxicity associated with combination therapy, clinical trials and studies focusing on sequential drug intake have been conducted (or ongoing) by different research groups and one of the observations from those trials is that toxicity at MTD is lower in sequential drug intake as compared with concomitant or simultaneous administration of anticancer agents.3–8 Thus, sequential regimen potentially lowers toxicity and provides a way to optimally deliver single-drug treatment and may improve patient’s quality of life.9 A review of trials which compares combination treatment with sequential regimens is provided by Miles et al,9 and an additional review of clinical trials and preclinical evidence in support of each strategy is given by Felici et al.10

To reduce toxicity, delay resistance to the drug combinations,11,12 and double the progression free survival of the agents,13,14 the authors from the research groups11,13–16 conducted clinical trials on administering the drug treatments sequentially or on alternating days. For instance, metastatic renal cell carcinoma (mRCC) treatment is focusing on agents that block tumors and vascular growth pathways.13 Sunitinib is directed at blocking the vascular endothelial growth factor receptors, whereas temsirolimus inhibits the mammalian target of rapamycin (mTOR). Sunitinib and temsirolimus are agents approved by the Food and Drug Administration (FDA) for mRCC treatment. Toxicity increases whenever such agents are taken together. The hope is that such agents can be safely administered sequentially at full dosages.13

Similar studies1,8,17 are experimental works in which the drug combinations were alternated, administered sequentially, or given simultaneously (at intervals) to ensure tolerable toxicity levels. In this study, we examined a mathematical framework for sequential drug treatment with anticancer agents and provide analytical insights into how effective such sequential treatment regimens are with a focus on GRN with switched (sequential) drug inputs. Specifically, we study how genetic regulatory systems respond to sequential (switched) drug inputs by treating the problem as a switched system with stable or weakly stable subsystems which can be addressed using the time-dependent switching mechanism. The design of the switching function guarantees the stability of the genetic regulatory system and the repression of the diseased genes with sequential (switched) drug inputs. For the proposed time-based switching approach, we provide simulation studies with proof-of-concept GRNs having switched drugs perturbations and a practical case of the proliferation and survival pathways with sequential drug inputs to show the effectiveness of the approach.

The model and problem formulation are discussed in detail in section “Problem Formulation.” Section “Switching Design and Stability Analysis for Genetic Regulatory Systems With Sequential (Switched) Drug Input” presents the time-based switching strategy and the stability analysis for GRNs. Simulation results are provided in section “Simulation Results.” Section “Discussions” provides further discussions and section “Conclusions” concludes the article.

Problem Formulation

Genetic regulatory networks can be modeled with rate equations that express the differences between production and degradation rates.18–20 The corresponding ordinary differential equation (ODE) model is defined as follows: where x⩾0 represents the expression level of the ith gene. g () is a nonlinear function denoting the rate of synthesis. γ corresponds to the degradation rate. For diseased genes not repressed, γ is reduced to γ⩽γ, which implies a weak or almost negligible negative feedback. When the gene loses self-regulation, this corresponds to the absence of negative feedback and the system will be unstable and the diseased gene expression level will grow very high. Whenever drugs are used as the control input for repressing the expression level of the target gene, say x, it is assumed that the drug supplies the negative feedback term −γ where γ is the drug effect factor.

Because multiple mutated genes are involved in the growth and development of cancer, multiple drugs are normally used simultaneously to attack the cancer cells, but the increased toxicity becomes a challenge. We therefore propose a switched systems mathematical model for taking the multiple drugs sequentially or alternately to inhibit those diseased genes while reducing toxicity, as corroborated by the clinical trials and previous studies.3–8 To illustrate the proposed sequential drug intake paradigm, consider a GRN with 2 diseased genes (x,x): The ODEs for the normal genes are omitted above. We assumed that the drug effects of the 2 drugs will not overlap, ie, the effect of drug 1 decays considerably when drug 2 is taken. This ensures that the toxicity of the 2 drugs will not add together such that the patient can tolerate. It should be noted that the systems become a switched system.

In various ODE models for biological pathways, simple linear approximations of g and g are used. Then, the state-space model for the pathways including the treated genes, equations (2) and (3), can be expressed as follows: where A and B are matrices of appropriate dimensions corresponding to the case that only the drug for x is taken. It is assumed that the matrix A is invertible. η is the drug effect factor related to the pharmacology model of the drug. σ denotes the switching logic. The switching rule between subsystems depends on drug administration.

As a simple conceptual introduction to the sequential drug intake strategy, consider a gene regulatory network, as shown in Figure 1, with a focus on gene 1 and gene 2 only. For simplicity and mathematical tractability, we examine ODE model with no cross talk among the genes as follows (i = 1,2): where x is the diseased gene, β is the constant synthesis rate of the genes, α is the gene degradation rates, and is the drug effect factor for drug i. In the matrices A, B, and η, we have

Figure 1

A gene regulatory network. The focus is on gene 1 and gene 2 with the effects of drugs 1 and 2 and neglecting the details in the dotted box for mathematical tractability.

For 2 diseased genes (x1,x2) with , the state-space trajectory is shown in Figure 2. The plot shows the different drug administration methods. For instance, with no drug intake, the expression levels of the diseased genes are too high signifying overexpression of the diseased genes. With only 1 drug intake, the other gene is overexpressed and its expression level is therefore too high. The sequential drug intake involves switching between the 2 drugs at intervals determined by the proposed time-driven drug switching logic. The state-space trajectory depicts that both the simultaneous drug intake and sequential drug intake significantly reduced the expression levels of the diseased genes. Although they are both effective against the diseased genes, toxicity is a concern with the simultaneous drug intake.3–8 Therefore, this study is focused on time-based switching control and stability analysis for the sequential drug intake approach due to the reduced toxicity associated with such drug intake methods.

Figure 2

State-space trajectory for the 2 diseased genes (x1, x2)with different drug intake methods where the drug effect factors are and for drugs 1 and 2, respectively.

Switching Design and Stability Analysis for Genetic Regulatory Systems With Sequential (Switched) Drug Input

This section aims to design the switching of the drug inputs so that the diseased genes are repressed. We refer to this as global asymptotic stability. Switched systems are made up of difference or differential equations with an associated rule that defines the switching strategy between them.21–24 Switched systems stability has been researched extensively, for instance, with stable subsystems,25–28 unstable subsystems,29–31 and a mixture of both.26,29,32 The switching strategy used in this study is referred to as the time-based switching approach in which all the subsystems or modes are stable or weakly stable based on the pathway models.33

Comment 1

Stability of the subsystems in the case of gene regulatory networks without (or with) drug input is based on the assumption that, for each subsystem, the genes may be overexpressed but the gene expression levels do not grow out of bounds as time t → ∞.

Time-based switching for GRNs with stable subsystems

For GRNs with sequential (switched) drug inputs in which all the subsystems are asymptotically stable, we modify the approach in the study by Geromel et al25 to analyze the switching design and stability of such GRNs. Geromel et al25 studied the stability of continuous time switched systems whose subsystems are stable by determining the minimum dwell time (DT) required for stability through a class of quadratic Lyapunov functions. In this case, the Lyapunov functions are not required to uniformly decrease at every instant of switching as a condition of stability. We denote piecewise Lyapunov functions as V(x). At each instant of switching, to bound the Lyapunov function increment, it is required that V (x)⩽µ V(x) where µ > 1 and i, j are subsystem’s indices before and after switching.

For switched systems with stable subsystems, DT constraints are based on the concept that at the switching instants, the likely increment of the Lyapunov function is compensated for by the decrease in the Lyapunov function within the DT. There is also the relaxed condition on the Lyapunov function that at each switching instant t, the sequence V(x(t)) for k = 0,…,∞, uniformly converges to 0.

Comment 2

The Lyapunov functions are employed in deriving the conditions that ensure a particular drug administration interval or the drug DT is effective in driving the overall system to the desired state, ie, driving the gene expression levels to the desired equilibrium point.

Consider the switched closed-loop genetic regulatory system,

The state is denoted by x ∈ R and w(t) denotes outward disturbances. σ(t) denotes the time-dependent switching logic that is based on whether the drug of interest is present or not and therefore chooses the appropriate subsystem’s sequences from the available N expressed as {A,C},i ∈ I[1, N]. We assume stable subsystems.

The following standard notations are used in this article: Set of real m × n matrices is R, S× denotes real, symmetric, n×n matrix and stands for positive definite matrices. I represents identity matrices of appropriate dimensions. The transpose of a matrix or vector is denoted as (’). For integers k1, k2, with k1

Stability analysis for time-based switching

The multiple quadratic Lyapunov functions are as follows: where and i is the active subsystem’s index. Each subsystem is associated with its own Lyapunov function.

Definition 1

The switching logic σ is defined to have a DT τ if t+1 − t⩾τ, ∀k where t, t+1 represent the successive switching instances.

Comment 3

The DT corresponds to the least time interval between 2 consecutive drug intakes which ensure that the gene expression levels eventually decay to the desired equilibrium point with the sequential drug administration.

The goal of the proposed drug switching strategy is to determine the minimum DT T* > 0 that guarantees the asymptotic stability of the equilibrium point of the gene regulatory system in equation (7). In other words, asymptotic stability is guaranteed if σ (t) is not changed for periods of time t ⩾ T*. The following proposition is modified from Geromel et al25 which provided the theorem that characterized an upper bound for T* as a possible solution to the problem.

Proposition 1

Assuming that for certain T > 0, there exists positive definite matrices , i ∈ I[1, N] such that25

Then, according to the DT switching approach σ(t) with t+1 − t ⩾T, the switched system (equation (7)) is globally asymptotically stable.

The proof is given in Appendix 1. An upper bound for the minimum DT T* is obtained from the optimal solution of the optimization problem25:

Comment 4

Equations (9) to (11) describe the conditions on the model parameters of the gene regulatory system and the least drug intake interval (or DT) that ensure the genes are not overexpressed as time t → ∞. This translates to the conditions under which the gene expression levels are regulated with a given sequential drug intake schedule.

Comment 5

An assumption behind the presented rationale of the mathematical framework is that a time factor (sequentially or alternating strategy) is crucial to allow the body to react properly to toxicity. This makes sense; a patient usually needs a few days to recover from cancer drugs. But the effects are often cumulative (for brain, heart, liver, kidney, skin, etc) depending on the drugs.

Simulation Results

We first examine simulations of simple GRNs with and without cross talk among the genes and with sequential drug perturbations using MATLAB/Simulink and the results obtained in section “Switching Design and Stability Analysis for Genetic Regulatory Systems With Sequential (Switched) Drug Input.” We then simulate the proliferation and survival pathways with sequential drug inputs using a sequential intake of lapatinib and temsirolimus which are FDA-approved drugs.

Case of 2 diseased genes (x1,x2) with sequential drug intake and absence of cross talk between the genes

We revisit the simple model in equation (5) with i = 1,2. We consider 2 diseased genes (x1,x2) with no cross talk between the genes for mathematical tractability. The following parameter values are used: β1 = β2 = 0.2, α1 = 0.2, α2 = 0.15, , and . Using the multiple quadratic Lyapunov functions v(x) = x, σ =1,2 and the analysis condition proposed in section “Time-based switching for GRNs with stable subsystems,” we have

Figure 3(A) to (D) depicts the switching of the drug input, dynamics of x1 and x2 with sequential drug intake, simultaneous drug intake, and no drug intake, respectively. Figure 3(E) is the state-space trajectory for x1 and x2 with different schedules of the 2 drugs including monotherapeutic cases where only 1 drug is administered. Figure 3(E) also shows the switching points (A, B, and C) of the drugs and the effects on the dynamics of x1 and x2 as we switch from one drug to the other in a sequential manner. It is observed that with no drug intake (Figure 3(D)), the diseased genes are overexpressed as indicated by the high expression levels compared with the other cases. With the sequential drug intake, Figure 3(B) depicts that the genetic regulatory system will be stabilized as the expression levels of the 2 diseased genes decay to the same equilibrium point as that of the simultaneous drug intake.

Figure 3

(A) The switching signal according to the analysis in section “Time-based switching for GRNs with stable subsystems.” (B) State dynamics of x1 and x2 with sequential drug intake according to the analysis in section “Time-based switching for GRNs with stable subsystems.” (C) State dynamics of x1 and x2 with simultaneous drug intake. (D) State dynamics of x1 and x2 with no drug intake. (E) State-space trajectory of the gene expression levels x1 versus x2 with different drug administrations. GRNs indicate genetic regulatory networks.

Case of 2 diseased genes (x1,x2) with sequential drug intake and with cross talk among the genes

With cross talk between the diseased genes, we consider a simple switched GRN where subsystem 1 (only drug 1 is taken) is defined as follows:

Similarly, when drug 2 is taken, subsystem 2 is defined as follows: where β1 and β2 denote the synthesis terms, γ1 and γ2 represent the degradation terms, and denotes the drug effect term. ψ1 > 0 and ψ2 > 0 are parameters of the positive feedback between the genes.

In state-space model form, we have where Setting

It can be observed that the subsystems have stable eigenvalues and are thus stable. We adopt the multiple quadratic Lyapunov functions , σ=1,2, and employing the analysis condition proposed in section “Time-based switching for GRNs with stable subsystems,” we have

From Figures 4 and 5, it is observed that with no drug intake, as shown in Figure 4(D), the diseased genes are overexpressed. With the sequential drug intake, Figure 4(A) and (B) depicts that the genetic regulatory system will be stabilized as the expression levels of the 2 diseased genes decay to the same equilibrium point as in the case of taking both drugs simultaneously. These observations are similar to the case when there is not any cross talk among the genes. However, unlike in Figure 3, it is also observed in Figure 4 that the value of this equilibrium point does not correspond to the ones that only 1 drug is taken. This is due to the positive feedback, and thus, the equilibrium points are different. In addition, unlike in Figure 3, due to the positive feedback, the switching is damped and is not very easy to observe in Figure 4.

Figure 4

(A) The switching signal according to the analysis in section “Time-based switching for GRNs with stable subsystems.” (B) State dynamics of x1 and x2 with sequential drug intake according to the analysis in section “Time-based switching for GRNs with stable subsystems.” (C) State dynamics of x1 and x2 with simultaneous drug intake. (D) State dynamics of x1 and x2 with no drug intake. (E) State dynamics of x1 and x2 with intake of drug 1 only. (F) State dynamics of x1 and x2 with intake of drug 2 only. The parameter values are ψ1 = 0.4, ψ2 = 0.4, γ1 = 0.4, γ2 = 1, , and β1 = β2 = 0.2. GRNs indicate genetic regulatory networks.

Figure 5

State-space trajectory of the gene expression levels x1 versus x2 with different drug administration schedules.

Case of the proliferation and survival pathway with sequential drug intake

The PI3K/AKT/mTOR pathway is a prototypical survival pathway which is essentially activated in several cancer types. The pathway is activated by various mechanisms some of which include mutating or amplifying the PI3K, losing the tumor suppressor (phosphatase and tensin homolog) functions, activating the growth factor receptors, amplifying or mutating AKT, exposures to carcinogens, and so on. After activation, signaling via AKT may propagate to different arrays of substrates, which includes mTOR, a vital regulator of protein translations. This pathway serves as an appealing drug target for cancer therapeutics as it functions as point of convergence for several growth stimuli. It is also responsible for regulating cellular processes contributing to cancer initiation and maintenance through its downstream substrates. In addition, activating the AKT/mTOR pathway aids resistances to several cancer treatment types, and it constitutes one of the poor prognostic factors for various cancer types.35

However, the RAS/RAF/MEK/ERK pathway is commonly associated with cell proliferations, prevention of apoptosis, and drug resistances. This pathway is employed by mitogens and growth factors in transmitting signals from the receptors for the regulation of gene expressions and prevention of programmed cell death or apoptosis. Some of this pathway components (eg, B-RAF, RAS) undergo mutation or overexpression in human cancer (eg, prostate and breast cancers).36

The aforementioned pathways are referred to as the proliferation and the survival pathways,33,34,37 which biologists presently understand, for instance, the Kegg collections of pathways (http://www.genome.jp/kegg/pathway.html) as well as National Institutes of Health BioCarta pathways collections (http://cgap.nci.nih.gov/Pathways/BioCarta_Pathways), are depicted in Figure 6 and Table 1. Table 1 shows the pathway dynamics for the proteins and complexes together with input from drugs lapatinib and temsirolimus.33

Figure 6

The proliferation and survival pathways with several drug inputs. Reproduced with permission from Li et al33 and Hua et al.34

Table 1

Notations and pathway dynamics (η1, η2, and η3 are the coefficients due to drug lapatinib acting on different proteins or complexes, and η6 is the coefficient due to drug temsirolimus acting on mammalian target of rapamycin33,37).

In this simulation, the drug Lapatinib is used to inhibit RAF/RAS pathway, and temsirolimus blocks mTOR. The 2 drugs were applied sequentially based on the switching logic derived analytically in section “Time-based switching for GRNs with stable subsystems.” Each drug was switched for roughly 50% of the simulation time duration. It can be observed from Figure 7(C) that the expression levels of the proteins and complexes given in Table 1 are stabilized and this translates to the regulation of the proliferation and survival pathway components and thereby regulating cell survival and proliferations. All the results here are in agreement with the analysis presented in “Time-based switching for GRNs with stable subsystems.” We simulated the pathway dynamics in Table 1 without drug intake and with both drugs taken simultaneously. Figure 7(A) corresponds to the case with no drug intake, and Figure 7(B) corresponds to the case with both drugs applied simultaneously. The expression levels of the genes in the drug intake case are significantly lower than the expression levels in the no drug case for both the sequential and simultaneous drug intake methods (see Figure 7(A) to (D)). For instance, the steady-state ratio of the expression levels of mTOR and FOS in the drug input case to the no drug case is approximately 40% (see Figure 8). This suggests that the drugs are effective when taken simultaneously or sequentially based on the reduction in the expression levels of the proteins and complexes. However, toxicity becomes a concern in the simultaneous intake case as already investigated in the previous studies.3–8 This suggests that the proposed sequential or switched drug inputs approach is a promising solution.

Figure 7

Proliferation and survival pathway simulations. (A) Pathway state dynamics according to Table 1 without drug intake. (B) Pathway state dynamics according to Table 1 with simultaneous drug intake. (C) Pathway dynamics with sequential drug intake based on analysis in section “Time-based switching for GRNs with stable subsystems,” (D) The Switching signal based on analysis in section “Time-based switching for GRNs with stable subsystems.” Parameter values are α1 =α2 = ⋯ = α13= α16= α19= 1, α14= α15= 0.6, and α17. = α18 = 0.8. GRNs indicate genetic regulatory networks.

Figure 8

The expression levels of mTOR and FOS with no drug and with sequential drug input.

Discussions

The dynamics of cancer cells can be described using continuous, discrete, or hybrid mathematical modeling framework. Continuous models are the seemingly appropriate candidates for modeling large-scale systems. They are able to capture large-scale behaviors of cancer growth and development at a reduced computational cost with the disadvantage of sacrificing the resolution of individual cells, especially when the cell properties vary over small spatiotemporal scales. Discrete models provide spatiotemporal representations of individual cells as well as cell-to-cell interactions. A main disadvantage of this modeling method is that the required computational cost is proportional to the number of cells being modeled. This limitation confines such modeling method to very small number of cells. Hybrid models combine the strengths of both continuous and discrete modeling approaches. They can also model randomness that may be inherent to the system being modeled. They are appealing for modeling GRNs under drug perturbation because biological systems are naturally nonlinear, have highly varied regulatory requirement, and possess a wide range of control strategies for meeting their needs. In the case of pathway dynamics, we employed the well-recognized ODE models which are widely used to model GRNs and pathways. An advantage of ODE models is that the mathematical analysis of the system structure is considerably simpler compared with other model types, and the solutions of ODE models are easier to simulate computationally with high efficiency. This means that such model structures can be made exceedingly complex before they become computationally infeasible.

Combination therapy is believed to be the standard of care in the prevention of gene mutations and drug resistance. Multiple drugs are used to attack cancer cells and these drugs synergistically disrupt distinct phases in the reproduction cycles of the cells.1 The benefits of this type of cancer therapy includes the improvement in patients’ compliance as a result of the reduction in number of administration, drug dosage decrement with associated reduction in toxic effects to healthy tissue, synergistic or additive effects of drug interactions, and delaying or overcoming multidrug resistances. These benefits are the motivations behind several combination therapy researches.1,17 They also motivated us to study combination therapy for GRNs from a sequential drug intake perspective. Because the perturbation of genes is usually responsible for oncogenesis, the theoretical analysis provided in this study presents a mathematical tool with a promising application to sequential cancer therapy which has been investigated as a way to reduce toxicity and improve the quality of life of the patient.

There are a number of challenges associated with the mode-ling framework proposed in this study. One is the assumption that the drug effects do not overlap which may not be the case in reality.1,17 Knowledge of the biological half-life of the drugs could help in scheduling the sequential drug administration such that the effect of one drug would have decayed significantly before the second drug is taken. Another way is to follow the methodology adopted in the study by Lewis13 where one of the drugs is taken for a certain period of time (eg, 4 weeks) followed by a long period of rest (2 weeks) and then the second drug is taken similarly to ensure that the drug effects do not overlap.

The state-space model adopted so far assumes that the regulatory systems are time-invariant. In reality, biological regulatory networks involve several different but interconnected phenomena, which are dynamic, nonlinear, stochastic, and may occur at different temporal and spatial scales.37 A way to partly tackle this challenge is by extending the current time-driven switching strategy to incorporate stochasticity in the analysis presented in section “Time-based switching for GRNs with stable subsystems.” The construction of a multiscale computational predictive model which is based on biological evidence and parameterized by biomedical data will undoubtedly be very useful. Advanced experimental technology and computational method should be mutually applied in a synergistic manner in addressing some of these challenges.

Conclusions

In this work, we have revisited combination therapy for cancer when toxicity is a concern. Instead of a patient taking all the drugs at the same time, the potential effects of sequential (switched) administration of drugs are examined. We believe that this is the first attempt to mathematically model and analyze such an approach to combination therapy for cancer. Specifically, the goal is to stabilize a GRN such that the oncogenes are repressed with sequential (switched) drug inputs. The analysis and design are based on multiple Lyapunov functions and linear matrix inequality. For linear approximations of the GRN, closed-form solution of the switching logic and thus the sequential (switched) administration of drugs are obtained. Simulation studies of proof-of-concept GRNs with switched/sequential drugs’ perturbations, and the proliferation and survival pathway with sequential drug inputs are provided to demonstrate the effectiveness of the method.

Experimental data have suggested that MTDs for combination therapy drugs are different from those of monotherapy in which the drugs are given sequentially or individually3–5,8 and toxicity of drug combinations are higher than single-drug treatments. The model presented in this work is based on some previous studies.3–8 It is shown that the sequential treatment schedule has lower drug-related toxicity when compared with the drugs taken simultaneously.13 Future studies will explore an extension of the current time-based switching approach to the stability analysis of tumor ODE models with sequential drug intake as well as the stability analysis for nonlinear gene regulatory systems. Such endeavor will be close to what obtains in wet lab experiments on tumor modeling with drug intake and it will also involve collaborations with experimental and clinical professionals.