Personalized Cancer Therapy Design: Robustness vs. Optimality
Intermittent Androgen Suppression (IAS) is a treatment strategy for delaying or even preventing time to relapse of advanced prostate cancer. IAS consists of alternating cycles of therapy (in the form of androgen suppression) and off-treatment periods. The level of prostate specific antigen (PSA) in a patient’s serum is frequently monitored to determine when the patient will be taken off therapy and when therapy will resume. In spite of extensive recent clinical experience with IAS, the design of an ideal protocol for any given patient remains one of the main challenges associated with effectively implementing this therapy. We use a threshold-based policy for optimal IAS therapy design that is parameterized by lower and upper PSA threshold values and is associated with a cost metric that combines clinically relevant measures of therapy success. We apply Infinitesimal Perturbation Analysis (IPA) to a Stochastic Hybrid Automaton (SHA) model of prostate cancer evolution under IAS and derive unbiased estimators of the cost metric gradient with respect to various model and therapy parameters. These estimators are subsequently used for system analysis. By evaluating sensitivity estimates with respect to several model parameters, we identify critical parameters and demonstrate that relaxing the optimality condition in favor of increased robustness to modeling errors provides an alternative objective to therapy design for at least some patients.
Several recent attempts have been made to develop mathematical models that explain the progression of cancer in patients undergoing therapy so as to improve (and possibly optimize) the effectiveness of such therapy. As an example, prostate cancer is known to be a multistep process, and patients who evolve into a state of metastatic disease are usually submitted to hormone therapy in the form of continuous androgen suppression (CAS) . The initial response to CAS is frequently positive, leading to a significant decrease in tumor size; unfortunately, most patients eventually develop resistance and relapse.
Intermittent Androgen Suppression (IAS) is an alternative treatment strategy for delaying or even preventing time to relapse of advanced prostate cancer patients. IAS consists of alternating cycles of therapy (in the form of androgen suppression) and off-treatment periods. The level of prostate specific antigen (PSA) in a patient’s serum is frequently monitored to determine when the patient will be taken off therapy and when therapy will resume. In spite of extensive recent clinical experience with IAS, the design of an ideal protocol for any given patient remains one of the main challenges associated with effectively implementing this therapy .
Various works have aimed at addressing this challenge, and we briefly review some of them. In  a model is proposed in which prostate tumors are composed of two subpopulations of cancer cells, one that is sensitive to androgen suppression and another that is not, without directly addressing the issue of IAS therapy design. The authors in  modeled the evolution of a prostate tumor under IAS using a hybrid dynamical system approach and applied numerical bifurcation analysis to study the effect of different therapy protocols on tumor growth and time to relapse. In  a nonlinear model is developed to explain the competition between different cancer cell subpopulations, while in  a model based on switched ordinary differential equations is proposed. The authors in  developed a piecewise affine system model and formulated the problem of personalized prostate cancer treatment as an optimal control problem. Patient classification is performed in  using a feedback control system to model the prostate tumor under IAS, and in  this work is extended by deriving conditions for patient relapse.
Most of the existing models provide insights into the dynamics of prostate cancer evolution under androgen deprivation, but fail to address the issue of therapy design. Furthermore, previous works that suggest optimal treatment schemes by classifying patients into groups have been based on more manageable, albeit less accurate, approaches to nonlinear hybrid dynamical systems. Addressing this limitation, the authors in  recently proposed a nonlinear hybrid automaton model and performed -reachability analysis to identify patient-specific treatment schemes. However, this model did not account for noise and fluctuations inherently associated with cell population dynamics and monitoring of clinical data. In contrast, in  a hybrid model of tumor growth under IAS therapy is developed that incorporated stochastic effects, but is not used for personalized therapy design.
A first attempt to define optimal personalized IAS therapy schemes by applying Infinitesimal Perturbation Analysis (IPA) to stochastic models of prostate cancer evolution was reported in . An IPA-driven gradient-based optimization algorithm was subsequently implemented in  to adaptively adjust controllable therapy settings so as to improve IAS therapy outcomes. The advantages of these IPA-based approaches stem from the fact that IPA efficiently yields sensitivities with respect to controllable parameters in a therapy (i.e., control policy), which is arguably the ultimate goal of personalized therapy design. More generally, however, IPA yields sensitivity estimates with respect to various model parameters from actual data, thus allowing critical parameters to be differentiated from others that are not.
In this paper we build upon the IPA-based methodology from  and  and focus on the importance of accurate modeling in conjunction with optimal therapy design. In particular, by evaluating sensitivity estimates with respect to several model parameters, we identify critical parameters and verify the extent to which the model from  is robust to them. From a practical perspective, the goal of this paper is to use IPA to explore the tradeoff between system optimality and robustness (or, equivalently, fragility), thus providing valuable insights on modeling and control of cancer progression. Assuming that an underlying, and most likely poorly understood, equilibrium of cancer cell subpopulation dynamics exists at suboptimal therapy settings, we verify that relaxing the optimality condition in favor of increased robustness to modeling errors provides an alternative objective to therapy design for at least some patients.
In Section 2 we present a Stochastic Hybrid Automaton (SHA) model of prostate cancer evolution, along with a threshold-based policy for optimal IAS therapy design. Section 3 reviews a general framework of IPA based on which we derive unbiased IPA estimators for system analysis. In Section 4 we evaluate sensitivity estimates with respect to several model parameters, identifying critical parameters and verifying the extent to which our SHA model is robust to them. We include final remarks in Section 5.
2 Problem Formulation
2.1 Stochastic Model of Prostate Cancer Evolution
We consider a system composed of a prostate tumor under IAS therapy, which is modeled as a Stochastic Hybrid Automaton (SHA). Details of the problem formulation are given in , but a condensed description of the SHA modeling framework is included here so as to make this paper as self-contained as possible. By adopting a standard SHA definition , a SHA model of prostate cancer evolution is defined in terms of the following:
A discrete state set , where (, respectively) is the on-treatment (off-treatment, respectively) operational mode of the system. IAS therapy is temporarily suspended when the size of the prostate tumor decreases by a predetermined desirable amount. The reduction in the size of the tumor is estimated in terms of the patient’s prostate specific antigen (PSA) level, a biomarker commonly used for monitoring the outcome of hormone therapy. In this context, therapy is suspended when a patient’s PSA level reaches a lower threshold value, and reinstated once the size of cancer cell populations has increased considerably, i.e., once the patient’s PSA level reaches an upper threshold value.
A state space , defined in terms of the biomarkers commonly monitored during IAS therapy, as well as “clock” state variables that measure the time spent by the system in each discrete state. We assume that prostate tumors are composed of two coexisting subpopulations of cancer cells, Hormone Sensitive Cells (HSCs) and Castration Resistant Cells (CRCs), and thus define a state vector with , such that is the total population of HSCs, is the total population of CRCs, and is the concentration of androgen in the serum. Prostate cancer cells secrete high levels of PSA, hence a common assumption is that the serum PSA concentration can be modeled as a linear combination of the cancer cell subpopulations. It is also frequently assumed that both HSCs and CRCs secrete PSA equivalently , and in this work we adopt these assumptions. Finally, we define variable , , where (, respectively) is the “clock” state variable corresponding to the time when the system is in state (, respectively), and is reset to zero every time a state transition occurs. Setting , the complete state vector is .
An admissible control set , such that the control is defined, at any time , as:
This is a simple form of hysteresis control to ensure that androgen deprivation will be suspended whenever a patient’s PSA level drops below a minimum threshold value, and that treatment will resume once the patient’s PSA level reaches a maximum threshold value. To this end, IAS therapy is viewed as a controlled process characterized by two parameters: , where is the lower threshold value of the patient’s PSA level, and is the upper threshold value of the patient’s PSA level, with . An illustrative representation of such threshold-based IAS therapy scheme is depicted in Fig. 1. Simulation driven by clinical data , was performed to generate the plot in Fig. 1, which shows a typical profile of PSA level variations along several treatment cycles.
An event set , where corresponds to the condition (i.e., ) and corresponds to the condition (i.e., ), where the notation indicates the time instant immediately preceding time .
System dynamics describing the evolution of continuous state variables over time, as well as the rules for discrete state transitions. The continuous (time-driven) dynamics capture the prostate cancer cell population dynamics, which are defined in terms of their proliferation, apoptosis, and conversion rates. As in , we incorporate stochastic effects into the deterministic model from  as follows:
where and are the HSC proliferation constant and CRC proliferation constant, respectively; and are the HSC apoptosis constant and CRC apoptosis constant, respectively; through are HSC proliferation and apoptosis exponential constants; is the HSC to CRC conversion constant; corresponds to the patient-specific androgen constant; is the androgen degradation constant; is the HSC basal degradation rate; and are the HSC basal production rate and androgen basal production rate, respectively. Finally, , , are stochastic processes which we allow to have arbitrary characteristics and only assume them to be piecewise continuous w.p. 1. The processes , , represent noise and fluctuations inherently associated with cell population dynamics, while reflects randomness associated with monitoring clinical data, more specifically, with monitoring the patient’s androgen level.
It is clear from (2)-(4) that and are dependent on , whose dynamics are affected by mode transitions. To make explicit the dependence of and on the discrete state (mode) , we let be the time of occurrence of the th event (of any type), and denote the state dynamics over any interevent interval as
We include as an argument to stress the dependence of the event times on the controllable parameters, but we will subsequently drop this for ease of notation as long as no confusion arises.
We thus start by assuming for . Solving (4) yields, for ,
It is then possible to define, for ,
where, for notational simplicity, we let
Next, let for , so that (4) implies that, for ,
Similarly as above, we define, for ,
It is then possible to rewrite (4) as follows:
Although we include as an argument in (7 ) and (9) to stress the dependence on the stochastic process, we will subsequently drop this for ease of notation as long as no confusion arises. Hence, substituting (7) and (9) into (2)-(3), yields
The discrete (event-driven) dynamics are dictated by the occurrence of events that cause state transitions. Based on the event set we have defined, the occurrence of results in a transition from to and the occurrence of results in a transition from to .
2.2 IAS Sensitivity Analysis
Recall that the main goal of this work is to perform sensitivity analysis of (2)-(6) in order to identify critical model parameters and verify the extent to which the SHA model of prostate cancer evolution is robust to them. Of note, several potentially critical parameters exist in the SHA model from . This work is a first step towards analyzing their relative importance, in which we select a subset of all model parameters in order to illustrate the applicability of our IPA-based methodology. The parameters we consider here are and (HSC proliferation constant and CRC proliferation constant, respectively), as well as and (HSC apoptosis constant and CRC apoptosis constant, respectively). These constants are intrinsically related to the cancer cell subpopulations’ net growth rate, whose value dictates how fast the PSA threshold values will be reached, and ultimately how soon treatment will be suspended or reinstated. As a result, correctly estimating the values of and , , is presumably crucial for the purposes of personalized IAS therapy design.
In this context, we define an extended parameter vector , where (, respectively) corresponds to the lower (upper, respectively) threshold value of the patient’s PSA level, (, respectively) corresponds to the HSC (CRC, respectively) proliferation constant, and (, respectively) corresponds to the HSC (CRC, respectively) apoptosis constant.
Within the SHA framework presented above, an IAS therapy can be viewed as a controlled process characterized by the parameter vector , as in (1), whose effect can be quantified in terms of performance metrics of the form . Of note, only the first two elements in vector are controllable, while the remaining parameters are not.
As in , here we make use of a sample function defined in terms of complementary measures of therapy success. In particular, we consider the most adequate IAS treatment schemes to be those that ensure PSA levels are kept as low as possible; reduce the frequency of on and off-treatment cycles. From a practical perspective, translates into the ability to successfully keep the size of cancer cell populations under control, which is directly influenced by the duration of the on and off-treatment periods. On the other hand, aims at reducing the duration of on-treatment periods, thus decreasing the exposure of patients to medication and their side effects, and consequently improving the patients’ quality of life throughout the treatment. Clearly there is a trade-off between keeping tumor growth under control and the cost associated with the corresponding IAS therapy. The latter is related to the duration of the therapy and could potentially include fixed set up costs incurred when therapy is reinstated. For simplicity, we disconsider fixed set up costs and take to be linearly proportional to the length of the on-treatment cycles. Hence, we define our sample function as the sum of the average PSA level and the average duration of an on-treatment cycle over a fixed time interval . We also take into account that it may be desirable to design a therapy scheme which favors over (or vice-versa) and thus associate weight with and with , where . Finally, to ensure that the trade-off between and is captured appropriately, we normalize our sample function: we divide by the value of the patient’s PSA level at the start of the first on-treatment cycle (), and normalize by .
Recall that the total population size of prostate cancer cells is assumed to reflect the serum PSA concentration, and that we have defined clock variables which measure the time elapsed in each of the treatment modes, so that our sample function can be written as
where and are given initial conditions. We can then define the overall performance metric as
We note that it is not possible to derive a closed-form expression of without imposing limitations on the processes , . Nevertheless, by assuming only that , , are piecewise continuous w.p. 1, we can successfully apply the IPA methodology developed for general SHS in  and obtain an estimate of by evaluating the sample gradient . We will assume that the derivatives exist w.p. 1 for all . It is also simple to verify that is Lipschitz continuous for in . We will further assume that , , are stationary random processes over and that no two events can occur at the same time w.p. 1. Under these conditions, it has been shown in  that is an unbiased estimator of , . Hence, our goal is to compute the sample gradient using data extracted from a sample path of the system (e.g., by simulating a sample path of our SHA model using clinical data), and use this value as an estimate of .
3 Infinitesimal Perturbation Analysis
where is a set of discrete states; is a continuous state space; is a finite set of events; is a set of admissible controls; is a vector field, ; is a discrete state transition function, ; is a set defining an invariant condition (when this condition is violated at some , a transition must occur); is a set defining a guard condition, (when this condition is satisfied at some , a transition is allowed to occur); is a reset function, ; is an initial discrete state; is an initial continuous state.
Consider a sample path of the system over and denote the time of occurrence of the th event (of any type) by , where corresponds to the control parameter of interest. Although we use the notation to stress the dependency of the event time on the control parameter, we will subsequently use to indicate the time of occurrence of the th event where no confusion arises. In order to further simplify notation, we shall denote the state and event time derivatives with respect to parameter as and , respectively, for . Additionally, considering that the system is at some discrete mode during an interval , we will denote its time-driven dynamics over such interval as . It is shown in  that the state derivative satisfies
with the following boundary condition:
when is continuous in at . Otherwise,
where is the reset function defined in (14).
Knowledge of is, therefore, needed in order to evaluate (16). Following the framework in , there are three types of events for a general stochastic hybrid system: (i) Exogenous event. This type of event causes a discrete state transition which is independent of parameter and, as a result, . (ii) Endogenous event. In this case, there exists a continuously differentiable function such that , which leads to
where . (iii) Induced event. Such an event is triggered by the occurrence of another event at time and the expression of depends on the event time derivative of the triggering event () (details can be found in ).
Thus, IPA captures how changes in affect the event times and the state of the system. Since interesting performance metrics are usually expressed in terms of and , IPA can ultimately be used to infer the effect that a perturbation in will have on such metrics. We end this overview by returning to our problem of personalized prostate cancer therapy design and thus defining the derivatives of the states and and event times with respect to , , , , as follows:
In what follows, we derive the IPA state and event time derivatives for the events identified in the SHA model of prostate cancer progression.
3.1 State and Event Time Derivatives
We proceed by analyzing the state evolution of our SHA model of prostate cancer progression considering each of the states ( and ) and events ( and ) therein defined.
1. The system is in state over interevent time interval . Using (15) for , we obtain, for ,
From (10), we have , , , and
It is thus simple to verify that solving (15) for yields, for ,
In particular, at :
where , , and are given from (23).
Similarly for , we have from (11) that , , , and
Combining the last four equations and solving for yields, for ,
where , , , .
In particular, at :