A new protocol for intermittent androgen suppression therapy of prostate cancer with unstable saddle-point dynamics

Journal of Theoretical Biology 350 (2014) 1–16

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

A new protocol for intermittent androgen suppression therapy of prostate cancer with unstable saddle-point dynamics Yasuyuki Suzuki a, Daichi Sakai a, Taishin Nomura a,n, Yoshito Hirata b, Kazuyuki Aihara b a b

Graduate School of Engineering Science, Osaka University, Toyonaka, Osaka 560-8531, Japan Institute of Industrial Science, The University of Tokyo, Komaba, Meguro-ku, Tokyo, Japan

H I G H L I G H T S








We propose a new protocol for intermittent hormonal therapy for prostate cancer. Tumor cell population dynamics are modeled by a hybrid dynamical system. Dynamics are determined by event-driven transition between two unstable subsystems. Saddle-point dynamics of the subsystems are deeply involved in the stabilization. Clever use of the stable manifold provides robust stability and successful therapy.

art ic l e i nf o

a b s t r a c t

Article history: Received 9 September 2013 Received in revised form 1 February 2014 Accepted 3 February 2014 Available online 10 February 2014

Intermittent androgen suppression (IAS) therapy is a class of hormonal treatment for prostate cancer, in which a drug-induced androgen deprivation can reduce the population of prostate cancer cells. In IAS therapy, drugs are administrated only in on-treatment periods that are separated intermittently by off-treatment periods. The presence of off-treatment periods may be beneficial for maintaining the sensitivity of the tumor cells to androgen deprivation. Thus, IAS can be superior to continuous androgen suppression (CAS) for delaying or possibly preventing relapse of a tumor. IAS therapy usually monitors the level of serum prostate-specific antigen (PSA), which is related to the population of tumor cells. Each on-treatment period begins when the PSA level is greater than an upper threshold; treatment results in a decrease in the PSA level. The on-treatment period is suspended when the PSA level falls below a lower threshold; the PSA level then rises again until the beginning of the next on-treatment period. To determine the transitions between on- and off-treatment periods, we propose a new IAS protocol that uses a model-based estimate of the state point in the phase space of the tumor dynamics. We show that the proposed protocol is effective if, in each of the on- and off-treatment periods, the tumor dynamics exhibits a saddle-point instability accompanied by a stable manifold. Mathematical analysis reveals that tumor dynamics can be controlled in a more effective and robust manner with the proposed protocol than with conventional IAS. We also discuss the clinical feasibility of the proposed protocol as an alternative to conventional IAS therapy. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Prostate cancer Intermittent androgen suppression Mathematical model Hybrid dynamical system Intermittent control

1. Introduction The number of patients with prostate cancer is increasing worldwide (Siegel et al., 2011). In many cases, patients with prostate cancer have few if any specific symptoms (Pentyala et al., 2000). Therefore, early detection is often achieved by the measurement of prostate-specific antigen (PSA) in the blood (Pentyala et al., 2000), and appropriate treatment can follow detection.

n

Corresponding author. Tel.: þ 81 6 6850 6532; fax: þ 81 6 6850 6534. E-mail address: [email protected] (T. Nomura).

http://dx.doi.org/10.1016/j.jtbi.2014.02.004 0022-5193 & 2014 Elsevier Ltd. All rights reserved.

It is known that androgen, the male hormone, enhances the growth of prostate cancer (Charles and Clarence, 1941; Tan et al., 2012). Thus, deprivation of androgen can be an efficient treatment for prostate cancer. Continuous androgen suppression (CAS) therapy is a class of hormonal treatment that exploits the fact that androgen deprivation, which is achieved by administration of pharmacological agents, can reduce the population of tumor cells in advanced prostate cancer. However, many patients suffer from a relapse of the cancer (Pollack et al., 1994; Yu et al., 1995; Feldman and Feldman, 2001; Heidenreich et al., 2008), which is clinically recognized by a re-increase in the concentration of PSA. Relapse is rooted in a mutation of the prostate cancer from androgen

2

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

dependent (AD) cancer cells to androgen independent (AI) cancer cells (Feldman and Feldman, 2001; Edwards and Bartlett, 2005). Prevention of relapse is one of the primary issues in the treatment of cancer (Sharifi et al., 2005; Oudard, 2013). Intermittent androgen suppression (IAS) is a strategy to delay or possibly prevent the relapse (Akakura et al., 1993; Gleave et al., 1998; Mackean et al., 1998; Crook et al., 1999; Mottet et al., 2012). In IAS therapy, drugs are administrated only in intermittent ontreatment periods, which are separated by off-treatment periods. Recent studies have revealed that the off-treatment periods in IAS therapy help to maintain the sensitivity of the tumor cells to androgen deprivation; i.e., the off-treatment periods inhibit the mutation from AD to AI cells. Ideta et al. (2008) were the first to mathematically describe the dynamics of a tumor being treated by IAS. They analyzed the population dynamics of AD cells, AI cells, and androgen, and demonstrated that IAS therapy was superior to CAS therapy. The model of Ideta et al. (2008) has been extended in various ways (Shimada and Aihara, 2008; Guo et al., 2008; Tao et al., 2009, 2010; Hirata et al., 2010; Jain et al., 2011; Portz et al., 2012; Jain and Friedman, 2013; Tao et al., in press). Among these, Hirata et al. (2010) performed nonlinear system identification on clinical time-series data of PSA in patients with prostate cancer under IAS therapy. They obtained dynamical system models that described the population dynamics of AD cells and two types of AI cells and that quantitatively reproduced the clinical data. In the studies of Ideta et al. (2008) and Hirata et al. (2010), the population dynamics of tumor cells with IAS therapy was modeled by hybrid dynamical systems. Specifically, tumor dynamics during ontreatment periods was modeled by one system, referred to here as the on-model, and that in off-treatment periods by another system, referred to as the off-model; thus, these two subsystems are alternated in governing the overall dynamics. Thus far, in conventional IAS therapy and thus also in the mathematical models of Ideta et al. (2008) and Hirata et al. (2010), each on-treatment period begins when the PSA level grows beyond a predetermined upper threshold, and results in a decrease in the PSA level. Treatment is suspended when the PSA level falls below a lower threshold; this is the onset of an off-treatment period during which the PSA level increases. Those two types of threshold-crossing events trigger the transitions between the on-model and the offmodel. We refer to this control protocol as the threshold method. Typically, the threshold method induces a sustained oscillation in the PSA level, in which the maximum and minimum levels are roughly determined by the upper and the lower thresholds, respectively. Mathematically, this oscillation corresponds to a limit cycle of the hybrid dynamical system model. In this way, IAS therapy with the threshold method achieves a bounded population of cancer cells. In this paper, we focus on the fact that the dynamics of each of the on-models and the off-models has a saddle-point instability for the typical cases explored in Ideta et al. (2008) and Hirata et al. (2010). In each on-treatment period of either IAS or CAS therapy, the effective degradation rate of the AD cells is responsible for the decrease in the PSA level during the early phase of the treatment, and this process is accompanied by a gradual increase in the number of AI cells. The corresponding dynamics is represented as follows: while moving along the stable manifold toward the saddle-point equilibrium located at the origin, the state point of the on-model is attracted by the unstable manifold.1 The increase in AI cells during the later phase of the treatment is responsible for a relapse; this decrease is represented by the movement of the state point away from the origin along the

1 In general, an unstable manifold of a saddle-point equilibrium is attractive. Thus, a state point approaches an unstable manifold as it moves away from such an equilibrium point.

unstable manifold. Similarly, in the absence of treatment (during each off-treatment period of IAS therapy), the population growth of AD cells is responsible for an increase in the PSA level; this is represented by the movement of the state point along the unstable manifold of the off-model. At the same time, the population of AI cells decreases; this is represented by the movement of the state point along the stable manifold of the off-model. In this way, the decrease and increase of the PSA level in the limit-cycle oscillation of IAS therapy, respectively, correspond to the movement of the state point along the stable manifold of the on-model and the unstable manifold of the off-model. Here, based on dynamics that was revealed in previous studies, we propose a new IAS-based tumor control protocol, which we will call the region-dividing method. In this method, we combine the unstable saddle-point dynamics of the on- and off-models in order to stabilize the overall dynamics of IAS therapy. We modify the IAS models of Ideta et al. (2008) and Hirata et al. (2010) in order to examine the proposed control protocol. In Section 2, we introduce the original IAS models. In Section 3, we illustrate how the new protocol can stabilize an equilibrium point located at the origin of the hybrid dynamical system for IAS therapy, and provide a theoretical basis for the region-dividing method. Then, in Section 4, we define modified IAS models that include the region-dividing method and are based on Ideta et al. (2008) and Hirata et al. (2010). In Section 5, we investigate the dynamics of the modified IAS models, and compare that with the dynamics of the original threshold method. In Section 6, we compare the stability of the models that use the region-dividing method with those that use the threshold method, when values of important parameters are changed. We also show that the parameter regions where tumor growth is successfully controlled by the region-dividing method are almost identical to the parameter regions where both the onand off-models have saddle-point equilibrium points. We then discuss, in Section 7, the clinical feasibility of using the proposed protocol as an alternative to conventional IAS therapy.

2. Models of conventional IAS therapy In this section, we briefly introduce the two mathematical models of conventional IAS therapy that were proposed by Ideta et al. (2008) and Hirata et al. (2010); we will refer to them as the ITTA model and HBA model, respectively. In each model, the dynamics of the IAS therapy are described by a hybrid dynamical system consisting of an on-model and an off-model. Transitions between the on-model and the off-model are triggered by threshold-crossing events in the serum PSA concentration, as monitored by a blood test. 2.1. ITTA model The ITTA model describes the androgen dynamics, the population dynamics of AD and mutated AI cells, and the serum concentration of PSA. The androgen dynamics in the presence of treatment is described by d aðtÞ ¼  γ aðtÞ; dt

ð1Þ

where a(t) represents the serum concentration of androgen, which decays to zero during the on-treatment period with a time constant of 1=γ . In the absence of treatment, the androgen concentration is regulated at a constant value of the normal androgen concentration in the untreated condition, a0, as follows: d aðtÞ ¼  γ ðaðtÞ a0 Þ: dt

ð2Þ

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

The population dynamics of AD and mutated AI cells are described as ! ! ! 0 v1;1 ðaðtÞÞ xAD ðtÞ d xAD ðtÞ ¼ ; ð3Þ v2;1 ðaðtÞÞ v2;2 ðaðtÞÞ xAI ðtÞ dt xAI ðtÞ where xAD(t) and xAI(t) represent the population of AD and AI cells, respectively. The net growth rates are determined by v1;1 and v2;2 , and the mutation rate by v2;1 . Since those rates are modeled as static functions of the serum androgen concentration, the population dynamics of the tumor cells are characterized differently during the on- and off-treatment periods. See Appendix A for details. The serum PSA concentration y is, for simplicity, defined as the linear weighted sum of xAD and xAI with the coefficients c1 and c2 as yðtÞ ¼ c1 xAD ðtÞ þ c2 xAI ðtÞ:

ð4Þ

In the ITTA model, the system that results when Eq. (3) is coupled with Eq. (1) defines the on-model, and when it is coupled with Eq. (2), it defines the off-model. In IAS therapy, the transition from the off-model to the on-model is triggered when the serum PSA concentration y grows beyond the upper threshold r on , and the transition from the on-model to the off-model is triggered when y falls below the lower threshold r off . Parameter values in the ITTA model are fixed according to Ideta et al. (2008), as shown in Table 1. A block diagram of the ITTA model is shown in Fig. 1(a).

2.2. HBA model The HBA model describes the population dynamics of AD cells, reversible AI (rAI) cells, irreversible AI (iAI) cells, and the serum concentration of PSA. In the presence of treatment, AD cells can change to either rAI cells or to iAI cells, but rAI cells can change only to iAI cells. It is noted, however, that in the absence of treatment, rAI cells can change to AD cells. In the HBA model, the

Table 1 Parameters for the original ITTA model and the modified ITTA model. Symbol

Description

Value/unit

a0 γ α1 β1 α2 β2 k1 k2 k3 k4 m1 d c1

Normal androgen concentration Rate constant of androgen dynamics Proliferation rate of AD cells Apoptosis rate of AD cells Proliferation rate of AI cells Apoptosis rate of AI cells Androgen sensitivity of AD cells Androgen sensitivity of AD cells Androgen sensitivity of AD cells Androgen sensitivity of AD cells Maximum mutation rate from AD to AI cells Net growth rate of AI cells Rate of contribution of AD cells to the serum PSA concentration Rate of contribution of AI cells to the serum PSA concentration Threshold for a transition from off to on model Threshold for a transition from on to off model Region-dividing parameter Region-dividing parameter Nonintervention period Parameter used to determine Kon Parameter used to determine Kon Parameter used to determine Koff Parameter used to determine Koff

30.0 nM 0.08 0.0204 days  1 0.0076 days  1 0.0242 days  1 0.0168 days  1 0.0 2.0 8.0 0.5 0.00005 1.0 1.0

c2 r on r off lon loff τ μon;1 μon;2 μoff ;1 μoff ;2

1.0 15.0 ng/ml 0.2 ng/ml 0.2 8.0 60 days 1 1 1 1

on-model is 0 1 0 on w1;1 xAD ðtÞ dB C B won ðtÞ x @ rAI A ¼ @ 2;1 dt won xiAI ðtÞ 3;1

3

0 won 2;2 won 3;2

0

10

xAD ðtÞ

1

0

xAD ðtÞ

1

B C B C 0 C A@ xrAI ðtÞ A  Won @ xrAI ðtÞ A; xiAI ðtÞ xiAI ðtÞ

won 3;3

ð5Þ and the off-model is 0 1 0 woff xAD ðtÞ 1;1 dB C B 0 @ xrAI ðtÞ A ¼ B @ dt xiAI ðtÞ 0

woff 1;2 woff 2;2 0

10

1 0 1 xAD ðtÞ xAD ðtÞ CB C B 0 C@ xrAI ðtÞ A  Woff @ xrAI ðtÞ C A; A xiAI ðtÞ xiAI ðtÞ woff 3;3 0

ð6Þ where xAD(t), xrAI(t), and xiAI(t) represent the populations of AD on cells, rAI cells, and iAI cells, respectively. The parameters won 1;1 , w2;2 , and won represent the growth rates during the on-treatment 3;3 periods of the populations of AD cells, rAI cells, and iAI cells, off off respectively, and the parameters woff 1;1 , w2;2 , and w3;3 , respectively, represent them during the off-treatment periods. The parameters on on won 2;1 , w3;1 , and w3;2 are the change rates during the on-treatment periods from AD to rAI cells, from AD to iAI cells, and from rAI to iAI cells, respectively. The parameter woff 1;2 is the change rate from rAI to AD cells during the off-treatment periods. The serum PSA concentration y is the linear weighted sum of xAD, xrAI, and xiAI, with the coefficients c1, c2, and c3, respectively, as follows: yðtÞ ¼ c1 xAD ðtÞ þ c2 xrAI ðtÞ þ c3 xiAI ðtÞ:

ð7Þ

As in the ITTA model, transitions between the on-model and the off-model are triggered by threshold-crossing events in the serum PSA concentration y. However, in the HBA model, once an on-treatment period has begun, it is maintained for at least 36 weeks, even if the serum PSA concentration drops below the lower threshold, according to their clinical references (Hirata et al., 2010; Bruchovsky et al., 2007). We refer to this 36-week period as the holding period. Therefore, in the HBA model, the transition from the on-model to the off-model is triggered when y falls below the lower threshold r off only if the holding period has been completed. The transition criterion from the off-model to the on-model is the same as it is for the ITTA model, i.e., it is triggered when y grows beyond the upper threshold r on . The parameter values in the HBA model are fixed according to Hirata et al. (2010), as shown in Table 2. A block diagram of the HBA model is shown in Fig. 1(b).

3. A new control protocol for IAS therapy In this section, we illustrate how, for the hybrid dynamical system for IAS therapy, a new protocol can stabilize the equilibrium point located at the origin. We then provide a theoretical basis for the region-dividing method. For simplicity, we consider a simplified version of the ITTA model; this simplifies our example and theoretical analysis while leaving unchanged the essential dynamics of the model. 3.1. The on-model and off-model of the simplified ITTA model The androgen concentration in the ITTA model, as shown in Eqs. (1) and (2), varies with the time constant 1=γ , which is about 12.5 days. This time constant is small enough in comparison with the typical time scale of the on- and off-treatment periods so that we can approximate the androgen concentrations in each of the on- and off-treatment periods by the corresponding steady-state values a ¼0 and a ¼ a0 , respectively. This approximation allows us to describe the population dynamics of prostate cancer cells, using

4

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

Fig. 1. Block diagrams of IAS therapy models with the threshold method and those of the proposed IAS therapy models with the region-dividing method. (a) Block diagram ITTA ITTA of the ITTA model (original model). Σ ITTA represent the androgen dynamics (Eqs. (1) and (2)), the population dynamics of AD and AI cells (Eq. (3)), 1;k ðk ¼ on; offÞ, Σ 2x , and Σ 2y HBA and the serum PSA concentration measurement (Eq. (4)). (b) Block diagram of the HBA model (original model). Σ HBA 2x;k ðk ¼ on; offÞ and Σ 2y represent the population dynamics ITTA of AD and AI cells (Eqs. (5) and (6)) and the serum PSA concentration measurement (Eq. (7)). (c) Block diagram of the modified ITTA model. Σ 1;k ðk ¼ on; offÞ represents the ITTA approximated androgen dynamics as a~ ¼ 0 for on-treatment periods and a~ ¼ a0 for off-treatment periods. Σ 2x represents the state observer (Eq. (47)). (d) Block diagram of HBA the modified HBA model. Σ 2x;k ðk ¼ on; offÞ represents the state observer (Eq. (51)).

the ITTA model, by the following set of two simple linear dynamical system models (Tanaka et al., 2008): ! ! xAD ðtÞ d xAD ðtÞ ¼ Ak ; ð8Þ xAI ðtÞ dt xAI ðtÞ

triangular matrices, their eigenvalues are easily obtained as

where Ak for k A fon; offg is the constant state matrix defined as ! v1;1 ðak Þ 0 Ak ¼ ; ð9Þ v2;1 ðak Þ v2;2 ðak Þ with aon ¼ 0 and aoff ¼ a0 . We refer to Eq. (8) as the simplified ITTA model. In order to characterize the tumor dynamics, we first observe the signs and values of the elements vi;j of the system matrix Ak for k A fon; offg. Recalling that v1;1 ðak Þ, v2;1 ðak Þ, and v2;2 ðak Þ represent, respectively, the net growth rate of AD cells, the mutation rate from AD cells to AI cells, and the net growth rate of AI cells, we have several inequalities satisfied by the elements of the system matrix. The first two inequalities are as follows: v2;1 ðaon Þ Z 0;

ð10Þ

v2;1 ðaoff Þ Z 0:

ð11Þ

These two inequalities hold because the mutation of AD cells to form AI cells always increases the population of AI cells. The next inequality that must be satisfied by the simplified ITTA model is v1;1 ðaon Þ o v1;1 ðaoff Þ;

ð12Þ

which represents the fact that the net proliferation rate of AD cells in the off-treatment period with a high androgen concentration is higher than that in the on-treatment period with a low androgen concentration. The signs of the elements v1;1 ðak Þ and v2;2 ðak Þ can be either positive or negative depending on the choice of values for the parameters used in the model. For each of k ¼ on and k ¼ off, the origin ðxAD ¼ xAI ¼ 0Þ is the equilibrium point of Eq. (8). Let λ Aon , and let λ

ð1Þ off

and λ

ð2Þ off

ð1Þ on

and λ

ð2Þ on

be the eigenvalues of

be those of Aoff . Since Aon and Aoff are

λð1Þ on ¼ v1;1 ðaon Þ;

ð13Þ

λð2Þ on ¼ v2;2 ðaon Þ;

ð14Þ

λð1Þ off ¼ v1;1 ðaoff Þ;

ð15Þ

λð2Þ off ¼ v2;2 ðaoff Þ:

ð16Þ

Since the elements vi;j of the system matrices are real numbers, we observe that Aon and Aoff each always possess two real eigenvalues. Moreover, since v1;1 ðak Þ a v2;2 ðak Þ in general, i.e., the net proliferation rates of AD cells and AI cells are not generally the same, the eigenvalues of Aon and Aoff are distinct. That is ð2Þ λð1Þ k a λk

ð17Þ

for k A fon; offg. Furthermore, for both the on-model and the offmodel, at least one of the two eigenvalues is positive; otherwise, the dynamics of the model would be stable, which is not the case for growing tumors. Denoting the eigenvectors corresponding to the eigenvalues ð2Þ ð1Þ ð2Þ ð1Þ ð2Þ ð1Þ ð2Þ λð1Þ on , λon , λoff , and λoff by x on , x on , x off , and x off , respectively, we

have xð1Þ on ¼

xð2Þ on ¼

xð1Þ ¼ off

¼ xð2Þ off

1

!

v2;1 ðaon Þ v1;1 ðaon Þ  v2;2 ðaon Þ

  0 1

;

;

ð19Þ 1

v2;1 ðaoff Þ v1;1 ðaoff Þ  v2;2 ðaoff Þ

  0 : 1

ð18Þ

! ;

ð20Þ

ð21Þ

For the parameter values presented in Table 1, each of Aon and Aoff has one positive and one negative eigenvalue. This means that,

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

5

Table 2 Parameters for the original HBA model and the modified HBA model. Symbol Description

Value/unit

won 1;1

Growth rate of AD cells during on-treatment periods  0.043844

won 2;1

Change rate from AD to rAI cells during 0.00130716 on-treatment periods Growth rate of rAI cells during on-treatment periods  0.004282

won 2;2

won 3;3

Change rate from AD to iAI cells during 5.8852  10  17 on-treatment periods Change rate from rAI to iAI cells during 9.01608  10  17 on-treatment periods Growth rate of iAI cells during on-treatment periods 0.00181

woff 1;1

Growth rate of AD cells during off-treatment periods 0.00223

woff 1;2

won 3;1 won 3;2

woff 2;2

Change rate from rAI to AD cells during 0.1 off-treatment periods Growth rate of rAI cells during off-treatment periods 0.00392

woff 3;3

Growth rate of iAI cells during off-treatment periods  0.2

c1

Rate of contribution of AD cells to the serum PSA concentration Rate of contribution of rAI cells to the serum PSA concentration Rate of contribution of iAI cells to the serum PSA concentration Threshold for a transition from off to on model Threshold for a transition from on to off model Region-dividing parameter Region-dividing parameter Nonintervention period Parameter used to determine Kon Parameter used to determine Kon Parameter used to determine Kon Parameter used to determine Koff Parameter used to determine Koff Parameter used to determine Koff

c2 c3 r on r off son soff τb μon;1 μon;2 μon;3 μoff ;1 μoff ;2 μoff ;3

1.0 1.0 1.0 10.0 ng/ml 4.0 ng/ml 7.0 1.0 60 days 1 1 1 1 1 1

for each case, there is an unstable saddle-point equilibrium accompanied by an unstable manifold and a stable manifold. Let us denote the stable and unstable manifolds of the simplified ITTA model in the presence of treatment (Eq. (8) with k ¼ on) by Eson and Euon , respectively (Fig. 2(a)). This situation corresponds to ð1Þ ð2Þ λon o 0 and λon 4 0. In this case, the stable manifold Eson is a line passing through the origin and in the direction of the eigenvector u of xð1Þ on , and the unstable manifold E on is a line passing through the origin and in the direction of the eigenvector of xð2Þ on that is identical with the xAI axis. Similarly, we denote the unstable and stable manifolds of the simplified ITTA model in the absence of treatment (Eq. (8) with k ¼ off) by Euoff and Esoff , respectively ð1Þ ð2Þ (Fig. 2(b)). In this case, λoff 40 and λoff o 0, and Euoff and Esoff are the lines passing through the origin and in the direction of xð1Þ and off xð2Þ , respectively. The stable manifold Esoff is identical with the off xAI axis. 3.2. The region-dividing method When the state point ðxAD ; xAI Þ is located close to Eson during an on-treatment period, it transiently approaches the equilibrium point along Eson , and then moves away from the equilibrium point along Euon (Fig. 2(a)). Similarly, when the state point is located close to Esoff during an off-treatment period, it transiently approaches the equilibrium point along Esoff , and then moves away from the equilibrium point along Euoff (Fig. 2(b)). In this way, the state point transiently approaches the equilibrium when it is near the stable manifold, either Eson or Esoff , depending on the treatment status. The region-dividing method cleverly uses these characteristics. When the state point is located in the neighborhood of Eson in the phase plane of the system, we select the on-model. On the other hand, when the state point is in the neighborhood of Esoff , we select the off-model. To achieve this state-dependent switching, we

Fig. 2. Phase-plane analysis of the simplified ITTA model. (a) Vector field, stable and unstable manifolds, and a sample trajectory of the simplified ITTA model in the presence of treatment (the on-model, a ¼0). (b) The same as (a) in the absence of treatment (the off-model, a ¼ a0 ). (c) Region-division used in the modified ITTA model. The phase plane is divided into three regions by two straight lines, Lon and Loff . Ron (the red region) and Roff (the blue region) represent the neighborhood of stable manifolds of the on-model and the off-model, respectively. The remaining region in the first quadrant is Rhold . (d) Typical trajectory of the simplified ITTA model. p0 is the initial state, p1 and p3 are the state points at the transitions from the on-model to the off-model, and p2 and p4 are the state points at the transitions from the off-model to the on-model. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

divide the xAD–xAI phase plane into three regions, Ron , Roff , and Rhold , by adding two straight boundary lines, Lon and Loff (Fig. 2(c)), defined as follows: Lon : xAI ¼ lon xAD ;

ð22Þ

Loff : xAI ¼ loff xAD ;

ð23Þ

where lon and loff are the parameters that determine the slopes of the boundary lines. The zigzag-shaped solid black curve in Fig. 2(d) represents a sample trajectory from an initial state point (p0) under the simplified ITTA model controlled by the region-dividing method. The segments from p0 to p1 and from p2 to p3 are the trajectories in the presence of treatment (the on-model), whereas those from p1 to p2 and from p3 to p4 are the trajectories in the absence of treatment (the off-model). We can see that by appropriately alternating between the on-model and the off-model, the state point converges to the origin ðxAD ; xAI Þ ¼ ð0; 0Þ, which corresponds to extermination of the cancer cells. In the first on-treatment period from p0 to p1, the state point moves according to the saddle-point flow of the on-model (Fig. 2(a)), where the population of AD cells decreases, but that of AI cells gradually increases. Then the state point enters the region Roff at p1 (Fig. 2(c) and (d)), and this triggers a transition from the on-model to the off-model. Since the state point is close to the stable manifold of the offmodel Esoff (Fig. 2(b)), for a period of time, it moves toward the origin, mostly in the region Rhold , and it then moves away from the origin along Euoff . When the state point enters the region Ron at p2 (Fig. 2(c) and (d)), this triggers a transition from the off-model to the on-model. The trajectory segments from p2 to p3 and from p3 to p4 are similar to those from p0 to p1 and from p1 to p2, but they

6

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

are located closer to the origin. By repeating these transitions, the state point converges to the origin, and thus the PSA level y goes to zero. 3.3. Conditions for stability of the simplified ITTA model with the region-dividing method The saddle-point dynamics of the on-model and the off-model are critical for the stability of the simplified ITTA model with the region-dividing method. We can show that the state point diverges if the equilibrium point of either the on-model or the offmodel is an unstable node. For example, the state point monotonically diverges in the direction of Euon if the equilibrium of the off-model is an unstable node, even if the on-model has a saddlepoint equilibrium. Another critical factor for stability is the geometrical configuration of the stable and unstable manifolds. In the simplified ITTA model, Eson almost coincides with Euoff , and Esoff exactly coincides with Euon . Indeed, if the state point is moving away from the origin along Euoff during an off-treatment period, it cannot be redirected toward the origin by turning the treatment on if the direction of Eson is much different from that of Euoff . Similarly, if the state point is moving away from the origin along Euon during an on-treatment period, it cannot be redirected toward the origin by turning the treatment off if the direction of Esoff is much different from that of Euon . In this way, the topological type of the equilibrium and the geometrical configurations of the manifolds of both the on-model and off-model determine whether the region-dividing method can stabilize the origin, i.e., whether it can successfully control relapse. Based on the observations described above, we claim that the equilibrium point (the origin) of the simplified ITTA model with the region-dividing method can be stabilized if the following three conditions are satisfied: (i) Both the on-model and the off-model have saddle-point equilibrium points. (ii) The direction of the stable manifold of the on-model is “close” to that of the unstable manifold of the off-model, and the direction of the unstable manifold of the on-model is “close” to that of the stable manifold of the off-model.2 (iii) The inequality loff 4 lon holds for the slopes of the two boundary lines that determine the model transition.

represent the stable and unstable manifolds (Eson and Euon ) if the on-model has a saddle-point equilibrium; this is the case for the parameter values in Table 1, as illustrated in Fig. 2. However, we do not use this fact in the current argument, because the equilibrium point could also be an unstable node for some other parameter values of the model. Instead, we show that, for stability of the simplified ITTA model with the region-dividing method, this must be a saddle-point equilibrium. We denote the slope of the line in the direction of xð1Þ on by kon , which is defined as kon ¼

v2;1 ðaon Þ ; v1;1 ðaon Þ  v2;2 ðaon Þ

using Eq. (18). The other boundary line of Don is the xAI axis with an infinite slope. Similarly, we define another region, also in the first quadrant, referred to as Doff . The boundary lines of the region Doff pass through the origin and are in the directions of the ð2Þ eigenvectors of xð1Þ and xoff . These lines represent the unstable and off u stable manifolds (Eoff and Esoff ) if the off-model has a saddle-point equilibrium, as in the case with the parameter values in Table 1. We denote the slope of the line in the direction of xð1Þ by kon , off which is defined as koff ¼

v2;1 ðaoff Þ ; v1;1 ðaoff Þ  v2;2 ðaoff Þ

2 Note that the latter half of condition (ii) is always satisfied because the unstable manifold of the on-model and the stable manifold of the off-model are always directed along the xAI axis regardless of the parameter values of the simplified ITTA model.

ð25Þ

using Eq. (20). The other boundary line of Doff is the xAI axis. Using the regions Don and Doff in the first quadrant of the phase plane, we observe that any trajectory governed persistently by the on-model cannot cross the boundary lines of the region Don , and similarly, any trajectory governed persistently by the off-model cannot cross the boundary lines of the region Doff . This is because the equilibrium points in neither the on-model nor the off-model are foci or improper nodes, since we have observed that for both the on-model and the off-model, they are either saddle points or unstable nodes. (See Appendix B for the analytical expression of solutions of the on-model and the off-model.) Thus, any state point located in the region Don remains within Don if it is persistently governed by the on-model. Similarly, any state point located in the region Doff remains within Doff if it is persistently governed by the off-model. Thus, for a transition from the on-model to the off-model to occur, the boundary line Loff must be located within the region Don . This condition can be written as kon o loff o 1;

To prove these claims, we first clarify the conditions required for transitions between the on-model and the off-model in the simplified ITTA model. Recall that the dynamics of the simplified ITTA model cannot be stabilized unless there are repeated transitions between the on-model and the off-model; this is because either the on-model or the off-model, when used alone, is unstable. Thus, one of the necessary conditions for the stability of the simplified ITTA model can be stated as follows: any trajectory must repeatedly transverse the boundary lines, Lon and Loff , and repeatedly transition between the on-model and the off-model. In the following discussion, it is sufficient to consider the dynamics of the simplified ITTA model only in the first quadrant, i.e., xAD Z0 and xAI Z 0, since xAD and xAI represent populations of cells and must therefore be nonnegative. In order to make our argument clear, let us define a region in the first quadrant of the phase plane, referred to as Don , that is bounded by two lines passing through the origin and in the ð2Þ directions of the eigenvectors xð1Þ on and x on . These two lines

ð24Þ

ð26Þ

where the infinity after the latter inequality represents the fact that one of the boundaries of Don is along the xAI axis. Moreover, immediately after a transition from the off-model to the on-model, the state point must be located within the region Don , since, after the transition, it cannot evolve into the region Don . Since the transition from the off-model to the on-model occurs on the boundary line Lon , this condition can be written as kon o lon o 1:

ð27Þ

Similarly, the condition that a transition from the off-model to the on-model can occur can be written as koff o lon o 1;

ð28Þ

and the condition that, immediately after the transition from the on-model to the off-model, the state point is located within the region Doff can be written as koff o loff o 1:

ð29Þ

Suppose that an initial state of the simplified ITTA model is located in the region Don when the first on-treatment period begins. The state point governed by the on-model can reach the boundary line Loff if Eq. (26) holds, and this triggers the first transition from the on-model to the off-model. Immediately after

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

the first transition, the state point is located in Doff if Eq. (29) holds. The state point is then governed by the off-model, and it can reach the boundary line Lon if Eq. (28) holds. This triggers the second transition from the off-model to the on-model. Immediately after the second transition, the state point is located in Don if Eq. (27) holds, and the state point is again governed by the onmodel. In this way, we have shown that the four inequalities (26)–(29) represent the necessary conditions for repeated transitions between the on-model and the off-model. As we have already discussed, those four inequalities are also the necessary conditions for the stability of the simplified ITTA model. It remains to derive further conditions necessary for the stability of the simplified ITTA model. Suppose inequalities (26)–(29) hold. Then it is easy to show that there are repeated transitions between the on-model and the off-model, and the duration of each on-treatment period t on and that of each off-treatment period t off can be obtained as follows: t on ¼

t off

log ðloff  kon Þ  log ðlon  kon Þ ð2Þ ð1Þ on  on

λ

λ

;

ð30Þ

log ðloff  koff Þ log ðlon  koff Þ ¼ : ð1Þ ð2Þ λoff  λoff

ð31Þ

See Appendix B for the derivations of t on and t off . We observe that Eqs. (30) and (31) are independent of the state points at the beginning of the on-treatment and off-treatment periods, which implies that t on and t off are constant. Note that inequalities (26)–(29) ensure that all arguments of the logarithmic functions are positive, which means that t on and t off are finite. Thus we conclude that the transitions between the on-model and the off-model occur with the period T ¼ t on þ t off . Considering the dynamics of the simplified ITTA model within a single cycle of period T, we have the following analytical expressions for the populations of AD and AI cells:   ð1Þ ð1Þ xAD ðTÞ ¼ exp λon t on þ λoff t off xAD ð0Þ; ð32Þ 

ð1Þ

ð1Þ



xAI ðTÞ ¼ exp λon t on þ λoff t off xAI ð0Þ:

ð33Þ

Moreover, the populations of AD cells at n and ðn þ 1Þ cycles of periodic transitions can be related to each other as xAD ððn þ 1ÞTÞ ¼ xAD ðnTÞe  1=ξ ;

ð34Þ

where ξ is the damping time constant given by

ξ¼



λð1Þ lon  kon on ð1Þ log loff  kon λð2Þ on  λon



T þ

λð1Þ off ð1Þ log λð2Þ  λoff off



lon  koff loff  koff

:

ð35Þ

Eqs. (30) and (35) imply that the necessary conditions for stabilizing the origin of the simplified ITTA model with the region-dividing method are as follows: ðIÞ t on 4 0;

ð36Þ

ðIIÞ t off 4 0;

ð37Þ

ð1Þ

ð1Þ

ðIIIÞ λon t on þ λoff t off o 0:

ð38Þ

Note that inequalities (26)–(29) ensure that values of t on and t off are finite, but they do not ensure that t on and t off are positive. ð1Þ ð1Þ For these three conditions to be satisfied, λon o 0 and/or λoff o 0 must hold because t on and t off are positive. Since we have observed ð1Þ ð1Þ v1;1 ðaon Þ o v1;1 ðaoff Þ in (12), which is equivalent to λon o λoff with Eqs. (13) and (15), we see that for condition (III) in (38) to be satisfied, it is necessary that ð1Þ λon o 0:

ð39Þ

7 ð1Þ

Since condition (III) in (38) requires that λon be negative, this results in ð2Þ λð1Þ on o0 o λon ;

ð40Þ ð1Þ on

ð2Þ on

since, by definition, either λ or λ must be positive. Thus, we conclude that the on-model must have a saddle-point equilibrium in order for conditions (I), (II), and (III) to be satisfied. This argument proves half of condition (i) of our claim. Inequality (40) implies that the denominator of Eq. (30) is positive. Thus, for condition (I) in (36) to be satisfied, the numerator of Eq. (30) must be positive. We thus have the inequality loff 4 lon ;

ð41Þ

using inequalities (26)–(29). This argument proves condition (iii) of our claim. Inequality (41) implies that the numerator of Eq. (31) is also positive, by which we have the inequality ð2Þ λð1Þ off 4 λoff ;

ð42Þ

in order that condition (II) in (37) is satisfied. Since, by definition, ð1Þ ð2Þ either λoff or λoff must be positive, we conclude that

λð1Þ off 4 0:

ð43Þ ð2Þ

If we can show that λoff o0, then we can conclude that the offmodel must have a saddle-point equilibrium point in order to satisfy conditions (I)–(III). This argument proves the remaining half of condition (i) of our claim. ð2Þ In order to show that λoff is negative, we substitute Eqs. (30) and (31) into inequality (38) for condition (III), which gives the following inequality: ð2Þ ð1Þ  λð1Þ =ðλð1Þ  λð2Þ Þ    λð1Þ on =ðλon  λon Þ loff koff off off off l kon o off : lon  koff lon kon

ð44Þ

In Appendix C, we show that, for inequality (44) to be satisfied, the following is necessary:

λð2Þ off o

ð2Þ λð1Þ off λon : ð1Þ λon

ð45Þ

We observe that the right-hand side of (45) is negative, since ð2Þ ð1Þ λð1Þ on o0, and λon 40 and λoff 40 by (40) and (43). Thus, inequality (44) implies that

λð2Þ off o 0:

ð46Þ

Thus, we conclude that the off-model must have a saddle-point equilibrium point in order for conditions (I)–(III) to be satisfied. This proves the remaining half of condition (i) of our claim. Condition (ii) remains to be proven. It claims that, for the simplified ITTA model to be stable, the direction of the stable manifold of the on-model must be “close” to that of the unstable manifold of the off-model. Condition (ii) can be rephrased to state that the values of kon and koff must be close to each other. Unfortunately, we are unable to provide a comprehensive proof of condition (ii) even for this simple model, although the validity of condition (ii) can be confirmed by numerical exploration of the simplified ITTA model. Thus, the necessity of condition (ii) for the stability of the simplified ITTA model is only a conjecture. Nevertheless, we have argued, at the end of Appendix C, that some inequalities required for condition (i) can be satisfied more easily, i.e., in a robust way against the choice of the parameter values of the model, if kon and koff are both small and close to zero. 3.4. Region-division dependency of stability and dynamics If we select values of the region-dividing parameters so that lon oloff holds, the region-dividing method may be able to successfully

8

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

ITTA model, the transitions between the on-model and the offmodel are triggered when the state point enters the neighborhood of a stable manifold. The region-dividing method requires the current location of the state point in the phase space of the system, not just the PSA value of y. However, since the state of prostate cancer cells cannot be measured directly in a clinical examination, we introduce a state observer that estimates the true state of ðxAD ; xAI Þ as the estimated state ðx~ AD ; x~ AI Þ from the measured PSA value of y (Fig. 1(c)). In this study, for the state observer, we employ the simplified ITTA model with a feedback term proportional to the difference between the actual serum PSA concentration y and the estimated serum PSA ~ The state observer is implemented as follows: concentration y. ! ! x~ AD ðtÞ d x~ AD ðtÞ ð47Þ ¼ Ak þKk ðy  y~ Þ; x~ AI ðtÞ dt x~ AI ðtÞ ~ ¼ c1 x~ AD ðtÞ þc2 x~ AI ðtÞ; yðtÞ

Fig. 3. (a) Relationship between the transition period T ¼ t on þ t off and the regiondividing parameters (lon and loff ) in the simplified version of the modified ITTA model. (b) Relationship between the damping time constant, Eq. (35), and the region-dividing parameters (lon and loff ). The equilibrium point of the simplified ITTA model is unstable if the set of ðlon ; loff Þ is located in the shaded triangular area defined by loff o lon , and it is stable otherwise. See the text.

control the simplified ITTA model. Convergence strongly depends on the values of the region-dividing parameters. Fig. 3(a) shows how the period T of the on–off transitions depends on the values of the regiondividing parameters lon and loff (Eqs. (30) and (31)). The closer are the values of lon and loff , the shorter is the period T. If the value of loff is set too close to lon but slightly larger than lon , then the region Rhold that separates the neighborhood of the stable manifold of the on-model and that of the off-model becomes narrow. This situation may cause a chattering-like motion of the state point, accompanied by very frequent transitions between the on-model and the off-model. On the other hand, if lon and loff are very different, then the period T becomes large, leading to less frequent transitions. The relationship between the damping time constant and the region-dividing parameters can also be clarified analytically, as we have argued above (Eq. (35)). Fig. 3(b) shows that the damping time constant becomes small (i.e., the convergence to the origin is rapid) as the values of lon and loff increase. However, there is little change in the time until convergence, so there is no need to tune the values of lon and loff . In summary, the values of the region-dividing parameters do not much affect the convergence of the simplified ITTA model, but they do affect the frequency of transitions between the on-model and the off-model. Indeed, the values of loff and lon in Table 1, which satisfy the condition lon o loff , do not result in frequent transitions.

4. Models of IAS therapy with the region-dividing method In this section, we modify the ITTA and HBA models by adding the region-dividing method; the modified ITTA model and the modified HBA model are shown in Fig. 1(c) and (d), respectively. 4.1. Modified ITTA model The modified ITTA model with the region-dividing method is defined using Eqs. (1)–(4). As illustrated above for the simplified

ð48Þ

where Kk for k A fon; offg is the observer gain matrix (see Appendix D for details). Note that the value of the simplified androgen concentra~ cannot be exactly the tion used in the state observer, denoted as a, same as the true value a in Eqs. (1) and (2), because the androgen dynamics used in the state observer has been linearized for simplicity. The difference between a~ and a induces an error between ðxAD ; xAI Þ and ðx~ AD ; x~ AI Þ, and the error can be large during the transient period after each model transition. To avoid the use of inaccurate estimates of ðx~ AD ; x~ AI Þ during the transient periods, after a transition, we introduce a dead time of τ days during which no additional transitions are triggered, even if the estimated state point newly enters the region Roff or Ron from the region Rhold . Thus, each transition from the off-model to the on-model takes place only when the estimated state ðx~ AD ; x~ AI Þ enters the region Ron and τ days have passed since the last transition. Similarly, each transition from the on-model to the off-model takes place only when the estimated state enters the region Roff and τ days have passed since the last transition from the off-model to the onmodel. The parameter values used for the modified ITTA model are shown in Table 1, and a block diagram of the modified ITTA model is shown in Fig. 1(c). 4.2. Modified HBA model The modified HBA model with the region-dividing method is defined using Eqs. (5)–(7). For the parameter values in Table 2, the system matrix of the on-model Won in Eq. (5) possesses two negative and one positive real eigenvalues. For the former, the corresponding eigenvectors span the two-dimensional stable manifold Eson , which is almost identical to the xAD–xrAI plane. For the latter, the corresponding eigenvector spans the xiAI-axis as the one-dimensional unstable manifold, Euon . The system matrix of the off-model Woff in Eq. (6) possesses two positive and one negative real eigenvalues. For the former, the corresponding eigenvectors span the two-dimensional unstable manifold Euoff , which is identical to the xAD–xrAI plane. For the latter, the corresponding eigenvector spans the xiAI-axis as the one-dimensional stable manifold, Esoff . Thus it is expected that the modified HBA model can be successfully controlled, because both the on-model and the off-model have saddle-point equilibrium points, and the direction of the stable manifold of one model is almost identical to the unstable manifold of the other model, i.e., Eson almost coincides with Euoff , and Esoff with Euon . We define a neighborhood region Ron of the stable manifold Eson of the on-model by its boundary with the conical surface Son , and a neighborhood region Roff of the stable manifold Esoff of the offmodel by its boundary with another conical surface Soff . More

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

specifically, Son and Soff are formulated as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Son : xiAI ¼ son ðxAD Þ2 þ ðxrAI Þ2 ; Soff : xiAI ¼ soff

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxAD Þ2 þ ðxrAI Þ2 ;

9

ð49Þ ð50Þ

where son and soff are the parameters that determine the slopes of the cones. We implement a state observer in the modified HBA model in order to obtain the state estimate ðx~ AD ; x~ rAI ; x~ iAI Þ for the true state ðxAD ; xrAI ; xiAI Þ from the PSA value of y, as follows: 0 1 0 1 x~ AD ðtÞ x~ AD ðtÞ dB ~ C B C ð51Þ @ x rAI ðtÞ A ¼ Wk @ x~ rAI ðtÞ A þKk ðy  y~ Þ; dt x~ iAI ðtÞ x~ iAI ðtÞ ~ ¼ c1 x~ AD ðtÞ þ c2 x~ rAI ðtÞ þ c3 x~ iAI ðtÞ; yðtÞ

ð52Þ

where Kk for k A fon; off g is the observer gain matrix. As in the modified ITTA model, each transition from the offmodel to the on-model takes place when the estimated state ðx~ AD ; x~ rAI ; x~ iAI Þ enters the region Ron and τ days have passed since the last transition. Similarly, each transition from the on-model to the off-model takes place when the estimated state point enters the region Roff and τ days have passed since the last transition. Note that the holding period in the original HBA model is not used in the modified HBA model because the dead time is equivalent to the holding time. The parameter values in the modified HBA model are shown in Table 2, and a block diagram of the modified HBA model is shown in Fig. 1(d).

5. Dynamics of the IAS therapy models In this section, we show that the region-dividing method is more effective than the threshold method at controlling the population of prostate cancer cells. To this end, we perform numerical simulations of the models using the forward Euler method with a time step of 1 day.3 We then compare the dynamics of the original ITTA model and HBA model with that of their modified models.

Fig. 4. Time series of the PSA concentration, treatment status (on/off), and trajectories for the original and the modified ITTA models. (a) PSA time series, y, in the ITTA model. The PSA time series in the CAS therapy model is superimposed (the gray dot-dashed curve). The dashed and dotted horizontal lines represent the upper and lower thresholds that determine the on–off transitions. (b) PSA time series for the modified ITTA model. The lower traces of (a) and (b) indicate times on and off of the treatment. The vertical gray bands represent the series of dead times. (c) and (d) are the phase-plane trajectories corresponding to the time series shown in the rectangular boxes in (a) and (b), respectively. In (d), the solid gray curve represents a trajectory of the true state ðxAD ; xAI Þ. The black solid and dot-dashed curves represent the corresponding trajectories of the estimated state ðx~ AD ; x~ AI Þ, where the latter curves represent the trajectories during the period of dead time. p0 is an initial state. p1 and p3 are the state points at the transitions from the on-model to the off-model. p2 and p4 are the state points at the transitions from the off-model to the on-model.

5.1. Dynamics of the ITTA model and the modified ITTA model Fig. 4 shows the typical dynamics of the ITTA model and the modified ITTA model. In Fig. 4(a), the time course of the PSA level, obtained from a simulation of a CAS therapy model, is also shown as a reference where the PSA level decreases at the beginning of the treatment, followed by an increase during a later phase of the treatment, which represents a relapse. The lower panels of Fig. 4(a) and (b) show transition patterns between the on-model and the off-model, and we can identify periodic cycles in the transitions. In Fig. 4(a), the y-value of the ITTA model shows a sustained oscillation, representing the limit cycle dynamics accompanied by the periodic model transitions in the steady state. In contrast, in Fig. 4(b), the y-value of the modified ITTA model with the region-dividing method shows a non-oscillatory and rapid convergence toward zero, corresponding to extermination of the tumor cells. However, despite the non-oscillatory dynamics in the modified ITTA model, the model showed periodic transitions. 3 In the simulations of the ITTA model and the modified ITTA model, the initial conditions are set as ðað0Þ; xAD ð0Þ; xAI ð0ÞÞ ¼ ð30:0; 15:0; 0:01Þ according to Ideta et al. (2008), and ðx~ AD ð0Þ; x~ AI ð0ÞÞ ¼ ð0:0; 0:0Þ. In the HBA model and the modified HBA model, the initial conditions are set as ðxAD ð0Þ; xrAI ð0Þ; xiAI ð0ÞÞ ¼ ð17:8395; 0:77803; 0:382485Þ according to Hirata et al. (2010), and ðx~ AD ð0Þ; x~ rAI ð0Þ; x~ iAI ð0ÞÞ ¼ ð0:0; 0:0; 0:0Þ.

The transition period of the modified ITTA model is shorter than that of the ITTA model with the threshold method. Indeed, the y-value of the modified ITTA model also oscillates synchronously with the transition cycle, although it is hard to see this in the figure because of the small amplitude of the oscillation and its rapid decay. Fig. 4(c) and (d) shows the steady-state limit cycle solution of the ITTA model and the convergence of the modified ITTA model, respectively, corresponding to Fig. 4(a) and (b). In Fig. 4(d), the true state point ðxAD ; xAI Þ exhibits a zigzag-shaped curve as in the simplified ITTA model (Fig. 2(d)). This demonstrates that the cancer cells are successfully exterminated by the region-dividing method with the modified ITTA model. It also shows that the dead time introduces a delay in the model transitions and that this is beneficial. That is, the estimated state ðx~ AD ; x~ AI Þ deviates from the true state in the periods of dead time after each model transition; by the end of this time, the estimation error has become small enough that we can use the estimated state for determining the subsequent transitions. More specifically, at the point p0, ðx~ AD ; x~ AI Þ is estimated to be almost identical to the true state ðxAD ; xAI Þ. In the first on-treatment period from p0 to p1, the estimated state point remains close to the true state point. When the estimated state

10

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

point enters the region Roff at p1, a transition from the on-model to the off-model is triggered. In the off-treatment period from p1 to p2, the true state point shows a hyperbolic trajectory. These dynamics is quite similar to that observed in the simplified ITTA model with the region-dividing method, as shown in Fig. 2(d). However, the corresponding estimated state point moves differently from the true state point. In particular, the estimated state point moves away from the true state point during the dead time (from p1 to pa, as shown by the black dot-dashed curve), but it then approaches the true state point (from pa to p2, the black solid curve). When the estimated state point enters the region Ron at p2, the dead time τ has already passed since the last transition at p1, and so the next transition from the off-model to the on-model is triggered. In the second on-treatment period from p2 to p3, the true state point approaches the origin for a period of time, and it then moves away from the origin. Once again, the corresponding estimated state point moves in a complicated manner away from the true state point during the dead time (from p2 to pb, the black dot-dashed curve), but then it approaches the true state point (from pb to p3, the black solid curve). Note that, in this case, the estimated state point enters the region Roff during the dead time. However, this event does not trigger a transition from the onmodel to the off-model, because the dead time τ has not passed yet. After the estimated state point becomes close to the true state, when the dead time τ has already passed since the last transition at p2, the estimated state enters again into the region Roff at p3. This triggers a transition from the on-model to the off-model. The segment from p3 to p4 and the successive segments repeat similar dynamics to that from p1 to p2 and from p2 to p3, and the state point becomes closer to the origin. The dynamics described here suggests that the dead time introduced in the modified ITTA model plays an indispensable role in avoiding chattering-like very frequent switches, and leads instead to an appropriate sequence of switches.

5.2. Dynamics of the HBA model and the modified HBA model Fig. 5 shows a typical time series for PSA, transition patterns of the drug administration status, and trajectories in the phase space, for the original and the modified HBA model. As in the original and the modified ITTA models, the original and the modified HBA models exhibit periodic transitions, as shown in Fig. 5(a) and (b). The periodicity of the modified HBA model is not apparent in the figure because there is such a long period between transitions. The transitions between p1 and p4 are repeated in the modified HBA model. The y-value for the HBA model (Fig. 5(a)) shows a sustained oscillation, which represents the limit cycle dynamics accompanied by periodic model transitions. In Fig. 5(b), the y-value for the modified HBA model shows a rapid convergence to zero. As in the modified ITTA model, the y-value for the modified HBA model oscillates synchronously with the transition cycle, although it is also hard to recognize from the figure because of the rapid decay and the long period of the transitions. Fig. 5(c) shows the transient dynamics and the steady-state limit cycle solution of the HBA model. These correspond to the sustained oscillation of y shown in Fig. 5(a). Fig. 5(d) shows the dynamics of the modified HBA model, in which the true state point ðxAD ; xrAI ; xiAI Þ and the estimated state ðx~ AD ; x~ rAI ; x~ iAI Þ have almost identical behavior. By repeatedly transitioning between the onmodel and the off-model, they converge to the origin, which leads to a successful treatment. The segments of the gray curve from p0 to p1 and from p3 to p4 are the trajectories in the presence of treatment (the on-model), whereas that from p1 to p3 is the trajectory in the absence of treatment (the off-model). In the first on-treatment period from p0 to p1, the state point moves along the

Fig. 5. Time series of the PSA concentration, treatment status (on/off), and trajectories for the original and the modified HBA models. (a) Time course of the PSA concentration y in the HBA model. The PSA time series in the CAS therapy model is superimposed (the gray dot-dashed curve). (b) PSA time series in the modified HBA model. See the caption of Fig. 4(a) and (b). (c) and (d) Phase-space trajectories corresponding to the time series shown in the rectangular boxes in (a) and (b), respectively. (d) Trajectories of the true state and the corresponding estimated state. See the caption of Fig. 4(d). p0 is an initial state. p1 and p4 are the state points at the transitions from the on-model to the off-model. p2 is the state point when it traverses the conical surface Son . p3 is the state point at the transition from the off-model to the on-model.

flow of the saddle-point vector field of the on-model. In this area, the populations of AD cells and rAI cells decrease, but that of the iAI cells gradually increases. In this period of time, the estimated state point closely follows the true state point, and then the estimated state point enters the region Roff , which triggers a transition from the on-model to the off-model. In the successive off-treatment period, since the true state point is located close to the stable manifold of the off-model Esoff (xiAI-axis), it moves toward the origin for a period of time, and then moves away from the origin along the unstable manifold of the off model Euoff (xAD–xrAI plane). The estimated state point follows this motion. Then the estimated state point enters the region Ron at p2. However, a transition from the off-model to the on-model is not triggered at this moment, because the dead time τ has not elapsed since the previous transition. When the true state point and the estimated state point arrive at p3, the dead time τ has passed, and the transition is triggered. Since the true state point is close to the stable manifold of the on-model Eson (the xAD–xrAI plane), it moves toward the origin for a period of time, and then moves away from the origin along the unstable manifold of the on-model Euon (the xiAI-axis). Then the true and estimated state points enter the region Roff at p4, triggering a transition from the on-model to the off-model. Similar model transitions from p1 to p4 are repeated, but with attenuated amplitude. In this way, the HBA model with the region-dividing method can successfully control the three-dimensional population dynamics of prostate cancer cells.

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

One may wonder why, in the modified HBA model but not in the modified ITTA model, the estimated state shown in Fig. 5(d) closely follows the true state. This happens because the state observer dynamics in the modified HBA model is exactly the same as the HBA model with no uncertainty, while that in the modified ITTA model employs simplified androgen dynamics. This means that the dead time τ in the modified HBA model can be much shorter than the value in Table 2. Indeed, if we reduce the dead time τ for the modified HBA model, the overshoot trajectory from p2 to p3 in Fig. 5(d) becomes small, and the region-dividing method can control the population dynamics more effectively (figures not shown).

6. Robustness of the control to variations in the parameters So far, we have considered the IAS therapy models with the typical parameter values shown in Tables 1 and 2. However, clinical practice has shown that although IAS therapy is effective for some patients, it is not effective for others (Bruchovsky et al., 2006, 2007). This indicates that the individual characteristics of the tumor cells influence the efficacy of the therapy. Thus, it is important to analyze how the efficacy of IAS therapy changes due to variations of the individual characteristics of tumor cells. To this end, we analyzed the robustness of IAS therapy to changes in several parameter values that may critically affect the stability of the models. Moreover, we clarified the relationship between the ability of the region-dividing method to control the cancer cells and the topological type of the equilibrium points of the IAS therapy models. These will be discussed below. Ideta et al. (2008) and Hirata et al. (2010) discussed the influence on the dynamics of their models of the maximum change rate from AD cells to AI cells in the presence of treatment and the net growth rate of AI cells in the absence of treatment. These parameters determine if it is possible for the threshold method to prevent a relapse. They showed that a relapse cannot be avoided if the growth rate of the AI cells in the absence of treatment is positive, i.e., if AI cells increase in the absence of treatment. In this case, the population of the AI cells increases both in the presence and in the absence of treatment, and thus the population of the AI cells inevitably diverges. Moreover, they showed that, even if the population of AI cells decreases in the absence of treatment, a relapse cannot be avoided if the mutation rate of AD to AI cells in the presence of treatment is large and the apoptosis rate of AI cells in the absence of treatment is small. These parameter values are also critical for the region-dividing method, since the topological type of the equilibrium points of the IAS therapy models may change depending on the values of these parameters. Therefore, we performed numerical simulations for various sets of parameter values, and we examined the stability of the various modified models. We judged that a model was stable if the value of y did not exceed 50.0 ng/ml for 10,000 days after the beginning of the treatment. 6.1. Stability of the original and the modified ITTA models As in Ideta et al. (2008), we examined how the maximum mutation rate from AD to AI cells, m1, and the strength of the decrease in the net growth rate of AI cells during the off-treatment period, d, affected the stability of the modified ITTA model (see Appendix A). As mentioned above, the larger the value of m1 is, the more actively the AD cells mutate into AI cells in the presence of treatment. The population of AI cells continues growing if d is small enough, whereas it is relatively easy to reduce the population of AI cells if d is large enough. In particular, if d 4 0:305, the equilibrium points of both the on-model and the off-model are saddle points. Otherwise, the equilibrium point of the off-model is

11

an unstable node. We set the ranges of m1 and d as m1 A ½0; 1:0  10  4  and d A ½0; 1:0, and analyzed the stability of the original and the modified ITTA models against changes in these parameter values. Fig. 6(a) and (b) shows the robustness of the original and the modified ITTA model, respectively, against changes in the values of m1 and d. The dynamics was bounded, namely asymptotically stable when the values of m1 and d were located within the blue region. Comparing Fig. 6(a) and (b), we can confirm that the region-dividing method results in stabilized dynamics in a much wider region, i.e., in a more robust way, than the threshold method. In other words, the modified ITTA model is stable even if m1 takes a value large enough that it makes the original ITTA model unstable. The robustness of the modified ITTA model against changes in d is almost the same as that of the original ITTA model. Both the original and the modified ITTA models are unstable if d is small, regardless of the value of m1. 6.2. Stability of the original and the modified HBA models on off We examined how the values of won 3;1 , w3;2 , and w3;3 altered the stability of the original and the modified HBA models. The larger the won 3;1 is, the more actively the AD cells change to iAI cells in the presence of treatment. Similarly, the larger the won 3;2 is, the more actively the rAI cells change to iAI cells in the presence of treatment. If woff 3;3 is positive, then the iAI cells increase in the absence of treatment; otherwise, they decrease. In particular, if woff 3;3 o 0, then the equilibrium points of both the on-model and the off-model are saddle points. Otherwise, the equilibrium point of the off-model is an unstable node. In this way, the parameter on off values of won 3;1 , w3;2 , and w3;3 are critical for the stability of both the original and the modified HBA models. We set biologically plauon sible ranges of those parameters, won 3;1 A ½0; 0:01, w3;2 A ½0; 0:1, and

woff 3;3 A ½  1; 1, and analyzed the resulting stability for both HBA models. Fig. 6(c) and (d) shows, respectively, the stability regions of the original and the modified HBA models in the parameter space of on off won 3;1 , w3;2 , and w3;3 . The state point is bounded or asymptotically on off stable if the set of parameter values ðwon 3;1 ; w3;2 ; w3;3 Þ is located within the blue region. Comparing Fig. 6(c) and (d), we can confirm that the region-dividing method successfully stabilizes the IAS dynamics in a much wider region, i.e., in a much more robust way, than does the threshold method.

6.3. Relationship between the control capability of the regiondividing method and the types of equilibria Based on our assessment of the necessary conditions for successful control of PSA with the region-dividing method, we analyzed the relationship between the robust control capability of the method and the topological type of the equilibrium points of the on- and off-models. Since the uncertainty involved in the state estimation and the presence of a dead time both affect the performance, we considered simplified cases in which the state of the prostate cancers could be directly and accurately measured, and thus neither the state observer nor the dead time τ was required. We first performed a stability analysis as we did in Sections 6.1 and 6.2, but this time for the modified ITTA and HBA models with the direct measurement of the state. In Fig. 6(b), the stable parameter region of the modified ITTA model with the direct measurement of the state is superimposed on the stability region of the modified ITTA model with the estimated state and the

12

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

Fig. 6. Stability regions of the four IAS models in the parameter spaces. In each panel, blue coloration indicates the stability region of the model. (a) ITTA model with the threshold method. (b) Modified ITTA model with the region-dividing method. (c) HBA model with the threshold method. (d) Modified HBA model with the region-dividing method. In each of (b) and (d), the black solid rectangle represents the parameter region in which the equilibrium point of the on-model and that of the off-model are both saddle points. The red dashed rectangle represents the stability region of the modified IAS model with the direct state measurement and with no dead time. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

inclusion of dead time. This shows that the stability regions of both cases are almost the same. We then obtained the parameter region in which both the on- and off-models have saddle-point equilibria. Comparing these two parameter regions, we confirmed that the stability region of the modified ITTA model is almost identical with (and included in) the parameter region with the saddle-point equilibria. This confirmation suggests that a necessary condition for successful stabilization of the model with the region-dividing method is that the equilibria of both the on- and off-models are saddle points and that the directions of the stable and unstable manifolds are close to each other. Moreover, the fact that the two regions are almost the same implies that this condition can also be considered to be approximately a necessary and sufficient condition for stability, providing that the region is divided appropriately. Fig. 6(d) shows the stability region of the modified HBA model with the direct state measurement and the parameter region with the saddle-point equilibrium for both of the on-model and the offmodel. We confirmed first that the stability region of the modified HBA model with the direct state measurement is slightly wider than that of the modified HBA model with the estimated state and the inclusion of the dead time. We then confirmed that the stability region of the modified HBA model with the direct state measurement is almost the same as the parameter region in which there are saddle-point equilibria. This confirmation means that the region-dividing method can stabilize the system as long as both the on-model and the off-model have saddle-point equilibria, providing that the directions of the stable and unstable manifolds are also close to each other; these are the same conditions as in the modified ITTA model. Note that the differences in the stability regions of the modified HBA model with and without the direct state measurement is due to the inclusion of the dead time τ. This can be understood by the fact that the estimated state point of the modified HBA model closely follows the true state point of the modified HBA model, as shown in Fig. 5(d). Thus, the times at which the model transitions are triggered by the estimated state are almost the same as those triggered by the true state. This coincidence means that the

robustness of the modified HBA model is further improved by the appropriate selection of the dead time τ.

7. Discussion In this study, we considered IAS therapy of prostate cancer, and we proposed the region-dividing method as an alternative to the conventional threshold method. Based on the theoretical analysis described in Section 3 and the original ITTA and HBA models that were proposed as mathematical models of IAS therapy with the threshold method, we constructed modified ITTA and HBA models that use the region-dividing method. By comparing the typical dynamics of the original and the modified ITTA models, and the original and the modified HBA models, we demonstrated that the region-dividing method is more effective than the threshold method for controlling the population dynamics of prostate cancer cells. Moreover, we compared the two control methods for robustness against changes in the parameter values that characterize the proliferation, degradation, and mutation rates of the cancer cells. We showed that the region-dividing method stabilized the system's dynamics in a much more robust way than did the threshold method. We clarified that the region-dividing method could lead to successful control of prostate cancer if the equilibrium points of both the on-model and the off-model are saddle points and if the directions (or the geometrical structure) of the stable manifold of one model (either on- or off-) and the unstable manifold of the other model (either off- or on-) are close to each other. Although the division of the regions should be made appropriately, the performance does not greatly depend on the selection of the regions. One may claim that the region-dividing method is identical with the sliding mode control (Utkin, 1977; Young et al., 1999), because the tumor dynamics controlled by the region-dividing method looks similar to typical dynamics of the sliding mode control. More specifically, in each method, the state point is confined to a subspace (a sliding surface in the sliding mode control and a stable manifold exhibited by the system with no

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

control signals in the region-dividing method) and it slides to the origin along the space. However, the region-dividing method differs essentially from the sliding mode control: the sliding mode control deals with a non-hybrid system that is not necessarily accompanied with a stable manifold when no control signal is applied (Utkin, 1977; Young et al., 1999), whereas the region-dividing method proposed in this study aims at stabilizing a hybrid system that consists of multiple (two in this study) sub-systems and at least one of the sub-systems exhibits saddle-point dynamics accompanied with a stable manifold. In the sliding mode control, it is required to define (or design) an appropriate sliding surface that includes the origin (the desired equilibrium state point), and then the vector field on the sliding surface is modified by an "optimal" control signal so that the dynamics on the sliding surface becomes asymptotically stable (Utkin, 1977; Young et al., 1999). Moreover, the controller should also be responsible for trapping the state point in a neighborhood of the sliding surface by repeatedly switching on and off another feedback control input according to the state point. Thus, the vector field of the original non-hybrid system and that of the controlled system might be largely different due to those two types of control signals. The designing those two types of control signals requires burdens of parameter tuning and also makes the sliding mode control less-robust due to possible instability of the tuned dynamics. On the other hand, the region-dividing method is applicable to a hybrid system whose sub-system(s) itself, with no modification of the vector field, possesses a stable manifold, meaning that there is neither need to design a sliding surface nor to generate optimal control signals to stabilize the dynamics on the stable manifold. Instead, the major issue to stabilize the hybrid system is just to design a neighborhood of the manifold that triggers switching among sub-systems. Conditions required for the neighborhood are simple. That is, a boundary of the neighborhood, and thus also the stable manifold, should be a global Poincaré section, which enables any state point to reenter the neighborhood of the stable manifold after a preceding transition from one sub-system to the other. This simple procedure for establishing the control in the region-diving method provides a basis of robust stability of the method. In this way, although systems controlled by both methods exhibit switching dynamics, two methods are not the same. Besides, in the sliding mode control, the chattering problems often arise due to strong and frequent input for trapping the state point in the neighborhood of sliding surface (Boiko and Fridman, 2005; Levant, 2010). In contrast, in the region-dividing method, the chattering problem rarely arises. See Figs. 3–5 to confirm moderately separated model transitions. Difference in the chattering characteristics between two methods can be understood by the discussion we made above. The sliding mode control is often considered as a stiff control, while the region-dividing method as a flexible control. See also references (Bottaro et al., 2005, 2008; Asai et al., 2009; Suzuki et al., 2012) that argue related issues on the intermittent control for hybrid systems with a stable manifold of saddle-point dynamics for stabilizing human movement. There are several issues to be addressed before employing the region-dividing method in clinical practice. The major issue is the estimation of the state from the measured values of the concentration of PSA. In the model simulations, it is easy to observe the PSA concentration, and the state observer can continuously estimate the full state of the system. In clinical practice, however, continuous or frequent monitoring of PSA requires a blood test and is laborious for patients, and it thus might lower the quality of life (QoL) of the patients. If we cannot use a sequence of frequent sampled PSA concentrations to appropriately estimate the state variables for the AD and AI cells, the region-dividing method cannot accurately control the tumor dynamics. Fortunately, there are several ways to resolve this difficulty; they include the use of an

13

advanced method for constructing a state observer that consists of the state estimation and modeling dynamics of the PSA concentration. Suzuki and Aihara (2013) have shown that the dynamics of the PSA concentration can be modeled from a sequence of discrete observations by using either a variational Bayesian method or a Gaussian process regression method for a nonlinear auto-regressive eXogeneous model. This may allow us to combine the estimated PSA values and the state observer defined in Section 3 to obtain accurate estimates of the number of AD and AI cells from a sequence of sparse samples of the PSA concentration. The second issue is how we construct a state observer that is appropriately characterized by patient-dependent parameters. If the estimated state of the tumor dynamics that is obtained by a state observer deviates greatly from the true state, the transitions between on-treatment and off-treatment cannot be performed appropriately. To address this problem, the model identification (parameter fitting) should be performed as accurately as possible. Hirata et al. (2010) showed that, for IAS therapy with the threshold method, the parameter values can be obtained using the PSA data for 1.5 cycles of treatment. Based on this report, we propose the following protocol, which is composed of two phases, for IAS therapy with the region-dividing method: (1) parameter estimation phase: obtain a state observer by using the parameter estimation method of Hirata et al. (2010); (2) treatment phase: use the state observer with the region-dividing method to determine the treatment. There is, however, a trade-off for clinical applications of the proposed method. On one hand, the parameter estimation phase should be completed as quickly as possible in order to reduce the inconvenience for the patients, but, on the other hand, it needs to take enough time to obtain an accurate state observer. One of the advantages of IAS therapy is its ability to improve the QoL of patients during the periods between treatments (Gleave et al., 1998; Crook et al., 1999; Mottet et al., 2012; Salonen et al., 2013; Mottet, 2013). Expected benefits include reductions of the adverse effects, such as hot flushes, decreased libido, and bone fractures. Several clinical studies of Phase II analysis showed that about 5–10 months of off-treatment were required for maximum recovery of QoL (Spry et al., 2006; Mearini et al., 2011; Ng et al., 2012). It is, of course, the patients’ preference that the offtreatment periods last for as long a time as possible. In this regard, the region-dividing method might be better than the threshold method, since an appropriate selection of the region-dividing parameter values might be able to increase the intervals between treatments, as shown in Fig. 3, while leaving the time until convergence almost unchanged. This suggests that we need to select the values of the region-dividing parameters with particular attention to increasing the intervals between treatments rather than increasing the rate of convergence. As shown in Fig. 3, if both lon and loff have moderate values and the value of lon is close to that of loff , the transition frequency is high, leading to short off-treatment periods that impose an enormous drain on the patients. This is not desirable in clinical practice. On the other hand, if lon is small and loff is large, the transition frequency can be low, leading to long off-treatment periods. However, it should be noted that if lon is small, xAD will transiently take a large value. Similarly, if loff is large, xAI will take a large value. These large values of xAD and xAI, even if transient, may not be desirable in clinical practice, since it cannot necessarily be guaranteed that the large populations of AD or AI cells can be eventually controlled and reduced in real situations. This is another trade-off between the frequency of intervention and the acceptable maximum population of tumor cells. Based on the above discussion, we propose another protocol, also in two phases: a treatment phase using the region-dividing method, followed by a long off-treatment phase. There is an

14

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

advantage in using the region-dividing method to provide a long off-treatment period, since the population of prostate cancer cells can be reduced to close to zero, as we demonstrated in this study. Thus, once the population of prostate cancer cells has been reduced to a small value, it might be a long time before it is necessary to again administer drugs, since it may take a long time for the cancer cells to regrow. After this long off-treatment period, once the PSA concentration exceeds a certain threshold, a new treatment phase will begin. By repeating these two phases, the population of the prostate cancer cells can be bounded. We expect that the long off-treatment period provided by this protocol will improve the QoL of patients. In this paper, we considered IAS therapy of prostate cancer. Similar intermittent treatments can be applied to other cancers, e.g., breast cancer (Traina et al., 2008), colorectal cancer (Scheithauer et al., 2002), and pancreatic cancer (Scheithauer et al., 2003). However, to the best of our knowledge, there are few mathematical studies on those cancers. Since applications of mathematical modeling and control engineering can improve clinical diagnosis and treatments, widespread and deeper analysis and the practical realization of such techniques are required.

The parameter a0 is the normal androgen concentration at the steady state in the untreated condition.

Appendix B. Dynamics of the simplified ITTA model B.1. Analytical solutions of the on-model and the off-model A solution of the simplified ITTA model in the presence of treatment (the on-model) can be described as follows: ð1Þ

xAD ðtÞ ¼ xAD ð0Þeλon t ;

ðB:1Þ ð1Þ

ð2Þ

ð2Þ

xAI ðtÞ ¼ xAD ð0Þkon ðeλon t  eλon t Þ þ xAI ð0Þeλon t ;

ðB:2Þ

where kon is defined in Eq. (24) in the main text. A solution of the simplified ITTA model in the absence of treatment (the off-model) can be described as follows: ð1Þ

xAD ðtÞ ¼ xAD ð0Þeλoff t ;

ðB:3Þ ð1Þ

ð2Þ

ð2Þ

xAI ðtÞ ¼ xAD ð0Þkoff ðeλoff t  eλoff t Þ þ xAI ð0Þeλoff t ;

ðB:4Þ

where koff is defined in Eq. (25) in the main text. Acknowledgments

B.2. Durations of t on and t off

This work was supported in part by JSPS Grants-aid 24800038 and 23300166. The work of Y.H. and K.A. was supported by the Aihara Innovative Mathematical Modelling Project, the Japan Society for the Promotion of Science (JSPS) through its ’Funding Program for World-Leading Innovative R&D on Science and Technology (FIRST Program)’, initiated by the Council for Science and Technology Policy (CSTP).

Here, we assume that inequalities (26)– (29) hold, and thus the transitions between the on-model and the off-model may occur repeatedly. Suppose that the initial state point ðxAD ð0Þ; xAI ð0ÞÞ is on the boundary line Lon and the treatment begins at t¼ 0. The first transition from the on-model to the off-model occurs when the state point reaches Loff at t ¼ t on . Using the equalities xAI ð0Þ ¼ lon xAD ð0Þ and xAI ðt on Þ ¼ loff xAD ðt on Þ, which are satisfied by the state points placed on the boundaries, and the analytical solution of the on-model given by Eq. (B.2), we have

Appendix A. Definitions of static functions of the androgen concentration in the ITTA model In the ITTA model, we define v1;1 ðaðtÞÞ, v2;1 ðaðtÞÞ, and v2;2 ðaðtÞÞ as follows: v1;1 ðaðtÞÞ ¼ α1 p1 ðaðtÞÞ  β1 q1 ðaðtÞÞ mðaðtÞÞ;

ðA:1Þ

v2;1 ðaðtÞÞ ¼ mðaðtÞÞ;

ðA:2Þ

v2;2 ðaðtÞÞ ¼ α2 p2 ðaðtÞÞ  β2 q2 ðaðtÞÞ;

ðA:3Þ

where a(t) represents the androgen concentration. α1 p1 and α2 p2 represent the proliferation rates of AD and AI cells, respectively. β1 q1 and β2 q2 represent the apoptosis rates of AD and AI cells, respectively. m is the mutation rate of AD into AI, which depends on the androgen concentration. m, p1, q1, p2, and q2 are also static functions of the androgen concentration, as follows:   aðtÞ mðaðtÞÞ ¼ m1 1  ; ðA:4Þ a0 p1 ðaðtÞÞ ¼ k1 þ ð1  k1 Þ

aðtÞ ; aðtÞ þ k2

ðA:5Þ

ð1Þ

ð2Þ

ð2Þ

p2 ðaðtÞÞ ¼ 1  d q2 ðaðtÞÞ ¼ 1;

aðtÞ ; a0

ðA:6Þ

ðB:5Þ

Solving this equation with respect to t on yields Eq. (30). Similarly, suppose that the state point ðxAD ðt on Þ; xAI ðt on ÞÞ is on the boundary line Loff and the off-treatment period begins at t ¼ t on . The second transition from the off-model to the on-model occurs when the state point reaches Lon at t ¼ t on þ t off . Using the equalities xAI ðt on Þ ¼ loff xAD ðt on Þ and xAI ðt on þ t off Þ ¼ lon xAD ðt on þ t off Þ, which are satisfied by the state points placed on the boundaries, and the analytical solution of the off-model given by Eq. (B.4), we have ð1Þ

ð2Þ

ð2Þ

ð1Þ

koff ðeλoff t off  eλoff toff Þ þloff eλoff t off ¼ lon eλoff t off :

ðB:6Þ

Solving this equation with respect to t off yields Eq. (31). Appendix C. Inequality (45) as the necessary condition for inequality (44) Here, we show that inequality (45) is the necessary condition for inequality (44) to be satisfied. We first claim that kon o 0 and koff 4 0. The negativity of kon is ð1Þ

aðtÞ ; q1 ðaðtÞÞ ¼ k3 þ ð1  k3 Þ aðtÞ þ k4

ð1Þ

kon ðeλon ton eλon ton Þ þlon eλon ton ¼ loff eλon ton :

ð2Þ

confirmed by λon o λon in inequality (40), which implies v1;1 ðaon Þ o v2;2 ðaon Þ and v2;1 ðaon Þ 4 0. The positivity of koff is confirmed by ð2Þ λð1Þ off 4 λoff in inequality (42), which implies v1;1 ðaoff Þ 4 v2;2 ðaoff Þ and

ðA:7Þ ðA:8Þ

where m1 represents a maximum mutation rate of AD cells, and d is related to the net growth rate of AI cells (see Ideta et al., 2008).

v2;1 ðaoff Þ 4 0. Thus inequality (44) can be rewritten as ð2Þ ð1Þ  λð1Þ =ðλð1Þ  λð2Þ Þ    λð1Þ on =ðλon  λon Þ loff jkoff j off off off l þ jkon j o off : lon jkoff j lon þ jkon j

ðC:1Þ

Using inequalities (26)–(29), we can confirm that the bases on both sides of (C.1) are positive. Moreover, both of them are greater

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

than unity because loff 4 lon in inequality (41), meaning that the numerator is greater than the denominator. We can also confirm that the exponents on both sides of (C.1) are positive, by using inequalities (40), (42), and (43). We then confirm that the inequality     loff  jkoff j l þ jkon j 4 off ðC:2Þ lon  jkoff j lon þ jkon j holds regardless of the values of jkon j and jkoff j, because this inequality is equivalent to the following trivial inequality: jkon j 4  jkoff j:

ðC:3Þ

Thus, inequality (C.2), together with the fact that the bases on both sides of (C.1) are greater than unity, implies that the inequality ð1Þ ð1Þ off o  ð2Þ on ð1Þ ð1Þ ð2Þ on  on off  off

λ

λ

λ

λ

λ

λ

ðC:4Þ

must hold for (C.1), and thus inequality (44) is satisfied. Rearranging (C.4), we have inequality (45), as claimed. At the end of this appendix, we argue that inequality (C.1) is satisfied more “easily” if jkon j and jkoff j are both small and close to zero. Note that this condition jkon j Cjkoff j C 0;

ðC:5Þ

means that the direction of the stable manifold of the on-model is close to that of the unstable manifold of the off-model and they almost coincide with the xAD-axis. Indeed, if (C.5) holds, the leftand right-hand sides of (C.2) become close to each other. In this case, inequality (C.1) can be satisfied even if the left- and righthand sides of (C.4) are close to each other, so long as (C.1) holds. This situation can be compared with other situations where the left- and right-hand sides of (C.2) are very different, i.e., the lefthand side of (C.2) is much larger than the right-hand side. In such cases, the right-hand side of (C.4) must be much larger than the left-hand side of (C.4). Thus, we can conclude that choosing a parameter that satisfies (C.5) is preferable for stabilizing the simplified ITTA model.

Appendix D. State observer Populations of the prostate cancer cells, xi for iA fAD; AI; rAI; iAIg, are estimated from an observed value of the serum PSA concentration y, as the estimated state x~ i , using a state observer. The state observer (Eq. (47)) used in the modified ITTA model (Fig. 1(c)) is defined as follows: ! ! x~ AD ðtÞ d x~ AD ðtÞ ~ ðD:1Þ ¼ Ak þ Kk ðy  yÞ; x~ AI ðtÞ dt x~ AI ðtÞ ! x~ AD ðtÞ ~ ¼C yðtÞ ; x~ AI ðtÞ

ðD:2Þ

where Kk for k A fon; offg is the observer gain matrix, and C is a coefficient matrix C ¼ ðc1 c2 Þ. Let e be an error vector between the actual state values, ðxAD ; xAI Þ, and the estimated state values, ðx~ AD ; x~ AI Þ. As shown in Fig. 2(d), the dynamics of the ITTA model is quite similar to that of the simplified ITTA model. Therefore, the dynamics of ðxAD ; xAI Þ can be approximated by Eq. (8). Then, from Eqs. (4), (8), (D.1) and (D.2), the dynamics of the error can be described by d e ¼ ðAk Kk CÞe: dt

ðD:3Þ

15

The error e may converge to zero for an appropriate Kk . In this study, we set Kk so that it satisfies the following equation: jλI  ðAk  Kk CÞj ¼ ðλ  μk;1 Þðλ  μk;2 Þ;

ðD:4Þ

where I is the identity matrix. μk;1 and μk;2 for k A fon; off g are the parameters representing the desired eigenvalues λ of Ak  Kk C. If the values of μk;1 and μk;2 are negative, then the error e converges to zero with the time constants of 1=jμk;1 j and 1=jμk;2 j. Because the estimation is based on the linearized androgen dynamics, x~ AD and x~ AI may differ from the actual values of xAD and xAI. Moreover, during each transient period after each model transition, it will take about 1=jμk;1 j or 1=jμk;2 j days for the error to be close to zero. A model transition triggered by the estimated state of x~ AD and x~ AI with an error may lead to failure of the control. For this reason, we introduce the dead time τ for the region-dividing method to avoid the use of x~ AD and x~ AI during the transient period after each model transition. See Tables 1 and 2 for the gain matrices Kk used for the modified ITTA model and the modified HBA model.

References Akakura, K., Bruchovsky, N., Goldenberg, S., Rennie, P., Buckley, A., Sullivan, L., 1993. Effects of intermittent androgen suppression on androgen-dependent tumors: apoptosis and serum prostate-specific antigen. Cancer 71, 2782–2790. Asai, Y., Tasaka, Y., Nomura, K., Nomura, T., Casadio, M., Morasso, P., 2009. A model of postural control in quiet standing: robust compensation of delay-induced instability using intermittent activation of feedback control. PLoS One 4 (7), e6169, http://dx.doi.org/10.1371/journal.pone.0006169. Boiko, I., Fridman, L., 2005. Analysis of chattering in continuous sliding-mode controllers. IEEE Trans. Autom. Control 50, 1442–1446. Bottaro, A., Casadio, M., Morasso, P., Sanguineti, V., 2005. Body sway during quiet standing: is it the residual chattering of an intermittent stabilization process? Human Mov. Sci. 24, 588–615. Bottaro, A., Yasutake, Y., Nomura, T., Casadio, M., Morasso, P., 2008. Bounded stability of the quiet standing posture: an intermittent control model. Human Mov. Sci. 27, 473–495. Bruchovsky, N., Klotz, L., Crook, J., Goldenberg, S., 2007. Locally advanced prostate cancer—biochemical results from a prospective phase II study of intermittent androgen suppression for men with evidence of prostate-specific antigen recurrence after radiotherapy. Cancer 109, 858–867. Bruchovsky, N., Klotz, L., Crook, J., Malone, S., Ludgate, C., Morris, W., Gleave, M., Goldenberg, S., 2006. 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: clinical parameters. Cancer 107, 389–395. Charles, H., Clarence, V.H., 1941. Studies on prostatic cancer. I. The effect of castration, of estrogen and of androgen injection on serum phosphatases in metastatic carcinoma of the prostate. Cancer Res. 1, 293–297. Crook, J., Szumacher, E., Malone, S., Huan, S., Segal, R., 1999. Intermittent androgen suppression in the management of prostate cancer. Urology 53, 530–534. Edwards, J., Bartlett, J., 2005. The androgen receptor and signal-transduction pathways in hormone-refractory prostate cancer. Part 2. Androgen-receptor cofactors and bypass pathways. BJU Int. 95, 1327–1335. Feldman, B., Feldman, D., 2001. The development of androgen-independent prostate cancer. Nat. Rev. Cancer 1, 34–45. Gleave, M., Goldenberg, S., Bruchovsky, N., Rennie, P., 1998. Intermittent androgen suppression for prostate cancer: rationale and clinical experience. Prostate Cancer Prostatic Dis. 1, 289–296. Guo, Q., Tao, Y., Aihara, K., 2008. Mathematical modeling of prostate tumor growth under intermittent androgen suppression with partial differential equations. Int. J. Bifurc. Chaos 18, 3789–3797. Heidenreich, A., Aus, G., Bolla, M., Joniau, S., Matveev, V., Schmid, H., Zattoni, F., 2008. EAU guidelines on prostate cancer. Eur. Urol. 53, 68–80. Hirata, Y., Bruchovsky, N., Aihara, K., 2010. Development of a mathematical model that predicts the outcome of hormone therapy for prostate cancer. J. Theor. Biol. 264, 517–527. Ideta, A., Tanaka, G., Takeuchi, T., Aihara, K., 2008. A mathematical model of intermittent androgen suppression for prostate cancer. J. Nonlinear Sci. 18, 593–614. Jain, H., Clinton, S., Bhinder, A., Friedman, A., 2011. Mathematical modeling of prostate cancer progression in response to androgen ablation therapy. Proc. Natl. Acad. Sci. U.S.A. 108, 19701–19706. Jain, H., Friedman, A., 2013. Modeling prostate cancer response to continuous versus intermittent androgen ablation therapy. Discrete Contin. Dyn. Syst. Ser. B 18, 945–967. Levant, A., 2010. Chattering analysis. IEEE Trans. Autom. Control 55, 1380–1389. Mackean, M., Planting, A., Twelves, C., Schellens, J., Allman, D., Osterwalder, B., Reigner, B., Griffin, T., Kaye, S., Verweij, J., 1998. Phase I and pharmacologic

16

Y. Suzuki et al. / Journal of Theoretical Biology 350 (2014) 1–16

study of intermittent twice-daily oral therapy with capecitabine in patients with advanced and/or metastatic cancer. J. Clin. Oncol. 16, 2977–2985. Mearini, L., Zucchi, A., Costantini, E., Bini, V., Porena, M., 2011. Intermittent androgen suppression in prostate cancer: testosterone levels and its implication. J. Sex. Med. 8, 1218–1227. Mottet, N., 2013. Intermittent androgen deprivation therapy in prostate cancer: is everything so clear?. Eur. Urol. 63, 121–122. Mottet, N., Van Damme, J., Loulidi, S., Russel, C., Leitenberger, A., Wolff, J., 2012. Intermittent hormonal therapy in the treatment of metastatic prostate cancer: a randomized trial. BJU Int. 110, 1262–1269. Ng, E., Woo, H., Turner, S., Leong, E., Jackson, M., Spry, N., 2012. The influence of testosterone suppression and recovery on sexual function in men with prostate cancer: observations from a prospective study in men undergoing intermittent androgen suppression. J. Urol. 187, 2162–2166. Oudard, S., 2013. Progress in emerging therapies for advanced prostate cancer. Cancer Treat. Rev. 39, 275–289. Pentyala, S., Lee, J., Hsieh, K., Waltzer, W., Trocchia, A., Musacchia, L., Rebecchi, M., Khan, S., 2000. Prostate cancer: a comprehensive review. Med. Oncol. 17, 85–105. Pollack, A., Zagars, G., Kavadi, V., 1994. Prostate specific antigen doubling time and disease relapse after radiotherapy for prostate cancer. Cancer 74, 670–678. Portz, T., Kuang, Y., Nagy, J., 2012. A clinical data validated mathematical model of prostate cancer growth under intermittent androgen suppression therapy. AIP Adv. 2, 011002. Salonen, A., Taari, K., Ala-Opas, M., Viitanen, J., Lundstedt, S., Tammela, T., 2013. Advanced prostate cancer treated with intermittent or continuous androgen deprivation in the randomised finnprostate study. VII. Quality of life and adverse effects. Eur. Urol. 63, 111–120. Scheithauer, W., Kornek, G., Raderer, M., Schüll, B., Schmid, K., Längle, F., Huber, H., 2002. Intermittent weekly high-dose capecitabine in combination with oxaliplatin: a phase I/II study in first-line treatment of patients with advanced colorectal cancer. Ann. Oncol. 13, 1583–1589. Scheithauer, W., Schüll, B., Ulrich-Pur, H., Schmid, K., Raderer, M., Haider, K., Kwasny, W., Depisch, D., Schneeweiss, B., Lang, F., Kornek, G., 2003. Biweekly high-dose gemcitabine alone or in combination with capecitabine in patients with metastatic pancreatic adenocarcinoma: a randomized phase II trial. Ann. Oncol. 14, 97–104. Sharifi, N., Gulley, J., Dahut, W., 2005. Androgen deprivation therapy for prostate cancer. J. Am. Med. Assoc. 294, 238–244. Shimada, T., Aihara, K., 2008. A nonlinear model with competition between prostate tumor cells and its application to intermittent androgen suppression therapy of prostate cancer. Math. Biosci. 214, 134–139.

Siegel, R., Ward, E., Brawley, O., Jemal, A., 2011. Cancer statistics, 2011: the impact of eliminating socioeconomic and racial disparities on premature cancer deaths. CA Cancer J. Clin. 61, 212–236. Spry, N., Kristjanson, L., Hooton, B., Hayden, L., Neerhut, G., Gurney, H., Corica, T., Korbel, E., Weinstein, S., McCaul, K., 2006. Adverse effects to quality of life arising from treatment can recover with intermittent androgen suppression in men with prostate cancer. Eur. J. Cancer 42, 1083–1092. Suzuki, T., Aihara, K., 2013. Nonlinear system identification for prostate cancer and optimality of intermittent androgen suppression therapy. Math. Biosci. 245, 40–48. Suzuki, Y., Nomura, T., Casadio, M., Morasso, P., 2012. Intermittent control with ankle, hip, and mixed strategies during quiet standing: a theoretical proposal based on a double inverted pendulum model. J. Theor. Biol. 310, 55–79. Tan, N., Margolis, D., McClure, T., Thomas, A., Finley, D., Reiter, R., Huang, J., Raman, S., 2012. Radical prostatectomy: value of prostate MRI in surgical planning. Abdom. Imaging 37, 664–674. Tanaka, G., Tsumoto, K., Tsuji, S., Aihara, K., 2008. Bifurcation analysis on a hybrid systems model of intermittent hormonal therapy for prostate cancer. Physica D 237, 2616–2627. Tao, Y., Guo, Q., Aihara, K., 2009. A model at the macroscopic scale of prostate tumor growth under intermittent androgen suppression. Math. Models Methods Appl. Sci. 19, 2177–2201. Tao, Y., Guo, Q., Aihara, K., 2010. A mathematical model of prostate tumor growth under hormone therapy with mutation inhibitor. J. Nonlinear Sci. 20, 219–240. Tao, Y., Guo, Q., Aihara, K. A partial differential equation model and its reduction to an ordinary differential equation model for prostate tumor growth under intermittent hormone therapy. J. Math. Biol. (in press), http://dx.doi.org/10. 1007/s00285-013-0718-y. Traina, T., Theodoulou, M., Feigin, K., Patil, S., Tan, K., Edwards, C., Dugan, U., Norton, L., Hudis, C., 2008. Phase I study of a novel capecitabine schedule based on the Norton–Simon mathematical model in patients with metastatic breast cancer. J. Clin. Oncol. 26, 1797–1802. Utkin, V.I., 1977. Variable structure systems with sliding modes. IEEE Trans. Autom. Control AC-22, 212–222. Young, K., Utkin, V., Özgüner, Ü., 1999. A control engineer's guide to sliding mode control. IEEE Trans. Control Syst. Technol. 7, 328–342. Yu, H., Diamandis, E., Prestigiacomo, A., Stamey, T., 1995. Ultrasensitive assay of prostatespecific antigen used for early detection of prostate cancer relapse and estimation of tumor-doubling time after radical prostatectomy. Clin. Chem. 41, 430–434.