Clinically relevant cancer chemotherapy dose scheduling via mixed-integer optimization

Computers and Chemical Engineering 33 (2009) 2042–2054

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Clinically relevant cancer chemotherapy dose scheduling via mixed-integer optimization John M. Harrold a,1 , Robert S. Parker a,b,∗ a b

University of Pittsburgh, Department of Chemical and Petroleum Engineering, United States University of Pittsburgh Cancer Institute, Molecular Therapeutics and Drug Discovery Program, United States

a r t i c l e

i n f o

Article history: Received 28 October 2008 Received in revised form 15 May 2009 Accepted 1 June 2009 Available online 11 June 2009 Keywords: Drug scheduling Cancer chemotherapy MILP MIQP Bilinear constraints Gompertz model

a b s t r a c t Cancer is a class of diseases characterized by an imbalance between cell proliferation and programmed cell death. Chemotherapy is commonly employed as a treatment by clinicians, who must deliver the agent on a schedule that balances treatment efficacy with the toxic side effects. Engineers have considered the development of drug administration schedules for simulated cancer patients constrained by pharmacokinetic (PK) and pharmacodynamic (PD) models. The results typically involve mathematically elegant solutions, although the clinical utility of such results is limited by the formulation of the problem as well as the level of abstraction. At issue is the common disconnect between solutions that are mathematically vs. clinically optimal. The focus of this work is to develop a methodology which can explicitly account for the constraints clinicians consider implicitly. To demonstrate the clinical relevance of this methodology two case studies were considered: a theoretical system from the literature and a preclinical mouse model. The problem formulation is accomplished in a mixed-integer programming framework that is capable of solving problems with complex objectives and constraints yielding results that are clinically relevant. © 2009 Elsevier Ltd. All rights reserved.

1. Introduction Cancer is the name given to a class of diseases characterized by an imbalance in cell proliferation and apoptosis, or programmed cell death (Slingerland & Tannock, 1998). If left untreated, cancer will typically metastasize and spread throughout the host organism leading to its death (Chambers & Hill, 1998). After a cancer diagnosis, accessible malignant tissues are removed by surgical methods (Kaufman & Chabner, 2001). Tumors that are not surgically resectable require an alternate treatment approach. Furthermore, it is assumed that once cancer has reached detectable sizes it has spread throughout the body, and a systemic form of treatment is needed (Clare, Nakhlis, & Panetta, 2000). Chemotherapy is the most commonly employed whole-body treatment; this systemic approach typically has effects on both healthy and diseased tissues. This creates a dichotomy for clinicians: how should treatment administration be designed such that it decreases the cancerous burden while maintaining tolerable levels of patient toxicity?

∗ Corresponding author at: 1249 Benedum Hall, Pittsburgh, PA 15261, United States. Tel.: +1 412 624 7364. E-mail address: [email protected] (R.S. Parker). 1 Current address: 548 Hochstetter Hall, Department of Pharmaceutical Sciences, State University of New York at Buffalo, Buffalo, NY 14269, United States. 0098-1354/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2009.06.005

Nominally, the optimal treatment schedule – here schedule is defined as the amount and frequency of drug delivered (Johnson, Monks, Hollingshead, & Sausville, 2001) – is the one found to be the most efficacious from the set evaluated during clinical trials. Preclinical trials in organisms with similar relevant physiology (e.g., mice) provide a starting point for evaluating dose schedules. New anticancer compounds that have been found to be efficacious in preclinical trials are then evaluated in clinical trials. During the early stages of clinical trials, patients that are nonresponsive to current treatments help to establish dose levels for new chemotherapeutic compounds in humans (Johnson et al., 2001). Phase I trials evaluate the toxicity of a compound on a given schedule, and Phase II trials establish the efficacy of the compound in humans. Several schedules may be considered and the schedule yielding the most promising results in a statistically controlled trial is considered optimal from a clinical perspective. New compounds or schedules found to be efficacious in Phase II trials are compared to the current standard of practice against specific forms of cancer in Phase III trials (Horstmann et al., 1991). For approved drugs, treatment is generally broken down into cycles. At the end of the first cycle toxicity is evaluated, and at the end of the second cycle efficacy is evaluated. This process is repeated, adjusting the schedule as needed, until an endpoint is met. Endpoints span the range from reduction of the cancer below measurable limits (termed complete clinical response, which may or may not be true cure) to disease progression. It is important to note that in the clinic, toxicity

J.M. Harrold, R.S. Parker / Computers and Chemical Engineering 33 (2009) 2042–2054

and efficacy drive treatment decisions (Cella, Peterman, Hudgens, Webster, & Socinski, 2003; Thigpen, 2003). The question addressed in this work is whether improved dosing decisions are possible based on measurements currently available to clinicians. In order to accomplish this, predictive models with clinically relevant measurements are needed to determine if the treatment is efficacious and to establish when toxicity will be encountered. This has been previously addressed by several authors considering macroscopic models of cancer (Afenya, 2001; Costa, Boldrini, & Bassanezi, 1995; Florian, Eiseman, & Parker, 2004; Kimmel & Gorlova, 2002; Ledzewicz & Schattler, 2002a, 2002b; Martin, 1992; Martin & Teo, 1994; Matveev & Savkin, 2002; Pereira, Pedreira, Pinho, Fernandes, & Sousa, 1990; Swan, 1986; Swierniak & Smieja, 2001). Swan utilized nonlinear macroscopic models to obtain homeostasis through continuous intravenous infusion (Swan, 1986). This method of treatment, continuous infusion, is prohibitively expensive (i.e., it requires the surgical implantation of an in-dwelling catheter) and can impact the patient quality of life. Panetta considered pulsed chemotherapy to characterize regimens which would or would not eliminate tumor burden in the absence of toxicity constraints (Panetta & Adam, 1995); the key drawback here is choice to neglect toxicity constraints, which inherently limit chemotherapeutic administration. Swierniak et al. considered a proliferation model which segregated the population into drug susceptible and drug resistant cells (Swierniak & Smieja, 2001). A periodic solution was developed with upper bounds placed on each dose and the total amount of drug administered over an interval. From a clinical perspective, the periodic solution is desirable to facilitate ease of implementation; however, an explicit (vs. implicit) toxicity formulation would more closely mirror clinical reality. The dynamic and algebraic constraints encountered in many chemotherapy treatment design problems motivates the use of optimal control theory as a solution engine. Afenya posed the problem as a minimum time problem by driving the tumor volume to a specified level as quickly as possible by finding the optimal time to switch between the maximum possible dose and no dose (Afenya., 2001). A more common problem found in the literature minimizes the tumor volume at a final time subject to toxicity constraints (Ledzewicz & Schattler, 2002a; Martin & Teo, 1994; Matveev & Savkin, 2002). These approaches provide an advantage, in that they can explicitly account for toxicity in the formulation. However, minimizing the tumor volume at a specified final time, while an intuitive objective, is not necessarily clinically relevant. In fact, focusing on the final state of the system completely ignores clinical practice and the importance of the dynamic trajectory of the patient. Clinicians are guided by toxicity and efficacy through regular patient feedback. Final treatment times are not known a priori. Rather, they are determined by disease response or progression. To aid clinicians in drug administration decisions, the aforementioned issues must be addressed.

2043

2.1. PK models Sets of linear ordinary differential equations (ODEs) are commonly used to describe drug PK. Compartmental modeling, for example, is a typical approach which treats drug concentration in tissues as a system of interconnected well-mixed compartments (D’Argenio & Schumitzky, 1997). The compartments are treated similar to continuous stirred tank reactors where the rate of drug leaving the compartment is proportional to the concentration of the drug in that compartment. However, it is also common to encounter nonlinearities such as saturating drug clearance; which can be included in the model if identified. Since model quality limits achievable performance (Morari & Zafiriou, 1989), it is desirable to construct the most detailed model possible. Physiologically based models (Fig. 1 (top)) require significant amounts of data to characterize, thereby leading to logistical and ethical concerns (e.g. taking tissue biopsies to evaluate drug concentration vs. time profiles). When constructed, such models are usually limited to animal descriptions (Florian et al., 2007). In contrast, macroscopic compartmental models (Fig. 1(bottom)) sacrifice detail for the ability to be informed from clinical measurements, which are typically less

2. Treatment design with macroscopic models Macroscopic pharmacokinetic (PK) and pharmacodynamic (PD) models are commonly used to characterize the quantities which are clinically measurable. PK models describe the distribution and metabolism of the drug, while PD models characterize the effects of the drug on the organism (Cobelli & Carson, 2001). PK models can include the distribution of drugs from the point of administration to the peripheral tissues, the elimination of the drug, and the amount of drug present at the site of action (Cobelli & Carson, 2001). The site of action under consideration here is the cancerous tissue and drug levels are coupled to the administration of the drug by the circulatory system.

Fig. 1. Two common PK models structures: physiologically based (top) and loworder compartmental (bottom). Doses can be administered intravenous (IV) or orally (PO).

2044

J.M. Harrold, R.S. Parker / Computers and Chemical Engineering 33 (2009) 2042–2054

invasive (e.g. plasma drug concentrations, tumor volumes established using imaging techniques, etc. (Lowe et al., 2000; Mukherji, Schmalfuss, Castelijns, & Mancuso, 2004)). The case studies below will introduce two drug- and system-specific PK models. 2.2. PD models While PK models describe drug distribution, PD models describe the effects of the drug. PD effects can be considered positive, such as reductions in tumor size (Martin, 1992), or negative, such as the suppression of the immune system (Barbolosi & Iliadis, 2001; Cui, 2002; Kloft, Wallin, Henningsson, Chatelut, & Karlsson, 2006; Latz, Karlsson, Rusthoven, Ghosh, & Johnson, 2006). PD models considered in this work are referred to as PK-driven models because the presence of the drug, as modeled by the PK model, is explicitly coupled to the PD effect. This modeling methodology is commonly employed in the literature (Cobelli & Carson, 2001; Ledzewicz & Schattler, 2002a; Martin & Teo, 1994) and it allows the rate of untreated cellular proliferation to be characterized separately from the effect of the drug. Cancer cells typically begin to proliferate in an exponential fashion, and in the absence of drugs this can be described as (Asachenkov, Marchuk, Mohler, & Zuev, 1994): N˙ e =

ln(2) Ne (t) e

(1)

The size of the cancerous mass, Ne , increases with a doubling time e . The size of the cancerous mass is measured experimentally as a volume, though this mass is often referred to in terms of the number of cells (106 cells ≈ 1 mm3 ) (Norton, 1988). Eventually the rate of proliferation slows as the number of cancer cells asymptotically approaches a plateau population. Cancerous masses reaching this stage of growth are typically described by the Gompertz equation (Norton, 1988):



N˙ g =





ln(g /N0 ) g 1 ln Ng (t) ln g Ng (t) ln(g /2N0 )



(2)

Here g is the pseudo-doubling time for the system, N0 is the initial tumor size for the mass under Gompertzian growth (Ng ), and g is the plateau population. The exponential model and Gompertzian tumor growth model are shown for comparison in Fig. 2. To account for drug effects, plasma drug concentration predictions can be obtained from PK models and used as approximations

for drug concentration at the tumor site. This is based on the rationale that tumors are well perfused and have permeable capillaries. The drug effect is then represented mathematically by adding a bilinear kill term to the tumor growth equation. The effect term is proportional to both the current size of the tumor and the concentration of the drug (Martin & Teo, 1994). The exponential and Gompertzian models are augmented in the following manner: N˙ e =

ln(2) Ne (t) − keff Ceff (t)Ne (t) e



N˙ g =



(3a)



ln(g /N0 ) g 1 ln Ng (t) ln g Ng (t) ln(g /2N0 )

 − keff Ceff (t)Ng (t)

(3b)

Here, Ceff (t) is the drug concentration at the site of the cancerous mass. This could be as simple as using the plasma drug concentration, as predicted from a PK model, or it might involve additional functionality specific to the problem, see Section 4. The proportionality constant for drug effectiveness is keff . 3. General problem statement Cancer chemotherapy administration design for simulated patients having linear PK and macroscopic PD descriptions with bilinear nonlinearities have been considered by several authors (Barbolosi & Iliadis, 2001; Ledzewicz & Schättler, 2002a; Martin, 1992). The most common approach has been to consider the fixed final time problem: min D(t)

s.t.

J = FJ (x(tf , D), (tf )) x˙ = Ft (x(t), D(t))

(4a) (4b)

Fi (x(t), D(t)) ≤ 0

(4c)

Fe (x(t), D(t)) = 0

(4d)

The functions FJ , Ft , Fi , and Fe represent arbitrary functions of their respective arguments which may or may not contain nonlinearities. The internal states and the drug administration profile are given by x(t) and D(t), respectively. The objective is to minimize some function, typically the number of cancerous cells or the volume of a tumor, at a final time, t = tf . This is accomplished while satisfying the dynamic constraints (4b) defined by the PK and PD of the drug (Martin, 1992; Ledzewicz & Schattler, 2002a). Also, inequality (4c) and equality (4d) constraints can be included to limit toxicity and account for logistical concerns. Considering the dynamic constraints, optimal control theory provides a convenient avenue for solving these types of problems. However, considering the objectives and constraints clinicians face, the above formulation (4) and its solution may not provide results likely to be implemented in the clinic. The remainder of this paper will focus on the analysis of two case studies. The first case study considers a theoretical system that has been thoroughly studied from an optimal control perspective (Martin & Teo, 1994). We develop a mixed-integer formulation to solve the optimal control problem, as well as a reformulation that provides a general structure for incorporating clinically relevant objectives. This latter methodology is then applied to a PK/PD model of a preclinical (i.e., animal) system. Before introducing the specific case studies, some theoretical underpinnings are used to reformulate the components of the cancer treatment problem. 4. Theory 4.1. PD variable transform

Fig. 2. Exponential (solid) and Gompertzian (dash-dot) tumor growth models for untreated cancer. To facilitate comparison, the normalized time, , is e for the exponential and g for the Gompertzian models, respectively.

Because of the bilinear term from the drug PD (either (3a) or (3b)) and the nonlinear nominal tumor growth in (3b), the

J.M. Harrold, R.S. Parker / Computers and Chemical Engineering 33 (2009) 2042–2054

2045

optimization in (4) is a nonlinear programming (NLP) problem. Nonlinear optimizations can lead to local minima, and efforts to eliminate nonlinearities can improve the likelihood of achieving the true optimum. The effect of the bilinear kill term can be included, but the explicit nonlinearity removed, by performing the following logarithmic transformation (Afenya., 2001): P(t) ˙ ˙ P(t) = ln(N(t)) ⇔ N(t) = eP(t) ⇔ N(t) = P(t)e

(5)

Notice that ln(N(t)) increases monotonically with N(t) such that the value of drug administration which minimizes N(tf ) will also minimize P(tf ). The transform in (5) can be applied to cancerous masses under both exponential and Gompertzian growth resulting in: ˙ P(t) =

ln(2) − keff Ceff (t) e



˙ P(t) =

(6a)



ln(g /N0 ) 1 ln (ln g − P(t)) − keff Ceff (t) g ln(g /2N0 )

(6b)

By logarithmically transforming the PD equations (3a) and (3b), cancer proliferating exponentially can be described by an ODE (6a) that is linear in effective drug concentration (Ceff (t)). Also PD models with bilinear kill terms in which cancer proliferates in a Gompertzian fashion can be reduced to an ODE (6b) which is linear in both effective drug concentration (Ceff (t)) and transformed tumor size (P(t)). 4.2. Constraint formulation

Fig. 3. Plasma PK profile for a drug with first order elimination illustrating the minimum therapeutic drug concentration, Cth (dashed line), the maximum tolerable drug concentration Cmax (dotted line), and the exposure (total exposure: area under the concentration vs. time curve, AUCexp ; effective exposure: shaded area, AUCeff ) over the dosing interval, t ∈ [t0 , tf ].

plasma concentration over the treatment interval:



tf

AUCexp =

C(t)dt ≤ Ccum

(7)

0

The constraints from the general problem statement (4) come from a variety of sources. The dynamic constraints (4b) represent the PK/PD models for drug distribution and host response. The inequality and equality constraints, (4c) and (4d), respectively, are derived primarily from three areas of concern: financial limitations, logistical restrictions, and toxicity constraints. Foregoing such restrictions, an algorithm to determine the optimal dosing regimen would deliver as much drug as possible as frequently as possible. While financial constraints inherently limit the available treatment methodologies, the present paper is not going to focus on public policy constraints and limitations induced by health plan providers. Instead, the present work focuses on toxicity and logistical treatment limitations, thereby formulating a clinically relevant best case outcome. This choice was made because the clinical best case, without the often-present financial limitations, must provide a treatment benefit vs. standard of practice; otherwise, there is little practical motivation to pursue algorithmic design of cancer chemotherapy under clinical constraints. Logistical constraints can be thought of as restrictions associated with managing treatment. For example, inconsistent drug administration times will lead to a patient who forgets to take his/her medication (Urquhart and De Klerk, 1998). Treatment requiring the direct attention of hospital staff, e.g., intravenous infusions, are fundamentally limited by the working schedules of the medical professionals. Similarly, drugs delivered in pill form are administered in discrete amounts, e.g., three pills every other day. Problem formulation and solution methodologies should be able to account for the discrete nature and limitations imposed by logistical constraints if the designs are to impact clinical practice. Toxicity constraints can be a function of reductions in the patient’s immune system performance, loss of body weight, or even pain experienced by the patient. It is common to characterize toxicity constraints in terms of the drug PK or measurable PD effects (Martin & Teo, 1994). Constraints related to drug PK are illustrated in Fig. 3. Placing a limit on exposure is a common toxicity constraint. Exposure, AUCexp is commonly calculated by integrating drug

Eq. (7) limits the toxicity of the drug to be less than an acceptable cumulative exposure, Ccum , the total area under the solid line in Fig. 3. Drugs may not become effective until a therapeutic plasma concentration is reached (Cth ). Once C(t) > Cth , the effective drug concentration (Ceff ) is that concentration above Cth . For such drugs, the effective exposure, AUCeff , is represented by the shaded region of the PK profile in Fig. 3. An effective drug concentration (Ceff ) can then be defined in terms of the therapeutic drug concentration (Cth ) in a piecewise fashion, as follows:



Ceff (t) =

0

C(t) ≤ Cth

C(t) − Cth

C(t) > Cth

(8)

Therefore, drug administration that does not increase the plasma concentration to at least Cth is ineffective but may still contribute to toxicity. The dose applied at time t3 shown in Fig. 3 contributes to the total exposure but not the effective exposure. Acute toxicity is reached when the drug plasma concentration exceeds some maximum, Cmax . This is a state constraint given by: C(t) ≤ Cmax

∀t ∈ [0, tf ]

(9)

In Fig. 3, this constraint is violated by the dose at t1 . It is also possible to bound the amount of drug that can be administered at any given point in time. This can be a result of discrete dosing (pills) or because the drug was found to be effective over a particular range of doses. In the latter case, a semicontinuous variable is encountered at each dosing opportunity, as follows: Dlb ≤ D(tq,i ) ≤ Dub or D(tq,i ) = 0 ∀i = 1, . . . , mq

(10)

Here Dlb and Dub are the lower and upper bounds of the continuous portion of the therapeutic dosing range. The upper bound is typically based on the maximum tolerated dose (MTD), either a direct consequence of (9) or the amount of drug that will produce chronic side effects. The lowest value of drug that has an observable effect defines the lower bound. Finally, tq,i is the ith dosing time, and tq is

2046

J.M. Harrold, R.S. Parker / Computers and Chemical Engineering 33 (2009) 2042–2054

the set of all possible times (q = 1, · · · , mq ) in which drugs may be administered during treatment.

of a binary variable:

4.2.1. Discrete formulations In order to cast the problem in the mixed-integer linear programming (MILP) framework, the constraints must be discretized. For the discrete formulations, k will denote the current time step with a system step size of h. Hence, k ∈ {1, . . . , mk }, and the final time step, mk , is defined as mk = tf /h. The subscriptd will be used to indicate the discrete form of continuous variables. For example, continuous states, x(t), and outputs such as plasma drug concentration, C(t), or logarithmically transformed tumor size, P(t), would have the discrete variable counterparts of xd (k), Cd (k) and Pd (k), respectively. The cumulative toxicity constraint (7) contains an integral term which can be approximated using the trapezoidal rule (Kreyszig, 1999), as follows:

When bd,u (q) is zero, D(q) must be zero. When the drug is to be administered, bd,u (q) must be one, and D(q) must lie between Dlb and Dub .



mk −1

h C (1) + h 2 d

Cd (k) +

h C (m ) ≤ Ccum 2 d k

Dlb bd,u (q) ≤ Dd (q) ≤ Dub bd,u (q) ∀q = 1, . . . , mq

(14)

5. Case study: intravenously administered drug with gompertzian proliferation The following theoretical system has been previously studied by Martin and Teo (1994) using optimal control theory, with solutions obtained via control vector parameterization: ˙ C(t) = D(t) − C(t)



N˙ g (t) =

(15a)





ln(g /N0 ) g 1 ln Ng (t) ln g Ng (t) ln(g /2N0 )

 − keff Ceff (t)Ng (t) (15b)

(11)

k=2

Discrete linear state inequality constraints, such as the acute toxicity constraint found in (9), are formulated by replacing the continuous variables with their discrete counterparts. Therefore (9) can be replaced with: Cd (k) ≤ Cmax

∀k = 1, . . . , mk

(12)

One issue which can arise from discretization is that the plasma concentration might be greater than Cmax between discretization points. This is not a problem for intravenously administered drugs described by first order kinetics. In this situation the peak associated with any given drug administration occurs at the start (bolus) or end (square-wave) of administration, which can be made to fall on a discretization point. More complicated PK models can lead to important dynamics occurring between discretization points, which can be most easily overcome by using suitably small values for h. The requirement that the effective drug concentration is nonzero only when the plasma concentration is greater than the therapeutic concentration creates a semicontinuous variable. This is written in a piecewise fashion in Eq. (8), but can be written as linear inequalities by introducing the binary variable bd,th (k) at each discrete time step, and requiring that the following inequalities be satisfied (Tyler and Morari, 1999):

Ceff (t) = (C(t) − Cth )H(C(t) − Cth )

(15c)

C(0) = C0

(15d)

Ng (0) = N0

(15e)

The PK of the drug is described by Eq. (15a) where the drug plasma concentration, C(t), increases with intravenous infusion of the drug, D(t), and decreases with first order elimination kinetics at a rate . The change in the number of cancer cells is described by Eq. (15b). The cancerous cells proliferate in a Gompertzian fashion described by g and g , as discussed in Section 2.2. The drug effect is proportional, with constant of proportionality keff , to the number of cancer cells, Ng (t), and the effective drug plasma concentration, Ceff (t), denoted as the drug concentration above the minimum therapeutic concentration, Cth . The initial drug concentration and number of cancer cells are given by C0 and N0 , respectively. 5.1. Constraint formulation

Cd (k) − Cth ≤ (Cmax − Cth )(1 − bd,th (k))

(13a)

0 ≤ Cd,eff (k)

(13b)

It was assumed that the drug could be administered weekly, and the selected final time was tf = 52 weeks (Martin & Teo, 1994). An upper bound was placed on the plasma drug concentration (9) and on the acute exposure (7). The final constraint placed on the system dealt with efficacy. Since it was undesirable for the status of a patient to decrease (e.g., tumor burden to increase), the number of cancer cells was not allowed to increase to a number larger than the initial condition:

(Cmax − Cth )(1 − bd,th (k)) ≥ Cd,eff (k)

(13c)

Ng (t) ≤ N0

Cth − Cd (k) ≤ Cth bd,th (k)

(13d)

0 ≤ Cd,eff (k) − (Cd (k) − Cth )

(13e)

Cth bd,th (k) ≥ Cd,eff (k) − (Cd (k) − Cth )

(13f)

Consider the situation when the drug plasma concentration is below the therapeutic value (Cd (k) ≤ Cth ). In this case, the left-hand side of (13d) will be greater than zero, which requires bd,th (k) to be one. This forces the left-hand side of (13c) to be zero, resulting in Cd,eff (k) being zero also. The left-hand side of (13a) will be negative, so this constraint will be satisfied. The right-hand side of constraints (13e) and (13f) will be positive on the range [0, Cth ], so these constraints will also be satisfied. A similar series of constraint satisfaction conditions can be shown for Cd (k) > Cth . The therapeutic range administration constraint, (10), can be transformed into a set of linear inequalities by a similar introduction

∀t

(16)

5.2. Optimal control problem The optimal control problem considered by Martin and Teo (1994) was given as: min Ng (tf ) D(t)

s.t.

(17)

(7), (9), (15), (16)

This equation had nonlinear growth and death terms as well as the discontinuity associated with the therapeutic drug concentration (Cth ). 5.3. MILP problem reformulation The first step in the reformulation involved the logarithmic transformation from Eq. (5). By performing this transformation on

J.M. Harrold, R.S. Parker / Computers and Chemical Engineering 33 (2009) 2042–2054

2047

the PD Eq. (15b), the transformed Eq. (6b) results, and the Ceff component of (6b) would be defined in the folloing manner:

Table 1 PK/PD parameters taken from the Martin study (Martin & Teo, 1994). Here [D] are the units of drug concentration.

Ceff (t) = (C(t) − Cth )H(C(t) − Cth )

Parameter

Value

g

5

months

g

1012

cells

N0

1010

cells

keff,1

2.7 ×10−2

1 days·[D]

keff,2

8.4 ×10−3

1 days·[D]

keff,3

1.5 ×10−3

1 days·[D]



.27

1 days

Cth

10

[D]

Cmax

50

[D]

Ccum

4.1 ×103

[D]·days

tf

364

days

(18)

The discontinuous effective drug concentration due to the presence of the threshold, Cth , was accounted for during discretization, where the dynamic constraints were converted into algebraic constraints at each time step. The plasma drug concentration was discretized at a time step h (Chen, 1984), resulting in a piecewise continuous function (19):

⎧ ⎨ − 1 (e−h − 1) Dd (q) + e−h C (k) when k = hq d  h (19) Cd (k + 1) = ⎩ −h e

Cd (k)

otherwise

The input is present along with a decay term at time steps that coincide with drug administration times. When no drug can be administered, the plasma concentration simply decays exponentially with a time constant 1/. The discontinuous effective drug concentration, Ceff (t), from Eq. (15c) was accounted for using the methodology discussed in Section 4.2.1. By applying the binary variable bd,th (k), and enforcing the constraints found in Eq. (13) at each time step, k, the discontinuous effective drug concentration was reduced to linear inequalities. The logarithmically transformed PD (6b) was discretized using Euler’s method (Kreyszig, 1999): Pd (k + 1) = Pd (k) + hFd (Pd (k), Cd,eff (k))

Units

(20)

Here Fd is the right-hand side of Eq. (6b), and Cd,eff are the discrete values of the continuous Ceff defined in (18). The continuous form of the acute toxicity constraint (9) was replaced with the discrete form (12), and cumulative toxicity was treated in the discrete domain using Eq. (11) in place of (7). The efficacy constraint (16) was incorporated as follows: Pd (k) ≤ ln(N0 );

∀k = 3, . . . , mk

(21)

It is important to note that the efficacy constraint was not enforced for the first two time steps. The nature of discrete systems without direct feedthrough dictates that the effects of an input change on states will not manifest until subsequent time steps. In this case, manipulated variable changes directly affect the plasma drug concentration. These manipulated variable changes then indirectly affect the number of cancer cells. Since the cells are continuously proliferating, the continuous-time efficacy constraint cannot be satisfied until the third time step. Based on the discretization results above, the continuous problem (17) was recast as the following MILP: min Pd (mk )

Dd (q)

(22)

s.t. (11), (12), (13), (19), (20), (21) 5.4. Results The parameters considered for this case study (Martin & Teo, 1994) are provided in Table 1, and a time step of h = 1 day was used. Three values of keff were considered representing different levels of drug efficacy (Martin & Teo, 1994). Highly effective, marginally effective, and ineffective drugs were represented by keff,1 , keff,2 and keff,3 , respectively. Solutions for each case were calculated and are shown in Fig. 4. The objective function values (N(tf )) for both the MILP and optimal control solutions are shown for each value of keff in Table 2. Problem formulations in GAMS contained 4106 equations, 1873 continuous variables and 364 discrete variables, and solution times in CPLEX were all less than a second on a dual AMD Athlon 1.8 GHz machine with one GB of RAM.

Fig. 4. Response to treatment. Top pane: tumor volume predictions for the three drug types: keff,1 (solid), keff,2 (dash-dot) and keff,3 (dotted). Remaining panes: dose schedule (bar/shaded region) and concentration (solid) for drugs keff,1 (2nd), keff,2 (3rd) and keff,3 (bottom); Cth (dashed) and Cmax (dotted) in panes 2–4.

The solutions to (22), for the three different efficacy levels given in Table 1, shown in Fig. 4 are qualitatively similar to those presented by Martin and Teo (1994). An initial dose was administered at the first time step in all instances. This was done to accommodate the efficacy constraint (21) and drives the number of cancer cells down initially. As the cancer population approached the initial Table 2 Objective function values for the MILP and optimal control (OC) solutions (Martin & Teo, 1994).

MILP OC

keff,1

keff,2

keff,3

<1 <1

977 1.2 × 103

1.9 × 109 1.8 × 109

2048

J.M. Harrold, R.S. Parker / Computers and Chemical Engineering 33 (2009) 2042–2054

number of cells, more drug was administered to satisfy the efficacy constraint. This was most evident for drug efficacies keff,3 and keff,2 . At the end of the treatment cycle, the drug was administered in large amounts to minimize tumor volume at tf such that the cumulative toxicity constraint (11) was met and the acute toxicity constraint (12) was not violated. When considering the highly effective drug, (keff,1 ), the optimal control and MILP solutions achieved essentially the same result. It should be noted that the absolute value predicted by the MILP solution was smaller than that predicted by the optimal control solution for the moderately effective drug (keff,2 ). This was likely due to numerical errors resulting from the discretization. However, for the moderately effective drug the MILP and optimal control solutions were clinically indistinguishable. There were noticeable differences between the two methodologies when considering the ineffective drug (keff,3 ). Each profile here fails to accommodate the objectives of treatment from a clinical perspective. Recall that clinicians use the first cycle of treatment to establish patient specific toxicity and the second cycle of treatment to titrate treatment in order to maximize efficacy. Asking a clinician to postpone treatment in order to minimize tumor volume at some arbitrary future time ignores real aspects of treatment (such as metastasis). It is important that the mathematical objectives reflect these clinical objectives, namely the rapid reduction and elimination of the tumor volume. 5.5. Objective reformulation: trajectory tracking Clinicians make use of periodic feedback to evaluate the efficacy of treatment and to adjust treatment as necessary to mitigate side effects. To approach the cancer chemotherapy dosing problem from a more practical perspective, the problem from problem (22) was reformulated as a receding horizon problem: min

Dd (q)

8  w=1

(Pd (w) − Td (w))2 + u

mq 

Dd (q)

(23)

q=1

s.t. (11), (13), (19), (20), (9), (21) This mixed integer quadratic programming problem (MIQP) minimizes the deviations between the transformed tumor volume, Pd (w), and a specified target, Td (w), every week w. The input penalty term, weighted by u , was added to penalize small drug doses. This formulation assumes the patient returns every eight weeks for evaluation and the prediction horizon for optimization purposes decreases by the number of treatment periods preceding the current period (a shrinking horizon formulation (Liotta, Georgakis, & El-Aasser, 1997; Thomas, Joseph, & Kardos, 1997)). The problem in (23) was modeled in GAMS and solved using CPLEX. This resulted in a series of 6 optimizations. Each optimization contained 1894 equations, 865 continuous variables, and 168 discrete variables, and the solution was found in less than a second. The results for the most efficacious drug, keff,1 , are shown in Fig. 5. The trajectory specified in Fig. 5 represents the desire of the clinician to decrease the tumor burden as quickly as possible to reduce the possibility of metastasis. Furthermore, the use of large doses, when drug was administered, reduces the chance of tumor cells developing resistance (as might be encountered for long-duration low-dose infusions). For the two cases presented here, u ∈ {0, 1}, a large amount of drug is administered during the first cycle. This rapidly eliminates a large number of the cancerous cells. With u = 0, the focus was on adhering to the trajectory and resulted in many small doses being administered. While dose levels were predicted which led to immediate plasma levels at or below the therapeutic level, Cth , these administrations did combine with subsequent administrations to increase the effective exposure.

Fig. 5. Receding horizon solution for keff,1 . Top pane: predicted tumor response profile, Pd , and desired response, Td , (×) for u = 0 (solid) and u = 1 (dotted). Other panes: drug administration levels (bar/shaded region) and plasma concentrations (solid) for u = 0 (middle pane) and u = 1 (bottom pane); Cth (dashed) and Cmax (dotted) in the middle and bottom panes.

Table 3 Total exposure, AUCexp , and effective exposure, AUCeff , for solutions to (23) and the given input weights, u .

AUCexp AUCeff AUCeff AUCexp

u = 0

u = 1

units

3438 1110

2069 1061

[d]·days [d]·days

0.32

0.54

—

While this result is mathematically optimal, it is suboptimal from a clinical perspective. A significant amount of drug was being administered that contributed to the cumulative exposure, while contributing little to the efficacy of the treatment. By increasing the penalty for drug administration to u = 1, the administration profile in the bottom pane of Fig. 5 was obtained. This tracked the trajectory well as indicated by the dotted line in the top pane of Fig. 5, with a greater fraction of the drug exposure above Cth . Table 3 contains the total exposure, AUCexp , the effective exposure, AUCeff , and the ratio of the effective to total exposure. By changing u from zero to one the fraction of the effective drug exposure increases from a third to over one half. 6. Case study: preclinical mouse model The system considered in the previous section was a theoretical construct created to represent a class of chemotherapy problems. This section will take the methodologies applied to the theoretical system and demonstrate that they can be extended to a nonlinear model of an accepted preclinical animal model. The system under study is the treatment of severe combined immunodeficient (SCID) mice bearing subcutaneous HT29 human colon carcinoma

J.M. Harrold, R.S. Parker / Computers and Chemical Engineering 33 (2009) 2042–2054

xenografts. Mice were treated with 9-nitrocamptothecin (9NC), a member of the camptothecin family that is readily metabolized into 9-aminocamptothecin (9AC). Both 9NC and 9AC exist in plasma in the active lactone and inactive carboxylate forms. The mechanism of action for the active forms of 9NC and 9AC is the inhibition of topoisomerase-I (Pommier, Pourquier, Urasaki, Wu, & Laco, 1999), an enzyme required for DNA synthesis.

The plasma disposition of the total amounts of 9NC and 9AC is given by the following nonlinear PK model (Harrold, 2005): (24a)

1 1 x˙ 2 (t) = Du (t) − x2 (t) u u

(24b)

1 1 x˙ 3 (t) = x4 (t) − x3 (t) l l

(24c)

x˙ 4 (t) =

1 1 (D (t) + x5 (t)) − x4 (t) l l l

(24d)

x˙ 5 (t) =

ˇ1 1 x3 (t) − x5 (t) r r

(24e)

x˙ 6 (t) =

Dlast D ˛2 (x1 (t − ) + ˛1 x3 (t − )) × H(t − ) − last x6 (t) nl nl (24f)

CNC = x1 (t)

Ku + Kl ˛1 x3 (t) Dlast

CAC = ˇ2 Ka (x1 (t) + ˛1 x3 (t)) + Ka x6 (t)

(24g) (24h)

Here Dlast is the value of the most recently administered nonzero dose in mg/kg,  represents a time delay, and the parameters ˛i and ˇi are complementary fractions. The concentrations of 9NC and 9AC are given by CNC and CAC , respectively. This model is characterized by two separate pathways which are dependent on the magnitude of the dose administered. Eqs. (24c) through (24f) define lower dose dynamics, and at higher doses, all of the states are active. This switching is accomplished by defining Dl (t) and Du (t) in terms of the dose mass, D(t), and a threshold, DT :



Dl (t) =

D(t)

D(t) ≤ DT

DT

D(t) > DT

 Du (t) =

0

D(t) ≤ DT

D(t) − DT

D(t) > DT

˙ N(t) =

ln(2) N(t) − keff Ceff (t)N(t) e (N(t))



e (N) =

e,f ,

N ≤ Nth

e,s ,

N > Nth

Ceff (t) = CNC (t) + CAC (t)

mm3 days days ml min−1 ng Camptothecin−1

(25b)

ng–9NC min min min min — — ml−1 ml−1 ml−1 min

(26a)

(26b) (26c)

If left untreated, the tumor size, N(t), will increase according to the first term on the right-hand side of Eq. (26a). This form of untreated tumor growth is similar to that suggested by Simeoni et al. (2004). Initially, cancer proliferates exponentially, and the rate of proliferation slows as the tumor mass increases. This behavior is approximated mathematically by defining two different regions of exponential growth (the latter rate slower than the initial rate). In (26), the tumor proliferates at a rate e,f until it reaches a tumor size of Nth , when the rate of proliferation slows to e,s . Gompertzian growth has been used to characterize this range of growth as well (Norton, 1988). However, the data required to resolve the plateau population in the Gompertz model is frequently unavailable in (pre)clinical studies. The bilinear kill term in equation (26a) utilizes the effective drug concentration (26c), which is defined as the sum of the total forms of both 9NC, predicted using Eq. (24g), and 9AC, predicted using Eq. (24h). The parameters used in (26) are given in Table 5. A second PD model was used to characterize the toxic effects of treatment by quantifying reductions in mouse body weight, as follows: B(t) = g(t) − e(t) + B0

(27a)

˙ g(t) =

kg 1 − g(t) g g

(27b)

˙ e(t) =

ke 1 AUClast − e(t) e e

(27c)

Btox (t) = B(t) −

Table 4 Parameters used in the nonlinear PK model, ((24) and (25)), of 9NC and 9AC total. 9.078 ×103 12.36 1.482 0.184 107.98 0.064 0.468 0.216 0.162 0.0491 49.20

477.36 9.56 13.49 3.49 × 10−5

Nth e,f e,s keff

(25a)

The parameters used in the PK model (24) and (25) are given in Table 4. The effects of 9NC were segregated into two pharmacodynamic (PD) models. The first PD model describes the efficacy of the drug,

DT u l r nl ˛1 ˛2 Knl Kl Ka 

Table 5 Parameters used for the efficacy model (26).

using a structure similar to the treatment for cancer under exponential growth (3a). In the present case, the rate of tumor volume change is defined as follows:

6.1. Discontinuous PD model driven by nonlinear PK

1 1 x2 (t) − x1 (t) x˙ 1 (t) = u u

2049

N(t) 1000

(27d)

Animals are provided with ample food and water, and changes in body weight are characterized by Eq. (27a). In the absence of treatment, the animal will grow according to Eq. (27b). During treatment, the body weight decrease is governed by Eq. (27c). The negative effects of the drug are described in terms of the exposure from the most recent dose administered:



t=end of day

AUClast =

(CNC (t) + CAC (t))dt

(28)

t=dosing time

The exposure is the cumulative area under the concentration vs. time curves (AUC) for 9NC and 9AC as predicted by (24). The corrected body weight, Btox (t), is calculated in Eq. (27d) by subtracting Table 6 Parameters for toxicity model (27). e ke g kg

2.85 ×105 8.00 ×10−3 3.45 ×104 11.69

min−1 g-mouse min−1 (ng(9NC +9AC) min)−1

2050

J.M. Harrold, R.S. Parker / Computers and Chemical Engineering 33 (2009) 2042–2054

⎧ ln 2 ⎪ ⎨  , P(t) < ln Nth

the mass of the tumor from the measured body weight, B(t). The mass of the tumor is calculated by assuming a density equal to that of water (Sickles, 1989). The parameters describing changes in body weight in the absence and presence of treatment are shown in Table 6.

Fs (P) =

6.2. Constraint formulation

Here, Eqs. (31a) and (31b) can be used as a substitute for (26a) and (26b), respectively. The logarithmically transformed PD Eq. (31) can then be discretized using a step size of h and Euler integration resulting in (32):

Constraints considered for this system relate to the logistics of treatment, the toxicity of 9NC and 9AC, and the need for efficacy of an administered dose. The laboratory personnel performing the efficacy experiments work Monday through Friday and administer drugs once per day. This defines the number and time of possible dosing opportunities, (q). The parameter Qmap (q) was used as a map between a dosing opportunity, q, and the corresponding timestep, k, in which it occurs. For example, the third dosing opportunity, q = 3, with a time step of one minute, occurs at k = 2881, so Qmap (3) = 2881. Acute drug toxicity was accounted for by placing an upper bound on the amount of drug that can be administered. The maximum dose, Dub = 1.0 mg/kg, was found to be tolerable for successive doses. The administration of 9NC causes diarrhea which may lead to fewer nutrients being absorbed from the gastrointestinal tract. These effects are modeled as shown above in terms of body weight reduction. Animals below a critical weight, Bmin , would be removed from study. In order to keep all simulated animals on study, a lower bound was placed on the corrected body weight: Btox (t) ≥ Bmin

(29)

Efficacy was accounted for in the constraints by placing a lower bound on the nonzero amount of drug which could be administered. The minimum effective dose, Dlb , was found in preclinical studies to be 0.44 mg/kg. The dosing decision variable was then constrained by Eq. (10) with the aforementioned bounds.

e,f

⎪ ⎩ ln 2 , P(t) > ln Nth e,s

Pd (k + 1) = Pd (k) + h(Fs,d (Pd (k)) − keff Ceff,d (k))

Fs,d (P(k)) =

min D(t)

mT 

N(Q (T ))

(30a)

T =1

s.t. Dynamics:(24), (25), (26), (27)

(30b)

Toxicity:(29)

(30c)

Dose:(10)

(30d)

The parameters, Q , represent the times at which the tumor size is being minimized. The set of minimization points, T , has a length of mT . The formulation found in (30) is a mixed integer nonlinear programming problems (MINLP). To reduce the complexity, and guarantee an optimal solution, the problem was reformulated as a MILP. 6.4. Problem reformulation 6.4.1. PD variable transform The effect of the bilinear kill term from (26a) can be included, but the explicit nonlinearity removed, by performing the transformation found in Eq. (5). The resulting efficacy equations are: ˙ P(t) = Fs (P(t)) − keff Ceff (t)

(31a)

⎧ ln 2 ⎪ ⎨  , Pd (k) < ln Nth e,f

⎪ ⎩ ln 2 , Pd (k) > ln Nth

(32a)

(32b)

e,s

The piecewise continuous function, Fs,d , was formulated mathematically by introducing the binary variable bP,d (k). This resulted in (32b) being replaced by the following Big-M constraints (Tyler & Morari, 1999): Pd (k) − ln Nth ≤ (ln(Nmax (k)) − ln(Nth )) × (1 − bP,d (k))

(33a)

(1 − bP,d (k))Lf ≤ Fs,d (k) −

ln(2) e,f

(33b)

(1 − bP,d (k))Uf ≥ Fs,d (k) −

ln(2) e,f

(33c)

ln Nth − Pd (k) ≤ ln(Nth )bP,d (k)

(33d)

bP,d (k)Ls ≤ Fs,d (k) −

ln(2) e,s

(33e)

bP,d (k)Us ≥ Fs,d (k) −

ln(2) e,s

(33f)

Lf = Uf =

ln(2) ln(2) − e,s e,f

(33g)

Ls = Us =

ln(2) ln(2) − e,f e,s

(33h)

6.3. Problem formulation The objective function here minimizes the tumor volume at specified points. As was demonstrated in Section 5.5, specifying a trajectory ensures efficacious treatment throughout the treatment cycle. In continuous time, the problem was formulated as the following optimization:

(31b)

The logical justification for these constraints is similar to that found above in (13). Here Nmax (k) is the value of N at any given step k in which no drug is administered (e.g., the maximum possible value the tumor volume can take at each time step). In the case where N ≥ Nth , the left-hand side of Eq. (33a) is positive, thereby forcing bP,d (k) to be zero. This forces left-hand side of Eqs. (33e) and (33f) to be zero and ensures that Fs,d is assigned to the slower growth rate. When N ≤ Nth , the left and side of Eq. (33d) is positive, and, therefore, bP,d (k) must be one. This forces the left-hand side of Eqs. (33b) and (33c) to be zero which results in Fs,d taking on the value of the faster growth rate. 6.4.2. Dosing bounds and parameterized PK The constraints from Eq. (10) were incorporated as discussed above in the development of Eq. (14) by introducing the binary variable bd,u (k). Nonlinearities such as those found in the PK model (24a)–(24f) can complicate optimizations. To eliminate these nonlinearities, it was observed that neither 9NC nor 9AC were observable 1 day after a dose of 9NC. While theoretically the range of non-zero doses is continuous (0.44–1.0 mg/kg), realistically there is a finite set of doses, l = 1, · · · , ml , which can be distinguished clinically from one another. The first value of l corresponds to Dlb and the last value of l, ml , corresponds to Dub . For each of these dose levels, a PK profile can be calculated a priori. The set of time steps with nonzero plasma concentration after a dose is given by z = 1, · · · , mz . The algorithm can then be given a set of profiles to choose from at each dosing opportunity, q = 1, · · · , mq , with the concentrations of 9NC

J.M. Harrold, R.S. Parker / Computers and Chemical Engineering 33 (2009) 2042–2054

2051

Table 7 Parameters used in the drug delivery profile optimization problem of case study II. Parameter

Value

Units

B0 Bmin N0 h tf QT

20 0.98 B0 40 5 8 [1, 2, · · · , 8]

g g mm3 min weeks weeks

centration profiles for dose levels between 0.44 and 1.0 mg/kg in increments of 0.01 mg/kg are shown in Fig. 6. This parameterization provides a near-continuous range of concentration profiles. 6.4.3. Effect of transform on body weight calculations The dynamics describing the body weight were also discretized using the same method used to discretize the tumor growth (32a): Bd (k) = gd (k) − ed (k) + B0



gd (k + 1) = gd (k) + h

ed (k + 1) = ed (k) + h

(35a)



kg 1 − g(k) g g

k

e

e

AUClast (k) −

(35b)



1 e(k) e

(35c)

Similar to the PK (34) in the PD effect model (32), the exposure areas were calculated a priori for each dose level in the set l. The total exposure for 9NC and 9AC at each dose level was stored in the parameter Qexp (l), and AUClast (k) was calculated as follows: AUClast (k) =

ml 

bdose,d (q, l), ∀k = Qmap (q), . . . , Qmap (q+1) − 1, ∀q

l=1

(36) The resulting corrected body weight calculation can now be written as: Btox,d (k) = Bd (k) − Fig. 6. Mesh-plot of concentration profiles for 9NC (top pane) and 9AC (bottom pane) for dose levels ranging from 0.44 to 1.0 mg/kg in 0.01 mg/kg increments.

and 9AC predetermined. This is achieved by introducing the binary variable bdose,d (q, l) at each dosing opportunity, q, for each possible dose level, l. The dynamic PK equations can then be replaced with the following constraints on dose days: CNC,d (k) =

ml 

ml 

(34a)

bdose,d (q, l)QAC (q, l),

(34b)

l=1

where k = Qmap (q) − 1 + z, ∀z = 1, . . . , mz 1=

ml 

To avoid the nonlinearity associated with the exponential term, a conservative estimate was made. The transformed tumor volume, Pd (k), was replaced with the value it would take if no drug were administered, QPnom (k): Btox,d (k) = Bd (k) −

eQPnom (k) 1000

(37)

The discretized bound on body weight is given by: bdose,d (q, l)QNC (q, l),

l=1

CAC,d (k) =

ePd (k) 1000

bdose,d (q, l), ∀q = 1, . . . , mq

(34c) (34d)

l=1

Ceff,d (k) = CAC,d (k) + CNC,d (k)

(34e)

Here QNC (z, l) and QAC (z, l) are matrices of precalculated concentrations of 9NC and 9AC, respectively. The constraint (34d) ensures that only one dose level is selected at any given dosing opportunity. Also, the dosing bounds from Eq. (10) were included explicitly in constructing QNC (z, l) and QAC (z, l). The possible nonzero con-

Btox,d (k) ≥ Bmin

(38)

Allowing the tumor size to increase at the nominal rate would yield large values at the end of each treatment cycle. This could then lead to conservative dose schedules being returned by the dosing algorithm. Two options could be used to address this problem: (i) a less conservative value for Bmin could be used; (ii) the duration of treatment cycles could be reduced to allow the clinician more frequent feedback—updating the current body weight and tumor volume at the end of each cycle. Finally, the second case study drug scheduling problem, (30), was reformulated as: min

Dd (q)

mT 

Pd (Qd(T ))

(39a)

T =1

s.t. (32a), (33), (34), (35), (37), (38)

(39b)

2052

J.M. Harrold, R.S. Parker / Computers and Chemical Engineering 33 (2009) 2042–2054

Table 8 Solution properties for (39) subject to: (i) (40a) or (ii) (40b). Dosing profile

Binary variables

Continuous variables

Solution time

(i) (ii)

18,448 16,208

115,217 112,977

11 h 22 min

Fig. 8. Response to dosing 1.0 mg/kg of 9NC QD × 5 × 2 every four weeks. Top pane: dosing profile ( ). Middle pane: normalized corrected body weight (—) and lower bound on corrected body weight (· · · ). Bottom pane: tumor volume (—).

Fig. 7. Optimal dosing profiles, normalized corrected body weights, and tumor volumes solving (39) with a desired tumor volume of zero mm3 at points q and subject to: (i) (40a) or (ii) (40b) Top pane: dosing profiles for (i) and (ii) given by ( ) and (×), respectively. Middle pane: corrected body weight for (i) and (ii) given by (—) and (– –), respectively, with the lower bound, Bmin , given by (· · · ). Bottom pane: tumor volume response to the dosing profiles (i) and (ii) given by (—) and (– –), respectively, with the trajectory of zero tumor volumes indicated by (×). Model predictions of normalized corrected body weight and tumor volumes come from continuous models, not the discrete approximation.

6.5. Results Optimization-specific parameters used in this case study are given in Table 7. Two different ranges of dose levels were considered: Dd (q) = 0, 0.44, 0.45, · · · , 0.99, 1.0 mg/kg

(40a)

Dd (q) = 0, 1.0 mg/kg

(40b)

The first, Eq. (40a), represents the near continuous range of dose possibilities between the upper and lower bounds of administration. Alternatively, a binary dosing option was considered (40b). The problem (39), for both sets of dose levels (40), was formulated in GAMS and solved using CPLEX. Solution characteristics are shown in Table 8. For an eight week treatment window, the problem using either (40a), or (40b) had mk = 16128 time steps and 193,575 equations. By allowing more dosing levels, the solution time was much longer than that found by allowing only the zero/max dosing options. The solution to (39) for dosing opportunities of (40a) and (40b) are shown in Fig. 7. The current standard of practice for this drug is to deliver 0.67 mg/kg of 9NC QD × 5 × 2 – Monday through Table 9 Performance results of different dosing profiles: (i) trajectory tracking using the near-continuous dosing profile, (ii) trajectory tracking allowing for 0 or 1.0 mg/kg of 9NC, and (iii) 1.0 mg/kg QD × 5 × 2 every four weeks.

 Btox (t) 

Dosing profile

max(N(t)) (mm3 )

min(N(t))

min

(i) (ii) (iii)

42.1 44.3 40.0

7.9 6.7 10.5

0.982 0.981 0.972

B0

mq

i=1

25.8 27.0 20.0

D(q) (mg/kg)

Friday for two weeks – repeated every four weeks. For comparison to (40b), the response to 1.0 mg/kg of 9NC QD × 5 ×2, every four weeks, is shown in Fig. 8. Table 9 summarizes the effects of the three different dosing profiles from Figs. 7 and 8. Both of the solutions to (39), (i) the finely discretized (40a) and (ii) the binary dosing options (40b), result in a systematic decrease in tumor burden. This is accomplished by applying more drug while allowing time for animal body weight to recover. As the concept of “more is better” (Thigpen, 2003) is generally employed in clinical dosing design, the ability to administer more drug without violating toxicity constraints is a close approximation to current clinical decision-making. However, toxicity is worse for the solution from Fig. 8 as the normalized corrected body weight is less than Bmin between days 10 and 14. In addition, the QD × 5 × 2 schedule is less effective in treating the cancer. The solutions shown in Fig. 7 are clinically indistinguishable when considered in terms of body weight and tumor volumes. While the solution to (39) with the finely discretized dosing levels is mathematically optimal when compared to the solution to the same optimization using only the binary dosing options, the dosing profile is also more complicated. From a (pre)clinical perspective, the less complicated solution is preferred because it decreases the potential for dosing error and increases the likelihood of patient adherence (Osterberg and Blaschke, 2005). The solution with the 0,1 dosing constraint has higher average tumor volumes but eventually converges (at approximately 28 days) to tumor volumes indistinguishable from the near-continuous dosing range. The finely discretized dosing problem delivered a total of 25.8 mg/kg of 9NC and has minimum and maximum tumor volumes of 7.9 and 42.1 mm3 , respectively. The solution to the more constrained problem delivered slightly more drug, 27 mg/kg of 9NC. The binary dosing option also had a wider range of tumor volumes, 6.7 and 44.2 mm3 for the minimum and maximum values, respectively. Both of the solutions had body weight values above the 98% level specified. The finely graded dosing option had a minimum normalized corrected body weight of 0.982, while the the binary dosing option had a minimum normalized corrected body weight of 0.9801. The standard of practice delivers 20 mg/kg over the eight week window with a minimum normalized corrected body weight of 0.972. The maximum tumor volume is lower than the two solutions previously discussed at 40 mm3 while the minimum tumor volume is higher at 10.5 mm3 .

J.M. Harrold, R.S. Parker / Computers and Chemical Engineering 33 (2009) 2042–2054

7. Summary The objective of this work was to pose and solve model-based cancer treatment problems that were clinically relevant. Dynamic models of drug administration, distribution, and effect, coupled with mixed-integer programming, can be used to recommend candidate treatment profiles for cancer patients. Based on a case study administering 9NC to SCID mice bearing human tumor xenografts, models of drug PK, PD, and toxicity were developed. The computational algorithm was formulated to penalize tumor volume across the horizon (contrary to the endpoint approach of optimal control), and algorithm results included treatment profiles that allow more drug to be dosed than the clinical standard of practice and provide smaller tumor volumes over the horizon studied, all while satisfying the toxicity constraint. The structure outlined herein provides an easily extensible tool for formulating and solving chemotherapy design problems with the potential for intermittent feedback, which is critical to clinical practice. The literature is rife with simulation examples where cancer models are successfully treated by delivery profiles developed using algorithms such as optimal control. So why are none of these treatment algorithms employed clinically? While the answer is anything but simple, a number of key shortcomings can be highlighted: • lack of acceptance (clinician and FDA), with expensive clinical trial needed as proof • questionable tailorability of unvalidated models to individual patients • convenient mathematical formulation vs. clinical practice It is this last item that is the first hurdle to clinical acceptance, because: (i) optimal control for cancer treatment is a poorly posed problem due to the lack of a final time in the clinical setting; (ii) mathematically convenient approaches typically employ a continuously variable dosing amount, while clinical practice employs pills or a limited number of dosing magnitudes (quantized variables); (iii) systemic toxins are not generally dosed continuously (over weeks), but are delivered in cycles with windows for recovery from the common toxicities that manifest as a result of treatment; etc. Of all the issues, toxicity is the key. Published formulations often constrain against toxicity, but that toxicity is generally an inferred or generalized toxicity handled by constraining input magnitude or rate rather than the measurable and specific toxic effect that may be specific to a particular drug-tumor pair. In the clinic, a model of toxicity on par in quality with the PK and efficacy PD models would provide a more significant benefit to treatment as it could provide guidance to the clinician prior to a patient undergoing a toxic event. Finally, mathematical and clinical optimality are not necessarily equivalent, and the eventual use of a treatment profile is in the clinical setting. In this case, a clinically relevant but mathematically suboptimal solution is superior to a mathematically optimal solution that is ignorant of clinical considerations. Dosing decisions are made within an existing framework (cycles with periodic feedback), and any attempt to augment this process must accommodate this framework. This makes close collaboration with clinicians the best, and perhaps only, path for rapid development and deployment of an algorithm such as the one outlined herein. Systems biology is the interface of basic science biology with systems engineering, thereby exploring the systems-type network interactions that play a significant role in fundamental response of biological systems. Problems like the one discussed herein are a result of systems medicine —the translational interface between medical practice (the application of biological knowledge to aid the human condition) and systems engineering. Commonly addressed case studies in cancer, diabetes, anesthesiology, and HIV are only the proverbial “tip of the iceberg”; by learning the biology, systems

2053

engineers can make valuable contributions to the design of novel treatments for disease. As translational scientists, engineers interested in medicine should embrace the challenges posed by clinical problems (e.g., uncertainty, variability, a variety of constraints, novel objectives, etc.), and in doing so construct practical and implementable systems medicine solutions through the development of novel tools in close collaboration with clinical practitioners. Acknowledgments This work would not have been possible without the input of our collaborators from the University of Pittsburgh Cancer Institute: Dr. Merrill J. Egorin, Dr. Julie L. Eiseman, and their laboratories. The work of Jeffry A. Florian, Jr., Ph.D., an alumnus from RSP’s lab, contributed significantly to our understanding of the limits imposed by toxicity constraints. Financial support for this research has been provided by the Whitaker Foundation (RG-01-0311), and the University of Pittsburgh. References Afenya, E. K. (2001). Recovery of normal hemopoiesis in disseminated cancer therapy—a model. Mathematical Biosciences, 172, 15–32. Asachenkov, A., Marchuk, G., Mohler, R., & Zuev, S. (1994). Disease dynamics. Boston: Birkhäuser. Barbolosi, D. (2001). Optimizing drug regimens in cancer chemotherapy: A simulation study using a PK-PD model. Computers in Biology and Medicine, 157–1732. Cella, D., Peterman, A., Hudgens, S., Webster, K., & Socinski, M. A. (2003). Measuring the side effects of taxane therapy in oncology. Cancer, 98(4), 822–831. Chambers, A. F., & Hill, R. P. (1998). Tumor progression and metastasis. In The basic science of oncology (3rd ed.). McGraw-Hill., 219–239, chapter 10. Chen, C. T., & Zuev, S. (1984). Linear system theory and design. Rinehart/Austin, TX: Holt/Winston. Clare, S. E., Nakhlis, F., & Panetta, J. C. (2000). Molecular biology of breast cancer metastasis. The use of mathematical models to determine relapse and to predict response to chemotherapy in breast cancer. Breast Cancer Research, 2(6), 430–435. Cobelli, C., & Carson, E. (2001). An introduction to modelling methodology. In Modelling methodology for physiology and medicine. San Diego, CA: Academic Press., pp. 1–44, chapter 1. Costa, M. I. S., Boldrini, J. L., & Bassanezi, R. C. (1995). Drug kinetics and drug resistance in optimal chemotherapy. Mathematical Biosciences, 125, 191–209. Cui, S. (2002). Analysis of a mathematical model for the growth of tumors under the action of external inhibitors. Journal of Mathematical Biology, 44(5), 395–426. D’Argenio, D. Z., & Schumitzky, A. (1997). ADAPT II users guide: Pharmacokinetic and pharmacodynamic systems analysis software. Biomedical simulations resource. Los Angeles, CA: University of Southern California. Florian, J. A., Eiseman, J. L., & Parker, R. S., Jr. (2004). A nonlinear model predictive control algorithm for breast cancer treatment. In Proc. DYCOPS 7 Boston, MA. Florian, J. A., Egorin, M. J., Zamboni, W. C., Eiseman, J. L., Lagattuta, T. F., Belani, C. P., Chatta, G. S., amd, H. I. S., Solit, D. B., & Parker, R. S., Jr. (2007). A physiologicallybased pharmacokinetic (PBPK) and pharmacodynamic model of docetaxel (Doc) and neutropenia in humans. In American society of clinical oncology annual meeting. Harrold, J. M. (2005). Model-based design of cancer chemotherapy treatment schedules. Ph.D. thesis. University of Pittsburgh, Pittsburgh, PA. Horstmann, E., McCabe, M. S., Grochow, L., Yamamoto, S., Rubinstein, L., Budd, T., Shoemaker, D., Emanuel, E. J., & Grady, C. (2005). Risks and benefits of phase 1 oncology trials, 1991 through 2002. The New England Journal of Medicine, 352(9), 895–932. Johnson, J. J., Monks, A., Hollingshead, M. G., & Sausville, E. A. (2001). Preclinical aspects of cancer drug discovery and development. In Cancer chemotherapy & biology (3rd ed.). Lippincott Williams & Wilkins., pp. 17–36, chapter 2. Kaufman, D. C., & Chabner, B. A. (2001). Clinical strategies for cancer treatment: The role of drugs. In Cancer chemotherapy & biology (3rd ed.). Lippincott, Williams, & Wilkins., pp. 1–16, chapter 1. Kimmel, M., & Gorlova, O. (2002). Stochastic models of progression of cancer and their use in controlling cancer-related mortality. In Proc. American control conf Anchorage, AK, May, (pp. 3443–3448). Kloft, C., Wallin, J., Henningsson, A., Chatelut, E., & Karlsson, M. O. (2006). Population pharmacokinetic and pharmacodynamic model for neutropenia with patient subgroup identification: Comparison across anticancer drugs. Clinical Cancer Research, 12(September 18), 5481–5490. Kreyszig, E. (1999). Advanced engineering mathematics (8th ed.). New York, NY: John Wiley & Sons. Latz, J. E., Karlsson, M. O., Rusthoven, J. J., Ghosh, A., & Johnson, R. D. (2006). A semimechanistic-physiologic population pharmacokinetic/pharmacodynamic

2054

J.M. Harrold, R.S. Parker / Computers and Chemical Engineering 33 (2009) 2042–2054

model for neutropenia following pemetrexed therapy. Cancer Chemotherapy and Pharmacology, 57, 412–426. Ledzewicz, U., & Schättler, H. (2002a). Analysis of a class of optimal control problems arising in cancer chemotherapy. In Proc. American control conf Anchorage, AK, (pp. 3460–3465). Ledzewicz, U., & Schättler, H. (2002b). Optimal bang–bang controls for a twocompartment model in cancer chemotherapy. Journal of Optimization Theory and Application, 114(3), 609–637. Liotta, V., Georgakis, C., & El-Aasser, M. S. (1997). Real-time estimation and control of particle size in semi-batch emulsion polymerization. In Proc. of the American control conf Albuquerque, NM, (pp. 1172–1176). Lowe, V. J., Boyd, J. H., Dunphy, F. R., Dunleavy, H. K., Collins, B. T., Martin, D., Stack, B. C., Hollenbeak, C., & Fletcher, J. W., Jr. (2000). Surveillance for recurrent head and neck cancer using positron emission tomography. Journal of Clinical Oncology, 18(3), 651–658. Martin, R., & Teo, K. L. (1994). Optimal control of drug administration in cancer chemotherapy. River Edge, NJ: World Scientific. Martin, R. B. (1992). Optimal control drug scheduling of cancer chemotherapy. Automatica, 28, 1113–1123. Matveev, A. S., & Savkin, A. V. (2002). Application of optimal control theory to analysis of cancer chemotherapy regimens. Systems & Control Letters, 46, 311–321. Morari, M., & Zafiriou, E. (1989). Robust process control. Englewood Cliffs, NJ: PrenticeHall. Mukherji, S. K., Schmalfuss, I. M., Castelijns, J., & Mancuso, A. A. (2004). Clinical applications of tumor volume measurements for predicting outcome in patients with squamous cell carcinoma of the upper aerodigestive tract. American Journal of Neuroradiology, 25, 1425–1432. Norton, L. (1988). A Gompertzian model of human breast cancer growth. Cancer Research, 48, 7067–7071. Osterberg, L., & Blaschke, T. (2005). Adherence to medication. New England Journal of Medicine, 353(5), 487–497.

Panetta, J. C., & Adam, J. (1995). A mathematical model of cycle-specific chemotherapy. Mathematical and Computer Modelling, 22(2), 67–82. Pereira, F. L., Pedreira, C. E., Pinho, M. R., Fernandes, M. H., & Sousa, J. B. (1990). An optimal control algorithm for multidrug cancer chemotherapy design. In Proc. IEEE EMBS ann. conf., vol. 12 (pp. 1021–1022). Pommier, Y., Pourquier, P., Urasaki, Y., Wu, J., & Laco, G. S. (1999). Topoisomerase I inhibitors: Selectivity and cellular resistance. Drug Resistance Updates, 2, 307–318. Sickles, E. (1989). Breast masses: Mammographic evaluation. Radiology, 173, 297–303. Simeoni, M., Magni, P., Cammia, C., De Nicolao, G., Croci, V., Pesenti, E., Germani, M., Poggesi, I., & Rocchetti, M. (2004). Predictive pharmacokineticpharmacodynamic modeling of tumor growth kinetics in xenograft models after administration of anticancer agents. Cancer Research, 64, 1094–1101. Slingerland, J. M., & Tannock, I. F. (1998). Cell proliferation and cell death. In The basic science of oncology (3rd ed.). McGraw-Hill., pp. 134–165, chapter 7. Swan, G. W. (1986). Cancer chemotherapy: Optimal control using the Verhulst–Pearl equation. Bulletin of Mathematical Biology, 48(3/4), 381–404. Swierniak, A., & Smieja, J. (2001). Cancer chemotherapy optimization under evolving drug resistance. Nonlinear Analysis, 47, 375–386. Thigpen, T. (2003). Maybe more is better. Journal of Clinical Oncology, 21(13), 2454–2456. Thomas, M. M., Joseph, B., & Kardos, J. L. (1997). Batch chemical process quality control applied to curing of composite materials. AIChE Journal, 43(October 10), 2535–2545. Tyler, M. L., & Morari, M. (1999). Propositional logic in control and monitoring problems. Automatica, 35, 565–582. Urquhart, J., & De Klerk, E. (1998). Contending paradigms for the interpretation of data on patient compliance with therapeutic drug regimens. Statistics in Medicine, 17, 251–267.