Adaptive non-linear control for cancer therapy through a Fokker-Planck observer.

In recent years, many efforts have been made to present optimal strategies for cancer therapy through the mathematical modelling of tumour-cell population dynamics and optimal control theory. In many cases, therapy effect is included in the drift term of the stochastic Gompertz model. By fitting the model with empirical data, the parameters of therapy function are estimated. The reported research works have not presented any algorithm to determine the optimal parameters of therapy function. In this study, a logarithmic therapy function is entered in the drift term of the Gompertz model. Using the proposed control algorithm, the therapy function parameters are predicted and adaptively adjusted. To control the growth of tumour-cell population, its moments must be manipulated. This study employs the probability density function (PDF) control approach because of its ability to control all the process moments. A Fokker-Planck-based non-linear stochastic observer will be used to determine the PDF of the process. A cost function based on the difference between a predefined desired PDF and PDF of tumour-cell population is defined. Using the proposed algorithm, the therapy function parameters are adjusted in such a manner that the cost function is minimised. The existence of an optimal therapy function is also proved. The numerical results are finally given to demonstrate the effectiveness of the proposed method.

1 Introduction

Since cancer has been one of the main reasons for deaths over the last decades [1], many attempts have been made to find an effective treatment. The mathematical modelling of tumour growth has interested many researchers. In practice, the drug dosage is obtained either from empirical data or clinical/preclinical assays. However, several factors, such as cost, time of disposal, and the difficulties of doing clinical tests, affect test results. Optimal control plays an important role in the validation and even the prediction of treatment strategy. Thereby, it decreases the cost of disease treatment. In the last decades, many attempts have been made to introduce effective treatment strategies through optimal control theory [2-7]. In [2, 3], a cost function based on measures obtained from the tumour‐cell population is presented, and an optimal dosage distribution of anti‐cancer drugs has been computed to minimise the tumour size. Bratus and Chumerina [4] considered the logistic model for tumour growth with the aim to determine an optimal treatment strategy. The suggested strategy includes the determination of the drug amount to minimise the number of tumour cells. Here, approximate and accurate solutions of the optimisation problem through the dynamic programming method have been obtained from the corresponding Hamilton–Jacobi–Bellman equation. El‐Gohary [5] considered optimal control and chaotic issues of tumour model with or without drugs, and studied the stability of its equilibrium states. A stability analysis of the model showed that it behaves chaotically for some values of model parameters. In the paper, optimal drug dosage values have been calculated by minimising the Hamiltonian function to control equilibrium states.

Genetic algorithm (GA) is an effective tool to seek an optimal or a near‐optimal solution to many optimisation problems. This algorithm has a higher performance compared with traditional optimisation methods, like the gradient approach [8]. GA can be considered as a multidirectional search method to handle optimisation problems and to find their solutions [8]. In many papers, GA uses medical images for tumour detection. Examples include diagnosis of pancreatic cancer [9], breast cancer [10], brain cancer [11] and liver cancer [12].

One of the important aims of chemotherapy is to kill a large number of cancer cells in a specific treatment period. This is why therapy scheduling is an important issue in chemotherapy. The authors in [13, 14] showed that GA is a useful tool to solve multidimensional and multi‐constraint cancer chemotherapy optimisation problems. The authors in [15, 16] presented a modified optimal control model for drug scheduling in the process of cancer chemotherapy. They defined a performance index by which an optimal drug scheduling can be determined by the use of adaptive elitist population‐based GA. The outcome of this research is consistent with the results of clinical treatments. Badakhshan and Khaloozadeh [17] repeated the stochastic method of [18] for optimal drug scheduling in chemotherapy by the use of GA.

The papers we reviewed and most of the other research works that used optimal control theory to determine an optimal treatment strategy considered a deterministic model for tumour growth. Such models do not give an appropriate description for tumour growth behaviour. In deterministic models, fluctuations or disturbances related to tumour dynamics are not taken into account. These fluctuations are unknown and immeasurable. In addition, there are sensible differences between clinical data and theoretical predictions due to environmental disturbances. To study the impact of environmental fluctuations on tumour growth, its model must be stochastic. One of the stochastic models that have been widely used in recent research works on tumour growth and cancer therapy is the Gompertz model [19-26]. In current research works, the treatment factor is considered implicitly within the drift term of the Gompertz model. Having fitted the model with empirical data, and by the use of statistical methods, such as linear regression [25] and the maximum likelihood method [20, 23, 25, 26], the parameters of the drift term are estimated. Consequently, the treatment strategy is determined.

In this paper, like many others [20, 21, 25, 26], the Gompertz model has been used to simulate the behaviour of tumour‐cell population. The therapy effect has been considered as an external variable to model. Since the Gompertz model is non‐linear, the model output (tumour‐cell population) will not be Gaussian. Hence, the control of first and second moments will not be sufficient. This is why the higher order moments of the output must be controlled. We know that all the moments of a stochastic process can be extracted from its probability density function (PDF), and PDF describes all the characteristics of a stochastic process. In fact, by controlling the PDF of tumour‐cell population, we can control all its moments. Moreover, to give a clear sense of decreasing tumour‐cell population, a Gaussian desired PDF is considered, in which both mean and variance values decrease exponentially by elapsing the time. The method of controlling the PDF of tumour‐cell population in this paper is similar to the method of [27, 28]. The differences between our paper and the two mentioned works are: the use of path integral method to solve Fokker–Planck equation, the definition of desired PDF for tumour‐cell population, the use of a time‐variant logarithmic function as control input instead of a constant function, and applying GA for optimisation purpose in each time window.

In this paper, we considered the therapy effect as an external control in which its parameters are optimised using a predictive control technique. The cost function of this optimisation is defined based on the difference between the PDF of tumour‐cell population and a desired PDF, where the PDF of the tumour‐cell population is computed by means of a Fokker–Planck non‐linear observer. The innovations of this paper are listed below:

In many research works [20, 23, 25, 26], the therapy effect has been included in the drift term of the stochastic Gompertz model. Then, by fitting the model with empirical data by the use of statistical methods like linear regression and maximum likelihood methods, the parameters of therapy function have been estimated. These works have not presented any algorithm to determine the optimal parameters of therapy function. In this paper, a logarithmic function is entered in the drift term of the Gompertz model to consider the therapy effect. The therapy function parameters are then predicted and adaptively adjusted by using the proposed control algorithm.

The growth of the tumour‐cell population is a stochastic process. To control it, the process moments must be monitored and manipulated. This paper employs the PDF control approach because of its ability to control all process moments. Such an approach has not been used in the reported research works to control the tumour‐cell population.

A non‐linear stochastic observer based on the Fokker–Planck equation is proposed to determine the PDF of tumour‐cell population at any instance of time.

The existence of an optimal therapy function is proved.

A desired PDF defined in this paper is theoretically meaningful. This PDF shows the decrease in tumour‐cell population with elapsing time. Also, an appropriate ‘desired PDF’ can be defined by adjusting the parameters of the PDF according to clinical observation.

The organisation of the paper is as follows. In Section 2, we will deal with the Gompertz model of tumour growth for which the therapy effect was added as an external control. A non‐linear Fokker–Planck observer is then designed for this stochastic model. In Section 3, a cost function has been defined in which the Fokker–Planck observer has been considered as a cost function constraint. In Section 4, we study the existence and uniqueness of our optimal solution. In Section 5, real‐coded GA is briefly described. In Section 6, the path integral method, by which the Fokker–Planck observer is numerically solved, is explained. In Section 7, a receding horizon model predictive control (RH‐MPC)‐based scheme is presented to determine the optimal parameters of in each time window. Here, we have used the real‐coded GA. Section 8 gives a numerical example to verify the performance and correctness of our approach. A desired PDF with Gaussian distribution has been designed. By the use of our suggested control algorithm, the optimal parameters of the therapy function are obtained. The paper is finally terminated with a short conclusion.

2 Modelling the system

2.1 Gompertz model

The Gompertz model is an umbrella for growth modelling of a variety of scientific issues like economic growth, biological growth, energy growth, and the expansion of greenhouse gases during climate change [20]. The stochastic Gompertz model can be obtained from the deterministic Gompertz model by adding noise. It is an appropriate model to describe tumour growth behaviour.

In many research works in which the Gompertz model has been used to describe tumour growth, the therapy effect has been entered either as a time‐dependent linear function [21, 22] or as a time‐dependent logarithmic function [20, 21, 23] into the drift term of the model. The model parameters are then estimated by the use of statistical methods, such as linear regression, maximum likelihood, and model fitting with empirical data. For instance, Moummou et al. [20] presented a Gompertz‐based stochastic model for tumour growth, in which the drift term is dependent on two time‐variant functions. The first function models the immunologic endogenous treatment factor. The second models the dynamic of external controllable treatment. In this work, the probabilistic characteristics of the model have been obtained by using its related It differential equation followed by the estimation of model parameters through the maximum likelihood method. Like [20, 25, 28], this paper has used the Gompertz model to describe tumour growth In (1), denotes the tumour‐cell population at time t, and represents the standard Wiener process. The parameters , and are positive constants representing birth rate, death rate, and stochastic oscillation width, respectively. Like [20, 25, 28], the time‐dependent function is used to insert the therapy effect into the model. Thus, (2) is obtained In (2), models the effect of therapy. It can be considered as either a constant function or a time‐dependent function. If is considered as drug concentration, , which means that the concentration of drug is constant during the therapy period. In constant treatment, tumour‐cell population converges into a non‐zero value. Such a therapy model was sought in this study by which the tumour‐cell population of (2) converges to zero. In [26], the function has been proposed to therapy model . Although the tumour‐cell population goes down to zero by this selection, this therapy scheduling is not acceptable to many clinicians. The authors in [19, 29] showed that to have a therapy scheduling with the capability of (i) asymptotically converging of tumour‐cell population to zero, and (ii) time‐increasing behaviour, a logarithmic therapy function is a very appropriate selection.

The present paper uses (3) to model the effect of therapy where e is the Euler number and and are parameters that must be adjusted. For a patient, this therapy is more tolerable than the one with linear therapy function.

2.2 Fokker–Planck equation for Gompertz model

The Fokker–Planck equation was introduced to describe the Brownian motion of particles [30]. This equation is, in fact, a parabolic partial differential equation of a given process model that indicates the evolution of the process PDF. It can also be considered as a non‐linear observer by which the PDF of a given process can be computed for a given initial PDF.

Suppose that denotes PDF of process at time t and in position . Then, the Fokker–Planck equation for (2) will be as follows [30]: In this paper, (4) is defined over , where is a bounded space.

3 Therapy effect prediction based on Fokker–Planck observer

3.1 Description of optimal control problem

The main aim of this paper is to introduce a methodology (based on MPC) to compute the optimal parameters of therapy function. To do this, the PDF of tumour‐cell population must be controlled. In this way, the entire stochastic characteristic of will be controlled. Here, the PDF must be computed at any instance of time. Since the stochastic model of tumour growth is known, a PDF observer is needed to compute its PDF.

We designed a Fokker–Planck observer for our stochastic process , that if solved the process PDF (deterministic process) will be obtained. Hence, the Fokker–Planck observer can be interpreted as a mapping vehicle from the stochastic domain to deterministic domain. To define our control problem, the time interval is divided into time windows . The size of each window is in which N is a positive integer and for . It is assumed that the control law for each time window is a logarithmic function (3) and the initial process PDF at is known in advance. The aim of control is to determine the value of u (by adjusting its parameters and ) for each time interval so that the process converges from the initial PDF at time towards a desired PDF at time . To do this, an optimal control problem is defined as follows. By solving this problem, control law u can be determined subject to where . In (5), the positive parameter denotes control weight. It should not be selected so small that the patient be unable to tolerate such harsh therapy and it should not be selected so large that the therapy be practically ineffective.

4 Existence and uniqueness of optimal solution

In this section, the existence and uniqueness of our optimal solution will be studied. The optimisation problem given by (5)–(8) shows a bilinear control problem in which the dependence of state y on u is not linear. Thus, this optimisation problem is non‐convex. The existence of an optimal solution can be proven by the use of [31], but because our optimal control is governed by a bilinear parabolic partial differential equation, the uniqueness of the optimal solution cannot be verified using [31]. Instead, the method of [32] is applied to prove the uniqueness of the optimal solution.

4.1 Setting of the problem

Like [32], we define , its dual space and their duality pairing . Let be identified with its dual . Then, we have the following Gelfand triple: It is assumed that all these embeddings are compact, continuous, and dense. We now define the following space: where . In the following equations, for simplicity, the time interval will not be written. Equations (6)–(8) can be rewritten as follows: where Operators A, B, and C are defined as follows:

, is a linear continuous operator and , where is equivalent to the norm defined in [32], and there exists a positive constant such that [31]

is a linear and continuous operator. Then we have

and is a linear and continuous operator. Then, there exists a positive constant such that [31]

is continuously differentiable on .

See Appendix 1. □

If and , then equation admits a unique solution y in . Also, we have and .

The mapping , where y is the solution of , is Frechet differentiable, and satisfies the following equation: where .

4.2 Existence of optimal solution

Assuming the cost function of (5) in which y is the unique solution of the following equation: Then there exists a pair such that is a solution of and minimises J in U.

The function is differentiable and we have in which p is the solution of the following adjoint equation:

4.3 Uniqueness of optimal solution

There are positive constants so that Then, by assuming Propositions 1 and 3, and adequately small initial PDF, the optimal solution is unique.

See Appendix 2.□

5 Real‐coded GA

GA is a search algorithm built on the basis of genetics and natural selection [33]. It has been widely used for optimisation. GA consists of three evolutionary operations: selection, crossover, and mutation. In this algorithm, a set of chromosomes is first selected randomly from a given search space to create a primary population. From this population, based on a predefined selection rule, the chromosomes that maximise a fitness function are selected as parents. In the crossover phase, a pair of chromosome exchanges meaningful information to generate two offspring. In the mutation phase, with random handling of a chromosome, a new chromosome is generated.

In traditional GA (binary‐coded GA), bit strings are used to represent chromosomes. This algorithm suffers from low convergence due to its long chromosome structure, especially in optimisation problems involving several parameters [34]. However, for engineering problems, real‐coded GA is more convenient than binary‐coded GA. The reason is that changing from real numbers to binary digits may lead to loss of accuracy [8]. This is why in this paper, we used real‐coded GA.

In the real‐coded GA used in this paper, represents a set of possible solutions to a given optimisation problem. Set is called ‘chromosome’, and for and is named ‘gene’. The search space of chromosome in real‐coded GA is defined as follows: It is also assumed that N indicates the number of chromosomes in the population space, and parameters and represent the rate of crossover and mutation processes, respectively.

6 Observer design

The PDF of tumour‐cell population must be available at any instance of time to control it. To obtain these PDFs, the Fokker–Planck equation is used as a non‐linear observer. The Fokker–Planck equation is obtained from the stochastic model of tumour growth. Solving this equation gives the instantaneous PDF. Various numerical methods have been developed to numerically solve the Fokker–Planck equation of which Monte Carlo [35], finite difference [36], spectral approximation [37, 38], and path integral [39, 40] methods are well known. Among these, the path integral method is simple, effective, and accurate. It computes the process PDF at any time given that its preceding PDF is known.

In this paper, the path integral method is used to solve the Fokker–Planck equation of stochastic model of tumour growth. The basis of the path integral method is the Chapman–Kolmogorov equation. This equation is written as follows for the Markov processes: Equation (13) can be considered as a recursive method to determine the process PDF for small values of . This means that if the conditional PDF (CPDF) and process PDF are known at time t, the process PDF at time can be easily computed for small values of . For Gaussian white noise, the process PDF can be obtained by the use of short‐time Gaussian approximation [38]. Using this approximation for (4), can be determined as follows:

6.1 Numerical solution

To numerically solve the Fokker–Planck equation through the path integral method, it must be discretised both in the time and space domains. To do this, like [27], a uniform spatial mesh is considered for spatial domain in which represents the mesh size and denotes the set of points located inside the mesh. It is defined as follows: The space‐time mesh is also defined as follows: where denotes the size of time step and is the number of time steps.

On grid , defines a grid function at at time t and . By this definition, the PDF of tumour‐cell population at time is defined as follows: where is the normalised CPDF for . Since the sum of columns of CPDF matrix must be one, will be normalised as follows: In which CPDF is obtained by the use of short‐time Gaussian approximation as follows:

7 Controller design

Our control aim in this paper is to find an optimal value for parameters of within such that the process of tumour growth converges to the predefined desired PDF at time when starting from the initial PDF at . To do this, an appropriate cost function in the form of (5) must be defined and minimised within each time window. To compute the PDF of tumour growth, and insert it into the selected cost function, the Fokker–Planck observer, defined in (4), is used.

Like [27], to determine optimal values for therapy function parameters, we have used an RH‐MPC scheme. So, the optimal problem in time interval is solved by the use of real‐coded GA. The value of PDF at time is then used as the initial PDF value to resolve the optimal problem within . This procedure is repeated until the last time window. The block diagram of this approach, to predict the optimal parameters of , is shown in Fig. 1.

Fig. 1

RH‐MPC based on Fokker–Planck observer to predict therapy function parameters

In the following, the RH‐MPC scheme defined in [27] is explained briefly. It is completed in six steps.

Step 1: Set the initial values – the initial PDF , for and desired PDF for .

Step 2: In time interval use subroutine 1 to solve to obtain an optimal value to therapy function u.

Step 3: Solve the Fokker–Planck equation by the use of path integral method to compute PDF for optimal control of u.

Step 4: Set the obtained PDF as the initial condition of Fokker–Planck observer in the next time window.

Step 5: If , set and go to step 2.

Step 6; Otherwise, end.

Subroutine 1: Real‐code GA

Step 1: Set the initial values of population size N, search space , maximum iteration and parameters , , , , and .

Step 2: Set .

Step 3: Generate randomly a primary population with N chromosomes in which all genes fall into .

Step 4: By using (3), compute the value of for each chromosome.

Step 5: Compute PDF for of each chromosome by the use of (17).

Step 6: Compute the cost function of each chromosome by the use of (5).

Step 7: Compute the selection probability of each chromosome by the use of the following equation:

In (20), parameters and represent the value of ‘pressure’ and ‘cost function for chromosome ’, respectively. Select chromosomes from the main population to carry out the crossover operation by the use of the Roulette wheel selection defined in [41]. The sign denotes the bracket operator and is the crossover rate.

Step 8: Create an offspring population by applying the arithmetic crossover operator to parents. Suppose that and denote the chromosomes of parents as

Then, offspring and are created by using the arithmetic crossover operation. The result is shown in the following equation:

In (22), is a random variable with uniform distribution in interval .

Step 9: Select randomly chromosomes from the primary population (note that ). Generate the mutated population by applying the Gaussian mutation operator to the selected chromosome. If denotes the selected gene for mutation operation, and shows the mutated gene, then (23) can be written as in which denotes the standard deviation of Gaussian distribution and shows a standard Gaussian distribution with average zero and variance 1 [33].

Step 10: Construct the new main population by replacing chromosomes with higher fitness values.

Step 11: If set and go to step 3 else go to the next step.

Step 12: End.

8 Simulation results

This section gives a numerical example to demonstrate the performance of or suggested control method. Consider a Gompertz stochastic model without any control input. It is described as Like [20], , , and are assigned. Moreover, is defined for tumour‐cell population. With these values, the Gompertz model is simulated in the Matlab environment.

Fig. 2 indicates ten sample paths of the process for the first 50 weeks without therapy. As the figure shows, without any therapy, the tumour‐cell population expands as time elapses. Fig. 3 indicates the initial PDF of tumour‐cell population with delta‐Dirac distribution. By the use of this PDF and numerical solving of Fokker–Planck equation (4), the evolution of the tumour‐cell population PDF is obtained. Fig. 4 depicts the PDF of tumour‐cell population for weeks 6, 12, 18, 24, and 30 when there is no control input.

Fig. 2

Ten sample paths of the process without therapy

Fig. 3

Initial PDF of tumour‐cell population

Fig. 4

PDF of tumour‐cell population for weeks 6, 12, 18, 24, and 30 without therapy

In the proposed control approach, we seek a control input by which the PDF of the tumour‐cell population converges to an expected ‘desired PDF’ by elapsing the time. It is assumed that the tumour‐cell grows continuously in the first 30 weeks during which the tumour‐cell population is maximised. During this period, no therapy has been carried out (i.e. therapy is started from week 30). The aim is the use of our proposed control algorithm to determine the optimal weekly therapy parameters so that in the next 15 weeks, the PDF of tumour‐cell population converges to the expected desired PDF.

The distribution of desired PDF is assumed to be Gaussian as in which and for .

The desired PDF for the treatment period (15 weeks) has been shown in Fig. 5. The values of its mean and deviation are selected such that, by elapsing the time, the width of the desired PDF and its mean are exponentially decreased. By tracking Fig. 5, all the moments of are decreased and controlled.

Fig. 5

Desired PDF graph

To obtain the optimal values of therapy parameters, the proposed algorithm is applied to the Gompertz model to ensure that in each time window (a week), the process PDF converges to its desired values. For this propose, and are selected. Moreover, for subroutine 1, it is assumed that . Search space is defined as follows: We set , , , , , and , respectively, for the parameters of subroutine 1. In addition it is assumed that , and .

Fig. 6 shows the desired PDF and the computed PDF graphs, and Fig. 7 indicates the optimal therapy function. As Fig. 6 shows, the PDF of tumour‐cell population converges to its desired values within 15 weeks. Fig. 8 indicates the optimal values of coefficients and calculated by real‐coded GA. Fig. 9 indicates ten sample paths of the process in the presence of optimal therapy. As the figure shows, by using the therapy after 30th weeks, the tumour‐cell population will be decreased. It can be seen that the mean and variance of the tumour‐cell population in the presence of therapy both exponentially decrease comparing with the tumour‐cell population without therapy (Fig. 2).

Fig. 6

Desired PDF and the computed PDF graphs

Fig. 7

Optimal therapy function graph

Fig. 8

Optimal values for coefficient and

Fig. 9

Ten sample paths of the process in the presence of therapy

9 Conclusion

The aim of this paper was to determine the parameters of therapy function by controlling the tumour‐cell population. Since the tumour‐cell population is a stochastic process, we were faced with a stochastic process control problem. We used the well‐known Gompertz model to describe the dynamics of tumour‐cell population. This model is non‐linear and the control of its first and second moments is insufficient for thorough control of a process. The higher order moments of the process must be controlled in this case. Since all the process moments can be extracted from its PDF, we used the process PDF as the control variable. We used a Fokker–Planck‐based non‐linear stochastic observer to determine the process PDF. The output of this observer was compared with a predefined ‘desired PDF’ and the result was minimised to predict and adjust the parameters of therapy function. To do this, a cost function was defined and a control algorithm was proposed to minimise that function. By using the proposed algorithm, the therapy function parameters were adjusted. The existence and uniqueness of an optimal therapy function was also proved. The numerical results are finally given to demonstrate the effectiveness of the proposed method.