OPTIMISATION MODELLING OF CANCER GROWTH

. Several computational models have been developed in the literature to describe the dynamics of the cell-cycle for the mammalian cell, in par- ticular for cancer cells, using both traditional and new techniques and yielding some positive results. In this paper, we discuss how to optimise model param- eters for these types of models and how this can serve to enhance numerical results. We pose the model parameter selection problem as an optimal param- eter selection problem on both a normal cell and cancer cell cycle of growth. Various possible objectives are discussed and we illustrate the process with some numerical results.


Introduction.
Cancer is now one of the biggest killers in western society, second only to cardiovascular disease, and it is still a largely unresolved issue in modern society [10]. Some of the key issues relating to cancer growth and its treatment can be addressed through the use of mathematical modeling, in particular computational modeling of cell growth dynamics.
At the subcellular level of growth, the cell cycle describes the activities that are present in a cell during its life. The activities occurring inside a cell and the signals it receives from its external environment are what determine its future. To be more precise, it has been recognized that there is a CPU of some sort in the nucleus, a cell-cycle clock, that analyzes the signals received by the cell and determines if the cell will continue through the cell cycle and undergo cell division, or become quiescent and proliferate at a later stage, or die [11].
The cell cycle is an intrinsic and complex process that essentially governs the fate of a single cell and encompasses four phases to reach division. The first of these phases is the G 1 or gap phase where cells make critical decisions about continuing in the growth cycle or entering into a non-active state, G 0 . A critical decision point, the "Restriction Point", is where the cell has to make a vital decision about its future. In other words, it has the option to continue through the cycle, go into G 0 or remain in G 1 . This is the only time at which the cell can make an autonomous decision about its fate. Research has suggested that this R point is where a cancer cell can form and it has also been proposed that this is true for many types of cancers [11].
The remainder of the phases are not the main focus of the model presented in this paper but will be explained briefly to show their importance and relevance during the normal cell cycle. The S or Synthesis phase is where DNA is replicated before entering into the second gap phase G 2 where the cell prepares for division. The final phase of the cell cycle is Mitosis or the M phase where the single cell divides and turns into two daughter cells.
The remainder of the paper is broken into 4 sections. The first section presents the model which we will use to illustrate the optimal parameter selection approach. We then give a brief introduction to the formulation of optimal parameter selection problems, including the notion of multiple characteristic time points. This formulation is then applied to the model we presented in Section 2 and results are discussed.
2. Model background. The model presented below looks at the G 1 /S transition phase and it is based on previous work by Alarcón, Byrne and Maini [1]. where x : average concentration of APC/Cdh1 complexes, m : maximum mass of a cell, z : average concentration of the protein p27, y : average concentration of complex formed from Cdk/Cyc complexes, u : average concentration of non-phosphorylating retinoblastoma protein (Rb).
The anaphase-promoting complex (APC) is a group of proteins that includes protein Cdh1. These are active during the late mitosis and G1 phases of the cell cycle. y represents a generic cyclin. The a's, b's, c's and d's all represent rate constants, J 3 and J 4 are Michaelis constants determined from experimental data, η is the growth rate parameter, B is a constant term relating to the oxygen tension, m * is the maximum mass of a cell and P represents oxygen tension.
3. Optimal parameter selection problem. Our motivation behind using optimal parameter selection techniques on the model by Alarcón et al [1] was the absence of any optimisation techniques used in their model development. The model parameter values that were used by Alarcón et al were based either on values found in literature or on estimates based on values determined by Dormann and Deutsch [2]. It is not clear, however, if these parameters, either estimates or literature based, were optimised to ensure a good match between the model and experimental data. It is also not clear if the estimates taken from Dormann and Deutsch [2] were optimised for model matching.
3.1. Generic problem. Much of the optimal parameter selection methodology we use here is based on the book by Teo, Goh and Wong [8]. We now give a brief description of a standard dynamic optimal parameter selection problem with accompanying canonical constraints based on this work. The general form encapsulates different constraints to allow the modelling of a range of real-life applications.
Consider the following dynamical system over a given time horizon [0,T ]: where is the vector of initial conditions for the state which may be dependent on the system parameter ζ,ẋ is the time derivative of the state x and f is a given function. We now define the following equality and inequality constraints in canonical form.
. . , N are both continuous real valued functions which are differentiable with respect to x and ζ. g i is also assumed to be a differentiable function on R and τ i is the characteristic time of the i th constraint. Furthermore, N e denotes the number of equality constraints while N is the total number of constraints. The aim is to minimise an objective functional subject to the dynamical system given by equations (6) and (7) and subject to the constraints (8). The objective functional is expressed in the form, where . . , s} is a compact and convex subset of R s . Here, a i and b i are real numbers for each i = 1, . . . , s. Also, Φ 0 and L satisfy the same assumptions as those used for the corresponding functions in (8). The above problem essentially constitutes a mathematical programming problem and can, in principle, be solved using non-linear programming techniques such as sequential quadratic programming (SQP). However, the presence of the dynamic constraints (6) means that the required gradients of the objective and constraint functionals need to be calculated in a round about manner [8].
For the problem of matching a dynamical model to experimental data, where x * i (t) are the observed output for i = 1 to 5, we can consider a range of objective functionals. One possible choice would be, The above objective function is assuming that all experimental states are observable over the entire time horizon. Clearly, in reality this is not the case.
A more realistic situation for experimental data is that data on the system state may only be available at a certain sequence of M observation times, τ i for i = 1, . . . , M . In this case, the objective should take the form where each h j , j = 1, . . . , m is some function of the model states that can be observed experimentally. The objective (11) does not fit into the general canonical form of (9) due to the presence of multiple characteristic times in the objective, τ i , i = 1, . . . , M . General forms of functionals of this type have been described in [7], [5] and [6]. Rory Martin [6] introduced the concept of Multiple Characteristic Time points to model treatment scheduling for cancer chemotherapy. Martin and Teo [7] presented a theoretical approach for dealing with multiple characteristic time points in objective and constraint functionals. The notation used by Martin and Teo is similar to what we have used above. In particular, they derived the gradients of such functionals and these were later incorporated into the optimal control software MISER3.3 [3]. The software MISER3.3 is used to generate the problem solution for our problem here. They also derived convergence results for solving control problems of this type by a sequence of numerical approximations.

First implementation.
To illustrate the parameter matching process, we first consider a model by Tyson and Novak [9] where we show how well our procedure works for the optimisation of parameters. The objective function that we chose to minimise is, where the system initial states are [0.5, 1.0, 1.8, 1.3, 0.7, 0.6] T . The terminal time is denoted by T and for our formulation has a value of 60, where time t is in minutes. x j (t) are the theoretical values and x * j (t) are the "experimental" values we obtained from a successful implementation of Tyson and Novak's model. The objective function (12) is subject to the dynamical system below, where the states are represented as , x 6 = m and the parameters are listed as follows, ζ 6 = k 4 , ζ 7 = k 5 , ζ 8 = k 5 , ζ 9 = k 6 , ζ 10 = k 7 , The Michaelis constants J 3 , J 4 , J 5 , J 7 and J 8 remain fixed in our formulation. The constant m * is the maximum size of an adult cell as it was for the Alarcòn model, and therefore will also remain fixed throughout the formulation.
[MAD] is a toggle and switches between 1 and a large number for the completion of chromosome alignment. A value of 1 for "on time" completion and a large number for a completion not on time. We will assume that in our formulation that chromosome alignment is "on time" and hence, [MAD] is fixed to a value of 1. The coefficient n is a Hill coefficient and has a value of 4, as published by Tyson and Novak [1]. The Hill coefficient represents the steepness of the relationship between concentration and response. As it is not closely related to the purpose of this research we will not give an in depth discussion relating to the Hill equation and its use in pharmacology. We use MISER3.3 to minimise the objective function (12) over the stated time horizon, using multiple time points τ 1 = 0, τ 2 = 10, τ 3 = 20, τ 4 = 30, τ 5 = 40, τ 6 = 50 and τ 7 = 60, where time is t minutes. The values returned by MISER3.3 are displayed in Table 1 and show either an exact match or very near match to those published by Tyson and Novak. A more detailed description of our results can be found in [4].

Second implementation.
We now give an objective function for the Alarcón et al model which we will use to simulate our numerical results.
where, for each state x = [x 1 , x 2 ] and ζ = [ζ 1 , ζ 2 , . . . , ζ 12 ] , we aim to minimise the objective functional (22) subject to equations (23) to (27) with initial conditions (28). Note that this once again constitutes an optimal parameter selection problem with multiple characteristic times, as outlined earlier in this section.
We chose to only use states x 1 and x 2 in the objective as these two states are the most affected during the G1 phase of a cancer cell cycle, in particular, the G1 phase up to the restriction point. The activity of these complex forming proteins is significant as the complexes can form cancerous mutations. The x 1 state is of particular interest, as the differential equation representing this state implements Michaelis-Menten type kinetics which can lead to varying results. We present the model from Alarcón et al [1] in canonical optimal parameter selection form. Here, the vector x = [x 1 , . . . , x 5 ] represents the state variables and vector ζ = [ζ 1 , . . . , ζ 12 ] represents the parameters we are optimising. The system is presented as follows:ẋ where the system initial states are, As in the model proposed by [1] the constant terms J 3 and J 4 represent Michaelis constants determined by biological observation. The constant variable m * still represents the maximum adult mass of a cell and constant P remains the oxygen tension. Descriptions of both the states and parameters used in our formulation are those given in Section 2. The states and parameters are listed below.
The parameter descriptions will relate to the parameters used in the Alarcón model previously described in Section 2, and are listed below, As we have no direct access to appropriate experimental data, we simply use the graphs in [1] to generate the 'experimental' state values x * j (τ i ), i = 1, 2, . . . , 11 and j = 1, 2. We chose to use 11 time points at which we 'measure' the states, x 1 and x 2 , leading to 11 multiple characteristic times. These are listed in the Table 2. Although the optimal parameter selection model we present is based on published data rather than experimental data, we can still produce a model that closely resembles data published by Alarcón et al. Note that the same technique could be used if actual experimental data was available. We will use the parameter values that were sent to us by Alarcón and associates [1] as our initial starting point for the optimal parameter selection technique (as seen in Table 3 under heading Alarcón). By optimising the parameters we are aiming to obtain a "better fit" that more accurately resembles the graphical data from Alarcón et al [1].
The optimised parameter values that we obtained using MISER3.3 are listed in Table 3 alongside the revised values given to us by Alarcón et al [1]. We were unable   Table 3 the only parameters which are close to the revised parameters by Alarcón et al are ζ 5 and ζ 6 . Some of the parameter values used by Alarcón et al [1] were taken from previously published work by Tyson and Novak [9]. We would expect that since the Tyson and Novak model was easily duplicated as seen in section 3.2 and that the use of their parameter values, in part, would be reasonable. After performing optimal parameter selection on the chosen parameters we are able to then put the parameters returned by MISER3.3 back into the system and plot the result.
The results we have obtained maybe due to one or more issues relating to this type of optimisation problem and specifically the problem we were analysing. Primarily, they are biological misinterpretation, incorrect parameter values, uncertainty in active states of proteins and machine accuracy and integration sensitivity. The biological misinterpretation has occurred with variable u, representing non-phosphorylated retinoblastoma protein, where its occurrence in the model is suppose to mimic the behaviour of non-phosphorylated Rb on Cdk activity [1]. Hypo-phosphorylation occurs during the G1 phase of the cell cycle up until the restriction point. After this, it enters into a hyper-phosphorylation state. To state that during G1 phase of the cell cycle Rb is non-phosphorylated indicates that Rb is in an unphosphorylated form which only occurs at the M/G1 transition of the cell cycle [11]. We now go on to discuss briefly the issues surrounding the machine accuracy and sensitivity with integration.
The equations (1) to (5), detailed above in Section 2, represent a stiff system of equations and therefore need to be dealt with carefully. This appears to be an indication of the sensitivity of the model to numerical precision of the machine used. The nature of the equations presents problems numerically due to the varying timescales associated with modelling biological phenomena.
The software MISER3.3 that we used to solve the optimal parameter selection problem encountered difficulties when using it's in-built integrator when attempting to solve the problem. For the dynamics presented by [1] the software would often switch to the in-built stiff integrator indicating that the system is indeed stiff.

4.
Conclusions. In this paper we have applied optimal parameter selection techniques to a known model in the field of cancer cell cycle growth. We have been able to show how using optimal parameter selection techniques we are able to achieve an optimal solution for the proposed problem. Our results showed that some of the published model parameters in Alarcon et al [1] were far from optimal when trying to replicate their numerical results. Despite this, we were unable to completely replicate their published results and we suggested a number of possible reasons for this discrepancy.
We have shown how dynamic cell cycle models can be effectively matched to experimental data by the use of an optimal parameter selection formulation and numerical optimisation tools. Existing papers on cell cycle models have not used this approach to date.