## Abstract

Hormone therapy in the form of androgen deprivation is a major treatment for advanced prostate cancer. However, if such therapy is overly prolonged, tumour cells may become resistant to this treatment and result in recurrent fatal disease. Long-term hormone deprivation also is associated with side effects poorly tolerated by patients. In contrast, intermittent hormone therapy with alternating on- and off-treatment periods is a possible clinical strategy to delay progression to hormone-refractory disease with the advantage of reduced side effects during the off-treatment periods. In this paper, we first overview previous studies on mathematical modelling of prostate tumour growth under intermittent hormone therapy. The model is categorized into a hybrid dynamical system because switching between on-treatment and off-treatment intervals is treated in addition to continuous dynamics of tumour growth. Next, we present an extended model of stochastic differential equations and examine how well the model is able to capture the characteristics of authentic serum prostate-specific antigen (PSA) data. We also highlight recent advances in time-series analysis and prediction of changes in serum PSA concentrations. Finally, we discuss practical issues to be considered towards establishment of mathematical model-based tailor-made medicine, which defines how to realize personalized hormone therapy for individual patients based on monitored serum PSA levels.

## 1. Introduction

Hormone therapy has been a major treatment for advanced or metastasized prostate cancer (Brawer 2006). The principle of hormone therapy is androgen suppression, which reduces the level of androgens and thereby inhibits the stimulation for growth of prostate cancer cells. Androgen deprivation, or castration, can be achieved by surgical removal of testes which are the primary source of androgens (Huggins & Hodges 1941) or administration of pharmacological agents such as luteinizing hormone-releasing hormone analogues. Androgen-receptor antagonists, or antiandrogens, are capable of inhibiting the binding of androgens to their receptors and obstructing the normal activation pathway. Combined usage of medical castration and antiandrogens, called maximal androgen blockade or combined androgen blockade, is often adopted for advanced prostate cancer (Chodak 2005).

A good biomarker for prostate cancer is serum prostate-specific antigen (PSA), which is a protein produced by cells in the prostate gland (Rao *et al.* 2008). When a rising serum PSA level exceeds a certain critical value, androgen suppression is initiated. In most cases, androgen deprivation is temporarily effective in reducing the serum PSA level, which implies tumour regression. Cancer cells sensitive to androgen suppression are called androgen-dependent (AD) cells. During several years of continuous androgen suppression (CAS), however, many patients will suffer a relapse, or recurrent tumour growth. The cause of a relapse is considered to be progression to androgen-independent (AI) prostate cancer, which is no longer responsive to androgen suppression and likely to grow even in androgen-poor conditions (Rennie *et al.* 1990; Bladou *et al.* 1996).

To delay progression to androgen independence, intermittent androgen suppression (IAS) was proposed first in animal models (Bruchovsky *et al.* 1990; Akakura *et al.* 1993; Sato *et al.* 1996) and since then has been explored in clinical studies to establish efficacy of the therapy (Akakura *et al.* 1993; Goldenberg 1995; Bruchovsky *et al.* 2000, 2001, 2006, 2007; Hurtado-Coll *et al.* 2002). The status of further phase-II and ongoing phase-III trials of IAS is the subject of several review articles (Gleave 1998; Rashid & Chaudhary 2004; Bhandari *et al.* 2005; Abrahamsson 2010). Figure 1 is a schematic illustration showing a continuation of cycles of on- and off-treatment periods under IAS with the serum PSA levels oscillating in a clinically feasible range.

To understand the mechanism of progression to androgen independence and examine the effects of CAS and IAS therapies, mathematical approaches to hormone therapy for prostate cancer have been developed recently. The present paper deals with mathematical modelling of prostate cancer growth and its application to clinical data analysis.

In an effort to assist in predicting the progression to hormone-refractory prostate cancer and a resulting relapse, a mathematical model was proposed that describes the growth of tumours consisting of AD and AI cells under CAS (Jackson 2004*a*). This model reproduces the experimentally observed three phases of an evolving tumour, i.e. the exponential pretreatment tumour growth, the transient androgen-sensitive period immediately after treatment and the eventual AI relapse of the tumour. However, the analysis of the model suggests that CAS can only be successful in preventing a relapse for a small range of biological parameters (Jackson 2004*b*).

Ideta *et al.* (2008) proposed a mathematical model describing prostate tumour growth under IAS based on the former model for CAS by considering intermittent treatment protocols. The model is formulated as a hybrid dynamical system, in which two different continuous systems for tumour growth in on-treatment and off-treatment periods are switched according to the IAS protocol. The model is capable of replicating three different kinds of qualitative behaviours depending on the parameter values, i.e. a relapse without undergoing an off-treatment period like that in CAS, a relapse after several cycles of on-treatment and off-treatment periods and prevention of a relapse repeating the cycles forever.

Phase transitions between the regimes of relapse and relapse prevention in the IAS therapy model were revealed from a mathematical viewpoint by numerical bifurcation analyses (Tanaka *et al.* 2008). It is notable that a grazing bifurcation of a limit cycle and a chaotic attractor, which cannot be observed in ordinary smooth dynamical systems, plays a crucial role in the phase transition (Tanaka *et al.* 2009).

By modifying the IAS therapy model, a competition between AD and AI cells in androgen-poor conditions was considered as a candidate mechanism to enhance the possibility of relapse prevention (Shimada & Aihara 2008). The IAS therapy model represented by switched ordinary differential equations was also developed into a switched partial differential equations version by considering the volume of the tumour (Guo *et al.* 2008; Tao *et al.* 2009).

Although these deterministic models are successful in qualitatively reproducing typical tumour growth under IAS therapy, in reality, noise and fluctuations are inevitably involved in the dynamics of cell populations and monitoring of the serum PSA levels. Therefore, it is important to investigate how the IAS therapy model is sensitive to stochastic noise. This is why the first purpose of this article is to develop a stochastic differential equation model with dynamical and observational noise under the protocol of intermittent hormone therapy used in recent clinical study (Bruchovsky *et al.* 2006, 2007).

Following the developments of the fundamental models of IAS therapy, we have proposed a method to analyse real time-series data of serum PSA concentrations based on a piecewise linear model (Hirata *et al.* 2010). To complement this study, we explain why the piecewise linear model, instead of conventional models, is necessary for characterizing actual observations of clinical PSA data. Moreover, we carry out additional clinical data analysis to further examine the validity of the method. This is the second purpose of this article.

In §2, we present an extended model of switched stochastic differential equations for IAS therapy based on the deterministic model (Ideta *et al.* 2008) and show how well it reproduces the properties of the real clinical data of the serum PSA level. In §3, using a piecewise linear model, we show time-series analyses of the serum PSA data and classification of patients based on time-series predictions. Finally, discussions are given in §4.

## 2. Stochastic model

### (a) Stochastic hybrid systems

The deterministic model (Ideta *et al.* 2008) is able to reproduce qualitative behaviour of prostate tumour growth under IAS and tractable in classifications of stationary states. However, some realistic factors should be taken into account to compare with the statistical results of clinical studies. Thus, we extend the former model as follows. First, stochastic noise is introduced to represent intrinsic fluctuations of the tumour dynamics and observation errors related to monitoring of the serum PSA levels. Second, the protocol of IAS is changed to that adopted in the clinical trials (Bruchovsky *et al.* 2006) with some simplification. The serum PSA level and the androgen level are monitored every 4 weeks. If the serum PSA level is less than *r*_{0}=4 ng ml^{−1} at 36 weeks after the initiation of therapy, the treatment is stopped at that time. Otherwise, the patient is not eligible for interruption of therapy and is withdrawn from the study. In the off-treatment period, therapy is resumed when the serum PSA level increases to *r*_{1}=10 ng ml^{−1}.

The extended model of stochastic differential equations is represented as follows:
2.1
2.2
2.3
2.4
and
2.5
where *x*, *y*, *z* and *v* denote the population of AD cells, the population of AI cells, the androgen concentration and the serum PSA concentration, respectively. The unit of time *t* is assumed to be given by ‘day’. The binary variable *u* denotes on-treatment (*u*=1) or off-treatment (*u*=0). The time from the start of each on-treatment term is denoted by Δ*t*_{on}. The stochastic noise is given by Gaussian white noise *ξ*_{i} (*i*=*x*, *y*, *z*, *v*) with 〈*ξ*_{i}(*t*)〉=0 and 〈*ξ*_{i}(*t*)*ξ*_{j}(*t*′)〉=*D*_{i}*δ*_{ij}*δ*(*t*−*t*′), where *D*_{i} is the noise intensity and 〈⋅〉 means time average.

The dynamics of tumour growth described in equations (2.1) and (2.2) depends on the androgen concentration *z*, which exponentially decreases during on-treatment periods and increases to the normal level *z*_{0} during off-treatment periods following equation (2.3). The growth and mutation rates will be discussed later.

For simplicity, the serum PSA level *v* is assumed to be representative of the sum of the populations of AD and AI cells as in equation (2.4). For the most part, serum PSA is correlated with tumour size, although the relationship may differ between groups of patients in different eras (Ochiai *et al.* 2007). Based on the serum PSA level *v*, interruption and restart of treatment are determined by the switching rule (2.5). This rule means that the treatment is resumed when an increasing PSA level reaches the upper threshold *r*_{1} and suspended when a decreasing PSA level falls below the lower threshold *r*_{0} at 36 weeks after the onset of each cycle. It should be noted that the patient is withdrawn from the study if the PSA level is not less than *r*_{0} at 36 weeks after resumption of treatment. The values of the clinically controllable parameters *r*_{1} and *r*_{0} are fixed in numerical simulations according to the protocol used in clinical trials (Bruchovsky *et al.* 2006).

### (b) Growth and mutation rates of cancer cells

The growth rate *G*_{x} of AD cells, the growth rate *G*_{y} of AI cells and the mutation rate *M*_{xy} from AD phenotype to AI phenotype, are, respectively, modelled as follows (Ideta *et al.* 2008):
2.6
2.7
and
2.8
On the right-hand side of equations (2.6) and (2.7), the first terms correspond to the proliferation rates and the second terms correspond to the apoptosis rates.

Table 1 lists the baseline parameter values for numerical simulations and comments on how they were determined.

Figure 2 shows the above functions of the growth and mutation rates of the tumour. The net growth rate of AD cells is positive with a high concentration of androgens but negative with a low concentration of androgens as shown in figure 2*a*.

The net growth rate of AI cells is assumed to bear a linear relationship with the androgen level for simplicity and the slope parameter *d* is introduced as shown in figure 2*b*, although the relationship can be nonlinear in general (Kokontis *et al.* 1994). It is unclear whether AI cells decrease in an androgen-rich condition, but there are some experimental results supporting this assumption (Kokontis *et al.* 1998, 2005; Chuu *et al.* 2005).

Figure 2*c* shows that the mutation rate from the AD cell phenotype to the AI cell phenotype is enhanced as androgens are depleted. Although only the irreversible mutational effect is considered in the present model for simplicity, reversible changes between AD and AI cell phenotypes and other pathways to AI are possible (Feldman & Feldman 2001; Scher *et al.* 2004; Dehm & Tindall 2005; Edwards & Bartlett 2005). This point will be discussed in §3.

### (c) Numerical simulations

The biochemical results from a prospective phase II study of IAS (Bruchovsky *et al.* 2006, 2007) were analysed for correlations to the onset of hormone-refractory disease. The analyses of time courses of the serum PSA and testosterone levels suggested some statistical trends. Testosterone is an androgen steroid hormone primarily secreted in the testes of males. We examine how well the stochastic model (2.1)–(2.5) is able to reproduce the features of the actual clinical data.

The numerical integration of the stochastic differential equations (2.1)–(2.5) was carried out using the Heun scheme with a time step of 0.005 (Kloeden & Platen 1992). The intensities of the dynamical and observational noise are assumed to be given as *D*_{x}=*D*_{y}=*D*_{z}=*D*_{d} and *D*_{v}=*D*_{o}, respectively. These parameter values are set at *D*_{d}=0.001 and *D*_{o}=0.0001, unless otherwise noted. The maximum observation time is fixed at 2500 days.

Figure 3 shows an example of the time evolutions of the serum PSA and testosterone levels. The therapy is terminated after 4.5 cycles of the on-treatment and off-treatment intervals, because the serum PSA level is not less than 4 ng ml^{−1} at 36 weeks after the beginning of the fifth cycle. The duration of the on-treatment period is 36 weeks for every cycle as specified by the protocol (2.5), while that of the off-treatment period gradually decreases as the cycle number increases. The decrease of the ratio of the off-treatment period in each succeeding cycle was also observed in the clinical data (Bruchovsky *et al.* 2006). The shortening of the off-treatment period may be correlated to the development of early AI. The time course of the serum PSA level also demonstrates that the serum PSA nadir level increases with more cycles of therapy. This tendency was noted in the clinical database of Bruchovsky *et al.* (2007).

The features of the serum PSA level depend upon the parameter *d*, which characterizes the growth rate of AI cells. As shown in figure 2*b*, the growth rate of AI cells decreases as the value of *d* increases. Therefore, both the average number of cycles and the average time to relapse are likely to increase as *d* itself increases (figure 4). In the deterministic case with *D*_{d}=0 and *D*_{o}=0, the serum PSA level evolves in the same manner for a fixed value of *d* without trial-to-trial variations. The low noise case with *D*_{d}=0.001 and *D*_{o}=0.0001 seems to be reasonable for comparison with the statistics of clinical data, because it produces variable time courses of the serum PSA levels without losing the tendency of the deterministic system. It is remarkable that the number of cycles is variable even for the same value of *d* owing to the low level of stochastic noise. This is caused by an interplay between the switching dynamics of the system and the stochastic noise. Noise has a great effect particularly on a trajectory near the boundary (*x*+*y*=*r*_{0}) for switching of dynamics: if a trajectory hits the boundary, then the next cycle follows; otherwise a trajectory diverges and no cycle follows. Therefore, this result suggests that even small stochastic noise can contribute to the diversity of clinical data in different patients. The other two noise levels in figure 4 are too high to mimic a time evolution of the serum PSA level in actual observations. In fact, while the maximum observation time is reached for *d*>0.6 as shown in figure 4*b*, the mean number of cycles is much less than that in the deterministic case as shown in figure 4*a*. Namely, the length of one cycle is very different from the actual data.

Figure 5 shows the effect of the initial testosterone level on the mean off-treatment time. In figure 5*a*,*b*, the mean length of the off-treatment periods decreases with each succeeding cycle. This trend can be found in the clinical data (Bruchovsky *et al.* 2006). However, the trend is not applicable to the case when the initial testosterone level is relatively high as shown in figure 5*c*. The reason for this is explained as follows. Since a higher initial testosterone level requires more time to be reduced by treatment, a longer delay ensues until the PSA level begins to fall. Therefore, the PSA nadir at 36 weeks after the start of therapy tends to be higher, even though it is below the lower threshold *r*_{0}. Consequently, in the following off-treatment period, the time needed for the PSA level to reach the upper threshold *r*_{1} is shorter. Thus, the length of the off-treatment period is closely related to the testosterone level at the start of therapy in each cycle.

Figure 6 shows correlations between the serum PSA nadir and the duration of the off-treatment period in each cycle. The length of the off-treatment period is longer for a lower nadir level of serum PSA. A high correlation is recognized in cycle 1 (correlation coefficient, cc=−0.61) as shown in figure 6*a*, while not in cycle 2 (cc=−0.38) and cycle 3 (cc=−0.26) as shown in figure 6*b*,*c*. In all cycles, the length of the off-treatment period is related inversely to the serum PSA nadir. This property is consistent with the result of the clinical data analysis (Bruchovsky *et al.* 2007).

In summary, the proposed stochastic model has generated variable time courses of the serum PSA and testosterone levels. The results demonstrate that stochastic noise can be responsible for the diversity of clinical PSA data in individual patients. The large effect of noise is caused by its interaction with switching properties of the system. Furthermore, we have demonstrated that the model is able to reproduce the characteristics observed in clinical trials, i.e. the shortening of the off-treatment period with more cycles and the inverse relationship between the serum PSA nadir level and the duration of the off-treatment period. Therefore, the stochastic model is considered to be more realistic than the former deterministic model.

## 3. Application to clinical data analysis

In this section, we apply our theoretical model to clinical data and show that it can approximate the behaviour of a tumour. To make predictions, a mathematical model should be accurate enough to reproduce growth and regression of prostate cancer under IAS. However, the variations of the serum PSA levels generated by system (2.1)–(2.5) appear to be different from the actual observations in figure 7. In system (2.1)–(2.5), the PSA level decreases sharply to a nadir during the on-treatment period as shown in figure 3. In reality, however, the PSA decline is biphasic with an initial steep slope followed by a more shallow one. Therefore, additional factors in system (2.1)–(2.5) need to be considered to reproduce the biphasic kinetics; moreover, redundant factors should be removed to limit the number of parameters to be estimated for modelling so as to avoid overfitting.

First, we observe that the androgen dynamics in equation (2.3) does not contribute much to the form of the reproduced PSA data sequence. Therefore, this property can be eliminated for the purpose of data-fitting by assuming that the androgen level instantaneously switches between the minimum level 0 and the normal level *z*_{0} at the moments of restart and suspension of treatment, i.e. in equation (2.3). At this stage, the model has been reduced to a two-dimensional system. Second, one additional dimension is needed to mimic the biphasic decline of the actual serum PSA level during the on-treatment period. We could reasonably increase the dimensions of the system to three by considering a population of AI cells that undergo reversible changes. This would be in line with the idea that there are reversible pathways based on epigenetic mechanisms that lead to AI (Rennie *et al.* 2005).

Following the above consideration, we proposed a modified model that contains three phenotypes of cancer cells as illustrated in figure 8 (Hirata *et al.* 2010). One phenotype corresponds to AD cells and the other two correspond to AI cells. During the on-treatment period, AD cells may change to AI cells of the first or the second phenotype. During the off-treatment period, AI cells of the first phenotype may return to AD cells, while AI cells of the second phenotype enter an irreversible state and cannot revert to AI cells of the first phenotype or AD cells. We also showed that a piecewise linear model that omits testosterone is adequate enough for the modelling of prostate cancer under IAS.

Thus, we consider the following mathematical model for clinical data analysis. Let *x*_{1}, *x*_{2} and *x*_{3} be the populations of AD cells, AI cells of the first phenotype and AI cells of the second phenotype, respectively. Letting be the contribution of cell population *x*_{j} to cell population *x*_{i} when the treatment is *m*, the dynamics under the on-treatment period is written as
3.1
Similarly, the dynamics under the off-treatment period is written as
3.2
The PSA level is given by *x*_{1}+*x*_{2}+*x*_{3} for simplicity.

The fitting of the above model to the serum PSA data faces some difficulties. Although the model needs to be able to reproduce a relapse to make a prognosis and/or personalize the treatment schedule, the relapsing part is not available for fitting the model. Therefore, simple phenomenological approaches will fail in this application.

To overcome this problem, three realistic constraints during the fitting are considered. The first constraint is that all non-diagonal parameters are non-negative. The second constraint is that each phenotype will change its volume by 20 per cent at most in a day and that at most 10 per cent of each phenotype may change to other phenotype within a day. The third constraint is that under CAS, the PSA level decreases at the beginning of the treatment but starts to increase within 5 years. We enforced these constraints using the penalty method (Price *et al.* 2005). The two subsystems (3.1) and (3.2) are switched according to the switching between on-treatment and off-treatment periods indicated in clinical data.

By using the above method, the actual data were well approximated by our mathematical model as shown in figure 9. Further, based on system (3.1) and (3.2), we can classify patients of clinical trials (Goldenberg 1995; Hurtado-Coll *et al.* 2002; Pether & Goldenberg 2004; Bruchovsky *et al.* 2006, 2007; Shaw *et al.* 2007) into three types.

The first criterion for classification is whether AI cells of the second phenotype will proliferate or not during an on-treatment period. If , the AI cells of the second phenotype can be diminished by treatment, and thus a relapse may be prevented; on the other hand, if *w*^{1}_{3,3}≥0, then the AI cells of the second phenotype will always increase and a relapse cannot be avoided. Therefore, if , we classify the patients to be type (i).

The second criterion is whether the relapse can be delayed or not. If and *w*^{1}_{2,2}≤*w*^{0}_{2,2} and , CAS is the method of choice to reduce PSA to the lowest possible level. If *w*^{1}_{1,1}>*w*^{0}_{1,1} or or *w*^{1}_{3,3}>*w*^{0}_{3,3}, IAS may delay the increase in PSA level, and thus IAS may be the better treatment strategy for these patients when compared with CAS. Therefore, if or *w*^{1}_{2,2}>*w*^{0}_{2,2} or , we classify patients as type (ii). Otherwise, they are classified as type (iii).

In summary, if patients are type (i), IAS may prevent a relapse of malignancy. If type (ii), IAS may delay a relapse. If type (iii), CAS is better than IAS. For details of the derivations of these criteria, see Hirata *et al.* (2010).

Using the model including irreversible and reversible changes and whole datasets, we have classified 90 patients into three types as shown in table 2. Most patients who experienced a relapse were classified as type (ii) or (iii). The majority of the patients without a relapse were classified as type (i). Therefore, the classification by the proposed model is consistent with judgements made by medical doctors of patient responses in the clinical setting.

## 4. Discussion

After reviewing mathematical modelling of IAS for prostate cancer, we have presented a stochastic model and shown that it is successful in reproducing statistical features of the time course of serum PSA in clinical trials (Bruchovsky *et al.* 2006, 2007), such as the shortening of the off-treatment periods with more on- and off-treatment cycles and the inverse relationship between the serum PSA nadir level and the length of the off-treatment interval. However, the stochastic model is not able to replicate the biphasic decrease of the serum PSA level. The model could be further extended by taking into account new findings on the possible pathways to AI in addition to irreversible changes from the AD cell phenotype to the AI cell phenotype. Another possible improvement is to model the growth of both normal and cancer cells.

For time-series analysis of clinical data of the serum PSA concentrations, we have presented a three-dimensional piecewise linear model. By fitting this model to actual data and estimating the model parameters, we have obtained criteria for classification of patients. The classification provides an objective guide for a prognosis and a better choice of treatments for individual patients. The results of predictions were almost consistent with judgements made by medical doctors.

Currently, we require the first 2.5 cycles to obtain a good fit such as the one shown in figure 9. Although it is desirable to estimate the parameter values from as short a time series as possible, we need at least the first 1.5 cycles to make predictions. If the time series does not contain the off-treatment period, we cannot obtain the information related to the dynamics under the off-treatment period. Even if the length of time series of the serum PSA level is long enough to determine a set of parameters, the estimation error can be large since the available time series of the PSA determinations would still be too short. To make a correct diagnosis, we have to overcome this estimation error. One of the ideas is to use past datasets of other patients as *a priori* knowledge for constructing the prior distribution for the parameters and initial conditions. Parameter estimation and overcoming the estimation error are two important topics to make our mathematical analysis standard practice.

Our final goal is to tailor and individualize hormone therapy for prostate cancer using a mathematical approach with hybrid systems modelling as described in this article. We infer that such an approach will be of diagnostic and prognostic benefit in a clinical setting. For this purpose, the optimal control method would help to provide insight into the best scheduling of hormone therapy (Suzuki *et al.* 2010).

## Acknowledgements

The authors thank Dr A. M. Ideta, Dr T. Shimada, Dr T. Suzuki, Dr S. Tsuji and Dr P. Rennie for valuable discussions and Dr J. Abersbach for processing of clinical data. This work was partially supported by IIS Evolving Research Fund of the University of Tokyo and the Japan Society for the Promotion of Science (JSPS) through its Funding Programme for World-Leading Innovative R&D on Science and Technology (FIRST Programme).

## Footnotes

One contribution of 12 to a Theme Issue ‘Theory of hybrid dynamical systems and its applications to biological and medical systems’.

- © 2010 The Royal Society