Bifurcation analysis on a hybrid systems model of intermittent hormonal therapy for prostate cancer

Physica D 237 (2008) 2616–2627 www.elsevier.com/locate/physd

Bifurcation analysis on a hybrid systems model of intermittent hormonal therapy for prostate cancer Gouhei Tanaka a,∗ , Kunichika Tsumoto a,b , Shigeki Tsuji a,b , Kazuyuki Aihara a,b a Institute of Industrial Science, University of Tokyo, Tokyo 153-8505, Japan b ERATO Aihara Complexity Modelling Project, JST, Tokyo 151-0065, Japan

Received 16 August 2007; received in revised form 30 January 2008; accepted 22 March 2008 Available online 4 April 2008 Communicated by A. Mikhailov

Abstract Hybrid systems are widely used to model dynamical phenomena that are characterized by interplay between continuous dynamics and discrete events. An example of biomedical application is modeling of disease progression of prostate cancer under intermittent hormonal therapy, where continuous tumor dynamics is switched by interruption and reinstitution of medication. In the present paper, we study a hybrid systems model representing intermittent androgen suppression (IAS) therapy for advanced prostate cancer. Intermittent medication with switching between ontreatment and off-treatment periods is intended to possibly prevent a prostatic tumor from developing into a hormone-refractory state and is anticipated as a possible strategy for delaying or hopefully averting a cancer relapse which most patients undergo as a result of long-term hormonal suppression. Clinical efficacy of IAS therapy for prostate cancer is still under investigation but at least worth considering in terms of reduction of side effects and economic costs during off-treatment periods. In the model of IAS therapy, it depends on some clinically controllable parameters whether a relapse of prostate cancer occurs or not. Therefore, we examine nonlinear dynamics and bifurcation structure of the model by exploiting a numerical method to clarify bifurcation sets in the hybrid system. Our results suggest that adjustment of the normal androgen level in combination with appropriate medication scheduling could enhance the possibility of relapse prevention. Moreover, a two-dimensional piecewise-linear system reduced from the original model highlights the origin of nonlinear phenomena specific to the hybrid system. c 2008 Elsevier B.V. All rights reserved.

PACS: 05.45.-a; 87.19.Xx; 87.53.Tf Keywords: Hybrid systems; Piecewise linear systems; Grazing bifurcations; Limit cycle; Chaos; Prostate cancer; Intermittent androgen suppression

1. Introduction Following great progresses of life science and nonlinear science in recent years, a mathematical approach is becoming a more promising methodology for advanced studies in biomedical science. In particular, much attention has been paid to mathematical and computational modeling of cancer dynamics involving nonlinear biological interactions [1,2]. In the context of prostate cancer, for example, dynamical systems have been extensively used to describe tumor growth [3,4] and temporal variations of biomarkers [5–7] in efforts to predict ∗ Corresponding author. Tel.: +81 3 5452 6693; fax: +81 3 5452 6694.

E-mail address: [email protected] (G. Tanaka). c 2008 Elsevier B.V. All rights reserved. 0167-2789/$ - see front matter doi:10.1016/j.physd.2008.03.044

medical conditions and help appropriate diagnoses. The main concern of these studies is to understand the mechanism of treatment-resistant tumor growth, or a cancer relapse, after a remission period induced by hormone deprivation therapy. A possible strategy to delay or prevent the progression to hormone-resistance caused by prolonged hormone suppression is to incorporate treatment interruption, which is known as intermittent hormonal therapy repeating cycles of on-treatment and off-treatment periods [8–15]. In this paper, we deal with nonlinear dynamics and bifurcations in a mathematical model of intermittent hormonal therapy for prostate cancer. Prostate cancer is a disease characterized by uncontrolled growth of cancer cells within the prostate gland in males. The prostate gland is dependent on hormonal secretion by

G. Tanaka et al. / Physica D 237 (2008) 2616–2627

the testes for growth and development. When production of the male hormone (androgen) decreases, the prostate begins to degenerate. Since this androgen-dependent (AD) nature of the prostatic cells is shared by unusual malignant tumor cells, namely AD cells of prostate cancer, androgen deprivation has been the mainstay for treating advanced or metastatic prostate cancer [16]. Most androgens synthesized and secreted from the testes can be suppressed by surgical or chemical castration. The influence of adrenal androgens remaining after castration can be eliminated by the addition of nonsteroidal androgen receptor antagonists, or antiandrogens. Combination of castration and antiandrogens, referred to as combined or maximum androgen blockade, can facilitate tumor regression [17]. The serum prostate-specific antigen (PSA), which is a widely used biomarker for prostate cancer, enables monitoring of the disease progression through blood tests. Androgen deprivation initiated in a state with a high PSA level promotes apoptotic death of malignant cells and leads to a significant decrease of the PSA value implying tumor regression. However, continuation of androgen suppression not only causes side effects from the treatment but also mostly results in a relapse by recurrent tumor growth after a remission period [18]. A relapse of prostate cancer is considered to be primarily due to the progression to an androgen-independent (AI) state where the tumor is no longer responsive to androgen deprivation [19]. Many efforts have been made to understand the properties of androgen receptors and signal-transduction pathways leading to a hormone-refractory prostate cancer [20–23]. Intermittent androgen suppression (IAS) is an alternative approach to maintain tumor sensitivity to androgen deprivation with the possibility to prolong a relapse [8–15]. Although clinical trials are still in progress for establishing the clinical efficacy of IAS therapy [11,14,15], it has already been confirmed that side effects and economic costs are at least reduced in comparison with continued androgen suppression. Following these advances in experimental and clinical studies of IAS, a mathematical model of IAS therapy for prostate cancer was presented [24–26]. IAS therapy aims at delaying prostate cancer relapse by keeping androgen dependency of the tumor through repeating on-treatment and off-treatment cycles. Correspondingly, the IAS therapy model alternates between two different kinds of dynamics. The androgen concentration tends to recover and maintain the normal level a0 (nmol/l) during off-treatment periods whereas it decays to almost zero during on-treatment periods. These different androgen levels make much difference in tumor dynamics regulated by the androgen concentration and in variations of the serum PSA concentration. Interruption and reinstitution of therapeutic administration for androgen suppression are determined based on changes of the serum PSA level in the model as well as in real treatments. It is assumed in our model that administration is suspended when a decreasing PSA level drops to a lower threshold value r0 (ng/ml) whereas it is resumed when an increasing PSA level reaches an upper threshold value r1 (ng/ml) [24–26]. If this dosing strategy is successful in relapse prevention, then the serum PSA level remains around a region between r0 and r1 . Otherwise,

2617

recurrent tumor growth eventually happens and thereby the serum PSA level continues to increase explosively. In the model study, occurrence and prevention of a relapse are characterized by divergent and nondivergent solutions, respectively. Thus, parameter conditions for relapse prevention can be revealed by investigating the transition between two qualitatively different regimes. In this paper, we examine bifurcation phenomena with respect to changes of clinically controllable parameters such as a0 , r0 , and r1 . For this purpose, we utilize shooting methods for locating bifurcation sets of a limit cycle in hybrid systems. Moreover, the origin of nonlinear and chaotic dynamics of the IAS therapy model is clarified through return plot analysis of a piecewise linear system derived from the original model by simplification. Hybrid dynamical systems involving both continuous and discrete variables have been pervasive especially in computer science and control engineering [27–29]. Analysis and control of nonlinear phenomena in hybrid dynamical systems are practically significant for industrial applications such as process control [30,31], mechanical systems [32], robotic control [33], power systems [34], and power electronics [35,36]. A number of methodologies have been developed for analysis and control design of hybrid (piecewise smooth) dynamical systems [37]. In nonlinear science, stability and bifurcations in hybrid systems have been the main focus since early studies [38,39], because interactions between a trajectory and borders for discrete events can bring about complex phenomena. Recently much attention has been paid to grazing phenomena [40] and border-collision bifurcations [41]. Return map analysis is often effective for understanding chaotic behavior in low-dimensional switched dynamical systems [42]. Numerical methods for locating bifurcation points of a periodic solution have been proposed by appropriately formulating conditional equations [43]. Similar methods have been developed to formulate grazing phenomena of periodic and non-periodic trajectories [44]. Tumor growth under intermittent therapy can also be viewed as a hybrid system composed of continuous tumor dynamics and discrete events with interruption and reinstitution of medication. In the model study, the clinical issue of how to schedule on–off cycles of administration for preventing an explosion of malignant cell populations can be reduced to a mathematical problem of how to adequately set adjustable parameter values for confining a solution in a finite region in continuous state space. Therefore, it is significant to understand nonlinear phenomena peculiar to hybrid systems, such as a grazing bifurcation which occurs when a trajectory tangentially hits a boundary set related to a discrete event [40]. Present research on the IAS therapy model not only supports a possible clinical advantage of IAS therapy but also highlights the complexity generated as a fundamental feature of hybrid systems. In Section 2, we will describe shooting methods to specify bifurcation sets of a limit cycle in a hybrid system by defining an appropriate Poincar´e map. Results of bifurcation analysis of the IAS therapy model will be shown in Section 3. The essential mechanism of nonlinear and complex behavior in the IAS therapy model will also be elucidated by analysis of a

2618

G. Tanaka et al. / Physica D 237 (2008) 2616–2627

simplified model. Finally, the last section will be devoted to conclusions. 2. Method 2.1. Hybrid systems In this section, we introduce a method to analyze bifurcations of a hybrid system consisting of two ordinary differential equations described as follows: Su :

dx = fu (x), dt

u = 0, 1,

(1)

Rn ,

where t ∈ R, x ∈ and u ∈ {0, 1} denote time, continuous state vector, and the binary discrete variable, respectively. For simplicity, we assume that the vector field fu (u = 0, 1) defined on Rn is a C ∞ -class map for any state variables and parameters. We consider a hybrid automaton [37] exhibiting transitions between subsystems S0 and S1 , which describes time evolution of the continuous state x and the discrete state u. We assume that switching from S1−u to Su (u = 0, 1) occurs when a trajectory in continuous state space hits the following codimension-one surfaces: Πu = {x ∈ R | qu (x) = 0}, n

u = 0, 1,

(2)

where qu is a scalar-valued function of the continuous state vector x. The solution of subsystem Su is represented as follows: x(t) ≡ ϕu (t, x0 ),

u = 0, 1,

be noted here that the global coordinate x0 ∈ Π0 ⊂ Rn is redundant to identify any point on the codimension-one Poincar´e section. Thus, we define a local coordinate w ∈ Σ ⊂ Rn−1 and introduce the projection map h and its inverse h −1 : h : Π0 → Σ ; x0 7→ w,

(7)

h −1 : Σ → Π0 ; w 7→ x0 ,

(8)

such that the Poincar´e map on the local coordinate is represented as follows: P : Σ → Σ;

w 7→ h ◦ T ◦ h −1 (w).

(9)

Similarly, a Poincar´e map for an m-folded limit cycle crossing the Poincar´e section for m times is given by P m = h ◦T m ◦h −1 .

2.3. Stability and local bifurcation

(3)

where x0 = x(0) is an initial condition of the continuous state. 2.2. Limit cycle and Poincar´e map

A limit cycle in continuous state space corresponds to a fixed point of the Poincar´e map P, which satisfies the following equation: Flc (w) ≡ P(w) − w = 0.

Let us begin with constructing a Poincar´e map in order to formulate the conditions of a limit cycle and its bifurcations in a hybrid system. We consider a situation where a trajectory of the hybrid system (1) exhibits a nearly periodic motion with repeating transitions between S0 and S1 alternately. Fig. 1 shows a schematic illustration of a part of a trajectory crossing the sections Π0 and Π1 alternately in continuous state space. A trajectory leaving from the initial point x0 on Π0 intersects with Π1 at x1 with traveling time τ0 under the dynamics in subsystem S0 . At the moment of intersection, the subsystem changes from S0 to S1 , and then the trajectory returns to Π0 with traveling time τ1 . We define maps from one section to the other one as follows: T0 : Π0 → Π1 ; x0 7→ x1 ≡ ϕ0 (τ0 , x0 ),

(4)

T1 : Π1 → Π0 ; x1 7→ ϕ1 (τ1 , x1 ).

(5)

By treating Π0 as a Poincar´e section, a Poincar´e map T : Π0 → Π0 is defined by a composite of two submaps (4) and (5) as follows [43]: T = T1 ◦ T0 ,

Fig. 1. Schematic illustration for construction of the Poincar´e map P on the local coordinate for a trajectory of the hybrid system (1), where the trajectory under consideration exhibits a nearly periodic motion crossing two sections of discrete events alternately.

(6)

where the total return time is given by τ = τ0 + τ1 . It should

(10)

In general this equation is not analytically solvable as usual for nonlinear systems; it can be numerically solved with respect to (n − 1) unknown variables w by using numerical methods such as Newton’s method. The Jacobian matrix of Flc required for the numerical computation is given by: ∂P − In−1 , (11) ∂w where In denotes the n × n identity matrix. See Appendix A for details on the calculation procedure of first-order derivatives of the Poincar´e map P. Stability of a fixed point of P, which is equivalent to that of a corresponding limit cycle, can be evaluated by characteristic (Floquet) multipliers of the Jacobian matrix of P at the fixed point, which are the solutions of the following characteristic equation:   ∂P χ (µ) ≡ det − µIn−1 = 0. (12) ∂w DFlc (w) =

If all the exponents are located inside the unit circle on the complex plane, then the fixed point is stable. Loss of stability of the fixed point occurs when one of the exponents crosses

2619

G. Tanaka et al. / Physica D 237 (2008) 2616–2627

Fig. 2. A grazing bifurcation point λ can be characterized by a trajectory tangent to the event-triggering section Π0 at x = x0 . A trajectory hits the section at λhit > λ whereas it misses to touch the section at λmiss < λ.

Fig. 3. Construction of maps T¯k and cumulative traveling time τ¯k for an mfolded limit cycle.

can be characterized by a trajectory tangent to the eventtriggering section. If such a qualitative change occurs for a periodic solution, then it is called periodic grazing [44]. At the right grazing bifurcation with respect to the section Π0 , the directional vector of the trajectory is orthogonal to the normal vector of the section at the point of tangency. Since the transversality condition (A.14) at Π0 (see Appendix A) does not hold for a tangency point of a grazing trajectory, we need to sort out the conditions of periodic grazing by regarding the return time of a limit cycle as an unknown variable. We consider a grazing bifurcation of an m-folded limit cycle corresponding to an m-periodic point of the Poincar´e map for a natural number m in the following. Let us suppose that a grazing trajectory starting from x0 on Π0 transversally intersects with Π0 and Π1 for (2m − 1)-times alternately and then finally becomes tangent to Π0 . For convenience, we define the following maps: T¯2k : Π0 → Π0 ; x0 7→ x2k ≡ (T1 ◦ T0 )k (x0 ), k = 0, 1, . . . , m − 1, (15) k+1 ¯ T2k+1 : Π0 → Π1 ; x0 7→ x2k+1 ≡ T0 ◦ (T1 ◦ T0 ) (x0 ), k = 0, 1, . . . , m − 1, (16)

the unit circle outward from the inside by changing values of system parameters, which is accompanied with a codimensionone local bifurcation. For specifying a local bifurcation of a fixed point, the fixed point condition (10) is combined with the eigenvalue condition as follows:   P(w) − w Flb (w, λ) ≡ = 0. (13) χ(µ)

as illustrated in Fig. 3. Correspondingly, the traveling time from xk to xk+1 is denoted by τk (xk ) for k = 0, 1, . . . , 2m − 1. Then, the cumulative traveling time from x0 to xk , denoted by τ¯k , satisfies the following formula:

These equations are simultaneously solved with respect to n unknown variables (w, λ) where λ is a bifurcation parameter. We fix the value of µ for a specific type of local bifurcation: µ = 1 for saddle-node bifurcation, µ = −1 for perioddoubling bifurcation, and µ = eiθ for Neimark–Sacker bifurcation. The Jacobian matrix of Flb needed to solve Eq. (13) by Newton’s method is written as follows: ∂ P ∂P  − In−1  ∂λ  . DFlb =  ∂w (14) ∂χ (µ)  ∂χ (µ) ∂w ∂λ The first- and second-order derivatives of P required to calculate components of the matrix (14) are described in Appendix A.

T¯ (τ, x0 ) = ϕ1 (τ − τ¯2m−1 (x0 ), x2m−1 ),

2.4. Grazing bifurcation A grazing bifurcation exhibited by hybrid systems is concerned with a sudden change in the interaction between a trajectory and an event-triggering boundary for switching of subsystems [44]. When a trajectory is transversal but almost tangent to the Poincar´e section Π0 as illustrated in Fig. 2, a slight change of a system parameter value can lead to loss of the intersection and induce a change of dynamics. This critical change in system behavior, called grazing bifurcation,

τ¯k+1 (x0 ) = τ¯k (x0 ) + τk (xk ),

k = 0, . . . , 2m − 1,

(17)

where τ¯0 (x0 ) ≡ 0. The above setting provides a map dependent on the total return time τ and the initial condition x0 as follows: (18)

where τ > τ¯2m−1 (x0 ). By using the map T¯ , the conditions of a grazing bifurcation of an m-folded periodic solution can be formulated as follows:  ¯  T (τ, x0 ) − x0  q (T¯ (τ, x ))  0 0  = 0. Fgb (x0 , τ, λ) ≡  (19)   ∂q0 · f1 (T¯ (τ, x0 )) ∂x The first equation means that a trajectory leaving from x0 on Π0 is periodic with return time τ . The second equation indicates that the trajectory is located on Π0 at time τ . The third equation represents the condition that the trajectory is tangent to section Π0 at time τ . The set of Eq. (19) are simultaneously solved with respect to (n + 2) unknown variables (x0 , τ, λ) where λ is a bifurcation parameter. The Jacobian matrix of Fgb required to numerically solve Eq. (19) is given by: DFgb 

∂ T¯ − In ∂x0 ∂q0 ∂ T¯ ∂x ∂x0  2  ∂ ϕ1 ∂ T¯ 0 · ∂x ∂t∂x ∂x0

    =    ∂q

∂ T¯ ∂τ ∂q0 ∂ T¯ ∂x ∂τ  2  ∂q0 ∂ ϕ1 ∂ T¯ · ∂x ∂t∂x ∂τ

∂ T¯ ∂λ ∂q ∂ T¯ ∂q0 + ∂x ∂λ ∂λ  2  ∂q0 ∂ ϕ1 ∂ T¯ ∂ 2 ϕ1 · + ∂x ∂t∂x ∂λ ∂t∂λ

     .   

(20)

2620

G. Tanaka et al. / Physica D 237 (2008) 2616–2627



(1 − k1 )a k1 + a + k2

Refer to Appendix B for details of the calculation procedure of the derivatives of T¯ .

g D (a) = α D

3. Results

g I (a) = α I (1 − ea) − β I ,   a m(a) = m 1 1 − . a0

3.1. Mathematical model of IAS therapy The numerical method for bifurcation analysis introduced in the previous section is used to investigate the mathematical model of IAS therapy for prostate cancer [24,25]. The purpose of this section is to elucidate how dynamical behavior is influenced by some control parameters with respect to clinical importance in prostate cancer therapy. As introduced in Section 1, at least two different kinds of prostatic cells are involved in androgen suppression therapy. AD cells are sensitive to androgen suppression and are likely to induce tumor regression by their apoptosis, while AI cells are insensitive to androgen deprivation and considered to be responsible for a relapse. The IAS therapy model represents growth of a tumor consisting of a mixed dynamical assembly of AD and AI cancer cells under intermittent administration. Population growths of both AD and AI cells are dependent on androgen concentration a (nmol/l) whose dynamics relies on the binary variable u indicating whether medication of maximum androgen blockade is administered (u = 1) or not (u = 0). Administration is reinstated when an increasing PSA concentration y (ng/ml) rises up to the upper threshold value r1 , whereas it is interrupted when a decreasing PSA level falls to the lower threshold value r0 . Due to these discrete events concerning administration, the IAS therapy model is formulated as a hybrid system [24,25] consisting of subsystems described as follows: dx D = (g D (a) − m(a))x D , dt dx I = m(a)x D + g I (a)x I , dt da = −γ (a − a0 (1 − u)), dt

(21)

where continuous state variables x D and x I represent populations of AD and AI cancer cells, respectively. The discrete variable u indicates on-administration periods for u = 1 or off-administration ones for u = 0. As in clinical practice, suspension and reinstitution of administration are conducted in the following way:  1 → 0 when y = r0 and dy/dt < 0 u: (22) 0 → 1 when y = r1 and dy/dt > 0. Serum PSA concentration y serves as a biomarker for growth of a prostatic tumor as follows: y = f (x D , x I ) = c D x D + c I x I where c D = c I = 1 for the sake of simplicity. Net growth rates of AD and AI cells are denoted by g D and g I , respectively, and the rate of mutation from AD cells to AI ones is indicated by m. These factors deeply affecting tumor dynamics are given by functions of androgen concentration as follows [24,25]:



− βD



(1 − k3 )a k3 + a + k4



, (23) (24) (25)

Here g D (a) is the difference between proliferation and apoptosis rates of AD cells, which is approximately equal to α D − β D in an androgen-rich environment. Values of α D and β D can be estimated from the experimental data [45] on cell proliferation and apoptosis rates of AD cells for hormonally untreated patients, respectively. In an androgendepleted environment, AD cells are not able to proliferate and their apoptosis rate can be estimated from real data [12– 15] of the serum PSA during androgen suppression. Plausible functions with these conditions can be realized by adjusting parameters k1 , k2 , k3 , and k4 . Growth rate of AI cells is also dependent on the androgen level [46] because their growth is still dependent on the androgen receptor, although details of their androgen dependence remain unclear. Therefore we assume that g I (a) is a linear function of the androgen level [24, 25] for the sake of simplicity, which is equal to α I − β I in an androgen-depleted environment. Values of α I and β I can be estimated from experimental data [45] for hormonally failing patients. The androgen dependence of the proliferation rate of AI cells is controlled by the parameter e. We also assume that the mutation rate m(a) is linearly increasing with a decrease of androgen level, because continuation of androgen suppression is considered to enhance the mutation from AD cells to AI ones. The IAS therapy model is regarded as a hybrid automaton with continuous state vector x = (x D , x I , a) and discrete state variable u. Event-triggering sections related to discrete events (22) are defined as follows: Π0 = {x | q0 (x) = f (x D , x I ) − r0 = 0}, Π1 = {x | q1 (x) = f (x D , x I ) − r1 = 0}.

(26) (27)

A Poincar´e map T with respect to the Poincar´e section Π0 can be derived as demonstrated in Section 2. By choosing the following transformation between local and global coordinates on Π0 : h : Π0 → Σ ; (x D , x I , a) 7→ (w1 , w2 ) = (x D , a), h −1 : Σ → Π0 ; (w1 , w2 ) 7→ (x D , x I , a) = (w1 , r0 − w1 , w2 ),

(28) (29)

we obtain the Poincar´e map P = h ◦ T ◦ h −1 on the local coordinate. Consequently, bifurcation analysis of a limit cycle in the IAS therapy model is reduced to that of a fixed or a periodic point of P. 3.2. Bifurcation phenomena In the previous study [24,25], we have investigated the IAS therapy model, mainly focusing on the effect of net growth rate of AI cells. On the other hand, in the present paper, we clarify the detailed influence of clinically controllable parameters such

G. Tanaka et al. / Physica D 237 (2008) 2616–2627

2621

as threshold values of the serum PSA, r0 and r1 , and the normal androgen level a0 . The other parameters are fixed as follows [24,25]: γ = 0.08,

m 1 = 0.00005,

α D = 0.0204,

β D = 0.0076,

β I = 0.0168, k1 = 0, k2 = 2,

k3 = 8,

e = 0.015, α I = 0.0242,

(30)

k4 = 0.5.

The values of γ , m 1 , and e are set so that a relapse occurs in numerical simulations around several years after the initiation of continuous androgen suppression, which are clinically plausible. We adopt growth rates of AD and AI cells in the bone metastasis case [45]. The values of α D and β D and those of α I and β I are estimated from experimental data [45] for hormonally untreated patients and hormonally failing patients, respectively. We set k1 = 0 because AD cells do not multiply without androgens. The value of k2 is set to make the growth rate curve of the AD cells plausible. The value of β D k3 , which is the apoptosis rate of AD cells in an androgen-depleted environment, is estimated by fitting an exponential function to the decreasing serum PSA concentration during androgen suppression in actual data for simplicity [12–15]. The value of k4 is set so that the evolution of AD cells changes from decrease to increase at around a = 5. See [24,25] for detailed discussions on parameter setting. Fig. 4 shows a bifurcation diagram for variations of a0 and r0 , indicating the parameter region of relapse prevention together with bifurcation sets of limit cycles. The colored region corresponds to relapse prevention, where a solution is not divergent but confined within a finite region. The diagram shows that the lower threshold value r0 to restart administration should be relatively small for preventing a relapse by the following reasons. The first reason is that the range of a0 for relapse prevention is wider for a smaller value of r0 as shown in Fig. 4. In this regard, however, a too small value of r0 resulting in a relapse should be avoided. The second reason is that if the lower threshold value r0 is near the upper one r1 , the frequency of the switching of on-treatment and offtreatment periods is too high to adequately evaluate progression of the disease and efficacy of IAS. The normal androgen level a0 also largely affects the tumor growth and the possibility of relapse prevention. In fact, when a0 is fixed to be less than a certain value, a relapse is inevitable for any choice of r0 value. This property suggests that the normal androgen level during off-treatment periods is related to the potency of IAS therapy. Although the mainstay for treatment of advanced prostate cancer is androgen suppression, some experimental studies suggested that AI tumors that have escaped from androgen deprivation therapy might be inhibited by high levels of androgen supplementation [47,48]. From a mathematical viewpoint, the reason why a diverging solution is unavoidable for a low value of a0 is that the net growth rate of AI cells is consistently positive independently of the androgen level as given by Eq. (24) and thereby the flow in continuous state space is always expanding. This result supports the possibility that androgen supplementation during off-treatment periods might

Fig. 4. Bifurcation diagram of the IAS therapy model (21)–(25) with r1 = 30. The parameter region of relapse prevention corresponds to that of nondivergent solutions. Period-doubling and grazing bifurcations of an m-folded limit cycle are denoted by PDm and Gm , respectively. Bifurcation phenomena along the horizontal arrowed line at r0 = 25 are shown in Fig. 5.

Fig. 5. 1-parameter bifurcation diagram of the IAS therapy model (21)–(25) with r0 = 25 and r1 = 30. For a fixed value of a0 , the intersecting points with the Poincar´e section Π0 are plotted. Stable and unstable limit cycles are indicated by solid and dashed lines, respectively. Period-doubling and grazing bifurcation points are denoted by PDm and Gm , respectively.

be beneficial for delaying or averting a relapse depending on the net growth rate of AI cells. Let us further examine detailed bifurcations of limit cycles in Fig. 4. Solid and dashed curves indicate period-doubling and grazing bifurcation sets, respectively. Fig. 5 shows a oneparameter bifurcation diagram along the horizontal arrowed line at r0 = 25 in Fig. 4. With decrease of a0 , a stable limit cycle undergoes successive period-doubling bifurcations and develops into chaotic solutions. Fig. 6(a)–(c) show a limit cycle, a period-doubled limit cycle, and a chaotic attractor, respectively. The limit cycle which loses its stability at the first period-doubling bifurcation PD1 finally disappears through a grazing bifurcation G1 as shown in Fig. 6(e). In a similar way, the period-doubled limit cycle generated at PD2 also exhibits a loss of stability and vanishes at another grazing bifurcation G2 as shown in Fig. 6(d). The chaotic attractor resulting from the period-doubling cascade is considered to disappear due to transient grazing [44] by phase space analysis. In fact, for a certain fixed parameter value, it depends on the initial condition on Π0 whether a trajectory returns to Π0 again or

2622

G. Tanaka et al. / Physica D 237 (2008) 2616–2627

a period-doubling one PDm seems to be nested with increase of m for variation of the parameter a0 . The mechanism of this structure is analogous to that of a similar bifurcation pattern with interplay between border collisions and a period-doubling cascade [49]. Next we fix the value of the normal androgen level a0 and focus on how to set threshold values of the serum PSA level to stop and restart administration. Bifurcation diagrams for different values of a0 are shown in Fig. 7. The size of the parameter region for relapse prevention becomes larger as a0 increases. We have more choices of threshold values for dosing by setting the normal androgen level as high as possible within a clinically feasible range. Boundaries of the parameter area for relapse prevention are approximately given by period-doubling bifurcation sets PD1 . Fig. 7 shows that the ratio between r0 and r1 is critical rather than the respective absolute values of r0 and r1 under the dosing strategy (22). This result stems from the scaling property of the IAS therapy model. Namely, Eqs. (21) and (22) are invariant under the following variable transformation: (x D , x I , r0 , r1 ) → ( px D , px I , pr0 , pr1 ),

Fig. 6. Oribital motions of the IAS therapy model (21)–(25) with r0 = 25 and r1 = 30: (a) A stable limit cycle with a0 = 36.4; (b) A stable 2-folded limit cycle with a0 = 35.2; (c) A chaotic solution with a0 = 34.9; (d) A grazing bifurcation of an unstable two-folded limit cycle with a0 = 33.5; (e) A grazing bifurcation of an unstable limit cycle with a0 = 30.5.

not. Qualitatively different motions are separated at a boundary set on the section, where a trajectory is grazing with respect to the section. Therefore, a contact between the edge of a chaotic attractor and the boundary set of initial conditions on the section can cause a disappearance of the chaotic attractor. We will characterize transient grazing with a simplified model in Section 3.3. The pair of a grazing bifurcation point Gm and

(31)

where p is any non-zero real value. It should be noted, however, that this scaling property seems not realistic for a large p. The transformation holds if the serum PSA concentration is given as a linear sum of the populations of the AD and AI cells, i.e. y = c D x D + c I x I . Thus it may be generally more appropriate to consider the serum PSA level y = f (x D , x I ) as a nonlinear function of the populations of cancer cells in more realistic modeling. 3.3. Origin of nonlinearity A simplified version of the IAS therapy model is investigated in order to clarify the origin of nonlinear phenomena and transient grazing in the hybrid system. Let us suppose, for the sake of simplicity, that androgen concentration instantaneously changes to steady values in response to switching between onadministration and off-administration periods, i.e. γ → ∞ in Eq. (21). Androgen concentration a is then considered to take alternative values 0 and a0 . Androgen concentration switches

Fig. 7. Bifurcation diagrams of the IAS therapy model (21)–(25). Parameter regions for relapse prevention correspond to those of nondivergent solutions. Perioddoubling and grazing bifurcations of a limit cycle are denoted by PD1 and G1 , respectively. (a) a0 = 28. (b) a0 = 30. (c) a0 = 32.

G. Tanaka et al. / Physica D 237 (2008) 2616–2627

2623

from a0 to 0 when an increasing PSA level y reaches the upper threshold value r1 , whereas it switches from 0 to a0 when a decreasing PSA level touches the lower threshold value r0 . Under this assumption, the IAS therapy model given by Eqs. (21) and (22) is reduced to a piecewise linear system described as follows: dx D = (g D (a) − m(a))x D , dt dx I = m(a)x D + g I (a)x I , dt  0 → a0 when y = r0 and dy/dt < 0 a: a0 → 0 when y = r1 and dy/dt > 0.

(32) (33) (34)

The simplified model can be viewed as a hybrid system consisting of the following two linear systems:      d xD v0 0 xD = , with a = a0 , (35) 0 w0 xI dt x I      d xD v1 0 xD = , with a = 0, (36) m 1 w1 xI dt x I where v0 = g D (a0 ), w0 = g I (a0 ), v1 = g D (0) − m 1 , and w1 = g I (0). Event-triggering sections are given by Eqs. (26) and (27) with the continuous state x = (x D , x I ). Analytical solutions of individual linear equations (35) and (36) can be easily obtained. We assume that a trajectory starting from (x0 , r0 − x0 ) on Π0 reaches Π1 with traveling time τ0 and returns to Π0 with traveling time τ1 . For such a trajectory, a one-dimensional Poincar´e map (return map) can be constructed with respect to an initial condition x0 [42] as follows: R(x0 ) = ev1 τ1 +v0 τ0 x0 ,

(37)

where τ0 and τ1 are the implicit functions of x0 , satisfying f (ev0 τ0 x0 , ew0 τ0 (r0 − x0 )) − r1 = 0, v1 τ1

f (e

v1 τ1

z(x0 ), α(e

−e

w1 τ1

(38) w1 τ1

)z(x0 ) + e

− r0 = 0,

(r1 − z(x0 ))) (39)

where α = m 1 /(v1 −w1 ) and z(x0 ) = ev0 τ0 x0 . Periodic motions can be well understood by observing a sequence of intersecting points between a trajectory and the section Π0 . Fig. 8 shows return plots calculated from Eqs. (37) to (39) for different parameter values of a0 . The fixed point of the return map corresponds to a limit cycle of the simplified model. In Fig. 8(a) and (b), the fixed point is globally stable in the domain of (0, r0 ) for the initial condition x0 . In Fig. 8(c), on the other hand, the return plot ends at a certain value of x0 , at which transient grazing occurs. Transient grazing is characterized by a trajectory of the system touching the section Π0 tangentially [44]. Namely, the returning trajectory loses an intersection with the section Π0 at the point of transient grazing. The location of a trajectory starting from the initial point (x0 , r0 −x0 ) at time τ > τ0 is denoted by (ϕ D (x0 , τ ), ϕ I (x0 , τ )) where ϕ D (x0 , τ ) = ev1 (τ −τ0 )+v0 τ0 x0 ,

Fig. 8. Return plots of the simplified model (32)–(34) with r0 = 5 and r1 = 30, which is represented by Eq. (37). (a) a0 = 30. (b) a0 = 25. (c) a0 = 20.

(40)

Fig. 9. A trajectory exhibiting transient grazing. The initial condition is given by the value of x0 at which the return plot ends in Fig. 8(c).

ϕ I (x0 , τ ) = α(ev1 (τ −τ0 ) − ew1 (τ −τ0 ) )ev0 τ0 x0 + ew1 (τ −τ0 )+w0 τ0 (r0 − x0 ).

(41)

Using these notations, we can formulate conditions for a trajectory that exhibits transient grazing with respect to Π0 as follows:   ϕ D (x0 , τ ) + ϕ I (x0 , τ ) − r0 Ftg (x0 , τ ) ≡ = 0, (42) (v1 + m 1 )ϕ D (x0 , τ ) + w1 ϕ I (x0 , τ ) which are solved with respect to unknown variables (x0 , τ ). Fig. 9 shows a grazing trajectory with initial condition x0 at the endpoint of the return plot in Fig. 8(c). Next we consider bifurcations of a limit cycle of the simplified model. Conditional equations for a period-doubling bifurcation of the fixed point of the return map with a bifurcation parameter λ are described as follows:   R(x0 ) − x0 F pd (x0 , λ) ≡ = 0, (43) R 0 (x0 ) + 1 which are solved with respect to unknown variables (x0 , λ). By combining the conditional equation (42) for transient grazing with a fixed point condition, on the other hand, the conditions for a grazing bifurcation of a limit cycle can be formulated as follows:   ϕ D (x0 , τ ) − x0 ,  ϕ D (x0 , τ ) + ϕ I (x0 , τ ) − r0 F pg (x0 , τ, λ) ≡  (v1 + m 1 )ϕ D (x0 , τ ) + w1 ϕ I (x0 , τ ) = 0,

(44)

which are solved with respect to unknown variables (x0 , τ, λ). Fig. 10 shows a bifurcation diagram of the simplified model (32)–(34) with r1 = 30, indicating the parameter region of relapse prevention. The return plots at the parameter values indicated by (a)–(c) of Fig. 10 are shown in Fig. 8(a)–(c), respectively. The period-doubling bifurcation denoted by PD1 runs along the boundary of the parameter region for relapse prevention. We can also confirm that the parameter set of the grazing bifurcation denoted by G1 exists very close to that

2624

G. Tanaka et al. / Physica D 237 (2008) 2616–2627

Fig. 10. Bifurcation diagram of the simplified model (32)-(34) with r1 = 30. The parameter region of relapse prevention corresponds to that of nondivergent solutions. Period-doubling and grazing bifurcations of a limit cycle are denoted by PD1 and G1 , respectively. The inset for enlargement shows that the grazing bifurcation set exists very close to the period-doubling one. The return plots in Fig. 8(a)–(c) are observed at the parameter values indicated by (a)–(c), respectively.

of the period-doubling bifurcation as shown in the inset of Fig. 10. From comparison between Figs. 4 and 10, we can see that the parameter dependency of the region of relapse prevention is quite similar between the IAS therapy model and its simplified version. Therefore, the dynamical property of the IAS therapy model is considered to be nearly equivalent to that of the simplified system consisting of two linear systems. It suggests that nonlinear phenomena of the IAS therapy model originate from the nonlinearity generated by switching between subsystems rather than from that of individual subsystems. 4. Conclusions We have studied nonlinear dynamics and bifurcations in the hybrid systems model representing intermittent androgen suppression therapy for advanced prostate cancer. First, we have introduced a numerical method to specify a limit cycle and its bifurcation sets in a general hybrid automaton. Stability and bifurcations of a limit cycle have been reduced to those of a fixed or a periodic point of the Poincar´e map. Then, the method has been applied to analysis of the IAS therapy model [24,25]. We have considered the effect of clinically controllable parameters on relapse prevention, where a solution is not divergent but confined in a finite region of continuous state space. It has been shown that the possibility of relapse prevention can be enhanced by an increase of the normal androgen level during off-medication periods. This result has suggested that androgen supplementation can be effective for relapse prevention depending on the net growth rate of a tumor. We have shown that the ratio between the upper threshold value of the serum PSA level to restart medication and the lower one to interrupt medication plays a significant role in the dosing strategy that we have adopted in this paper. From bifurcation analysis, we have found that the model generates

complex chaotic solutions concerned with the nested structure of a pair of period-doubling and grazing bifurcations. Such a chaotic attractor should also be effective for IAS therapy because it remains within a bounded region of state space like stable limit cycles, although it may be difficult to clinically observe nonlinear behavior just by several on–off cycles in real treatment. Furthermore, we have investigated a piecewise linear system reduced from the IAS therapy model. Return plot analysis has been helpful for understanding the nonlinearity generated by the hybrid system and transient grazing. It has also been clarified that the bifurcation structure of the IAS therapy model is approximately inherited by the simplified model. Hybrid or switched dynamical systems would be a useful framework in modeling intermittent therapy which is worth considering not only for prostate cancer but also for other diseases. Treatment interruption has been practiced based on empirical or heuristic factors so far. Currently, systematic understanding of disease progression under intermittent therapy is becoming more and more important. In this regard, medication control based on a mathematical model of intermittent therapy could open a new frontier for pursuit of better treatment in biomedical science. For developing this attempt further, modeling of biomedical systems should be carefully advanced in collaboration with medical doctors and biomedical researchers as an interdisciplinary research. Acknowledgments The authors are grateful to Prof. Nicholas Bruchovsky of The Prostate Center at Vancouver General Hospital, Dr. Takashi Shimada of University of Tokyo, and Dr. Yoshito Hirata of ERATO Aihara Complexity Modelling Project of JST for valuable comments and fruitful discussions. Appendix A. Derivatives of the Poincar´e map P We need to numerically compute first- and the secondorder derivatives of the Poincar´e map P with respect to the initial condition w and a system parameter λ in order to solve conditional equations for a limit cycle or its bifurcation by Newton’s method. From Eq. (9), the derivatives of P are represented by those of T as follows: ∂P ∂h ∂ T ∂h −1 = , ∂w ∂x ∂x0 ∂w   ∂P ∂h ∂ T ∂ T ∂h −1 = + , ∂λ ∂x ∂λ ∂x0 ∂λ  2 ∂2 P ∂h ∂ 2 T ∂h −1 = , ∂x ∂x20 ∂w ∂w2 ∂2 P ∂h = ∂w∂λ ∂x

∂ 2 T ∂h −1 ∂ 2 T ∂h −1 ∂h −1 + 2 ∂x0 ∂λ ∂w ∂x0 ∂w ∂λ

(A.1) (A.2) (A.3) ! .

(A.4)

2625

G. Tanaka et al. / Physica D 237 (2008) 2616–2627

From Eq. (6), the derivatives of T are represented by those of the submaps T0 and T1 as follows: ∂T ∂ T1 = ∂x0 ∂x1

 ∂ T0 , ∂x0   ∂ T1 ∂ T1 ∂ T0 ∂T + = , ∂λ ∂x1 ∂λ ∂λ !   ∂ T1 ∂ 2 T0 ∂2T ∂ 2 T1 ∂ T0 2 , = + ∂x1 ∂x20 ∂x20 ∂x21 ∂x0      ∂ 2 T1 ∂ T0 ∂ T0 ∂ T1 ∂ 2 T0 ∂2T + = ∂x0 ∂λ ∂x1 ∂x0 ∂λ ∂λ ∂x21 ∂x0   2 ∂ T1 ∂ T0 . + ∂x1 ∂λ ∂x0 

(A.5) (A.6) (A.7)

(A.8)

From Eqs. (4) and (5), the derivatives of each submap Tu (u = 0, 1) are described as follows: ∂ Tu ∂τu ∂ϕu = + fu , ∂xu ∂xu ∂xu ∂τu ∂ϕu ∂ Tu = + fu , ∂λ ∂λ ∂λ ∂ 2 τu ∂ 2 Tu ∂ 2 ϕu ∂ 2 ϕu ∂τu ∂ 2 ϕu ∂τu = + + + f , u ∂xu ∂t ∂xu ∂t∂xu ∂xu ∂x2u ∂x2u ∂x2u ∂ 2 Tu ∂ 2 ϕu ∂ 2 ϕu ∂τu ∂ϕu ∂τu = + + ∂xu ∂λ ∂xu ∂λ ∂xu ∂t ∂λ ∂t∂λ ∂xu ∂ 2 τu . + fu ∂xu ∂λ

(A.9) (A.10) (A.11)

(A.12)

(A.13)

where u 0 = 1 − u. If a trajectory transversally intersects with the section, ∂qu 0 · fu 6= 0. ∂xu

(A.14)

Under this assumption, by differentiating both sides of Eq. (A.13) with respect to the initial condition xu , the derivatives of τu are explicitly described as follows: ∂τu 1 = − ∂q 0 u ∂xu fu



∂τu 1 = − ∂q 0 u ∂λ fu



∂x

∂qu 0 ∂ϕu ∂x ∂xu



,

(A.15)

 ∂qu 0 ∂qu 0 ∂ϕu + , (A.16) ∂x ∂λ ∂λ ∂x   1 ∂qu 0 ∂ 2 ϕu ∂ 2 ϕu ∂τu ∂ 2 ϕu ∂τu ∂ 2 τu , + = − ∂q 0 + u ∂x ∂xu ∂t ∂xu ∂t∂xu ∂xu ∂x2u ∂x2u ∂x fu (A.17)

(A.18)

The derivatives of ϕu can be calculated from first- and secondorder variational equations given as follows:     ∂fu ∂ϕu d ∂ϕu = , (A.19) dt ∂xu ∂xu ∂xu     d ∂ϕu ∂fu ∂ϕu ∂fu = + , (A.20) dt ∂λ ∂xu ∂λ ∂λ       ∂fu ∂ 2 ϕu ∂ 2 fu ∂ϕu 2 d ∂ 2 ϕu = + , (A.21) dt ∂x2u ∂xu ∂x2u ∂x2u ∂xu  2     2   ∂ ϕu ∂fu ∂ 2 fu ∂ϕu ∂ ϕu ∂ϕu d = + dt ∂xu ∂λ ∂xu ∂xu ∂λ ∂λ ∂x2u ∂xu   2 ∂ fu ∂ϕu + . (A.22) ∂xu ∂λ ∂xu The fourth-order Runge–Kutta method is employed for numerical integration of the above equations with the following initial conditions: ∂ϕu ∂ϕu = I , = On , n ∂xu t=0 ∂λ t=0 (A.23) ∂ 2 ϕu ∂ 2 ϕu = On , = On , ∂x ∂λ ∂x2 u

A trajectory leaving from the initial condition xu on Πu (u = 0, 1) is assumed to reach the other section Πu 0 with traveling time τu under the system Su as follows: qu 0 (ϕu (τu (xu ), xu )) = 0,

∂ 2 τu 1 ∂qu 0 = − ∂q 0 u ∂xu ∂λ f ∂x ∂x 2u  ∂ ϕu ∂ 2 ϕu ∂τu ∂ 2 ϕu ∂τu + + . × ∂xu ∂λ ∂xu ∂t ∂λ ∂t∂λ ∂xu

t=0

u

t=0

where In and On indicate n × n identity and null matrices, respectively. Appendix B. Derivatives of the map T¯ From Eq. (18), the derivatives of T¯ with respect to the initial condition x0 , the return time τ , and the bifurcation parameter λ, are calculated as follows: ∂ϕ1 ∂ T¯2m−1 ∂ τ¯2m−1 ∂ T¯ = − f1 , ∂x0 ∂x1 ∂x0 ∂x0 ∂ T¯ = f1 , ∂τ ∂ T¯ ∂ϕ1 ∂ T¯2m−1 ∂ τ¯2m−1 ∂ϕ1 = − f1 + . ∂λ ∂x1 ∂λ ∂λ ∂λ

(B.1) (B.2) (B.3)

In a similar way, derivatives of τ¯2m−1 can be sequentially calculated by using the following recurrence formula: ∂ τ¯k ∂ τ¯k−1 ∂τk ∂ T¯k−1 = + , k = 1, . . . , 2m − 1, (B.4) ∂x0 ∂x0 ∂xk ∂x0 ∂ τ¯k ∂ τ¯k−1 ∂τk ∂ T¯k−1 ∂τk = + + , k = 1, . . . , 2m − 1, ∂λ ∂λ ∂xk ∂λ ∂λ (B.5) which are derived from Eq. (17). From Eqs. (15) and (16), derivatives of T¯2m−1 can be calculated by using recurrence

2626

G. Tanaka et al. / Physica D 237 (2008) 2616–2627

equations for derivatives of T¯2k and T¯2k+1 given as follows:  ¯  ∂ T1 ∂ T2k−1 ∂ T¯2k = , k = 0, 1, . . . , m − 1, (B.6) ∂x0 ∂x2k−1 ∂x0   ∂ T¯2k−1 ∂ T0 ∂ T¯2k ∂ T1 + = , k = 0, 1, . . . , m − 1, ∂λ ∂x2k−1 ∂λ ∂λ (B.7)  ¯  ¯ ∂ T0 ∂ T2k ∂ T2k+1 , k = 0, 1, . . . , m − 1, (B.8) = ∂x0 ∂x2k ∂x0   ∂ T¯2k+1 ∂ T0 ∂ T¯2k ∂ T1 = + , k = 0, 1, . . . , m − 1. ∂λ ∂x2k ∂λ ∂λ (B.9) The derivatives of T0 and T1 can be calculated by Eqs. (A.9) and (A.10). References [1] H.M. Byrne, T. Alarcon, M.R. Owen, S.D. Webb, P.K. Maini, Modelling aspects of cancer dynamics: A review, Phil. Trans. R. Soc. Lond. A 364 (2006) 1563–1578. [2] A.L. Garner, Y.Y. Lau, D.W. Jordan, M.D. Uhler, R.M. Gilgenbach, Implications of a simple mathematical model to cancer cell population dynamics, Cell Prolif. 39 (1) (2006) 15–28. [3] T.L. Jackson, A mathematical model of prostate tumor growth and androgen-independent relapse, Disc. Cont. Dyn. Syst.-Series B 4 (1) (2004) 187–201. [4] T.L. Jackson, A mathematical investigation of the multiple pathways to recurrent prostate cancer: Comparision with experimental data, Neoplasia 6 (6) (2004) 697–704. [5] K.R. Swanson, L.D. True, D.W. Lin, K.R. Buhler, R. Vessella, J.D. Murray, A quantitative model for the dynamics of serum prostate-specific antigen as a marker for cancerous growth, Am. J. Pathol. 158 (6) (2001) 2195–2199. [6] P.W.A. Dayananda, J.T. Kemper, M.M. Shvartsman, A stochastic model for prostate-specific antigen levels, Math. Biosci. 190 (2004) 113–126. [7] D. Veestraeten, An alternative approach to modelling relapse in cancer with an application to adenocarcinoma of the prostate, Math. Biosci. 199 (2006) 38–54. [8] K. Akakura, N. Bruchovsky, S.L. Goldenberg, P.S. Rennie, A.R. Buckley, L.D. Sullivan, Effects of intermittent androgen suppression on androgendependent tumors: Apoptosis and serum prostate-specific antigen, Cancer 71 (9) (1993) 2782–2790. [9] A.M. Evens, T.M. Lestingi, J.D. Bitran, Intermittent androgen suppression as a treatment for prostate cancer: A review, The Oncologist 3 (1998) 419–423. [10] M.S. Bhandari, J. Crook, M. Hussain, Should intermittent androgen deprivation be used in routine clinical practice? J. Clin. Oncol. 23 (32) (2005) 8212–8218. [11] N. Mottet, C. Lucas, E. Sene, C. Avances, L. Maubach, J.M. Wolff, Intermittent androgen castration: A biological reality during intermittent treatment in metastatic prostate cancer? Urol. Int. 75 (2005) 204–208. [12] N. Bruchovsky, L.H. Klotz, M. Sadar, J.M. Crook, D. Hoffart, L. Godwin, M. Warkentin, M.E. Gleave, S.L. Goldenberg, Intermittent androgen suppression for prostate cancer: Canadian prospective trial and related observations, Mol. Urol. 4 (3) (2000) 191–199. [13] N. Bruchovsky, S.L. Goldenberg, N.R. Mawji, M.D. Sadar, Evolving aspects of intermittent androgen blockade for prostate cancer: Diagnosis and treatment of early tumor progression and maintenance of remission. In: Andrology in the 21st Century, Proceedings of the VIIth International Congress of Andrology, 2001, pp. 609–623.

[14] N. Bruchovsky, L. Klotz, J. Crook, S. Malone, C. Ludgate, W.J. Morris, M.E. Gleave, S.L. Goldenberg, Final results of the canadian prospective phase II trial of intermittent androgen suppression for men in biochemical recurrence after radiotherapy for locally advanced prostate cancer, Cancer 107 (2) (2006) 389–395. [15] N. Bruchovsky, L. Klotz, J. Crook, S. Larry, S.L. Goldenberg, Locally advanced prostate cancer – biochemical results from a prospective phase II study of intermittent androgen suppression for men with evidence of PSA relapse after radiotherapy, Cancer 109 (5) (2007) 858–867. [16] D.G. McLeod, Hormonal therapy: Historical perspective to future directions, Urology 61 (2003) 3–7. [17] L. Klotz, Combined androgen blockade in prostate cancer: Meta-analyses and associated issues, BJU Int. 87 (2001) 806–813. [18] E.D. Crawford, M. Rosenblum, A.M. Ziada, P.H. Lange, Overview: Hormone refractory prostate cancer, Urology 54 (Suppl. 1) (1999) 1–7. [19] M. Diaz, S.G. Patterson, Management of androgen-independent prostate cancera, Cancer Control 11 (6) (2004) 364–373. [20] B.J. Feldman, D. Feldman, The development of androgen-independent prostate cancer, Nat. Rev. Cancer 1 (2001) 34–45. [21] S.M. Dehm, D.J. Tindall, Regulation of androgen receptor signaling in prostate cancer, Expert Rev. Anticancer Ther. 5 (1) (2005) 63–74. [22] J. Edwards, J. Bartlett, The androgen receptor and signal-transduction pathways in hormone-refractory prostate cancer. Part 1: Modifications to the androgen receptor, BJU Int. 95 (2005) 1320–1326. [23] M.-E. Taplin, Drug insight: Role of the androgen receptor in the development and progression of prostate cancer, Nat. Clin. Pract. Oncol. 4 (4) (2007) 236–244. [24] A.M. Ideta, G. Tanaka, T. Takeuchi, K. Aihara, A mathematical model of intermittent androgen suppression remedy for prostate cancer, Technical Report, Department of Mathematical Informatics, Graduate School of Information Science and Technology, University of Tokyo, 2006-32, May 2006. [25] A.M. Ideta, G. Tanaka, T. Takeuchi, K. Aihara, A mathmatical model of intermittent androgen suppression for prostate cancer, J. Nonlinear Sci. (in press). [26] K. Aihara, G. Tanaka, T. Suzuki, Y. Hirata, A hybrid systems approach to hormonal therapy of prostate cancer and its nonlinear dynamics, in: Noise and Fluctuations: 19th International Conference on Noise and Fluctuations, Sep 2007, pp. 479–482. [27] A.J. van der Schaft, J.M. Schumacher, An Introduction to Hybrid Dynamical Systems, in: Lecture Notes in Control and Information Sciences, vol. 251, Springer, 2000. [28] A.V. Savkin, R.J. Evans, Hybrid Dynamical Systems: Controller and Sensor Switching Problems, Springer-Verlag, 2001. [29] D. Liberzon, Switching in Systems and Control, Birkhauser, Boston, 2003. [30] H.E. Garcia, A. Ray, R.M. Edwards, A reconfigurable hybrid systems and its application to power plant control, IEEE Trans. Control Syst. Tech. 3 (2) (1995) 157–170. [31] B. Lennartson, M. Tittus, B. Egardt, S. Pettersson, Hybrid systems in process control, Contr. Syst. Mag. 16 (1996) 45–56. [32] B. Brogliato, Nonsmooth Mechanics: Models, Dynamics, and Control, 2nd edition, Springer-Verlag, New York, 1999. [33] M. Spong, M. Vidyasagar, Robot Dynamics and Control, Wiley, New York, 1989. [34] I. Hiskens, Power system modeling for inverse problems, IEEE Trans. Circuits Syst.-I 51 (2004) 539–551. [35] S. Banerjee, G.C. Verghese (Eds.), Nonlinear Phenomena in Power Electronics: Attractors, Bifurcations, Chaos, and Nonlinear Control, IEEE Press, New York, 2001. [36] M. di Bernardo, C.K. Tse, Chaos in power electronics: An overview, in: Chaos in Circuits and Systems, World Scientific, Singapore, 2002, pp. 317–340. [37] P.J. Antsaklis, X.D. Koutsoukos, Hybrid systems: Review and recent

G. Tanaka et al. / Physica D 237 (2008) 2616–2627

[38] [39]

[40]

[41]

[42] [43]

progress, in: T. Samad, G. Balas (Eds.), Software-Enabled Control: Information Technology for Dynamical Systems, IEEE Press, 2003, pp. 272–298. S.D. Johnson, Simple hybrid systems, Int. J. Bifurcat. Chaos 4 (6) (1994) 1655–1665. J. Guckenheimer, S. Johnson, Planar hybrid systems, in: Hybrid Systems II, in: Lecture Notes in Computer Science, vol. 999, 1995, pp. 202–225. M. di Bernardo, C.J. Budd, A.R. Champneys, Grazing and bordercollision in piecewise-smooth systems, Phys. Rev. Lett. 86 (12) (2001) 2553–2556. H.E. Nusse, E. Ott, J.A. Yorke, Border-collision bifurcations: An explanation for observed bifurcation phenomena, Phys. Rev. E 49 (1994) 1073–1077. T. Saito, S. Nakagawa, Chaos from a hysteresis and switched circuit, Philos. Trans. R. Soc. Lond. A 353 (1995) 47–57. T. Kousaka, T. Ueta, H. Kawakami, Bifurcation of switched nonlinear dynamical systems, IEEE Trans. Circ. Syst.-II 46 (7) (1999) 878–885.

2627

[44] V. Donde, I.A. Hiskens, Shooting methods for locating grazing phenomena in hybrid systems, Int. J. Bifurcat. Chaos 16 (3) (2006) 671–692. [45] R.R. Berges, J. Vukanovic, J.I. Epstein, M. CarMichel, L. Cisek, D.E. Johnson, R.W. Veltri, P.C. Walsh, J.T. Isaacs, Implication of cell kinetic changes during the progression of human prostatic cancer, Clin. Cancer Res. 1 (1995) 473–480. [46] J. Kokontis, K. Takakura, N. Hay, S. Liao, Increased androgen receptor activity and altered c-myc expression in prostate cancer cells after long-term androgen deprivation, Cancer Res. 54 (1994) 1566–1573. [47] Y. Umekita, R.A. Hiipakka, J.M. Kokontis, S. Liao, Human prostate tumor growth in athymic mice: Inhibition by androgens and stimulation by finasteride, Proc. Natl. Acad. Sci. 93 (1996) 11802–11807. [48] R.T. Prehn, On the prevention and therapy of prostate cancer by androgen administration, Cancer Res. 59 (1999) 4161–4164. [49] Y. Ma, H. Kawakami, Bifurcation analysis of switched dynamical systems with periodically moving borders, IEEE Trans. Circ. Syst.-I 51 (6) (2004) 1184–1193.