System identification and parameter estimation in mathematical medicine: examples demonstrated for prostate cancer

Abstract

We review our studies on how to identify the most appropriate models of diseases, and how to determine their parameters in a quantitative manner given a short time series of biomarkers, using intermittent androgen deprivation therapy of prostate cancer as an example. Recently, it has become possible to estimate the specific parameters of individual patients within a reasonable time by employing the information concerning other previous patients as a prior. We discuss the importance of using multiple mathematical methods simultaneously to achieve a solid diagnosis and prognosis in the future practice of personalized medicine.

Graphical abstract

Keywords

mathematical medicine / dynamical model / parameter estimation

INTRODUCTION

Developments inmeasurement techniques have enabled us to model diseases quantitatively and mathematically. For example, several mathematical models [ 1– 5] have been proposed for intermittent androgen deprivation therapyforprostate cancer [ 6– 10]. When we employsuch a mathematical dynamical model in a clinical setting, there are certain problems that commonlyappear, even if the target diseases are different. For example, we need to estimate a set of parameters characterizing each patient from a short time series of biomarkers [ 11, 12]. Further problems include how to optimize a treatment schedule given such a personalized set of parameters [ 13– 15].

In this manuscript, we focus on the first problem of mathematical medicine, i.e., how to estimate a set of parameters for quantitatively characterizing each patient. Because themathematical modeling of intermittent androgen deprivation therapyforprostate cancerwith the sensitive biomarker PSA (Prostate-Specific Antigen) is well developed, we employit as an example.

PHYSIOLOGY OF PROSTATE CANCER

Androgen deprivation therapy (ADT) is often used to treat prostate cancer. ADT suppresses the growth of prostate cancer by lowering the androgen level in patients. Although ADT often showsremarkable effects in tumor regression initially, prostate cancer acquiresthe ability to grow without androgen in long-term ADT treatment,resulting in eventual relapse.

Intermittent ADT [ 6– 10] has been introduced in order to overcome this problem and sustain the hormonal sensitivity in prostate cancer. Intermittent ADT is an approach that isbased on the premise that cycles of androgen deprivation followed by re-exposure may delay “androgen independence”, reduce treatment morbidity, and improve quality of life. In intermittent ADT, patients continue ADT until the PSA levelsdeclinesufficiently, andthenceaseADT. After the PSA levels elevate again to upper threshold values, ADTis resumed. From a mathematical viewpoint,oneproblem concerningintermittent ADTis that the protocol of intermittent ADT is dependent onthe decisions of individualmedical doctors, and it is not clearwhether the appliedprotocol is optimal or not for each patient.

In order to personalize intermittent ADT, several mathematical models of intermittent ADThave been proposed [ 1– 5]. Among these models, Ideta et al. [ 1] were the first to mentionthe mathematical possibility that intermittent ADT can preventthe relapse of prostate cancer if the population of androgen independent cancer cells shrinksafterhormonal therapy is stopped. However, this model is unable toexplain the bi-phasic decline of PSA in periods that hormonal therapy is carried out. In order to describe the bi-phasic decline during on-treatment periods, Hirata et al. [ 3] have proposed a three-dimensional piecewise-linear model. This modeldescribes the dynamics of prostate cancer effectively and quantitatively,and it will be introduced mathematically in the next section.

CONSIDEREDMODEL

In general, a dynamical system describing disease dynamics can be written in the following form:

(1) {\displaystyle \frac{dx}{dt}}=f_m(x),

where x describes a state of the target disease, m corresponds to a regimen of treatment, and f_m describes the dynamics of the disease given a regimen m of treatment.

Such a system encompasses intermittent ADT of prostate cancer [ 6– 10]. Here, we describe the model of Hirata et al. [ 3]. The model describes thedynamical behavior of prostate cancer, and thesedynamics have been comprehensivelystudied [ 11, 12, 14– 22].To define the model, we first define the followingthree variables: x_{\mathrm{1}} corresponds to the non-dimensionalized number of androgen dependent cancer cells, and x_{\mathrm{2}} and x_{\mathrm{3}} to thenon-dimensionalized numbersof androgen independent cancer cells. When hormone therapy is active( m=\mathrm{1} ),the cancer population specified by x_{\mathrm{1}} may change to that specified by x_{\mathrm{2}} or x_{\mathrm{3}} , and the cancer population specified by x_{\mathrm{2}} may change tothat specified by x_{\mathrm{3}} (see Figure 1). When the hormone therapy is suspended ( m=\mathrm{0} ), androgen independent cells, expressed by x_{\mathrm{2}} , may change back to androgen dependent cells, expressed by x_{\mathrm{1}} . However,androgen independent cells expressed by x_{\mathrm{3}} ,cannot change back to cancer populationspecified by either x_{\mathrm{1}} or x_{\mathrm{2}} (see Figure 1).That is, the androgen independent cellsexpressed by x_{\mathrm{2}} are cells generated through reversible changes, and theandrogen independent cells expressed by x_{\mathrm{3}} are cells generated through irreversible changeslike mutation. Mathematically, the dynamics among these variables is written as

(2) {\displaystyle \frac{d}{dt}}\left(\begin{array}{l}{x}_{\mathrm{1}}\\ x_{\mathrm{2}}\\ x_{\mathrm{3}}\end{array}\right)=\left(\begin{array}{lll}{w}_{\mathrm{1},\mathrm{1}}^{\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ w_{\mathrm{2},\mathrm{1}}^{\mathrm{1}}& w_{\mathrm{2},\mathrm{2}}^{\mathrm{1}}& \mathrm{0}\\ w_{\mathrm{3},\mathrm{1}}^{\mathrm{1}}& w_{\mathrm{3},\mathrm{2}}^{\mathrm{1}}& w_{\mathrm{3},\mathrm{3}}^{\mathrm{1}}\end{array}\right)\left(\begin{array}{l}{x}_{\mathrm{1}}\\ x_{\mathrm{2}}\\ x_{\mathrm{3}}\end{array}\right),

while the hormone therapy is active(m= 1), and

(3) {\displaystyle \frac{d}{dt}}\left(\begin{array}{l}{x}_{\mathrm{1}}\\ x_{\mathrm{2}}\\ x_{\mathrm{3}}\end{array}\right)=\left(\begin{array}{lll}{w}_{\mathrm{1},\mathrm{1}}^{\mathrm{0}}& w_{\mathrm{1},\mathrm{2}}^{\mathrm{0}}& \mathrm{0}\\ \mathrm{0}& w_{\mathrm{2},\mathrm{2}}^{\mathrm{0}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& w_{\mathrm{3},\mathrm{3}}^{\mathrm{0}}\end{array}\right)\left(\begin{array}{l}{x}_{\mathrm{1}}\\ x_{\mathrm{2}}\\ x_{\mathrm{3}}\end{array}\right),

while the hormone therapy is inactive (m= 0). The level of prostate specific antigen (PSA) is now assumed to be p(t)=x_{\mathrm{1}}(t)+x_{\mathrm{2}}(t)+x_{\mathrm{3}}(t) , for simplicity. The unit of t isoneday.

In the paper of Hirata et al. [ 3], a method of fitting with constraints was also proposed to achievephysiological validity of the model. In particular, the method focusedon the non-negativities of x_{\mathrm{1}} , x_{\mathrm{2}} , and x_{\mathrm{3}} ; moderateamounts of changes within a day; and the reproducibility of relapse when hormone therapyis continued. To achievethe non-negativity of x_{\mathrm{1}} , x_{\mathrm{2}} , and x_{\mathrm{3}} , we assumethat w_{i,j}^m\ge \mathrm{0} when i\ne j and a\in \{\mathrm{0},\mathrm{1}\} . To ensure moderate amountsof changes within one day, we assume that -\mathrm{0.2}\le w_{i,i}^m\le \mathrm{0.2} , w_{i,j}^m\le \mathrm{0.1} when i\ne j , and -\mathrm{0.2}\le {\sum }_i{w}_{i,j}^m\le \mathrm{0.2} for j\in \{\mathrm{1},\mathrm{2},\mathrm{3}\} and m\in \{\mathrm{0},\mathrm{1}\} . For the reproducibility of relapse, we assumethat w_{\mathrm{3},\mathrm{3}}^{\mathrm{1}}\ge \mathrm{0} , \overline{p}(\mathrm{360})\le \mathrm{2} , and \overline{p}(\mathrm{360}\times \mathrm{5})\ge \mathrm{10} , where \overline{p}(t) is the level of PSA at time in the case that we keep the hormone therapyactive ( m=\mathrm{1} ) for all time. We use the penalty method to enforce these conditions. If we would like to realize x\ge \mathrm{0} , then we minimize the following quantity h(x) :

(4)

The function h(x)=\mathrm{0} when x\ge \mathrm{0} is satisfied. Thus, bycombining thiswith the fitting error represented by the L_{\mathrm{1}} norm, the cost function to determine the parameters and initial conditions q=(w,x(\mathrm{0})) can be written as follows:

(5) \underset{q}{\min }\{{\sum }_{i=\mathrm{1}}^I|o_i-p(t_i)|+C(q)\},

where \{(t_i,o_i)|i=\mathrm{1},\mathrm{2},\cdots ,I\} aretheobserved PSA values o_i after t_i days, and

(6) \begin{array}{l}C(q)={\displaystyle \sum _{m,i\ne j}h(w_{i,j}^m)}+{\displaystyle \sum _{i,m}h(w_{i,i}^m+\mathrm{0.2})+{\displaystyle \sum _{i,m}h(-w_{i,i}^m+\mathrm{0.2})+{\displaystyle \sum _{m,i\ne j}h(\mathrm{0.1}-w_{i,j}^m)}}}\\ +{\sum }_{m,i\ne j}h({\sum }_i{w}_{i,j}^m+\mathrm{0.2})+{\sum }_{m,i\ne j}h({\sum }_i{w}_{i,j}^m+\mathrm{0.2})+h(\mathrm{2}-\overline{p}(\mathrm{360}))+h(\overline{p}(\mathrm{360}\times \mathrm{5})-\mathrm{10}).\end{array}

The last two terms of Equation (6) correspond toconstraints realizing the assumptions thatthe PSA value will once go down and reach the value less than 2 ng/ml one year after the start of the hormone therapy, while it will become greater than 10 ng/ml after 5 years if we continue the hormone therapy. Theseassumptions came from general common observations made by urologists [ 3].

Asummary of recent methods [ 11, 12, 17, 19, 20] for the system identification and parameter estimation is presentedin Table 1. The methods are roughly classified into three classes: parametric estimation, semi-parametric estimation, and non-parametric estimation. In the following sections, we introduce these methods one at a time.

PARAMETRIC ESTIMATION

Two parameter-estimation methods have beenproposed for estimating a set of parameters for dynamical diseases. The first method is to use cross entropy [ 19], and the second method is to use the Bayesian theorem with the data ofpast patients [ 12].

The method of cross entropy [ 19]involveschoosing a set of parameters by recursively minimizing the Kullback-Leibler divergence [ 23] between the previous truncated distribution and the current distribution without the truncationsuch that the cost function is sufficiently small.The method of cross entropy can be formulated as follows. First, we define a cost function F(q) as

(7) F(q)=\sqrt{{\sum }_i{(o_i-p(t_i))}^{\mathrm{2}}}+C(q).

At each step of iteration, we randomly draw N sets of parameters q_{\mathrm{1}} , q_{\mathrm{2}} ,…, q_N such that each i^th parameter value follows the Gaussian distribution with mean {\mu }_i^{(k-\mathrm{1})} and standard deviation {\mu }_i^{{\mathrm{2}}^{(k-\mathrm{1})}} , where k represents the iteration number. Then, we obtain the r-quantile {\gamma }^{(k)} among

\{F(q_n)|n=\mathrm{1},\mathrm{2},\dots ,N\}

. After obtaining the optimizer

q*=\underset{q\in \{q_n|n=\mathrm{1},\mathrm{2},\dots ,N\}}{\arg \min }F(q)

, we move each q_n slightly closer tothe optimizer q*, namely {\widehat{q}}_n=\delta q*+(\mathrm{1}-\delta )q_n . Then, we choose mean {\widehat{\mu }}_i^{(k)} and deviation {\widehat{\sigma }}_i^{{\mathrm{2}}^{(k)}} such that we can minimize the Kullback-Leibler divergence between the previous truncated distribution and the current distribution.Equivalently,we maximize the following quantity in terms of the parameters \left\{{\widehat{\mu }}_i^{(k)}\right\} and \left\{{\widehat{\sigma }}_i^{{\mathrm{2}}^{(k)}}\right\} of the current distribution, expressed by the product of univariate Gaussian distributions:

(8)

Where \varsigma is a weight, and f_k follows the product of univariate Gaussian distributions of mean {\widehat{\mu }}_i^{(k)} and deviation {\widehat{\sigma }}_i^{{\mathrm{2}}^{(k)}} . By differentiating Equation (8) in terms of {\widehat{\mu }}_i^{(k)} and {\widehat{\sigma }}_i^{{\mathrm{2}}^{(k)}} , we obtain

(9)

and

(10)

Then, we adjust {\mu }_i^{(k-\mathrm{1})} and {\sigma }_i^{{\mathrm{2}}^{(k-\mathrm{1})}} towards {\widehat{\mu }}_i^{(k)} and {\widehat{\sigma }}_i^{{\mathrm{2}}^{(k)}} , as

(11) {\mu }_i^{(k)}=\alpha {\widehat{\mu }}_i^{(k)}+(\mathrm{1}-\alpha ){\mu }_i^{(k-\mathrm{1})},

(12) {\sigma }_i^{{\mathrm{2}}^{(k)}}={\displaystyle \frac{\beta }{k}}{\widehat{\sigma }}_i^{{\mathrm{2}}^{(k)}}+\left(\mathrm{1}-{\displaystyle \frac{\beta }{k}}\right){\sigma }_i^{{\mathrm{2}}^{(k-\mathrm{1})}},

and completea set of routines for an iteration. Here,we set \alpha \in (\mathrm{0.5},\mathrm{1}) and \beta \in (\mathrm{0.8},\mathrm{1}) . We stop the iterations when eitherthe cost function reaches our intended lower bound or the number of iterations reaches a certain number.

This method was demonstrated by Guo et al. [ 19]. The parameter values can be estimated using the first two and half cycles of intermittent ADT. The prediction intervals for the PSA values were also estimated.

The second method is based on the Bayesian theorem, with a prior that is constructed based on the data ofpast patients [ 12]. Suppose that we maximize the posterior probability Q(q|o) of a set q of parameters, given a dataset o. Using the Bayesian theorem, the posterior probability can be written as follows:

(13) Q(q|o)\propto Q(o|q)Q(q).

By taking the logarithm, Equation (13) can berewritten as

(14) \ln Q(q|o)=\ln Q(o|q)+\ln Q(q)+V,

where V is constant. We now model the prior distribution Q(q) byusing the multivariate Gaussian distribution estimated from past patients. Letting \overline{q} and W denotethe mean and the covariance matrix, and enforcing the constraints of Equation (6) for the physiological appropriatenesss, the problem of estimating a set of parameters for the current patient can be written as

(15) \underset{q}{\min }\left\{{\sum }_{i=\mathrm{1}}^I{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\sigma }^{\mathrm{2}}}}{(o_i-p(t_i))}^{\mathrm{2}}+{(q-\overline{q})}^T{W}^{-\mathrm{1}}(q-\overline{q})+C(q)\right\},

where \sigma is the standard deviation for the observations of thePSA values.In the paper of Hirata et al. [ 12], we employed the entire datasets of 36other patients, who were not used later in evaluating the fitting accuracyfor constructing the prior distribution.

In the paper of Hirata et al. [ 12], it was shown that the estimated parameter values with the first one and half cycles exhibitedcorrelations with those obtained usingthe entiredatasets.

SEMI-PARAMETRIC ESTIMATION

By a semi-parametric method,we mean a method that combines the parametric model described in Equations (2) and (3) with a non-parametric approach such that we can quantify the uncertainty in terms of the estimation of the parameters in Equations (2) and (3). Two methods have been proposed for estimating a set of parameters byusing the semi-parametric estimation. The first method is that of bootstrapping [ 17]. The second method is calledthe temporal expert advicemethod [ 11].

The bootstrapping method [ 17] is quite simple,and takes into account the uncertainty owing to a few measurements of PSA.Weresample a given dataset with replacements 100 times. For each resampled dataset, we obtain a set of parameters resulting from Equation (5), and represent the uncertainty as the distribution ofthe 100 estimated sets of parameters.The problem withthe bootstrapping method is its huge computational cost.

In the temporal expert advice method [ 11], we prepare experts by using the sets of parameters for past patients and fitting only the initial conditions with the data of thetarget patient. Then, the initial value obtained using the set of parameters for each past patient provides a prediction for the future. Thus, we use the first few time points to evaluate whether the past patients have similar dynamical disease behavior asthe current patient or not.Then,depending on the similarities, we either take the weighted average over the prediction or wetake an ensemble of predictions to obtain theprediction intervals. This method has been tested in Morino et al. [ 11] on the datasets of patients who hadpreviously undergonethe first radical treatments and were waiting for the start of their next additional treatment.

NON-PARAMETRIC ESTIMATION

There havealso beentwo non-parametric-estimation methods proposed in the context of themathematical modeling of prostate cancer under intermittent ADT. The first method is thevariational Bayes [ 20]method,andthe second is theGaussian process [ 20]method.

In the variational Bayesmethod [ 20], we use a piecewise affine model. Consider aphase space X that includes (\mathrm{2}m+\mathrm{1}) dimensional vectors x_t=[y_{t-\mathrm{1}},y_{t-\mathrm{2}},\cdots ,y_{t-m},u_t,u_{t-\mathrm{1}},\cdots ,u_{t-m}], where y_t corresponds to the logarithm of the PSA value at time t, and u_t to the treatment option at time t. Suppose that a partition of thephase space X is given, satisfying X=U_i{X}_i . Then, a piecewise affine model to be usedin this caseis defined as

(16) f(x_t)=w_i^T\left[\begin{array}{l}{x}_t\\ \mathrm{1}\end{array}\right],

for x_t\in X_i . We would like to approximate y_t by f(x_t) . We denote D=\{(x_t,y_t)\} .

When we applythe Gaussian mixture for modeling x_t , it can be written as

(17) Q(x|\Psi ,M)={\sum }_i{\varphi }_{i}N(x|{\rho }_i,S_i^{-\mathrm{1}}),

Where {\varphi }_i denotesthe weight ofthe ith Gaussian distribution,N represents the G aussian distribution, and {\rho }_i and S_i^{-\mathrm{1}} denote the mean and the covariance matrix for the ith Gaussian distribution. We also use a latent variable z_i^n , where if the nth data point belongs to the ith component, then z_i^n=\mathrm{1} , otherwise z_i^n=\mathrm{0} . Then, each x_t is modeled as

(18) Q(x_t|z_i^n=\mathrm{1},\Psi )=N(x_n|{\rho }_i,S_i^{-\mathrm{1}}),

and each y_t is modeled as

(19) Q(y_t|x_t,z_i^n=\mathrm{1},\Theta )=N(y_t|f(x_t),{\beta }_i^{-\mathrm{1}}).

Therefore, the totalprobability that we have data D and latent variables Z=\{z_i^n\} for the Gaussian mixture M, a model of input \Psi , and a model of output \Theta , is given by

(20) Q(D,Z|M,\Psi ,\Theta )={\Pi }_n{\Pi }_i\{{\varphi }_{i}N(x_t|{\rho }_i,S_i^{-\mathrm{1}})N(y_t|f(x_t),{\beta }_i^{-\mathrm{1}})\}.

The estimates of the parameters are obtained as the mean with respect to the posterior distribution constructed by the Bayes theorem.However,the computation of the posterior is a difficult task in the Gaussian mixture model. Thus, in the variational Bayes method we approximate the posterior as a computationally tractable form and obtain the parameters M, \Psi , and \Theta . This is thevariational Bayesmethod for non-parametric estimation.

The Gaussian processmethod [ 20] can be summarized as follows. Suppose that \tilde{f}(x) obeys the zero-mean Gaussian distribution.That is, for any finite input points (x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_k) , the joint distribution of \tilde{f}(x_{\mathrm{1}}),\tilde{f}(x_{\mathrm{2}}),\cdots ,\tilde{f}(x_k) is zero-mean multivariate Gaussian. Then,we can define a covariance function as

(21) k(x,x^{\prime })=E_{\tilde{f}-GP}[\tilde{f}(x)\tilde{f}(x^{\prime })],

where E_{\tilde{f}-GP} represents the operation oftaking the expectation in terms of the Gaussian process \tilde{f} .Let f={[\tilde{f}(x_{\mathrm{1}}),\tilde{f}(x_{\mathrm{2}}),\cdots ,\tilde{f}(x_n)]}^T , and f*=\tilde{f}(x) . Then, we can define the joint distribution Q({[f^{T}f*]}^T) as

(22) N\left({[f^{T}f*]}^T|\mathrm{0},\left[\begin{array}{ll}{K}_{\mathrm{11}}& K_{\mathrm{12}}\\ K_{\mathrm{12}}^T& K\mathrm{22}\end{array}\right]\right),

where we define

(23) K_{\mathrm{11}}={(k(x_n,x_{n^{\prime }}))}_{n=\mathrm{1},n^{\prime }=\mathrm{1}}^{N,N},

(24) K_{\mathrm{12}}={[k(x_{\mathrm{1}},x),k(x_{\mathrm{2}},x),\cdots ,k(x_n,x)]}^T,

and

(25) K_{\mathrm{22}}=k(x,x).

We employ the distribution introduced above as a prior distribution of [f^T,f*] . We alsomodel the likelihood of f as a Gaussian distribution, as follows:

(26) Q(y|f,{\sigma }^{\mathrm{2}}I).

Because the posterior distribution of f* also follows the Gaussian distribution, the mean \overline{f*} of f* becomes

(27) \overline{f*}=K_{\mathrm{12}}^T{(K_{\mathrm{11}}+{\sigma }^{\mathrm{2}}I)}^{-\mathrm{1}}y.

Writing the mean in a different way, we have

(28) \overline{f*}={\Sigma }_{n}k(x,x_n){\alpha }_n,

where \alpha ={(K_{\mathrm{11}}+{\sigma }^{\mathrm{2}}I)}^{-\mathrm{1}}y . The variance of f* is given by

(29) K_{\mathrm{22}}-K_{\mathrm{12}}^T{(K_{\mathrm{11}}+{\sigma }^{\mathrm{2}}I)}^{-\mathrm{1}}{K}_{\mathrm{12}}.

Suzuki and Aihara [ 20] discussed the factthat the variational Bayes method is slightlymore effectivefor the long-term predictions, while the Gaussian process is better for short-term predictions.

DISCUSSIONS

Here,we return to Table 1. The methods can be discussed fromthree different aspects. The first is whetherwe can provide the parameter distribution and prediction interval. The methods of the bootstrapping [ 17], the cross entropy [ 19], the variational Bayes [ 20], and the Gaussian process [ 20] can provide the parameter distribution, and thus the prediction interval as well. The method of the temporal expert advice [ 11] can provide the prediction interval by weighting predictionsaccording to each expert with Gaussian distributions, although thismethod cannot provide the parameter distribution itself.

The second aspectis whether we use constraints for theparameters. By using constraints for theparameters, we can reproduce relapse under continuous ADT. Although the parametric and semi-parametric estimations can be combined with constraints, and thus reproduce relapse of prostate cancer under continuous ADT, this is not the case fornon-parametric estimation.

The third aspectis whetherwe can use the data ofpast patients simultaneously with the proposed methods. The methods of the temporal expert advice [ 11] and the Bayes methodwiththe prior of past patients [ 12] can be combined with the information of past patients, while the other methods currentlycannot. This property is directly linked with how short thetime series of the current patient can be. Namely, if we usethe data ofpast patients, we can use a shorter time series to estimate theset of parameters for the current patient.

An important point in the actual clinical setting is whether we can provide multiple indices to make diagnosis and prognosis. Because the methods for parameter estimation discussed here can be used in conjunctionwith other mathematical models, and other mathematical models for intermittent androgen deprivation therapy of prostate cancer have already been proposed [ 1, 2, 4, 5], the use ofmultiple methods with various models will provide useful information. This is afuture direction thatwe should explore in order to make furtherprogress in mathematical medicine.

References

[1] Ideta, A. M., Tanaka, G., Takeuchi, T. and Aihara,K. (2008) A mathematical model of intermittent androgen suppression for prostate cancer. J. Nonlinear Sci., 18, 593–614

[2] Shimada,T. and 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

[3] Hirata,Y., Bruchovsky,N. and Aihara,K. (2010) Development of a mathematical model that predicts the outcome of hormone therapy for prostate cancer. J. Theor. Biol., 264, 517–527

[4] Jain,H. V., Clinton,S. K., Bhinder,A. and 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

[5] Portz, T., Kuang,Y. and Nagy,J. D. (2012) A clinical data validated mathematical model of prostate cancer growth under intermittent androgen suppression therapy. AIP Adv., 2, 011002

[6] Akakura, K., Bruchovsky, N., Goldenberg, S. L., Rennie, P. S., Buckley, A. R. and Sullivan, L. D. (1993) Effects of intermittent androgen suppression on androgen-dependent tumors. Apoptosis and serum prostate-specific antigen. Cancer, 71, 2782–2790

[7] Bruchovsky, N., Klotz, L., Crook, J., Malone, S., Ludgate, C., Morris, W. J., Gleave, M. E. and Goldenberg, S. L. (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

[8] Bruchovsky,N., Klotz, L., Crook, J. and Goldenberg, S. L. (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

[9] Crook, J. M., O’Callaghan, C. J., Duncan, G., Dearnaley, D. P., Higano, C. S., Horwitz, E. M., Frymire, E., Malone, S., Chin, J., Nabid, A., (2012) Intermittent androgen suppression for rising PSA level after radiotherapy. N. Engl. J. Med., 367, 895–903

[10] Hussain, M., Tangen, C. M., Berry, D. L., Higano, C. S., Crawford, E. D., Liu, G., Wilding, G., Prescott, S., Kanaga Sundaram, S., Small, E. J., (2013) Intermittent versus continuous androgen deprivation in prostate cancer. N. Engl. J. Med., 368, 1314–1325

[11] Morino, K., Hirata, Y., Tomioka, R., Kashima, H., Yamanishi, K., Hayashi, N., Egawa, S. and Aihara, K. (2015) Predicting disease progression from short biomarker series using expert advice algorithm. Sci. Rep., 5, 8953

[12] Hirata, Y., Morino, K., Akakura, K., Higano, C. S., Bruchovsky, N., Gambol, T., Hall, S., Tanaka, G. and Aihara, K. (2015) Intermittent androgen suppression: Estimating parameters for individual patients based on initial PSA data in response to androgen deprivation therapy. PLoS ONE, 10, e0130372

[13] Suzuki, T., Bruchovsky, N. and Aihara, K. (2010) Piecewise affine systems modelling for optimizing hormone therapy of prostate cancer. Philos Trans A Math Phys Eng Sci, 368, 5045–5059

[14] Hirata, Y., di Bernardo, M., Bruchovsky, N. and Aihara, K. (2010) Hybrid optimal scheduling for intermittent androgen suppression of prostate cancer. Chaos, 20, 045125

[15] Hirata, Y., Azuma, S. and Aihara, K. (2014) Model predictive control for optimally scheduling intermittent androgen suppression of prostate cancer. Methods, 67, 278–281

[16] Tanaka, G., Hirata, Y., Goldenberg, S. L., Bruchovsky, N. and Aihara, K. (2010) Mathematical modelling of prostate cancer growth and its application to hormone therapy. Philos Trans A Math Phys Eng Sci, 368, 5029–5044

[17] Kuramae, H., Hirata, Y., Bruchovsky, N., Aihara, K. and Suzuki, H. (2011) Nonlinear systems identification by combining regression with bootstrap resampling. Chaos, 21, 043121

[18] Hirata, Y., Akakura, K., Higano, C. S., Bruchovsky, N. and Aihara, K. (2012) Quantitative mathematical modeling of PSA dynamics of prostate cancer patients treated with intermittent androgen suppression. J. Mol. Cell Biol., 4, 127–132

[19] Guo, Q., Lu, Z., Hirata, Y. and Aihara, K. (2013) Parameter estimation and optimal scheduling algorithm for a mathematical model of intermittent androgen suppression therapy for prostate cancer. Chaos, 23, 043125

[20] Suzuki, T. and Aihara, K. (2013) Nonlinear system identification for prostate cancer and optimality of intermittent androgen suppression therapy. Math. Biosci., 245, 40–48

[21] Suzuki, Y., Sakai, D., Nomura, T., Hirata, Y. and Aihara, K. (2014) A new protocol for intermittent androgen suppression therapy of prostate cancer with unstable saddle-point dynamics. J. Theor. Biol., 350, 1–16

[22] Tao, Y., Guo, Q. and Aihara, K. (2015) A partial differential equation model and its reduction to an ordinary differential equation model for prostate tumor growth under intermittenthormone therapy. J. Math. Biol., 69, 817–838

[23] Cover, T. M. and Thomas, J. A. (1991) Elements of Information Theory. New York: Wiley Interscience Publication

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg