A measurement-based control design approach for efficient cancer chemotherapy

Information Sciences 333 (2016) 108–125

Contents lists available at ScienceDirect

Information Sciences journal homepage: www.elsevier.com/locate/ins

A measurement-based control design approach for efficient cancer chemotherapy Sofiane Khadraoui a,∗, Fouzi Harrou b, Hazem N. Nounou c, Mohamed N. Nounou b, Aniruddha Datta d, Shankar P. Bhattacharyya d a

Department of Electrical and Computer Engineering, University of Sharjah, United Arab Emirates Department of Chemical Engineering, Texas A&M University at Qatar, Doha, Qatar c Department of Electrical and Computer Engineering, Texas A&M University at Qatar, Doha, Qatar d Department of Electrical and Computer Engineering, Texas A&M University, College Station, TX 77843, USA b

a r t i c l e

i n f o

Article history: Received 13 August 2014 Revised 27 September 2015 Accepted 13 November 2015 Available online 27 November 2015 Keywords: Measurement-based cancer control Cancer chemotherapy Killing cancerous cells PID controller design Frequency response data

∗

a b s t r a c t This paper presents a new measurement-based control design method for cancer chemotherapy. Cancer chemotherapy aims basically at simultaneously eradicating or significantly reducing the number of cancer cells and maintaining tolerable levels of drug concentration and toxicity. To achieve such aim, drugs are often injected into the patient’s body according to a drug schedule specifying the drug dose and delivery time. Several strategies for planning cancer chemotherapy have been developed in the literature, where evolutionary algorithms have been applied to find optimal drug schedules of cancer treatment under constraints on some key treatment parameters such as drug concentration and toxic side effects. In such methods, the amount of drug doses, delivered in the body at each time during the treatment, does not depend on the current drug concentration, toxicity level, and/or number of cancer cells. Successful design of chemotherapy drug scheduling requires the availability of an accurate mathematical model that perfectly predicts the number of cancerous cells and describes effects of treatment. Several models with either complex or simple structures are available in the literature. Complex-structure models are proposed to deeply understand interactions between cancer and normal cells that affect the performance of the cancer chemotherapy. Nevertheless, such complex models are based on a high-order set of differential equations which can be difficult to solve. Simple-structure models, which are often obtained on the basis of some simplifying assumptions, can be viewed only as an approximation of the cancer system. Hence, designing chemotherapy drug schedules on the basis of simplified models may result in unsuccessful cancer treatment. Unlike conventional control strategies for cancer chemotherapy, our attempt in this paper is to address the problem of designing a control system for cancer treatment using a set of frequency-domain data. Hence, a two-degree-offreedom PID (proportional-integral-derivative) control scheme is proposed to control cancer growth. These PID controllers are designed to simultaneously provide the optimal amount of drug doses to be delivered into the patient’s body according to the current drug concentration and toxicity level, and maintain the drug concentration and toxicity levels within their pre-specified ranges. The proposed cancer control technique is validated through a first simulation example. Another example to control biological systems is also presented to show the

Corresponding author. Tel.: +97433147830. E-mail address: [email protected], sofi[email protected] (S. Khadraoui).

http://dx.doi.org/10.1016/j.ins.2015.11.026 0020-0255/© 2015 Elsevier Inc. All rights reserved.

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

109

feasibility of the proposed method. Simulation results obtained have demonstrated the capability of the proposed control scheme to address cancer chemotherapy problems. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Clinical chemotherapy is one of the most common and effective treatment methods for cancer that has acquired a growing interest and importance over the past years. Cancer chemotherapy aims at reducing as much as possible the number of cancerous cells after a fixed period of treatment or even entirely eradicating them. The killing procedure of cancerous cells is based on injecting chemotherapy drug into the patient’s body. However, it is well known that chemotherapy drugs are toxic in nature, which may damage normal cells and therefore affect the life of patients. Hence, drug dosage should be appropriately selected to induce cancerous cell death without affecting normal cells. In clinical applications of cancer chemotherapy, it is often required to find suitable amounts of chemotherapy drug doses so that the maximum drug concentration level in the body and the cumulative toxicity of drug remain below predefined thresholds over the whole period of treatment. Nevertheless, it is worth mentioning that finding optimum drug doses that can simultaneously kill cancerous cells and prevent toxic side effects is not a trivial task. Various approaches to solve the problem of optimizing cancer chemotherapy drug doses have been developed in the literature. The first cancer chemotherapy that takes into account constraints on the drug resistance of the patient and toxicity level was proposed by Martin in [15,16]. The authors in [11,12] used evolutionary computation and computer modeling for optimizing and automating the chemotherapy drug delivery. A direct search optimization procedure to solve chemotherapy drug scheduling problems has been presented in [3]. This solution has been improved by Luus et al. [14], where the drug scheduling problem is solved through direct search optimization using random numbers and search region contraction. Evolutionary algorithms have been widely used to optimize the drug dosage for cancer treatment. The Paladin-distributed evolutionary algorithm has been used to solve the problem of optimal control of drug scheduling in cancer chemotherapy [22]. The authors in [5] proposed a control approach based on an adaptive stochastic algorithm that optimizes the chemotherapy drug scheduling. The authors in [23] solved a multi-drug scheduling optimization problem using an adaptive elitist genetic algorithm combined with a local search algorithm. In [6], an adaptive neural network control methodology is introduced to optimize cancer chemotherapeutic treatment. The above-mentioned control strategies developed for cancer chemotherapy are viewed as open-loop control methods, where the authors attempted to find optimal chemotherapy drug schedules (drug dose and time of delivery) of cancer treatment. In such cancer chemotherapy strategies, no controller is used, meaning that any deviation of the drug concentration and toxicity levels, during the period of treatment, cannot automatically be corrected, since chemotherapy drug schedules are determined off-line. To overcome this problem, the concept of closed-loop control should be used in cancer chemotherapy treatment. Using a closed-loop control scheme, the authors in [1] developed a control technique for designing PID (proportional, integral and derivative) and I-PD (modified PID) controllers for cancer chemotherapy. In this latter work [1], multi-objective optimization techniques are used to tune parameters of the controller that can provide optimal drug dosage and achieve cancer chemotherapy purposes. All the aforementioned control approaches developed for cancer chemotherapy rely on a mathematical description of the tumor growth and effects of treatment. Mathematical models of cancer chemotherapy developed in the literature are generally defined by a set of ordinary differential equations describing interactions between drug dosage and cancerous cell killing, along with the effects of chemotherapy. Nevertheless, such ordinary differential equations can be viewed only as an approximation of cancer chemotherapy systems because of the complex nature of such cancer systems and lack of complete knowledge about their dynamics. Hence, the use of a restricted information on the system for cancer chemotherapy in the absence of a reliable and accurate mathematical description may result in degradation of the performance and outcomes of the cancer treatment. In other words, the success of chemotherapy in cancer treatment is strongly limited by the quality and reliability of the mathematical model in describing the dynamics of the tumor growth and toxic effects during treatment cycle. To overcome this issue, we propose in this paper a control strategy based only on a set of measured data for efficient cancer chemotherapeutic treatment. The main motivation behind introducing the concept of measurement-based control design is that the modeling process is avoided since all relevant information about the system is included in the set of collected data [8,9,20]. Our proposed measurementbased control technique is based on a two degree-of-freedom control scheme, where two PID controllers (see [2]) are tuned to control the drug concentration and toxicity levels and hence enhancing the cancer treatment. The PID controller design problem is formulated as a constrained optimization problem which is solved using Genetic Algorithms (GA). The main advantage of such a proposed measurement-based control technique with respect to the cancer control methods outlined beforehand is that no mathematical model is required here for the design of chemotherapy drug scheduling. Moreover, both PID controllers designed through the proposed technique are able to simultaneously achieve constraints on drug concentration level, toxicity level and drug doses injected into the patient’s body. In addition, the cancer control scheme suggested here does not need specific control policies for the drug dosage as proposed in the aforementioned cancer chemotherapy strategies. For example, the control policy in [1] is based on two patterns of desired drug concentration. In the first pattern, which is called Repeat & Continuous, it is suggested to: (1) set the desired drug concentration to a maximum value of 50 for the first two days, (2) reduce it by 5% for the next two days, (3) set it to zero for the next couple days, and finally maintain it at 40 over the remaining period of treatment. In

110

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

Cancer treatment system

Cell killing

Drug concentration

drug dose

number of cancerous cells drug concentration level

Toxicity

toxicity level

Fig. 1. Cancer chemotherapy system.

Dr (t)

+ -

PID controller-d d

u(t)

Drug concentration

D(t)

Fig. 2. Diagram of the feedback closed-loop control of drug concentration.

the second pattern (called Repeat), the desired drug concentration is switched each day between the maximum allowable value of 50 and zero, over the cycle of treatment. The authors in [1] did not explain how they determined such patterns. The paper is organized as follows. The proposed control method for cancer chemotherapy is presented in Section 2, where the tuning of two measurement-based PID controllers for the drug concentration and toxicity levels is discussed. In Section 3, a simulation example is presented to illustrate and demonstrate the feasibility and efficacy of the proposed cancer control method. Simulation results are discussed and compared with those obtained for other existing cancer chemotherapy techniques. Finally, some concluding remarks are given in Section 4. 2. Proposed control design method This section aims at developing a measurement-based control design approach that can be used to enhance cancer treatment process. Here, we consider the cancer chemotherapy system which involves one input variable (drug dose) and three output variables (number of cancerous cells, drug concentration, and toxicity level) as shown in Fig. 1. 2.1. Feedback control scheme Successful cancer chemotherapy can be achieved by simultaneously maximizing cancerous cell killing and preventing toxic effects on normal cells. To this end, the drug concentration level should be greater than a given threshold β , which is typically 10, in order to eliminate or significantly reduce the number of cancerous cells at the end of treatment. However, such a drug concentration level should not exceed 50 over the whole period of treatment in order to protect normal cells from toxic side effects as discussed in [12]. A drug concentration level between the lower and upper limits specified above can be achieved by manipulating the drug dose injected into the patient’s body. Towards this end, the feedback closed-loop control scheme presented in Fig. 2 is proposed to maintain the drug concentration level, D(t), within its desired range D(t) ∈ [10, 50] at the end of treatment. In such a figure, the proportional-integral-derivative (PID) controller is used to select the drug dose, u(t), to be delivered into the patient’s body. Such a PID controller-d (d is used to indicate the drug concentration) involves three constant gains to be computed: the proportional Kp , integral Ki and derivative Kd gains. At each time of the treatment cycle, the drug concentration level, D(t), is compared with the reference level, Dr (t), to generate an error signal, d (t ) = Dr (t ) − D(t ). The reference level, Dr (t), is the desired set-point value of the drug concentration to be maintained during the treatment process, which must be within the range [10, 50]. The PID controller-d takes the current error signal and provides the suitable drug dose, u(t), to be injected into the patient’s body. The drug dose signal, u(t), can be expressed in terms of three terms as follows:

u(t ) = Kp d (t ) + Ki



0

tf

d (t )dt + Kd

d  (t ), dt d

(1)

where tf is the final time of treatment which is generally about t f = 84 days. The first term of the right-hand-side of (1) provides a drug dose directly proportional to the error, the second term is the accumulation of past errors, and the last term gives a prediction of future errors based on the rate of change of past errors. At the beginning of treatment, the error is positive large, meaning that the drug concentration level, D(t), is lower than the desired set-point drug concentration, Dr (t). In such a case, the cancerous cell killing is minimum and not much efficient. The PID controller-d attempts to reduce the error by increasing the drug dose infused to the body, and therefore improve the cell killing. After a certain time of drug dose injection which depends on the PID controller gain values, the drug concentration level, D(t), reaches its desired set-point value, Dr (t), meaning that the cancerous cell killing is more efficient. At that time where the error tends to zero, the PID controller-d attempts to keep the drug dose constant in order to avoid that the drug concentration level exceeds or falls below its desired value along the

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

outer loop Tr (t)

+ -

t

111

Cell killing

N(t)

Toxicity

T(t)

inner loop PID controller-t

+

Dr(t) -

d

Drug PID controller-d u(t) concentration

D(t)

Fig. 3. Cancer chemotherapy control scheme.

remaining time of treatment. To allow the PID controller to successfully perform duties outlined above, its parameters should be adequately tuned. Some effective tuning methods of PID controllers can be found in the literature [4,19,21,26]. The feedback closed-loop diagram presented in Fig. 2 allows us to maintain the drug concentration level within the pre-defined range [10, 50], and therefore reduce the number of cancerous cells significantly. However, it is very worth noting that a non-proper drug dosage may result in toxic side effects due to the fact that the cumulative drug toxicity is quantified by the integral of drug concentrations over the whole cycle of treatment. To this end, the maximum level of toxicity inside the patient’s body should be constrained in order to avoid any unexpected toxic side effects. According to medical knowledge and experience, the maximum level of toxicity should not exceed 100 during the entire cycle of treatment [13]. Hence, a second PID controller can be used to maintain the level of toxicity below the threshold of 100. To be able to simultaneously maintain the drug concentration and toxicity levels within their corresponding pre-specified ranges, we propose the feedback control diagram shown in Fig. 3. In such a figure, the feedback control scheme of Fig. 2 used to control the drug concentration level D(t) is considered as the inner loop, while the outer loop is used to control the toxicity level, T(t). The system linking the toxicity output, T(t), to the desired set-point of the drug concentration, Dr (t), which consists of the inner loop in cascade with the toxicity part, is considered to be the controlled system for the outer loop. It can be seen from the control diagram shown in Fig. 3 that the desired set-point of the drug concentration, Dr (t), is the system control input (output of the PID controller-t). The inner loop is designed to maintain the drug concentration level, D(t), at the desired set-point value, Dr (t), provided by the PID controller-t, where t stands for toxicity. However, it should be mentioned here that the desired drug concentration level must be below 50 over the entire period of treatment as outlined above. Hence, this point should be considered for the design of the PID controller-t. The controller output (reference level of the drug concentration), Dr (t), can be stated as follows:

Dr (t ) = Kp t (t ) + Ki



0

tf

t (t )dt + Kd

d t (t ), dt

(2)

where K p , Ki and Kd are respectively the proportional, integral, and derivative gains of the PID controller-t to be computed. 2.2. Measurement-based PID controller design We aim here to design the above PID controllers that can be used to control the cancer chemotherapy system. The design of such controllers is based on a finite set of measured data without going through the use of any mathematical model. 2.2.1. Design of PID controller-d Consider the control configuration shown in Fig. 2, where the drug concentration system is controlled by a PID controller-d. The PID controller-d is defined in Eq. (1) with unknown parameters Kp , Ki and Kd . In this contribution, we suppose that an appropriate model describing the drug concentration behavior is unknown, but a set of input/output measurements (drug dose, u(t)/drug concentration, D(t)) is assumed to be available over a range of distinct frequencies {ω1 , ω2 , . . . , ωN }. Based on such input/output measurements, the frequency response of the drug concentration system, denoted by Gd (jω), over the set of frequencies ω ∈ {ω1 , ω2 , . . . , ωN } can easily be obtained using Fourier analysis. The natural question in our control design problem is: based on the frequency response data of the drug concentration system, Gd (jω), how can one derive the PID controller-d that meets some desired closed-loop performance measures? Indeed, such desired closed-loop performance measures are generally given in terms of overshoot, settling time, steady-state error, as well as constrained control input, as follows: • drug concentration response should not exceed the reference level, Dr (t), over the whole period of treatment. This means that closed-loop behavior should be without overshoot. Nevertheless, closed-loop behavior with small overshoot can be allowed in some practical situations; • the time that the drug concentration response takes to reach and stay within 5% of its final value, should be about a given value ts5% . Such a settling time, ts5% , gives us an indication about the cancer treatment speed;

112

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

• the final value of the drug concentration response should be very close to the desired set-point value (reference level), Dr (t), meaning that the steady-state error is null or very close to zero (high precision level); • for a maximum value of the desired drug concentration level Dr (t ) = 50, the drug dose (control input) provided by the PID controller-d should be effective in treatment, but should not exceed a given threshold of drug dosage u, i.e., u(t ) ≤ u. The upper bound u should be selected to be relatively low to ensure effective treatment. The objective here is to find suitable values of the PID controller-d gains for which the above requirements hold. In other words, parameters of such a PID controller-d should be tuned so that the closed-loop system of Fig. 2 best matches a reference system that describes the desired performance measures outlined above. The first three requirements on the drug concentration response can be defined by the following first-order reference model, Hd∗ (s ), s is the Laplace variable:

Hd∗ (s ) =

D (s ) K , = Dr (s ) 1 + τs

(3)

where K and τ are respectively the static gain and time constant. The static gain is generally selected as K = 1 in order to ensure zero steady-state error, meaning that the drug concentration is equal or very close to its desired set-point value at the steadyts state regime. The time constant, which depends on the desired settling time, is defined as τ = 5% . The transfer function, Hd∗ (s ), 3 describes a desired closed-loop response with zero overshoot. It is important to note that the above first-order transfer function, Hd∗ (s ), is selected due to its simple structure, however, a reference model, Hd∗ (s ), with high-order can also be used. The desired frequency response of the drug concentration can easily be derived from the above transfer function, Hd∗ (s ), as follows:

Hd∗ (s = jω ) =

K , 1 + jτ ω

∀ω ∈ {ω1 , ω2 , . . . , ωN }.

(4)

To meet the first three requirements regarding the drug concentration response, the PID controller parameters should be computed such that, for each frequency ω ∈ {ω1 , ω2 , . . . , ωN }, the frequency response of the closed-loop system given in Fig. 2 (denoted by Hdcl ( jω )) is very close to the desired frequency response, Hd∗ ( jω ). The frequency response Hdcl ( jω ) can be obtained as follows: let Cd (s) be the transfer function of the PID controller-d which can easily be derived from (1) as follows:

Cd (s ) = Kp +

Ki + Kd s. s

(5)

Its frequency response, Cd (s = jω ), can be expressed for all frequencies ω ∈ {ω1 , ω2 , . . . , ωN } as follows:

Cd ( jω ) = Kp − j

Ki

ω

+ jKd ω.

(6)

The frequency response, Hdcl ( jω ), of the closed-loop system depicted in Fig. 2 (i.e., the inner loop of Fig. 3), can be derived using the frequency response data of the drug concentration system, Gd (jω), and the controller frequency response, Cd (jω), as follows:

Hdcl ( jω ) =

Cd ( jω )Gd ( jω ) . 1 + Cd ( jω )Gd ( jω )

(7)

Minimizing the following cost function allows us to find parameter values of the PID controller, Cd (s), for which the first three requirements hold:

J = min

Kp ,Ki ,Kd

N    Hd∗ ( jωk ) − Hdcl ( jωk )2 .

(8)

k=1

It should be mentioned that the PID controller, Cd (s), with parameter values resulting from the minimization of the above objective function, J, does not necessarily ensure the constraint on the drug dose delivered to the patient’s body (the fourth requirement outlined above). To this end, the constraint on the drug dose should be considered along with the above minimization problem. From the closed-loop architecture presented in Fig. 2, the drug dose, u, can be expressed as

u = Cd Sd Dr ,

(9)

where Sd is the sensitivity function defined as

Sd =

1 . 1 + Cd Gd

(10)

According to the fourth requirement, to maintain u ≤ u for a maximum desired drug concentration level of Dr = 50, we need, according to Eq. (9), to ensure that:

sup |Cd ( jω )Sd ( jω )| ≤ ω

u . 50

(11)

This inequality constraint related to the amplitude of the control input, u, should also be considered for the controller design. Hence, our control design problem can be stated as finding suitable values of the PID controller-d gains for which the closedloop frequency response, Hdcl ( jω ), fits the desired frequency response, Ht∗ ( jω ), over the range of frequencies {ω1 , ω2 , . . . , ωN },

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

113

and the above inequality constraint is satisfied. This problem can mathematically be formulated as finding gains Kp , Ki and Kd that minimize the following objective function, over the range of frequencies ωk , k = 1, 2, . . . , N, subject to one inequality constraint:

J = min

Kp ,Ki ,Kd

N    Hd∗ ( jωk ) − Hdcl ( jωk )2 , k=1

s.t. sup |Cd ( jωk )Sd ( jωk )| ≤ ωk

u . 50

(12)

The above control design problem is a constrained nonconvex nonlinear optimization problem which can be solved using evolutionary algorithms, such as the well-known Genetic Algorithm described below in Section 2.3. Once suitable values of gains Kp , Ki , and Kd are obtained, it is important to check the stability of the closed-loop system of Fig. 2 for the PID controller with such gains. Since the plant model is assumed here to be unknown, a possible way to check closed-loop stability is via the Nyquist criterion using frequency response data of open-loop system, Cd (jωk )Gd (jωk ), k = 1, 2, . . . , N. According to the Nyquist criterion, the closed-loop system given in Fig. 2 (i.e., the inner loop of Fig. 3) is stable if and only if the Nyquist plot of the open-loop frequency response, Cd (jωk )Gd (jωk ), traversed in the sense of growing frequencies (ω1 < ω2 < . . . < ωN ), leaves the critical point (−1, 0 j ) on the left-hand side [10,18]. It should be mentioned that Nyquist plot can be obtained by plotting the imaginary part of the frequency response against its real part for all frequencies ω1 < ω2 < . . . < ωN . 2.2.2. Design of PID controller-t The inner loop with the PID controller-d designed after solving the above problem is capable of maintaining the drug concentration at its desired level, which must be within the range [10, 50]. Nevertheless, it has been pointed out beforehand that maintaining the drug concentration level within [10, 50] enables us to efficiently kill cancerous cells, but not necessarily to ensure a toxicity level below the pre-defined threshold of 100. Hence, the control of toxicity level is of a great importance to successfully achieve the overall cancer chemotherapy purposes. To this end, we aim here to design the PID controller-t used for the outer loop presented in Fig. 3. It is important to note that the same control design procedure outlined above can be used here to design such a PID controller-t using a finite set of frequency-domain data. It can be seen from Fig. 3 that the overall system to be controlled here consists of the inner loop in cascade with the toxicity part, i.e., the system linking the toxicity output, T(t), to the desired drug concentration level, Dr (t). It should be mentioned that the frequency response data of the inner loop, Hdcl ( jωk ) (for ωk ∈ {ω1 , ω2 , . . . , ωN }), is available since the PID controller-d has been previously computed. The frequency response of the toxicity part, denoted Pt (jωk ), can be obtained from a set of drug concentration level/toxicity level measurements. Hence, the frequency response data, Gt (jωk ), of the system that relates the toxicity level T(t) to the desired set-point of the drug concentration, Dr (t), can be obtained as follows:

Gt ( jωk ) = Hdcl ( jωk )Pt ( jωk ).

(13)

It should be mentioned here that the frequency response, Gt (jωk ), can also be obtained directly from a set of desired concentration Dr (t)/toxicity level T(t) measurements carried out on the controlled system (inner loop system+toxicity part). Once the frequency response data, Gt (jωk ), is available, the PID controller-t can be designed to meet the following requirements for any desired set-point value of the toxicity, Tr (t) ≤ 100: • no overshoot is allowed for the toxicity level response of the closed-loop system linking the toxicity level, T(t), to the desired toxicity Tr (t); • settling time of the outer loop should be around a given value ts5% . Such a response time, ts5% , should be selected such that the response of the inner loop is faster than that of the outer loop. A suitable choice can be ts5% = 2ts5% ; • the final value of the toxicity response should be very close to the desired level of toxicity, Tr (t) (zero steady-state error); • for a maximum value of the desired level of toxicity Tr (t ) = 100, the reference level of the drug concentration Dr (t) (control input) provided by the PID controller-t should not exceed the upper bound Dr (t) ≤ 50. Similar to the case of the drug concentration control, a first-order reference model that describes the above first three requirements can be defined as:

Ht∗ (s ) =

K , 1 + τs

where K = 1 and τ =

Ht∗ ( jω ) =

(14)

ts5% . The frequency response of the above reference model can be expressed as follows: 3

K . 1 + jτ ω

(15)

The frequency response of the outer loop shown in Fig. 3, can be expressed in terms of the frequency response data, Gt (jω), K

and the controller frequency response, Ct ( jω ) = K p − ωi + jKd ω, as follows:

Htcl ( jω ) =

Ct ( jω )Gt ( jω ) . 1 + Ct ( jω )Gt ( jω )

(16)

114

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

Gains K p , Ki and Kd of the PID controller-t, Ct (s), for which the first three requirements outlined above are met, can be obtained by minimizing the following objective function J:

J = min

Kp ,Ki ,Kd

N    Ht∗ ( jωk ) − Htcl ( jωk )2 .

(17)

k=1

It is very worth mentioning again that the minimization of this objective function, J, may result in gain values, K p , Ki and Kd , for which the fourth requirement on the drug concentration set-point (which should not exceed 50 for a maximum desired level of toxicity of 100) is not satisfied. Hence, it is required to account, in addition to the first three requirements, for this condition in the design process. From Fig. 3, the relationship linking the toxicity set-point, Tr , to the drug concentration set-point, Dr , can be expressed as follows:

Dr = Ct St Tr ,

(18)

1 1+Ct Gt

where St = . As outlined above, it is required that Dr ≤ 50 for a maximum desired toxicity level Tr = 100. To this end, the PID controllert should be designed to achieve, in addition to the first three requirements, the following inequality constraint for ω ∈ {ω1 , ω2 , . . . , ωN }:

sup |Ct ( jω )St ( jω )| ≤ ω

50 . 100

(19)

Similar to the first case dealing with the control of the drug concentration, the design problem here consists in finding suitable values of the PID controller parameters K p , Ki , and Kd that minimize the following constrained nonconvex nonlinear optimization problem, over the range of frequencies ωk , k = 1, 2, . . . , N:

J = min

Kp ,Ki ,Kd

N    Ht∗ ( jωk ) − Htcl ( jωk )2 , k=1

s.t. sup |Ct ( jωk )St ( jωk )| ≤ 0.5. ωk

(20)

This constrained nonconvex nonlinear optimization problem can be solved in a similar way using Genetic Algorithm. After obtaining parameters of the controller, Ct (s), the stability of the outer loop of Fig. 3 should be analyzed through the Nyquist criterion using the open-loop frequency response data Ct (jωk )Gt (jωk ), k = 1, 2, . . . , N. 2.3. Genetic Algorithm Solving both design problems (12) and (20) outlined above, which are constrained nonlinear programming problems, requires the use of global optimization methods. Genetic Algorithm (GA), which is a numerical global optimization technique, is used here to solve our control design problems (12) and (20). GA is a stochastic optimization algorithm originally inspired by the mechanism of natural selection and evolutionary genetics [7,24]. The basic idea of this evolutionary optimization algorithm has been developed by John Holland in the late 1960s and early 1970s [7]. GA is based on an initial population of chromosomes, where each chromosome is viewed as a candidate solution (here values of the PID controller parameters [Kp , Ki , Kd ] or [K p , Ki , Kd ]) to the problem. The size of this initial population (called first generation) is typically selected between 20 and 100. In GA, an encoding of variables of the optimization problem is used rather than using variables themselves. In other words, each candidate solution (chromosome) is encoded as a string of binary digits (sequence of 0s and 1s). All chromosomes in the first generation are evaluated using the objective function and constraints, to ultimately assign a fitness score to each chromosome. The fitness score determines the solution’s quality (level of optimality of each chromosome as a solution to the problem at hand). In the next step, GA attempts to create a new population from the current population using three main operators: selection, crossover, and mutation. The main objective of the selection operator (reproduction operator) is to identify, from the current population, chromosomes with best fitness function (good solutions), eliminate worst solutions, make multiple copies of better solutions and place them in a mating pool (having the same size as the original population). It should be mentioned that the selection operator does not attempt to generate new solutions from the current solutions, but just to make multiple copies of good solutions. New solutions can be created by crossover and mutation operators. The crossover operator aims at taking some pairs of chromosomes (called parent’s chromosomes) at random from the mating pool, and exchanging some portions of strings (substrings) of each pair to produce a new pair of offspring chromosomes. After the crossover operation, parent chromosomes in the mating pool are replaced by their offspring chromosomes, while keeping the same number of chromosomes in the mating pool. Next, the mutation operator is applied to the resulting offspring chromosomes. The mutation operation takes each chromosome from the new mating pool and attempts to flip a bit of the string with a low probability (such as 0.001) from 0 to 1, or vice versa. The main motivation behind this mutation phase is to avoid premature convergence of the algorithm to a local optimum. After the application of these three operators, a new generation of offspring chromosomes with the same size as the original population is obtained. The same procedure is repeated using the new generation of chromosomes, and after a user-specified number of iterations, the algorithm stops and provides the solution. The accuracy of the obtained solution depends greatly on the maximum number of iterations as well as on the size of population selected.

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

115

Table 1 Parameter values of the cancer chemotherapy system. Parameter

Value

λ κ β γ η

1.5 × 10−4 per day 9.9 × 10−3 per day per drug unit 10 drug unit 0.27 per day 0.4 per day

Remark 1. It is important to note that GA, which is used as an effective numerical global optimization algorithm, is capable of providing, after a finite number of iterations, only an approximate numerical solution, rather than an accurate optimal solution. Remark 2. It is worth noting that the study given in this paper relies on PID controller structures due to their simplicity (low dimensionality) and their wide use in practical applications. However, other different controller structures can also be used to control cancer growth. Nevertheless, it should be mentioned that the efficacy of the cancerous cell killing depends on the controller structure selected. 3. Application of the proposed control method 3.1. Case study 1 This section is mainly devoted to a simulation study in which we show the efficacy of the proposed measurement-based cancer chemotherapy control approach. Over the few past years, more attention has been given to the study of tumor growth behavior. One important relevant work is the one by Martin [17], where he developed drug scheduling models for cancer chemotherapy purposes. According to Martin’s work, the effects of cancer chemotherapy are described by the following differential equations:

d X (t ) = −λX (t ) + κ (D(t ) − β )H (D(t ) − β ) dt

(21a)

d D(t ) = u(t ) − γ D(t ) dt

(21b)

d T (t ) = D(t ) − ηT (t ) dt

(21c)

where X(t), D(t), T(t), and u(t) (t is the time variable) are the transformed variable, drug concentration level, toxicity level, and drug dose, respectively. λ, κ , β , and η are the growth speed of cancerous cells, cells killed per unit time per unit drug concentration, threshold level of drug concentration, and the elimination rate constant of the drug toxicity. H is the Heaviside unit function defined as follows:



H (D(t ) − β ) =

1 0

if D(t ) ≥ β if D(t ) ≤ β

(22)

Over the whole period of the cancer chemotherapy, the transformed variable X(t) is inversely related to the number of cancerous cells per unit time, N(t), as follows:

N (t ) = 1012 × exp(−X (t )).

(23)

Values of the parameters in the above equations are summarized in Table 1. Initial values of the number of cancerous cells, drug concentration level, toxicity level, and drug dose are respectively N (0 ) = 1010 , D(0 ) = 0, T (0 ) = 0, and u(0 ) = 0. Remark 3. It is very worth noting that the differential equations outlined above are assumed to be unknown for our control design purposes. They are used here only to generate frequency-domain data over a finite set of distinct frequencies. Indeed, the design process of both PID controllers does not require the availability of the plant model (differential equations). In this simulation example, we design both PID controllers Cd (s) and Ct (s) to achieve the desired closed-loop performance specifications outlined beforehand in Sections 2.2.1 and 2.2.2, where ts5% = 10 days and ts5% = 20 days and the upper bound of drug dosage considered here is u = 20 for the maximum desired drug concentration level Dr (t ) = 50. The final time of treatment is selected to be t f = 84 days. The frequency range of interest considered in this simulation study ranges from 0 to 0.5Hz, in which 120 distinct frequencies {ω1 , ω2 , . . . , ωN=120 } have been selected. Let us first design the PID controller, Cd (s), for the drug concentration. In our case, equations that perfectly describe the cancer behavior are assumed to be unavailable. Hence, the frequency response data, Gd (jωk ), should be computed experimentally from a set of drug dose/drug concentration measurements as outlined in the previous section. However, since the objective in this numerical example is only to validate the proposed cancer control method, equations given in (21) are used here to generate the frequency response data, but they are not used for

116

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125 Table 2 Results obtained for different couples of population sizes and number of iterations. Couple

Cd (s)

J

Kp (5, 1) (5, 5) (5, 10) (5, 20) (5, 50) (20, 1) (20, 5) (20, 10) (20, 20) (20, 50) (50, 1) (50, 5) (50, 10) (50, 20) (50, 50)

Ki

0.314 0.315 0.31 0.3 0.3 0.376 0.271 0.301 0.299 0.3 0.302 0.301 0.3 0.3 0.3

0 0.098 0.08 0.0805 0.081 0.0095 0.0818 0.0816 0.0809 0.081 0.074 0.0814 0.0812 0.081 0.081

Ct (s) Kp

5.587 0.08 0.004 1.8 × 10−7 4.55 × 10−8 2.319 0.037 1.64 × 10−4 2.5 × 10−8 1.19 × 10−10 0.0144 10−4 8.22 × 10−7 1.64 × 10−8 6.6 × 10−11

0.431 0.219 0.27 0.3 0.315 0.407 0.336 0.332 0.328 0.327 0.357 0.338 0.329 0.331 0.327

J Ki 1.88 × 10 0.101 0.082 0.054 0.0547 0.002 0.049 0.0538 0.055 0.054 0.028 0.061 0.0534 0.054 0.054

Kd −7

0.105 0.416 0.388 0.345 0.339 0.22 0.258 0.323 0.328 0.328 0.23 0.31 0.326 0.325 0.328

8.382 1.328 0.5 0.107 0.098 4.96 0.251 0.098 0.0958 0.094 0.592 0.14 0.0948 0.0944 0.094

the control design process as mentioned above in Remark 3. Using Eq. (21b), frequency response data of the drug concentration system, Gd (jωk ), can easily be derived over the considered range of frequencies {ω1 , ω2 , . . . , ω120 } as follows:

G d ( j ωk ) =

1

γ + j ωk

,

k = 1, 2, . . . , 120.

(24)

The above frequency response data, Gd (jωk ), is substituted in the closed-loop and sensitivity frequency responses, Hdcl ( jωk ) and Sd (jωk ), used in the constrained optimization problem (12). The resulting constrained optimization problem is then solved using Genetic Algorithm to find the optimal values of the PID controller-d parameters Kp , Ki , and Kd . In this example, an initial population with size 20 and maximum number of iterations of about 20 are selected. The selection of such values will be justified later. After applying Genetic Algorithm, the following PID controller gains are obtained: Kp ≈ 0.3, Ki ≈ 0.081, and Kd = 0. Now, we attempt to tune parameters of the PID controller-t, Ct (s), used for the outer loop. As discussed in the previous section, the plant frequency response data of the toxicity system, Gt (jωk ), can be obtained from a set of input/output measurements. Here, it is defined over the considered range of frequencies {ω1 , ω2 , . . . , ω120 } as follows:





Kp − j ωKi 1  k , Gt ( jωk ) = η + jωk Kp + γ + j ωk − ωKi k

(25)

where Kp ≈ 0.3 and Ki ≈ 0.081 are gains of the PID controller-d obtained above. The frequency response data, Gt (jωk ), is substituted in the constrained optimization problem (20) to find gains of the PID controller-t, Ct (s). The same size of population and number of iterations are used again here. Applying Genetic Algorithm to the optimization problem at hand results in suitable values of the PID controller-t parameters K p = 0.328, Ki = 0.055, and Kd = 0.328. According to Remark 1, gain values of both PID controllers (Cd (s) and Ct (s)) obtained above are only approximate solutions for the design problems. Since the quality of the solution depends greatly on the size of population selected, as well as on the maximum number of iterations, we attempt in the sequel to observe the pattern of both PID controllers’ gains and analyze their impact on values of the objective functions (performance of cancer treatment). To accomplish this, both PID controllers are re-designed for different couples of population size and number of iterations. Here, three population sizes, pop_size = 5, pop_size = 20 and pop_size = 50, are selected, whereas the number of iterations, iter, is selected to be {1, 5, 10, 20, 50}, where their combinations result in 15 couples (pop_size, iter). Gains of PID controllers (Cd (s) and Ct (s)) and minimum values of the cost functions, J and J, obtained for each couple selected (pop_size, iter) are given in Table 2. It should be mentioned that the values of the derivative gain of the drug concentration PID controller, Cd (s), is neglected for all cases, i.e., Kd ≈ 0. It can be observed from such results that increasing the population size and number of iterations provides minimum objective functions, J and J. Controller parameter values that provide the minimum objective functions (good fitting between the computed and desired frequency responses) allow us to meet the desired performance measures. Indeed, a high quality fit between closedloop frequency responses (Hdcl ( jωk ) and Hd∗ ( jωk ), as well as Htcl ( jωk ) and Ht∗ ( jωk )) means a suitable accuracy of the control procedure. It can be verified from results given in Table 2 that the gains, Kp and Ki , of the PID controller, Cd (s), converge to the approximate optimal solution for pop_size = 5 and iter = 20. It can also be verified that the convergence of gains K p , Ki , Kd to their approximate optimal values requires a population size of about 20, and more iterations ite ≥ 20. Hence, in this particular example, a population size of about 20 and number of iterations of about 20 are quite enough to find suitable parameter values of both PID controllers for an effective cancer chemotherapy treatment.

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

117

500 400 300 Imaginary axis

200 100 0

+

-100

critical point (-1,0j)

-200 -300 -400 -500 -1.2

-1

-0.8

-0.6 -0.4 Real axis

-0.2

0

0.2

Fig. 4. Nyquist plot of the open-loop system (forward path system of the inner loop), Cd (jωk )Gd (jωk ).

3.2. Stability and performance analysis Now we attempt to evaluate the effectiveness of the PID controllers, obtained for pop_size = 20 and iter = 20, in achieving the closed-loop requirements in terms of overshoot, settling time, steady-state error, and limits on the drugs and toxicity level. As mentioned in the previous section, before the performance evaluation phase, we need to check whether or not the obtained controllers stabilize the entire closed-loop system in Fig. 3. This stability check can be performed using the Nyquist criterion since the plant model is assumed here to be unknown. Let us first start by analyzing the stability of the inner loop of Fig. 3 using the open-loop frequency response, Cd (jωk )Gd (jωk ) (k = 1, 2, . . . , 120). The Nyquist plot of the forward path system of the inner loop, over the set of frequencies, {ω1 , ω2 , . . . , ω120 }, is shown in Fig. 4. It can be seen from Fig. 4 that the open-loop Nyquist contour, traversed in the sense of growing frequencies, {ω1 < ω2 < . . . < ω120 }, does not encircle the critical point (−1, j0 ). Hence, according to results in [10,18], the designed controller, Cd (s), stabilizes the inner loop of Fig. 3. The same procedure is utilized here to check whether or not the PID controller-t, Ct (s), ensures the stability of the outer loop of Fig. 3. In such a case, the open-loop frequency response, Ct (jωk )Gt (jωk ) (k = 1, 2, . . . , 120), shown in Fig. 5, will be used for stability analysis. According to such a figure, it is clear that the Nyquist plot, traversed in the sense of growing frequencies, {ω1 < ω2 < . . . < ω120 }, leaves the critical point (−1, j0 ) on the left-hand side, meaning that the entire closed-loop system of Fig. 3 is stable. At this stage, it has been verified that the inner and outer loops of Fig. 3 are stable for the PID controllers designed, Cd (s) and Ct (s). Thus, these obtained controllers can be implemented to evaluate the resulting responses. To evaluate the performance of the closed-loop system shown in Fig. 2 (or the inner loop in Fig. 3), the PID controller-d, Cd (s), with gain values, K p = 0.3, Ki = 0.081, and Kd = 0, is first placed in the closed-loop scheme, then a step reference input of a maximum amplitude Dr (t ) = 50 is applied. The simulated response of the drug concentration level obtained is plotted in Fig. 6(a). It can be seen from Fig. 6(a) that the drug concentration response obtained is without any overshoot and with very small steady state error. It has been found that the settling time is about ts5% ≈ 9.9 days. Fig. 6(b) presents the drug dose to be injected into the patient’s body. It is clear from Fig. 6(b) that the drug dose resulted from the application of the maximum reference input, Dr (t ) = 50, ranges from 13.5 and 15 (below the maximum allowable drug dose u = 20) over the whole period of treatment 84 days. Such simulation results have demonstrated the efficacy of the computed PID controller-d to achieve the desired performance for the drug concentration system. In order to demonstrate the efficacy of the proposed method, we attempt next to evaluate the performance of the cancer chemotherapy control scheme shown in Fig. 3. To this end, both PID controllers, Cd (s) and Ct (s), obtained for pop_size = 20 and iter = 20, are first placed in the closed-loop scheme shown in Fig. 3, then a simulation test is performed for 84 days by applying a step reference input of maximum amplitude Tr (t ) = 100. Simulation results obtained, such as toxicity level T(t), drug concentration level D(t), desired drug concentration Dr (t), and chemotherapy drug dose u(t), are presented in Fig 7. It can be seen from Fig. 7(a) that the toxicity response rapidly rises in the first two weeks of cancer treatment, then gradually increases to its steady-state level over a period of three weeks. Over the remaining period of treatment, the toxicity finally approaches its desired value Tr (t ) = 100. Such a toxicity response does not exhibit any overshoot or steady-state error. The settling time obtained here is about ts5% ≈ 20.5 days. Fig. 7(b) shows the chemotherapy drug concentration which increases, in a similar manner as observed in toxicity response, to reach a level of 40 in approximately six to seven weeks. It is clear from Fig. 7(b) that the drug concentration level, over the whole period of treatment, does not exceed its maximum allowed value D(t) ≤ 50. The drug concentration reference signal, Dr (t), provided by the computed controller-t, Ct (s), is presented in Fig. 7(c). It can be viewed from Fig. 7(c) that the drug concentration reference, over the entire period of treatment, varies in the range

118

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

300 Imaginary axis

200 100 0

−100 −200 −300 -1.2

-1

-0.8 -0.6 -0.4 -0.2 Real axis (a)

0

0.2

-0.2

0

15

Imaginary axis

10 5 0

+ critical point (-1,0j)

-5 -10 -15 -1.2

-1

-0.8

-0.6 -0.4 Real axis (b)

Fig. 5. (a) Nyquist plot of the open-loop system (forward path system of the outer loop), Ct (jωk )Gt (jωk ). (b) Zoom on the part of the Nyquist plot around the critical point.

Drug concentration [drug unit]

60 50 40

D r (t)

30 20

D(t)

10 0

0

10

20

30

40

50

60

70

80 90 Time [days]

50

60

70

80 90 Time [days]

(a) 15.5

Drug dose [drug unit/day]

15 14.5

u(t)

14 13.5 13

0

10

20

30

40 (b)

Fig. 6. Simulation results. (a) Drug concentration step response, (b) response of the drug dose injected into the patient’s body.

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

Drug concentration [drug unit]

Toxicity [drug unit . day]

50

100

40

80

Tr (t)

60

30 D(t)

20

40 20 0

119

T(t) 0

20

10 40

0

60 80 Time [days]

(a)

0

20

40 (b)

Drug concentration ref [drug unit]

11

42

60 80 Time [days]

Drug dose [drug unit/day]

40 10.5

38

u(t)

D r (t)

36

10

34 32

9.5 0

20

40 (c)

60 80 Time [days]

0

20

40 (d)

60 80 Time [days]

Fig. 7. Simulation results. (a) Toxicity response T(t), (b) drug concentration D(t), (c) drug concentration reference Dr (t), (d) drug dose u(t).

Table 3 Performance measures obtained. Drug dose

Drug concentration

Toxicity

Average

Max

Average

Max

Average

Max

10.668

10.8

37.55

40

90.96

100

9

x 10 number of cancerous cells 12 10 8 6 4

N(t) 2 0 0

10

20

30

40

50

60

70 80 90 Time [day]

Fig. 8. Number of the cancerous cells, N(t), remaining at the tumor side.

of 32–40 allowing to achieve a drug concentration level lower than 50. The chemotherapy drug doses are shown in Fig. 7(d) over the 84 days of treatment, where it is shown that the drug doses range from 9.928 to 10.8. It is important to note that the maximum drug dose, u(t ) = 10.8, obtained for a maximum desired drug concentration Dr (t ) = 40, remains much lower than 20 the allowable upper bound of drug dose u = 40. 50 = 16. In this simulation test, the average and maximum values of drug doses, drug concentration level and toxicity level during the whole period of treatment are summarized in Table 3. Based of the above simulation results and Eqs. (21a) and (23), the number of the remaining cancerous cells over the period of treatment is obtained and drawn in Fig. 8. According to Fig. 8, it is clear that drug doses injected into the patient’s body significantly reduce the initial number of the cancerous cells (N0 = 1010 ) over the first three weeks of treatment, where the

120

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125 Table 4 Results obtained for η = 0.36 and η = 0.44.

η

0.36 0.44

Drug concentration

Toxicity

Number of

Average

Max

Average

Max

Remaining cells

34.137 40.7

36 3.997

91.584 89.846

99.999 99.992

18.336 0.0785

number of cells remaining is nearly 1.55 × 108 . The drug doses continue killing cancerous cells till the end of treatment, where one single cancerous cell remains at the last day of treatment, N (84 ) = 1.078. According to such results, it is clear that the cancer chemotherapy control scheme proposed in this paper has the ability to successfully kill cancerous cells and avoid any unexpected issues of toxicity. It has been shown that values of the drug concentration, drug dose as well as toxicity level are below their maximum limiting values predefined in inequality constraints of the optimization problems. 3.3. Sensitivity analysis of the cancer chemotherapy control scheme At this stage, we have shown that the implemented PID controllers, Cd (s) and Ct (s), have successfully achieved the desired closed-loop performance measures. Nevertheless, due to the fact that GA algorithm can provide only approximate numerical values of both PID controller parameters, it is important to analyze the robustness of the proposed cancer chemotherapy control scheme with respect to the plant parameter variations. We aim here to perform this analysis using different values of the following four parameters, λ, κ , η, and γ , given in Table 1. To this end, each parameter will be perturbed by 10% of its nominal value to evaluate the effectiveness of PID controllers, Cd (s) and Ct (s), obtained above, for the cancer treatment. For a variation of the parameter λ within the range [1.35 × 10−4 , 1.65 × 10−4 ], it has been verified that the drug concentration, toxicity level, as well as the number of cancerous cells at the end of the treatment remain approximately similar to those obtained for the nominal value of λ. This variation in the parameter λ does not have a significant effect on the treatment. Now, we attempt to analyze, in a similar way, the effect of the variation of the parameter κ (cells killed per unit time per unit drug concentration) on the performance of the proposed cancer treatment scheme. The parameter κ is varied here within the range [8.9 × 10−3 , 10.9 × 10−3 ], while the remaining parameters are set at their nominal values. The variation of κ over this range results in a small number of cancerous cells at the end of treatment, which is between 0.1177 and 11.546. The drug concentration and toxicity level are not affected by the parameter κ . The variation of the parameter κ has an effect on the speed of cancerous cell killing. Indeed, the number of cancerous cells remaining after the first three weeks of treatment are 1.832 × 108 , 1.55 × 108 , and 0.75 × 108 for the minimum, nominal, and maximum values of κ , respectively. The effect of the rate of toxicity elimination, η, on the cancer treatment is also analyzed. The variation of this parameter, η, within the range [0.36, 0.44] affects the performance of the cancer treatment. Results obtained for the upper and lower bounds of η are summarized in Table 4. For both values of η, the maximum drug concentration and toxicity level obtained do not exceed their limit values. The parameter η has also an effect on the speed of cancerous cell killing, where large values of η help us to quickly kill cancer cells. When η is decreased from 0.4 to 0.36, it has been observed that the number of cancerous cells remaining at the end of treatment has been increased from 1.078 to 18.336, which is acceptable from a practical point of view. For η = 0.44, PID controllers have almost killed cancerous cells. We have also analyzed the effect of varying the drug decay parameter, γ , within the range [0.243, 0.297] on the proposed cancer chemotherapy treatment. Over the range of variation of γ ∈ [0.243, 0.297], it has been verified that the drug concentration, toxicity level, as well as number of cancerous cells are approximately similar to those obtained for nominal values. Such results obtained above demonstrate the robustness of the proposed control scheme for cancer chemotherapy treatment when the plant parameters are perturbed by 10% of their nominal values. Nevertheless, it very worth mentioning that a very large variation of the plant parameters from their nominal values may lead to unsatisfactory performance results. 3.4. Comparative study 3.4.1. Comparison with some existing methods Now, we attempt to conduct a comparative performance analysis of the proposed cancer chemotherapy control method with different available works. Such a comparative performance analysis considered here is given in terms of the number of the cancerous cells, N(tf ), remaining at the end of treatment, i.e., at t = t f = 84 day. Table 5 presents the comparative performance analysis of the proposed approach with some relevant works in the field [1,3,5,12,14,17]. Cancer chemotherapy control methodologies presented in [3,5,12,14,17] are viewed as open-loop control strategies, where an optimal amount of drug doses to be injected into the patient’s body is computed so that a minimum number of cancerous cells is achieved under constraints on the drug resistance of the patient and toxicity level of the drug. In these strategies based on the concept of open-loop control (called also non-feedback control scheme), no controller is used. Such approaches are limited by the fact that, at each time during the treatment, the amount of drug doses delivered in the patient’s body does not depend on the current drug concentration, toxicity level, and number of cancer cells. When disturbances act on the cancer system, a large drift in the output responses (such as drug concentration and toxicity responses) may occur. In such a case, the open-loop control system cannot correct such deviations to

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

121

Table 5 Number of cells remaining for different methods. Method

Numbers of cells, N(tf )

Martin and Teo [17] Bojkov et al. [3] Luus et al. [14] Carrasco and Banga [5] Liang et al. [12] Algoul et al. [1] Proposed method

4.878 × 104 3.31 × 104 2.57 × 104 1.534 × 104 1.76 × 103 15 1.078

maintain the drug concentration and toxicity levels within their ranges. In cancer chemotherapy control procedures proposed in [1], a closed-loop control scheme (feedback control system) is utilized, where the drug concentration is fed back for comparison with its set-point value. This closed-loop cancer chemotherapy control system involves one controller, where its parameters are adjusted, at each time during the treatment, using a multi-objective Genetic Algorithm so that the number of remaining cancerous cells is continuously reduced and the drug concentration and toxicity level remain within their specified ranges. In our work, two control loops (inner and outer loops) are used for cancer chemotherapy control, where each loop involves one controller. The controller in the inner loop is designed to automatically provide, at each time during the period of treatment, suitable drug doses based on the difference between the actual drug concentration and its set-point. The controller in the outer loop is designed to provide, at each time during the treatment, suitable drug concentration set-point values according to the difference between the actual toxicity level and its desired set-point value. Closed-loop cancer chemotherapy control is advantageous in the sense that any deviation, during the cancer treatment, of the output responses from their set-points can automatically be corrected by the controller. The above-mentioned cancer chemotherapy control strategies are different in nature, but their one common aspect is the use of Martin’s model (21). According to results presented in Table 5, the first five drug scheduling techniques have reduced the number of cancerous cells to a value in the range of 1.76 × 103 –4.878 × 104 . The drug scheduling technique developed in [1] has significantly reduced the number of cancerous cells at the end of treatment to 15. From a practical point of view, totally killing cancerous cells or reducing as much as possible their size is of a great importance. It has been shown in this paper that the proposed technique has almost eliminated cancerous cells, where only one single cell remains at the end of treatment. Hence, the proposed cancer chemotherapy control strategy has outperformed all results reported in [1,3,5,12,14,17]. 3.4.2. Comparison of results obtained for other controller structures In the previous section, a PID structure has been used for both controllers, Cd (s) and Ct (s). In order to show that the performance of the cancer chemotherapy depends on the controller structure selected (as outlined in Remark 2), we attempt here to re-apply the proposed methodology by assuming new different structures for both controllers, Cd (s) and Ct (s) as follows:

a1 s + a0 a2 s + 1 b1 s + b0 . Ct (s ) = b2 s + 1

Cd (s ) =

(26)

The tuning method presented above results in the following controller parameter values: a0 = 6.447, a1 = 21.606, a2 = 73.80, b0 = 4.422, b1 = 32.354 and b2 = 96.494. Both controllers above, with such parameter values, have been placed in the feedback control system of Fig. 3 to evaluate their ability in cancer chemotherapy. Simulation results obtained are shown in Fig. 9 with solid red line. For comparison purposes, we have plotted in the same figure results (responses with solid blue line) obtained above for PID controllers. It is clear from Fig. 9 that quite different responses are obtained when other controller structures are selected. According to Fig. 9(a), the toxicity response obtained for controllers in (26) does not reach the desired level of toxicity reference, Tr (t ) = 100. This steady-state error is due to the fact that the controller structure, Ct (s), in (26) does not involve an integral action. Moreover, it can be verified that the drug dose and drug concentration level (obtained for controllers given in (26)) remain below their allowable maximum values, over the whole period of treatment. Table 6 summarizes the average and maximum values of drug doses, drug concentration level and toxicity level during the whole period of treatment. Simulation results obtained above are used to compute the number of remaining cancerous cells over the period of treatment, which is drawn in Fig. 10 with solid red line. The number of cancerous cells obtained for PID controllers and presented in Fig. 8

Table 6 Performance measures obtained for controllers given in (26). Drug dose

Drug concentration

Toxicity

Average

Max

Average

Max

Average

Max

9.642

10.588

33.9

36.54

82.09

91.34

122

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

Drug concentration [drug unit]

Toxicity [drug unit . day]

50

100

40

80 Tr (t)

30

60 40

20

20 0

10 0

20

40 (a)

0

60 80 Time [days]

0

11

40

10.5

38

10

36

9.5

34

9 20

40 (c)

60 80 Time [days]

Drug dose [drug unit/day]

42

0

40 (b)

Drug concentration ref [drug unit]

32

20

8.5 0

60 80 Time [days]

20

40 (d)

60 80 Time [days]

Fig. 9. Comparison of simulation results. Responses obtained for PID controllers plotted with solid blue line. Responses obtained for first-order controllers plotted with dashed black line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

9

x 10 number of cancerous cells 12 10 8 6

with first-order controller (26)

4 2

with PID

0 0

10

20

30

40

50

60

70 80 90 Time [day]

Fig. 10. Number of the remaining cancerous cells, N(t).

are plotted here again in Fig. 10 for comparison purposes. According to Fig. 10, it can be seen that both responses are very close to each other. Nevertheless, although controllers designed above (26) significantly reduce the initial number, N0 = 1010 , of the cancerous cells, 22.36 cancerous cells remain at the end of treatment. Hence, it is clear, in this particular example, that the selection of PID controller structures allows us to successfully kill cancerous cells.

3.5. Case study 2 Due to the lack of clinical facilities for cancer chemotherapy, we were unable to experimentally validate the presented method. Nevertheless, we attempt here to show the feasibility of the proposed method to control a complex biological system: Purine Metabolism Pathway. Purine metabolism, which is an essential process of the production of nucleotides and building blocks for synthesis of DNA and RNA, is crucial in living cells. Imbalances in the enzymic pattern of purine metabolism greatly affect the transformation and/or tumor progression. Hence, regulation of the purine metabolism is important for cancer therapy. We are interested here in showing the feasibility of our proposed method to control the following purine metabolism pathway model

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

123

presented in [25]: 0.65 0.7 −0.165 −0.067 0.036 −0.237 x18 − 31.3335x1.27 x4 x 6 x8 x˙ 1 = 0.898x1−0.03 x4−0.04 x8−0.04 x17 1 x2 0.085 0.075 −0.0092 x13 x15 x18 ,

x˙ 2 = 1.629x10.72 x−0.28 x0.319 x7 − 0.03x8−0.024 x0.144 x−0.06 − 5.584x0.454 x−0.156 4 13 2 2 18 4 x−0.012 x80.126 x−0.11 , 7 18 −0.05 x˙ 3 = 3.5932x20.4 x−0.24 x0.2 − 66544.7x0.99 x−0.95 , 8 x18 3 4 4 0.988 x˙ 4 = 0.082082x10.00025 x30.004 x4−0.00418 x0.0023 x0.00037 x11 − 614.27x40.054 x5−0.0045 5 6 0.00087 −0.00028 x80.128 x9−0.00003 x10 x18 ,

x˙ 5 = 7.2067x40.2 x−0.6 − 9.006x50.368 , 5 , x˙ 6 = 0.29x50.9 − 220.14x10.495 x4−0.79 x0.748 6 x˙ 7 = 1.2823x20.15 x7−0.09 x8−0.03 − 0.3738x40.12 x0.16 7 , x0.996 x0.001 − 408.11x2−0.000057 x0.0498 x˙ 8 = 0.048842x10.0034 x40.00014 x70.00019 x−0.0034 11 15 4 8 −0.000029 −0.001 x−0.00029 x80.133 x9−0.00009 x10 x18 , 7 0.017 0.98 x˙ 9 = 0.00229x40.002 x−0.0058 x10 x12 − 3.3x0.431 x0.324 , 9 10 9 0.986 x˙ 10 = 0.00151x80.0058 x9−0.017 x−0.0058 x12 − 2.2655x0.413 x0.34 9 10 , 10

x˙ 11 = 1024.195x40.05 x0.13 − 0.115389x11 , 8 0.33 x˙ 12 = 5.5062x90.42 x10 − 0.0032545x12 ,

− 8.2827x0.817 x2−0.66 x0.528 , x˙ 13 = 0.1033x20.43 x40.409 x90.04 x−0.19 1 13 18 0.343 0.236 x˙ 14 = 0.7234x13 x15 − 0.9325x0.569 , 14 0.02 −0.33 x˙ 15 = 0.2615x0.881 x10 x18 − 136.6x0.925 x−0.925 x0.438 , 8 10 15 8 0.55 x˙ 16 = 0.949x14 − 0.00008744x2.21 16 ,

(27)

where xi , i = 1, 2, . . . , 16 are the dependent variables, while x17 and x18 are the independent variables. Each variable involved in the above purine metabolism pathway model is defined in [25, p. 347]. It has been outlined in [25] that parameters of this purine metabolism pathway model have been identified using experimental information such as steady-state concentrations, steady-state fluxes, and kinetic characteristics (Michaelis constants and maximal rates). Since such a purine metabolism pathway model has been obtained from experimental information, we use it here for data generation for control design purposes. A key metabolite variable in purine biosynthesis is the PRPP (PhosphoRibosylPyroPhosphate) which should be maintained at a desired level. The PRPP variable, x1 , is affected by variations in the Ribose-5-Phosphate (x17 in the model). Hence, the PRPP variable x1 is used here as an output variable, while the derivative of the R5P variable (i.e. u(t ) = x˙ 17 ) is considered as a control variable. Initial values of all variables in (27) are as follows: x1 (0 ) = 5, x2 (0 ) = 100, x3 (0 ) = 0.2, x4 (0 ) = 2500, x5 (0 ) = 4, x6 (0 ) = 1, x7 (0 ) = 25, x8 (0 ) = 400, x9 (0 ) = 6, x10 (0 ) = 3, x11 (0 ) = 28600, x12 (0 ) = 5160, x13 (0 ) = 10, x14 (0 ) = 5, x15 (0 ) = 5, x16 (0 ) = 100, x17 (0 ) = 18. x18 (t ) = 1400 is kept here fixed ∀t. Our purpose, in this example, is to design a PID controller, using the methodology presented in the previous section, able to provide suitable control signal, u(t), that brings the PRPP variable, x1 , from its initial value, x1 (0 ) = 5, to the target value xd1 = 10. In addition to bringing the PRPP to the desired value (very small steady state error), its response should be with small overshoot and settling time of about 20min. Moreover, the control input, u(t), should not exceed 3. It should be mentioned that the above purine metabolism pathway model will be used here only to generate data for control design purposes, meaning that the control design process relies only on a set of data without assuming the availability of the model. Based on a set of data obtained from Eqs. (27), the application of the proposed control method results in the following PID controller parameter values:

Kp = 0.5806, Ki = 0, Kd = 0.0147.

(28)

Following the same procedure for stability analysis outlined above, it has been verified that the designed PID controller with parameters obtained above stabilizes the closed-loop system. In order to evaluate the performance of the obtained PID controller, it is placed in the closed-loop system as shown in Fig. 11. In this example, a white noise, n(t), has been introduced and added to the control input, u(t), in order to check the ability of the controller to reject the effect of such a noise. Fig. 12 shows simulation result obtained when a set-point, xd1 = 10, is applied. According to such a figure, it is clear that the response of PRPP, x1 (t), is without overshoot, with settling time of about 20.6 min, and it is maintained at the desired level of xd1 = 10 at the steady-state. It can also be seen from such a figure that the control input, u(t), provided by the PID controller, and

124

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

n(t)

x1d(t)

+ (t)

-

+ PID controller u(t)

+ ^ u(t)

Purine metabolism pathway

x1 (t)

Fig. 11. Closed-loop control system.

12

Output x1(t)

10 8 6 4

0

10

20

30

40 50 Time [min]

20

30

40 50 Time [min]

20

30

40 50 Time [min]

Control input u(t) 3 2 1 0 -1 0

10 ^ u(t)=u(t)+n(t)

4 3 2 1 0 -1

0

10

Fig. 12. Time responses obtained.

the control input, uˆ (t ), applied to the purine metabolism pathway model, are below the limit value 3. Hence, the controller has played its role in meeting the desired performance measures in presence of noise. 4. Conclusion In this paper, the control design problem for cancer chemotherapeutic treatment is investigated. We have proposed a measurement-based control design method that has been used to efficiently kill cancerous cells while avoiding toxic side effects. Such a proposed control method utilizes a set of measured data to design two PID controllers capable to maintain tolerable levels of drug dose, drug concentration and toxicity. It has been shown that parameters of each PID controller can be obtained by solving a constrained optimization problem using Genetic Algorithms. The key advantage of the proposed cancer control technique with respect to the existing methods is that no mathematical model is required for the design process, and its ability to automatically correct any deviation of the drug concentration response and/or toxicity response. A simulation example has been presented to illustrate and validate the proposed cancer control method. Simulation results obtained have effectively demonstrated the efficacy of the proposed control approach in cancer chemotherapeutic treatment. The robustness of the proposed control scheme for cancer chemotherapy treatment with respect to the plant parameter variation has also been shown. Moreover, a comparative performance analysis has shown that our results outperform those reported in other existing cancer chemotherapy techniques. To show the feasibility of the presented method, another example to control a complex biological system (purine metabolism pathway) has also been presented.

S. Khadraoui et al. / Information Sciences 333 (2016) 108–125

125

Acknowledgments This work was made possible by NPRP grant NPRP09-1153-2-450 from the Qatar National Research Fund (a member of Qatar Foundation). The statements made herein are solely the responsibility of the authors. References [1] S. Algoul, M.S. Alam, M.A. Hossain, M.A.A. Majumder, Multi-objective optimal chemotherapy control model for cancer treatment, Med. Biol. Eng. Comput. 49 (2011) 51–65. ˚ [2] K. Aström, T. Hagglund, PID Controllers: Theory, Design, and Tuning, second ed., Instrument Society of America, 1995. [3] B. Bojkov, R. Hansel, R. Luus, Application of direct search optimization to optimal control problems, Hung. J. Ind. Chem. 21 (1993) 177–185. [4] A. Besançon-Voda, Iterative auto-calibration of digital controllers, Methodol. Appl. Control Eng. Pract. 6 (1998) 345–358. [5] E. Carrasco, J.R. Banga, Dynamic optimization of batch reactors using adaptive stochastic algorithms, Ind. Eng. Chem. Res. 36 (1997) 2252–2261. [6] A. Floares, C. Floares, M. Cucu, L. Lazar, Adaptive neural networks control of drug dosage regimens in cancer chemotherapy, in: Proceedings of the International Joint Conference on Neural Networks, Portland, Oregon, 2003, pp. 154–159. [7] J. Holland, Adaptation in Natural and Artificial System, University of Michigan Press, Ann Arbor, 1975. [8] L.H. Keel, S.P. Bhattacharyya, Controller synthesis free of analytical models: Three term controllers, IEEE Trans. Automatic Control 53 (2008) 1353–1369. [9] S. Khadraoui, H. Nounou, M. Nounou, A. Datta, S.P. Bhattacharyya, A model-free design of reduced-order controllers and application to a DC servomotor, Automatica 50 (2014) 2142–2149. [10] L.H. Keel, S.P. Bhattacharyya, A Bode plot characterization of all stabilizing controllers, IEEE Trans. Automatic Control 55 (2010) 2650–2654. [11] Y. Liang, K.S. Leung, T.S.K. Mok, A novel evolutionary drug scheduling model in cancer chemotherapy, IEEE Trans. Inf. Technol. Biomed. 10 (2006) 237–245. [12] Y. Liang, K. Leung, T. Mok, Evolutionary drug scheduling models with different toxicity metabolism in cancer chemotherapy, Appl. Soft Comput. 8 (2008) 140–149. [13] Y. Liang, K.S. Leung, S.K. Mok, Evolutionary drug scheduling models for cancer chemotherapy, in: Proceedings of the International Conference GECCO, Seattle, 2004, pp. 1126–1137. [14] R. Luus, F. Hartig, F. Keil, Optimal drug scheduling of cancer chemotherapy by direct search optimization, Hung. J. Ind. Chem. 23 (1995) 55–58. [15] R.B. Martin, Optimal control of drug administration in cancer chemotherapy Ph.D. thesis, School of Computer & Information Sciences, University of Western Australia, Perth, Australia, 1991. [16] R.B. Martin, Optimal control drug scheduling of cancer chemotherapy, Automatica 28 (1992) 1113–1123. [17] R. Martin, K.L. Teo, Optimal Control of Drug Administration in Chemotherapy Tumour Growth, World Scientific, Singapore, 1994, pp. 1–10. [18] N. Mohsenizadeh, S. Darbha, L.H. Keel, S.P. Bhattacharyya, Model-free synthesis of fixed structure stabilizing controllers using the rate of change of phase, in: IFAC Conference on Advances in PID Control, Brescia, Italy, vol. 2, 2012, pp. 745–750. [19] R.E. Precup, S. Preitl, PI and PID controllers tuning for integral-type servo systems to ensure robust stability and controller robustness, Electrical Eng. 88 (2006) 149–156. [20] J.C. Spall, J.A. Cristion, Model-free control of nonlinear stochastic systems with discrete-time measurements, IEEE Trans. Automatic Control 43 (1998) 1198– 1210. [21] R.R. Sumar, A.A.R. Coelho, L. dos S. Coelho, Computational intelligence approach to PID controller design using the universal model, Inf. Sci. 180 (2010) 3980–3991. [22] K. Tan, E.F. Khor, J. Cai, C.M. Heng, T.H. Lee, Automating the drug scheduling of cancer chemotherapy via: evolutionary computation, Artif. Intell. Med. 1 (2002) 908–913. [23] S. Tes, Y. Liang, K.S. Leung, K. Lee, T.S. Mok, A memetic algorithm for multiple-drug cancer chemotherapy scheduling optimization, IEEE Trans. Syst. Man Cybern. 37 (2007) 84–91. [24] S. Venkatraman, G.G. Yen, A generic framework for constrained optimization using genetic algorithms, IEEE Trans. Evol. Comput. 9 (4) (2005) 424–435. [25] E.O. Voit, Computational analysis of biochemical system, in: A Practical Guide for Biochemists and Molecular Biologists, Cambridge Univ. Press, 2000. [26] D. Wang, P. Bao, Robust impulse control of uncertain singular systems by decentralized output feedback, IEEE Trans. Automatic Control 45 (2000) 795–800.