WO2010092379A1

WO2010092379A1 - A method, system and computer program for developing an interrogatable dynamic model of a breast cancer and burkitt's lymphoma

Info

Publication number: WO2010092379A1
Application number: PCT/GB2010/050203
Authority: WO
Inventors: Robert George Clyde
Original assignee: The University Of Abertay Dundee
Priority date: 2009-02-11
Filing date: 2010-02-10
Publication date: 2010-08-19
Also published as: EP2396740A1; GB0902292D0

Abstract

A method of developing an interrogatable model of a breast cancer in order to determine characteristics of the breast cancer which may aid in the diagnosis and treatment thereof; the method comprising the steps of : - defining a generic model of a cell having a conserved dynamic state, which model includes: - one or more mathematical equations which substantially represent the dynamic state of the cell; - modifying the or each of the mathematical equations of the generic model to represent a one or more changes caused by the breast cancer to the molecular interactions; - determining an initial condition of the generic model by an iterative process until the initial conditions meet a predetermined performance requirement; - solving the one or more equations to produce a plurality of solutions which represent the time-dependent concentrations of one or more selected molecular species; - testing the solutions to determine the best solution which fits a predetermined requirement at least partially established from qualitative clinical data of a patient.

Description

A METHOD, SYSTEM AND COMPUTER PROGRAM FOR DEVELOPING AN INTERROGATABLE DYNAMIC MODEL OF A BREAST CANCER AND BURKITT'S LYMPHOMA

Field of the Invention This invention relates to a method, system and computer program for developing an interrogatable dynamic model of a disease of a cell cycle; and, more particularly a method, system and computer program for developing an interrogatable dynamic model of a breast cancer and Burkitt's lymphoma.

Background to the Invention

The inherent inter-disciplinarity of modelling a disease such as cancer presents a number of unique challenges. Indeed, one of the most intractable problems is that of overcoming the barriers between traditionally non-overlapping disciplines (e.g. mathematics and biochemistry) to fruitfully marry the specialist skills and knowledge therein.

To attain a degree of veracity, it is clear that a model of a disease must be founded on its underlying biology. Cancer is a genetic disease manifested in failures of a cell's normal regulatory processes. Thus, a model of a cancer must incorporate cell cycle regulatory processes such as growth and division, growth restriction, survival and programmed cell death. However, the body of scientific literature on these subjects is vast, complex and sometimes ambiguous. Furthermore, whilst a large number of papers have been published on the subject of cell cycle regulation, these papers rarely present information in a form suitable for mathematical modelling. For example, whilst qualitative data (e.g. from blot analysis) is readily available, quantitative data (e.g. time course data) is more scarce. In an effort to circumvent these problems, a number of data-mining techniques have been applied to static biological data to search for patterns indicative of cancer. However, since these data-mining techniques do not attempt to analyse or characterise the dynamic state of a cell, they provide little insight into the underlying disease process of a cancer. In view of this limitation, a number of mathematical models have been proposed that dynamically model a small number of cancers. However, these models typically employ a number of assumptions (e.g. constant kinase concentration) that render the models unrealistic and misleading. Similarly, most, if not all, of the prior art models are non- interrogatable. In other words, the prior art models do not readily allow users to alter input (or other internal) variables (e.g. different drug dosages etc.) to enable the effect of changing these variables to be studied. This seriously limits the use of the prior art models for assessing potential treatment strategies.

As an aside, it will be realised that real biological systems tend to be extremely complex, featuring extensive connectivity between different regulatory components and numerous feedback loops. Both of these attributes provide redundancy and stability to biological systems. In contrast, many of the above-mentioned prior art models are oversimplified, omitting the above connectivities and feedback loops. As a result, these models are fragile and incapable of expressing the inherent robustness of real biological systems. Furthermore, prior art models often incorporate only a small part of a cell's regulatory apparatus. Thus, even the marginally interrogatable models are generally incapable of considering the effects of pleiotropy. This limits the use of these models for developing combinatorial multi-drug treatment regimes that simultaneously target a number of points in a cell's regulatory cycle. Summary of the Invention

According to a first aspect of the invention there is provided a method of developing an interrogatable model of a breast cancer in order to determine characteristics of the breast cancer which may aid in the diagnosis and treatment thereof; the method comprising the steps of:

- defining a generic model of a cell having a conserved dynamic state, which model includes: o one or more molecular species, the or each of which undergo molecular interactions within the cell with one or more other molecular species; o a set of parameters associated with the cell; o one or more mathematical equations which substantially represent the dynamic state of the cell, the or each equations comprising: - a one or more phase variables at least partially corresponding to the or each molecular species; ■ parameters to allow variations thereto to be investigated; and o one or more input functions which can be used to generate the variations;

- modifying the or each of the mathematical equations of the generic model to represent a one or more changes caused by the breast cancer to the molecular interactions;

- determining an initial condition of the generic model by an iterative process until the initial conditions meet a predetermined performance requirement;

- solving the one or more equations to produce a plurality of solutions which represent the time-dependent concentrations of one or more selected molecular species based on one or more input functions to determine the effect of the input function on the dynamic state of the cell; - testing the solutions to determine the best solution which fits a predetermined requirement at least partially established from qualitative clinical data of a patient; and

- selecting the best solution to demonstrate the effect of the input functions on the dynamic state of the cell with respect to the breast cancer.

According to a second aspect of the invention there is provided A system of developing an interrogatable model of a breast cancer in order to determine characteristics of the breast cancer which may aid in the diagnosis and treatment thereof; the system comprising:

- a model generating module adapted to define a generic model of a cell having a conserved dynamic state, which model includes: o one or more molecular species, the or each of which undergo molecular interactions within the cell with one or more other molecular species; o a set of parameters associated with the cell; o one or more mathematical equations which substantially represent the dynamic state of the cell, the or each equations comprising:

^■ a one or more phase variables at least partially corresponding to the or each molecular species;

^■ parameters to allow variations thereto to be investigated; and o one or more input functions which can be used to generated the variations;

- a modification module adapted to modify the or each mathematical equation of the generic model to represent a one or more changes caused by the breast cancer to the molecular interactions; - a determination module adapted to determine an initial condition of the generic model by an iterative process until the initial conditions meet a predetermined performance requirement;

- a calculating module adapted to solve the one or more equations to produce a plurality of solutions which represent the time-dependent concentrations of one or more selected molecular species based on one or more input functions to determine the effect of the input function on the dynamic state of the cell;

- a testing module adapted to test the solutions to determine the best solution which fits a predetermined requirement at least partially established from qualitative clinical data of a patient; and

- a selection module adapted to select the best solution to demonstrate the effect of the input functions on the dynamic state of the cell with respect to the breast cancer.

According to a third aspect of the invention there is provided a method of developing an interrogatable model of Burkitt's Lymphoma in order to determine characteristics of the Burkitt's Lymphoma which may aid in the diagnosis and treatment thereof; the method comprising the steps of: - defining a generic model of a cell having a conserved dynamic state, which model includes: o one or more molecular species, the or each of which undergo molecular interactions within the cell with one or more other molecular species; o a set of parameters associated with the cell; o one or more mathematical equations which substantially represent the dynamic state of the cell, the or each equations comprising:

^■ a one or more phase variables at least partially corresponding to the or each molecular species; ^■ parameters to allow variations thereto to be investigated; and o one or more input functions which can be used to generate the variations; - modifying the or each of the mathematical equations of the generic model to represent a one or more changes caused by the Burkitt's Lymphoma to the molecular interactions;

- solving the one or more equations to produce a plurality of solutions which represent the time-dependent concentrations of one or more selected molecular species based on one or more input functions to determine the effect of the input function on the dynamic state of the cell; - testing the solutions to determine the best solution which fits a predetermined requirement; and

- selecting the best solution to demonstrate the effect of the input functions on the dynamic state of the cell with respect to the Burkitt's Lymphoma.

According to a fourth aspect of the invention there is provided a system of developing an interrogatable model of a Burkitt's Lymphoma in order to determine characteristics of the Burkitt's Lymphoma which may aid in the diagnosis and treatment thereof; the system comprising: - a model generating module adapted to define a generic model of a cell having a conserved dynamic state, which model includes: o one or more molecular species, the or each of which undergo molecular interactions within the cell with one or more other molecular species; o a set of parameters associated with the cell; o one or more mathematical equations which substantially represent the dynamic state of the cell, the or each equations comprising:

^■ parameters to allow variations thereto to be investigated; and o one or more input functions which can be used to generated the variations; - a modification module adapted to modify the or each mathematical equation of the generic model to represent a one or more changes caused by the Burkitt's Lymphoma to the molecular interactions;

- a determination module adapted to determine an initial condition of the generic model by an iterative process until the initial conditions meet a predetermined performance requirement;

- a testing module adapted to test the solutions to determine the best solution which fits a predetermined requirement; and a selection module adapted to select the best solution to demonstrate the effect of the input functions on the dynamic state of the cell with respect to the Burkitt's Lymphoma.

In contrast with the simplifying assumptions and abstractions of prior art cancer models, the method of the preferred embodiment provides a mechanism for developing holistic models that enables more systems- level simulation than the prior art. Since the models developed by the method of the preferred embodiment are based on readily measurable variables, the models are inherently more realistic than those of the prior art. Furthermore, in building models based on qualitative and quantitative data, the method of the preferred embodiment is particularly well-suited to the biological scientific literature. Thus, in further contrast with the prior art models, the models developed by the method of the preferred embodiment are capable of making full use of the information contained in the literature.

In further contrast with the prior art, the method of the preferred embodiment develops interrogatable models from which it is possible to determine the interventions needed to achieve a particular desired state or outcome. More particularly, by establishing a mechanistic link between biopsy pathology and drug therapies, the models developed by the method of the preferred embodiment can be used to predict the viability of apoptotic pathways; and design optimal drug treatment regimes for a given pathology.

Whilst the following discussions are focussed on the application of the method of the preferred embodiment to modelling cancer, the skilled person will realise that the method of the preferred embodiment is not limited to this application. In particular, the method of the preferred embodiment is also applicable to modelling other diseases of a cell cycle.

Brief Description of the Drawings

An embodiment of the invention will now be described, by way of example only, with reference to the accompanying drawings in which:

Figure 1 is a flowchart providing an overview of the method of the preferred embodiment; Figure 2 is a block diagram of a network architecture in a generic model developed by the method of Figure 1 ;

Figure 3A is a schematic of a trophic factor pathway in a regulatory control pathways module of the architecture of Figure 2; Figure 3B is a schematic of an anti-growth pathway in the regulatory control pathways module of the architecture of Figure 2;

Figure 3C is a schematic of an apoptotic pathway in the regulatory control pathways module of the architecture of Figure 2;

Figure 3D is a schematic of the initiation of a growth cycle pathway comprising expression of an Myc transcription factor and activation of the Ras GTPase in the regulatory control pathways module of the architecture of Figure 2;

Figure 3E is a schematic of a RAS-α cascade pathway in the growth and division module of the architecture of Figure 2; FIGURE 3F is a schematic of the Cyclin D pathway in a growth and division cycle module of the architecture of Figure 2;

Figure 4 is a flowchart of a first algorithm for solving ordinary differential equations (ODEs) in the generic model of Figure 2;

Figure 5 is a flowchart of a second algorithm for solving ODEs in the generic model of Figure 2;

Figure 6A is a graph of the concentration profiles of a plurality of proteins in a healthy cell, simulated by an extension on the generic model of Figure 2;

Figure 6B is a graph of the concentration profiles of a plurality of proteins in a Burkitt's Lymphoma (BL) diseased but pre-malignant cell, simulated by an extension on the generic model of Figure 2;

Figure 6C is a graph of the concentration profiles of a plurality of proteins in a pre-malignant cell in which the BL has progressed further than in Figure 6B, wherein the concentration profiles are simulated by an extension on the generic model of Figure 2; Figure 6D is a graph of the concentration profiles of a plurality of proteins in a cell in which P73-V is over-expressed, wherein the concentration profiles are simulated by a further extension on the generic model of Figure 2; Figure 6E is a graph of the concentration profiles of a plurality of proteins in a cell in which P53 is lost, wherein the concentration profiles are simulated by a further extension on the generic model of Figure 2;

Figure 6F and 6G are graphs of the concentration profiles of a plurality of proteins in a healthy cell under the influence of mitogenic growth factor signalling at 90% and 80% respectively of that allowed for in the generic model of Figure 2;

Figure 6H is a graph of the concentration profile of a plurality of cyclins in a fully malignant cell, wherein the concentration profiles are simulated by an extension on the generic model of Figure 2; Figure 7 A to 7E are schematics of amendments to the pathways of

Figures 3A to Figure 3F to include breast cancer drugs;

Figure 8 is a schematic of interactions between drugs used in the treatment of ER- HER2+ breast cancer;

Figure 9a is a graph of the concentration profiles of a plurality of proteins in a healthy breast cell as simulated by the network of Figure 3A to 3F amended in accordance with Figures 7A to 7E and 8;

Figure 9b is a graph of the concentration profiles of a plurality of proteins in a cancerous breast cell with loss of PTen and gene amplification of Her2, as simulated by the network of Figure 3A to 3F amended in accordance with Figures 7A to 7E and 8;

Figure 9c is a graph of the concentration profiles of the proteins of Figure 9b in a breast cancer cell subject to the drug regime determined using the network of Figures 7A to 7E and 8;

Figure 10 is a schematic of interactions between drugs used in the treatment of ER+ HER2+ breast cancer; Figure 1 1 a is a graph of the concentration profiles of a plurality of proteins in a healthy breast cell as simulated by the network of Figure 3A to 3F amended in accordance with Figures 7 A to 7E and 10;

Figure 1 1 b is a graph of the concentration profiles of a plurality of proteins in a cancerous breast cell with the deletion of PTen and over- expression of Her2 simulated with the same network as that of Figure 1 1 a;

Figure 1 1 c is a graph of the concentration profiles of the proteins of Figure 1 1 b in a breast cancer cell subject to the drug regime determined using the network of Figures 7 A to 7E and 10; Figure 1 1 d is a graph of the concentration profiles of the proteins of

Figure 1 1 a in a healthy cell subjected to the drug regime of Figure 11 c;

Figure 12 is a schematic of interactions between drugs used in the treatment of ER+ HER2- breast cancer;

Figure 13a is a graph of the concentration profiles of a plurality of proteins in a healthy breast cell as simulated by the network of Figure 3A to 3F amended in accordance with Figures 7A to 7E and 12;

Figure 13b is a graph of the concentration profiles of a plurality of proteins in a cancerous breast cell with the deletion of PTen and over- expression of Her2 simulated with the same network as that of Figure 13a; Figure 13c is a graph of the concentration profiles of the proteins of a healthy breast cell of Figure 13a showing activation of the pro- translational proteins S6-P and elF-4G/A/B/E by mTOR-P;

Figure 13d is a graph of the target concentration profiles of the pro- translational proteins S6-P and elF-4G/A/B/E reliant on mTOR-P; Figure 13e is a graph of the concentration profiles of the proteins of

Figure 13b in a breast cancer cell subject to the drug regime determined using the network of Figures 7A to 7E and 12; and

Figure 14 is a block diagram of a system of developing an interrogtable model of breast cancer and Burkitt's lymphoma. Detailed Description

1. Introduction

The description will comprise a broad overview of the method of the preferred embodiment, followed by a more detailed analysis in which each of the steps of the method is discussed in turn. Two examples of applications of the method will then be described. The first example relates to the modelling of Burkitt's lymphoma and cellular regulatory aberrations caused thereby. The second example relates to the modelling of ER+Her2+ breast cancer. Beyond demonstrating how the models developed by the method of the preferred embodiment are extendable to embrace different forms and types of cancers, the second example also demonstrates how the models may be used and interrogated to design optimal drug treatment regimes for a cancer. As an aside, for brevity, the method of the preferred embodiment will henceforth be known as the 'modelling method.

2. Overview of the Modelling Method

Referring to Figure 1 , the modelling method comprises the steps of: - creating (10) a generic model based on the molecular biology of a cell of non-specific type; treating the cell as a conserved dynamical system and describing (12) its dynamics by reference to a one or more non-linear ordinary differential equations (ODEs) whose variables represent the concentrations of a one or more molecular species involved in the cell's critical regulatory processes; solving (14) the ODEs; and optionally implementing (16) a one or more calibration procedures to convert the generic model into a cell-specific form. The generic model must be of sufficient complexity to provide a comprehensive and thorough description of the molecular interactions, redundancies and feedback loops that occur in real biological systems. In the present example, two hundred and eighteen variables that represent the molecular pathways implicated in the most prevalent forms of cancer are incorporated in the generic model. This ensures that all the known key molecular pathways are included; thereby enabling analysis of all the major processes involved in cell cycle regulation. Nonetheless, it will be appreciated that the modelling method is not limited to the above- mentioned number of variables. On the contrary, the modelling method is adaptable to embrace as many variables as required to describe a given disease state.

In the present example, the generic model is autonomous, insofar as it represents a cell's behaviour solely by relevant, molecular concentrations.

However, the modelling method also embraces a non-autonomous model that allows for external control thereof. In particular, the non-autonomous model allows for:

- control and simulation of regulatory processes (through external input data relating to ligand concentration); simulation of externally and internally directed DNA damage; and

- identification of appropriate intervention processes (e.g. by assessing the impact of various drug regimes).

Since an initial condition of the generic model is unknown, the state of the generic model at any time is unknown. Thus, the values of the model parameters that control molecular concentration kinetics are also unknown. The modelling method solves this problem by providing a number of algorithms that link solutions to the dynamic system with parameter sets that evolve according to preset criteria chosen to (quantitatively or qualitatively) represent increasingly acceptable solutions.

3. Detailed Description of the Modelling Method

3.1 Creating a Generic Model

Referring to Figure 2, the generic model comprises a network whose architecture embraces a large number of factors relating to a cell and its regulatory control pathways. The network architecture also embraces molecular interactions and intervention procedures etc. More particularly, the network architecture 20 comprises a mitogenic signalling module 22 and an intervention procedures module 26, both of whose outputs are coupled to a cell regulatory control pathways module 24.

The mitogenic signalling module 22 defines an extracellular matrix. Accordingly, the mitogenic signalling module 22 comprises the concentrations of mitogenic signalling factors for survival, anti-growth, extrinsic apoptosis; and growth and division. In use, the survival, anti- growth and extrinsic apoptopic mitogenic signalling factors respectively activate : a plurality of survival pathways (including the phosphoinositide 3-kinase (PI-3K) and the serine/threonine AKT pathway); a plurality of anti-growth pathways (through the SMAD pathway); and an extrinsic apoptopic pathway (e.g. through the Myc and the Ras pathway) of the cell regulatory control pathways module 24.

The intervention procedures module 26 comprises variables representing external damage, internal damage and drug inputs to the modelled cell. These variables activate an intrinsic apoptotic pathway of the cell regulatory control pathways module 24 (through Ataxia Telangiectasia Mutated (ATM) protein, the checkpoint kinase ChK1 , the phosphatise Cdc25A and tumour suppressor protein P53).

The outputs from the cell regulatory control pathways module 24 are coupled to a transcription factors module 28, which, in use, calculates the concentrations of a variety of transcription factors such as TFX, Myc, TCF, P53 and E2F. The transcription factor concentrations are in turn coupled to a gene expression module 30, which, in use, calculates the concentrations of P15, P16, P21 , P27, SKP2, Cyclin dependent kinases Cdk4, Cdk2, Cdk1 , Cyclin D, Cyclin E, Cyclin A, Cyclin B, E2F, Hdm2 and P14.

The calculated concentrations of transcription factors (from the transcription factors module 28); Myc and Ras proteins (from the cell regulatory pathways module 24) and other proteins (from the gene expression module 30) are transmitted, in use, to a growth and division module 34. From this information, the growth and division module 34 calculates the values of variables representing the extent of Cdk4, Cdk2, Cdk1 , Cdc25A, Cdc25B/C, Cdc20/Cdh1 and Rb pathway activation. The calculated concentrations of the transcription factors and other proteins are also fed back to the cell regulatory pathways module 24.

Figures 3A to 3D illustrate some of the interactions between the mitogenic signalling module and the regulatory control pathways module of the generic model. More particularly, referring to Figures 3A and 3B, acquiescent cells are subject to a survival signal 40 and an anti-growth signal 42. Thus, in the generic model, these signals have a default state of ON. This situation is maintained by ligands in the extra-cellular matrix including trophic factors (TF) and anti-growth factors such as TGF-β. Referring to Figure 3A, trophic factors bind to trophic factor receptors (TFR). This causes the activation of phosphatidylinositol 3'-kinase (PI-3K), which phosphorylates phosphatidylinositol 4,5 biphosphate (PIP2) to phosphatidylinositol 3,4,5 triphosphate PIP3. PIP3 binds to 3-phosphoinositide-dependent protein kinase 1 (PDK1 ) and protein kinase B (PKB). The PIP3:PDK1 complex activates PKB, resulting in the dissociation of PIP3, so that a quantity of active PKB (PKB-a) is made available. This protein ensures cell survival by phosphorylating the cytosolic protein BAD, to hold in check the production of caspases.

Referring to Figure 3B, anti-growth factors (e.g. TGF-β) bind to R1 1 receptors, forming a complex (TGF-B:R1 1 ) which binds to and activates R1 receptors. The resulting trimer (TGF-B: R1 1 : R1 -a) phosphorylates SMAD3 to produce SMAD3-P, which binds with SMAD4 (with assistance of IMP-b) to form a complex (SMAD-3P:SMAD4: IMP-B) which dissociates leaving the dimer SMAD3-P:SMAD-4. The dimer binds with a nuclear transcription factor (TFX); and the resulting trimer initiates expression of the inhibitor P27 (not shown) to halt any progression through the cell's growth and division cycle.

Non-acquiescent cells can undergo programmed cell death or enter into a growth and division cycle. Referring to Figure 3C, programmed cell death (or apoptosis) can be initiated by a pathway involving an increase in the caspases.

Referring to Figure 3D, the growth and division cycle is activated when there is an increase in available growth factors in the extra-cellular matrix, sufficient to bind with and activate the receptor tyrosine kinases (RTKA and RTKB) to activate three separate pathways, namely: (1 ) additional phosphorylation of PIP2 resulting in increased availability of PKB (which influences gene expression of SKP2; phosphorylation of GSK-3β; and secondary phosphorylation of MYC);

(2) initiation of the expression of the MYC transcription factor; and

(3) binding of the GRB2:SOS dimer with active RTKA to form a trimer whose phosphorylation of the RAS-α protein results in an activation cascade (Figure 3E) which leads to activation of the SRF-a:TCF-a group of transcription factors; and gene expression of many of the protein species involved in the cell-cycle including Cyclin D and P21. Referring to Figure 3F, as Cyclin D (CycD) becomes available, the cell can progress to a G1 phase of a cell cycle, wherein Cyclin D binds with the Cyclin Dependent Kinase (Cdk4) resulting in the formation of the CyclinD:Cdk4 dimer (DK4). The extent of this process is tightly controlled by the inhibitors P15, P16, P21 and P27, which bind to DK4 restricting its activity.

3.2 Describing the Dynamics of the Generic Model

In developing the non-linear ODEs that describe the dynamics of the generic model, it is assumed that the cell behaves as a chemical reactor whose mass varies only under the influence of rate definable sources and sinks. The generic model may also include a growth function to take into account any growth within the system. The sources for the model comprise mRNA and protein species including extra-cellular ligands, whose rates of creation may be modelled, for example, as direct linear and exponential synthesis; or synthesis regulated by strictly increasing continuous sigmoidal functions controlled by other proteins. The sinks relate to degradation of mRNAs and proteins (which may be subject to influences from other related proteins) and are generally modelled as negative exponentials. In order to permit external control of the dynamic system, mitogenic signalling is modelled by a set of continuous sigmoidal functions. Intervention processes including radiation input, simulated internal DNA damage and drug input are similarly modelled. The modelling method is not limited to these intervention processes; and may comprise additional or alternative intervention processes.

The ODEs employed in the modelling method are constructed from a one or more molecular reactions, substantially each of whose rates may be established according to a relevant kinetic equation. The kinetic equations are based on generic expressions for:

^■ linear synthesis (V = k, where k is a constant);

^■ exponential synthesis (v = k_: Jk₂[dmRNA] wherein k-i and k₂ are constants and Jk₂[dmRNA] represents the accumulated level of relevant mRNA);

^■ sigmoidal synthesis (v=k.f(S).y_n, where f(S) is the value of a sigmoidal function and y_n is the level of a protein); and

^■ exponential degradation (V=-k.y_a.f(y_b) where k is a constant, y_a is the level of the protein being degraded and f(y_b) is a function of a relevant protein)

However, it will be appreciated that the modelling method is not limited to these generic expressions; and may include any number of additional or alternative generic expressions.

Take for example, the sigmoidal synthesis expression, the sigmoidal functions used therein employ variations on the equation

Y y = , This equation is strictly increasing and continuous on the open interval, ]-°°, +°°[. The shape and position of the sigmoidal function

can be controlled using wherein a, b and c

respectively define the location of a point of inflection on the y-axis; the sharpness of the sigmoidal switch; and the magnitude of the function on the y-axis. A variety of increasing or decreasing sigmoidal functions may be created using predefined parameters for switching location, sharpness and magnitude. These functions may be combined to form on-off switches. By using first derivatives of these functions in the generic model, external control is achievable by parameter input alone. Since a numerical integrator is used to solve the ODEs, the sigmoidal functions must be mathematically well-behaved and exist only within the pre-defined limits of ]0 c[. In the case of mitogenic signalling, x represents time and lies in the closed interval, [0, T], where T is the total period of integration. In the case of transcription, where strictly increasing functions are employed, x represents the quantity of a relevant transcription factor and lies in the semi-closed interval, [0, +°°[.

The above-mentioned generic expressions are supplemented by the Law of Mass Action (wherein rate is based on the quantities of protein available) and Michaelian kinetics (wherein rate varies hyperbolically according to the quantity of substrate).

In the present example, the parameters of the kinetic equations may be classified as: • multipliers relevant to the steady state of the system;

• magnitude of mitogenic input signal;

• sharpness of application/ withdrawal of mitogenic input signal;

• switch-on/off time for mitogenic input signal; • rate of production/degradation of mRNA;

• rate of synthesis/degradation of proteins;

• constant level of certain proteins;

• association/disassociation constants; • Vmax for Michaelis Menten phosphorylation/dephosphorylation;

• Michaelis constant for Michaelis Menten phosphorylation/ dephosphorylation;

• inhibition and efficiency constants; • threshold levels for transcription;

• sharpness of transcription switch-on; and

• magnitude factor applicable to transcription.

However, it will be appreciated that the modelling method is not limited to the above classifications. Instead, the parameters of the modelling method may be classified according to any relevant criteria.

Referring to the first of the above classes, phase variables can be considered to be members of a steady-state group if they share an associated non-zero reaction rate. If the steady state is achieved with parameter set {S}, the steady state will also be achieved with parameter set α{S} (where α is a multiplier). Thus, some variation can be achieved in the parameters restricted by the steady state requirements. This enables inter steady state group behaviour to be matched to observable data without disturbing any previously identified steady state condition. It should be noted that steady state occurs through association and disassociation processes between the members of the class and as a result of post-translational modification processes.

Many proteins are always present in a cell; and can be treated as constants in the generic model. Similarly, the parameters of the magnitude of mitogenic input signal class represent the level of transcription factor which would enable transcription at a level of, for example, 50% of maximum DNA saturation. As an aside, the parameters, association, disassociation, Vmax etc. should not be interpreted as having values constrained by their literal descriptions. Instead, these parameters permit characterisation of the dynamics of a modelled cell in line with experimentally observed data.

3.3 Solving the Ordinary Differential Equations of the Generic Model

The first step in solving the ODEs of the generic model is to establish its initial condition. The initial condition may be established in an iterative manner, wherein an estimate of the initial condition is formed and the acceptability of the dependent parameter subset checked by equating the autonomous parts of the modelled system to zero. The initial condition may then be revised until both it and the parameter subset are consistent with the timescale assumed in the model for completion of a cell cycle. In the present example, the timescale was chosen to encompass most of the cell population doubling times from available experimental data. However, it will be appreciated that other criteria may be used for assessing an initial condition and parameter subset.

In a non-cycling state, it is assumed that all molecular species involved in a cell's regulatory processes are constant. Thus, in establishing an initial condition, the survival and anti-growth mitogenic signals are set to an ON state and all other mitogenic and intervention signals set to an OFF state. With the non-autonomous part of the generic model fixed in this way, the expression, dx/dt = 0 may be solved for an assumed initial condition (xθ). If the ODEs are expressed as dy/dt = f (y, t, c) (wherein y, t and c are respectively a set of variables, time and a parameter set), the expression can be solved initially for c, where t is defined for the non-autonomous part of the generic model, and y is a set of assumed values. Two algorithms have been developed for performing this process.

Algorithm 1

Algorithm 1 involves linking a solver for the ODEs with a genetic algorithm (GA). The GA evolves progressively better values of c to yield a (quantitatively or qualitatively) pre-defined value of y. More particularly, the problem is to find a parameter set c which permits an acceptable solution to the dynamic system, dy/dt = f (y, t, c). Let Ce 9t^pxq be the population of chromosomes, wherein p and q respectively represent the number of chromosomes in the population and the number of genes (or parameters encoded) on a chromosome. Further, let ae 9t^pxl be a fitness vector of the population.

Referring to Figure 4, in a first step, ranges for the parameters are generated 40. In one approach a 'best guess' value may be generated for each parameter and a range defined above and below these best guess values. Alternatively, a range may be established using minimum and maximum values initially defined for each parameter. Chromosomes are initially seeded 42 by populating genes therein with values randomly (or otherwise) selected from the ranges of the associated parameters.

Since many of the ODEs of the generic model are extremely nonlinear, they tend to be 'stiff. In other words, the solutions of the ODEs are extremely unstable. Thus, error can dominate the calculations unless integration is performed with very short time steps. However, this would result in an unmanageably long computational time. The modelling method solves this problem by using a variable step integrator, which allows an integration time step to be reduced so that error is kept within predetermined limits, whilst detecting the high proportion of ODEs that are not integratable or take too long to be integrated and whose corresponding chromosomes must be discarded quickly.

The ODE expressed by each chromosome is solved 44 over a predefined integration period. The results thereof form a matrix Re 9t^πxm wherein: n-1 is the number of time steps used in the integration process; the first column comprises the cumulative time at each integration step; and the remaining (m-1 ) columns comprise the value of each phase variable at the corresponding time step.

The results are tested 46 to determine how well they represent the system to be modelled. A first fitness function tests the results by measuring the distance between initially defined target values for any number of chosen variables and those generated by solving a given chromosome. More particularly, define a target matrix TM e 9t^πxm , whose first column comprises cumulative time and remaining columns comprise the values of a chosen subset of phase variables at the respective times in the first column. If a matrix R₁ is formed from results matrix R using the same subset of phase variables as TM, then the fitness of a given result

P -Oy (TM(n,m)- R₁(n,m))² can described as Fitness = — , wherein P and Q

100 are variable scaling constants used to present the fitness results as a percentage. A second fitness function tests the results in terms of demonstrated oscillations (wherein target oscillations may be pre-defined in terms of amplitude and period) therein. More particularly, a set of maxima M_n for any phase variable (y_n) can be obtained by analysing a corresponding column r, of the results matrix R. If the number of oscillations sought over a timescale T is given by OT_n then the fitness for a

100 ^* χp(P -|θT_n - length(M_n)|) given result is given by Fitness = . This

formula can be extended to address the requirement for periodic oscillations of given amplitude (a).

A third fitness function tests the results in terms of the sequencing of events. More particularly, whilst scientific literature generally provides little quantitative protein information, it often provides qualitative descriptions (i.e. biochemical activity) of an order of events. Thus, the third fitness function tests the results in terms of the sequencing of maxima and non-zero minima against a pre-defined qualitative description. More particularly, define an ordered subset se 9l^lxπ whose elements are phase variables listed in order of occurrence of their first maxima over a pre-defined time. A second ordered subset, R₂ e 9l^lxπ may be derived from the results matrix R, wherein R₂ lists the corresponding phase variables in the time order of their first maxima. Thus the fitness function

100*fp -∑{((find(R₂ == s(n))-n))) can be expressed as Fitness = —

Q

The modelling method is not limited to the above fitness functions. In particular, other fitness functions can also be employed. Similarly, the fitness functions can also be used in combination. For example, to model the behaviour of a group of proteins where experimental data shows each protein achieving a maximal or minimal value once over a given time, the second fitness function would be applied first, followed by the third and the first.

Having evaluated the chromosomes, some of them are retained 48 for future participation in the GA. To this end, the modelling method employs three selection methods. The first method selects only the "fittest" chromosomes for retention. In the second method, the probability of a chromosome being retained is proportional to the relationship between its fitness and those of its competitors. The third method comprises an elitism-enabled augmentation of the second method. Here, a chosen number of the fittest chromosomes are retained regardless of their proportionality relationships.

Thereafter, a crossover and mutation process 50 may occur. Crossover involves combining two of the retained chromosomes, so that certain aspects of each chromosome are used to create a new chromosome (i.e. an Offspring'). Following the creation of an offspring, random mutations can be implemented, by selecting a frequency and defining a range of possible mutations. In this way, new chromosomes are created which retain many of the characteristics of the original retained chromosomes whilst comprising some new and potentially advantageous features. A revised population (of similar size to the original) is then created 52, by including a number of retained chromosomes together with a number of evolved chromosomes. The new population is used to restart the process of solving ODEs. After a number of repetitions of the solution/evolution process the initial population may be reset (e.g. by taking the fittest chromosome from the current population as a centroid for a new range for the parameters). After all the above processes are completed, the solutions are extracted from the surviving chromosomes.

Algorithm 1 may also include a tool for plotting a fitness landscape, wherein the tool produces a graphical representation of any two chosen parameters in terms of the resulting fitness value. The algorithm may also include a phase plane analysis tool, which produces a graphical representation of the trajectory of any two phase variables over time.

Algorithm 2

Algorithm 2 is a solution simulation algorithm. The algorithm includes an ODE solver, which is adapted to allow examination of the dynamic system under conditions of parameter change. Referring to Figure 5, Algorithm 2 comprises similar elements to those of Algorithm 1. However, Algorithm 2 also comprises a parameter slider mechanism 54, which enables the effects of changes in parameter values to be studied for any combination of selected parameters. The mechanism works by taking a pre-defined set of values for a given parameter and finding a solution to the system for each value. In simulations each solution is retained to visualise the effects of parameter change.

Algorithm 3

Algorithm 3 is a qualitative data optimizer algorithm. This algorithm is used to match a potential solution to the dynamic system with any available qualitative data. In accordance with most experiments carried out in molecular cell biology, the results are shown in the form of blots, wherein the size and intensity of a particular blot is taken as an indication of the amount of a respective molecular species. In this way, trends and relationships between particular proteins can be derived. The blots are not quantifiable, particularly when a blot is completely 'black'. However, for the purposes of calibrating a model, trends in intensity over time can be used through relationships between proteins, expressed in the timing of their respective maxima and relevant gradients.

Consider a target matrix TM₁ e 9l^mxπ whose first column comprises a cumulative time. Further, let each of the remaining columns of TM₁ comprise the sign (+1 , -1 or 0) of a rate of change of a protein at the times listed in the first column. Usually TM₁ will be derived from the results of an experiment using image analysis software. A matrix R₂ may be constructed from the results matrix R, wherein the first column of R₂ denotes the closest integration time steps to those of the first column of TM₁ and the remaining columns list the sign (+1 , -1 or 0) or the rate of change of the set of proteins in TM₁. A fitness function can be constructed as follows Fitness = P(lOO - ∑ ^TM₁ (n,m)- R₂ (n,m|))/ 100

Since it is important that models are stable and not overly sensitive to small parameter changes, an essential part of the modelling method is testing the stability and conditioning of models generated thereby. The dynamic systems that represent the complex networks involved in cell regulation are highly non-linear. However, the nature of equilibrium in these situations can be investigated through linear approximations. In a complete dynamic system, the relative conditioning of a particular solution can be examined by considering the relative changes in a parameter compared with the relative changes in the results.

3.4 Calibration of the Generic Model

The generic model is not specific to a particular cell type. To make it suitable for specific applications, a series of calibration procedures may be required. Each calibration procedure may involve adjustment of the network topology (to include new molecular species, reactions, increased complexity of reactions, etc.) or, more frequently, kinetic parameters/ reaction kinetics (to include diffusion processes and other spatial considerations). The calibration process is not limited to the adjustment of the above-mentioned model features. In particular, other features of the model could be adapted in the calibration process.

Application 1 : Use of Generic Model to Study Burkitt's Lymphoma

Burkitt's Lymphoma (BL) is an aggressive malignancy of highly proliferating partially differentiated B-cells. The formation of tumours occurs despite a high apoptotic rate which fails to compensate for incessant cell proliferation. BL typically displays three characteristic chromosomal translocations which render a Myc gene susceptible to point mutations. BL can be modelled by adding to the generic model a few:

^■ variables (e.g. mRNA:P73V, P73F, P73V, Puma and Bcl-6),

^■ reaction equations (e.g. for the opening of Bax and the synthesis of P53 with or without Bcl-16 etc.), ■ differential equations, initial conditions and parameters (e.g. rate of production of mRNA:P73V by P73F).

A healthy cell that is not influenced by the mitogenic stimulation of growth factors or apoptosis-promoting ligands is in a completely stable state. This state is maintained by the mitogenic influences of atrophic factors and anti-growth factors only. None of the regulatory proteins (apart from inhibitors, which are necessary to initiate and progress a cell cycle) are present. In contrast, a healthy cell subject to mitogenic stimulation by growth factors, progresses through a growth and division cycle. This process is simulated by the generic model, by setting the SURVIVAL (Figure 3A) and the GROWTH variable to an ON state and setting the ANTI-GROWTH variable to an OFF state. Referring to Figure 6A, the output from the generic model shows that the level and timing of the stable form of Myc (Myc-P), is highly regulated, which causes the level of an E2F transcription factor to be controlled. The E2F transcription factor is necessary for the synthesis of sufficient quantities of Cyclin E and Cyclin A for ordered entry into an S-Phase of a cell cycle.

A cell in a diseased but pre-malignant state has suffered at least one chromosomal translocation. This causes the cell to continually express Myc despite the absence of mitogenic stimulation by growth factors. The level of Myc in the diseased cell, is not higher than it would normally be in a healthy cell. However, in contrast with a healthy cell, the level of Myc in the diseased cell is maintained at that level. This results in a sustained level of E2F which causes a rise in Cyclin E:Cdk2 consistent with progression through the cell cycle into its late G1 phase. However, the low level of active Cyclin A:Cdk2 in the diseased cell is insufficient to drive the cell through the S-Phase. In particular, the cell cycle does not complete, since the sustained level of E2F causes activation of the intrinsic apoptotic pathway whose continued viability prevents a potentially malignant transformation.

The generic model simulates this process by setting both the SURVIVAL and ANTI_GROWTH variables to an ON state and setting the GROWTH variable to an OFF state. Continual expression of Myc is modelled by setting a rate equation variable v80 to a value of 37.5; which is the equivalent rate of Myc protein synthesis occurring in the generic model when the GROWTH signal is set at ON. Finally, an event function is provided which results in cell death when Caspase 3 reaches 16,000. Referring to Figure 6B, the results from the simulation showed that the maximum level of free E2F transcription factor is less than that of a healthy cell. However, in contrast with the peaked profile of E2F in a healthy cell (in Figure 6A), the level of E2F is sustained.

In a cell which has progressed to a further stage in the disease but is still in a pre-malignant state, Myc is not only continuously expressed but also stabilized as a result of a mutation in Exon 1 or as a result of a mutation at Threonine 58 in Exon 2. This condition is simulated using the signalling for the previous case and the revised rate equation for v80. Total Myc is modelled by amending the rate of degradation of unphosphorylated Myc to that applicable for Myc-P. Referring to Figure 6C, results of the simulation suggest that the over-expression of Myc has not caused a deterioration of the situation in terms of free E2F. Furthermore, the intrinsic apoptotic pathway is still able to initiate and progress programmed cell death. Although the total E2F available is of the order 27,000, the amount of free E2F is still low, confirming that inappropriate proliferation does not occur with the Rb pathway intact

It has been observed that in some cases of BL the intrinsic apoptotic pathway remained viable (in respect of the high apoptotic rates associated with the disease and evident from the histopathology). In contrast, the model shows that even where Myc is continuously expressed and stabilised, the intrinsic apoptotic pathway will still result in cell death thereby preventing malignancy. Thus, the model suggests that development of full malignancy in BL requires at least one further defect over and above chromosomal translocations affecting Myc and the stabilization thereof.

First, mutations in the gene encoding P73 are considered. P73 is a homologue of P53 and can act either as a positive or negative regulator of apoptotic cell death. Thus, a number of positive feedback loops and inhibitory controls need to be considered in the generic model. It is also necessary to slightly extend the generic model to encompass 'new' proteins. More particularly, the reactions and phase variables can be integrated into the generic model by extending the model to include appropriate rate equations, differential equations, initial conditions and parameters.

In order to investigate the effects of changes in P73 and P53, the simulation is repeated, with parameters relevant to P73-V (a splice variant version of P73) set at levels consistent with over-expression of P73-V. This situation is simulated by reducing the rate of degradation of mRNA:P73-V in the generic model by 80%. Referring to Figure 6D, the results from the simulation suggest that the intrinsic apoptotic pathway can be disabled. More particularly, P73-V is increased (from a base of 100 to approximately 260), which leads to a fall in P73F and P53. The reduction in both of these proteins effectively disables both arms of the intrinsic apoptotic pathway, (i.e. through P53 and through P73F, Puma and Bax). This removes any protection for the cell against inappropriate proliferation.

In 30% of BL cases, loss of heterozygosity (LOH) at q17 can result in the loss of P53, possibly causing the intrinsic apoptotic pathway to fail. However, the simulation results suggest that the pathway still operates, but apoptosis is regulated in a P53 independent manner through P73, Puma, Cyt-c Release and the Caspases cascade. Referring to Figure 6E, in particular, the intrinsic apoptotic pathway is still viable (even with a P53 knockout and provided there are no mutations in the P73 splice variants), thereby possibly still preventing malignant transformation. Finally the model is used to simulate the incidence of uncontrolled and inappropriate proliferation in fully malignant cells. Referring to Figures 6F and 6G, simulations of a healthy cell under the influence of growth factor, indicate that the growth factor present in the cell's microenvironment must be between 80% and 90% of that used in the generic model. This is sufficient to drive the cell through the growth and division cycle; and complete mitosis. However, simulation of a fully malignant cell (in which all of the chromosomal aberrations and mutations leading to the cancer are present and the mitogenic growth factor signal is set at 0) indicate that the cell is still capable of executing the cell cycle. Referring to Figure 6H, more particularly, the cell is still capable of producing sufficient Cyclin E and Cyclin A to complete S-Phase, mitosis and cell division. The cell has become proliferatively autonomous and will continue to proliferate in an uncontrolled fashion with the intrinsic apoptotic pathways disabled.

Application 2: Use of Generic Model to Study Breast Cancer

Breast cancer can occur when any of the mechanisms that stimulate proliferation in breast epithelia become disordered. In particular, over- expression of estrogen receptors can bring about inappropriate proliferation in response to low quantities of EGF and NRG. Thus, breast cancers are usually classified according to the apparent expression levels of estrogen nuclear receptors and Her2 growth factor membrane receptors. For example, an ER⁺Her2^" breast cancer is one caused by a fault in the endocrine pathways resulting in over-expression of estrogen receptors with normal expression of the growth factor receptor Her2. In addition, to breast cancers brought about through over expression of receptors, it should be noted that of the order of 10% to 20% of cases arise from hereditary causes usually due to mutations in certain genes. For example, this may occur where one copy of a gene is defective, perhaps PTen or BRCA1/2. When the second copy is lost through mutation then the relevant pathways may become defective and cancer may result.

Cancer diagnosis is usually confirmed by a pathology study of a biopsy sample of an apparently diseased tissue, wherein levels of estrogen receptors and Her2 membrane receptors are assessed by a scoring system which considers both the number of affected cells in a sample and the level of expression within individual cells. Samples are usually allocated a score of between 1 and 8, wherein a score of 1 is taken to be normal. Estrogen Receptor status (ER) (one of two breast cancer biomarkers used to aid prognosis and treatment) is also routinely measured by immunohistochemistry.

It is possible to extend the generic model to incorporate molecular pathways relevant to breast cancer and functions for modelling the molecular behaviour of various drugs. However, despite these extensions, two major challenges remain, namely:

• establishing a link between molecular pathology data, acquired at a single point in time, from a multi-cellular tissue sample (wherein the cells are all at different stages in their cell cycles) and the single cell dynamic models produced by the modelling method; and

• interpreting the results from the models produced by the modelling method in terms of a clinical outcome.

Referring to the first challenge, it is difficult to obtain the quantitative time-dependent data needed for direct calibration of the models produced by the modelling method. Thus, it is necessary to redefine the parameterisation of the models in terms of a given pathology. In particular, in a first step, a dynamic model dy/dt = f (y, c) (y being a set of phase variables and c being a parameter set) is defined for normal breast epithelial cells. Thus, the behaviour of a variable y over a given time

T interval T, can be expressed as E = Jydt . In a second step, an equivalent

0 dynamic system for a model (dY/dt = f (Y, C)) of a breast cancer cell for a

T T given pathology, is constructed such that Jydt = A J Ydt , where A is

0 0 defined as the ratio of the average intensity scoring on cancer proliferating cells to the average scoring on proliferating normal cells. The dynamic system is then solved for C by incorporating the above equation into a customized fitness function.

A further difficulty is that there may be no data available in terms of drug activity in relation to specific molecular targets. For example, a drug may be available which is designed and known to inhibit a particular protein, but there are no measurements which relate its dosage to its activity with respect to the target. Dosage is usually related firstly to toxicity tolerances and secondly to the success of the drug in cancer treatment, e.g. tumour shrinkage. This problem can be solved by calibrating the generic model with respect to drug regime and clinical outcome, after which the generic model can be developed to a point where its accuracy and precision of prediction becomes sufficient for direct application. In summary, the modelling method may be amended by: amending the generic model to incorporate: molecular pathways implicated in breast cancer; and ■ functions for modelling proposed drug inputs; - resetting the parameters of the generic model to describe the regulation of a 'normal' breast epithelial cell; using molecular pathology results to revise the parameters of the ODEs; redefining the parameters to describe any 'abnormal' pathology; identifying useful drugs; and adjusting Algorithm 1 of the modelling method to support variable drug parameters and thereby enable the design of a suitable drug regimen.

For brevity, the above-amended modelling method and the amended generic model will be henceforth referred to as the revised modelling method and revised generic model respectively. In a similar fashion, mathematical models produced by the revised modelling method, will be referred to henceforth as revised mathematical models. To be useful, the revised mathematical models must relate to a set of proteins expected to vary under drug intervention, wherein the proteins can be taken to represent a cell in both its healthy and cancerous state. For the sake of brevity, these proteins will be referred to henceforth as indicator proteins.

The network used in the revised generic model is detailed in Figures 7 A to 7E. In particular, comparing the network architecture of Figures 3A to 3F with that of Figures 7A and 7D, it can be seen that the network architecture of the original generic model has been amended to include:

- growth factor receptors Her1 , Her2 and Her3 shown separately (Figure 7A) (note: heregulin (HRG/NRG) binds to Her3 and epidermal growth factor (EGF) binds to Her1 , both Her3 and Her1 recruit Her2 as a co-activator and the resultant complexes activate the Pl-3k pathway and the Ras pathway respectively);

- estrogens binding to nuclear receptors ERa and ERβ while gestagens bind to PR (Figure 7B - the resultant complexes then form a self regulatory sub-network where Est:ERα is restricted by Est:ERβ and Gest:PR); - inclusion in the DNA damage response pathway of DNA repair proteins BRCA1/2 which are activated by ATM and phosphorylate and activate Chk1 (Figure 7C);

- extension of the AKT pathway to incorporate activation of the translation promoting proteins S6 and eif-4G/A/B/C (Figures 7D and 7E- note: S6-P is activated through PKB-a (usually referred to as AKT), m- TOR and P70. eif-4G/A/B/C is activated through m-TOR and eif- 4E:Phas1 :4E-AP1 -P. Also, in this revised network, PTen is treated as a variable).

Whilst a large number of drugs are used to treat breast cancer, the present example restricts itself to five, namely, Trastuzumab (Herceptin), Letrozole, Tamoxifen, RAD001 and Doxorubicin. Accordingly, referring to Figures 7A to 7E, the network architecture of the original generic model has been further amended to include:

- RAD001 activity (Figure 7D) enabled by binding with FKB12 which then binds with and inhibits mTor;

- tamoxifen binding with and restricting the activity of the nuclear receptors ERa and ERβ (Figure 7B); - letrozole restricting the level of estrogen production (Figure 7B);

- trastuzumab (also known as herceptin) binding (Figure 7A) with Her2 (the growth factor receptor co-activator) and thereby inhibiting the level of active growth factor receptors; and

- doxorubicin causing double strand breaks (DSBs) in DNA through its ability to disrupt Topol 1 and therefore activating ATM (Figure 7C).

The revised generic model takes a simple approach to modelling the effect of the drug inputs. More particularly, the model employs a single compartment paradigm, wherein absorption is modelled by a double sigmoidal function. This allows the input time and the absorption rate to be pre-set. More particularly, in the present example, the drugs are assumed to be absorbed over 30 minutes. Elimination is taken as a negative exponential with a specific half life. Each drug dose is permitted to fluctuate in the closed interval [0 10⁵ (maximum toxicity tolerance)] and the genetic algorithm combined with the ODE solver is run with all parameters fixed, apart from those related to drug input. The rate equations, differential equations, initial condition and the single parameter set, used throughout all of the simulations have been updated to reflect the amendments to the generic model and the modelling method.

In order to show how the amended modelling method and amended generic model can be used, three separate examples of breast cancer subclasses are considered. These are an EFTER⁺, an ER⁺ER⁺ and an ER⁺ER^".

Example 1 : ER^"Her2⁺ Breast Cancer

In an ERΗer2⁺ cancer anti-endocrine therapy would normally have little effect. This form of breast cancer is characterised by a number of pathology problems, namely:

- a major defect in the PI-3K pathway due to the loss of PTen; - significant over expression of Her2 which would indicate sustained activity in terms of the regulatory proteins involved in the cell cycle;

- a high proliferation rate which has resulted in tumour formation and may be due to the defects in PTen and Her2 only, although other factors may also be involved.

Accordingly, there is little point in administering an aromatase inhibitor such as Letrozole as the estrogens and gestagens are not over expressed. This is not the case with Tamoxifen since it is capable of inducing TGF-β and down regulating Her2 even though this drug is normally used in anti-endocrine therapy. Trastuzumab is certainly necessary to cope with the Her2 gene amplification and a further reason for using this drug is that a reduction in Her2 may result in a reduction in the active Her2:NRG:Her3:PI-3K complex which in turn will inhibit the conversion of PIP2 to PIP3 and hence reduce AKT. In the case of RAD001 there is certainly a good case for using an mTOR inhibitor. High AKT levels induced by the lack of PTen increase the pro-translation protein S6-P which is activated following a phosphorylation cascade involving mTOR-P and P70-P. AKT also activates the pro-translation complex elF-4G/A/B/E and the intervention of RAD001 should also counteract this. Doxorubicin would appear to be necessary to deal with the existing cancer cells. Since BRCA1/2 is normal and there is no reason from the pathology to suppose that P53 is affected. It may be useful to use an anthacycline such as Doxorubicin particularly since Her2 gene amplification would include over expression of Topol 1 a and thus an increase in sensitivity to this type of drug.

Thus, all the drugs available will be used in this example with the exception of Letrozole. Accordingly, the network of Figures 7A to 7E are adjusted to include the interactions between the drugs themselves as summarized in Figure 8, which indicates that certain equations will need to be rewritten to allow for the induction of TGF-β and inhibition of Her2.

The indicator proteins for the present example are mTOR-P, Her2, PIP3 and TGF-β. Referring to Figure 9a, using a growth factor stimulus of 10% and SURVIVAL and ANTI-GROWTH signals set to ON, the indicator proteins behave as expected in a healthy non-cancerous cell, with PIP3 increasing as the PI-3K pathway is activated and Her2 reducing as the receptor binds to the ligands. The cancerous state is simulated by cancelling PTen and changing the parameters for the Her2 expression so that Her2 is increased by a factor of 4. The result is shown in Figure 8c, wherein it can be seen that mTOR-P and Her2 are significantly increased, but TGF-β is unaltered.

Accordingly, the therapeutic aim is to determine a drug regime using Tamoxifen, Trastuzumab and RAD001 which will convert the behaviour of the indicator proteins from the cancerous state shown in Figure 9b to the normal healthy state shown in Figure 9a. Since the cell under consideration is epithelial, wherein TGF-β acts as an anti-growth factor which can indirectly induce P27 (and thereby prevent further proliferation), it would also be beneficial if TGF-β could be increased.

In the present example, the model is uncalibrated. Thus, the results produced using the parameter set will be of limited use and not sufficiently accurate to make the model clinically viable. However, with appropriate calibration this will be ameliorated. Furthermore, as previously discussed, the amended generic model relates to the dynamics of a single cell. However, drug treatment of cancer addresses the problem of diseased tissue comprising multiple cells at various stages of the cell cycle. This can present a problem of determining when drug input should be applied, but the model can be adapted to take this into consideration. This is especially the case with RAD001 where at certain stages the drug is ineffective. In the following simulations, all of the drugs are input early on in the cycle.

Figure 9c shows a possible solution which was obtained using an initial population of 128 parameter sets evolving over 10 generations. mTOR-P and Her2 are reduced to levels to be much closer to that of the healthy cell target. Also, TGF-β has been increased by some 50% which is beneficial. Less ideal is the fact that the drugs have been unable to reduce Pl P3. This can be corrected or improved by applying other parameter sets and running a further iteration to the model. The parameters of the algorithm can be changed but it is clear from Figure 9c that the results are converging and this is dependent on the application of the following drug regime.

Tamoxifen - 16954 units

Trastuzumab - 18278 units

RAD001 - 20108 units

Also, Doxorubicin would normally be used at a dose close to toxicity tolerance for the purpose of cancer cell elimination.

Example 2: ER⁺HER2⁺ Breast Cancer

In an ER⁺Her2⁺ cancer there are defects in the growth receptor pathways and also in the endocrine pathways. Four of the drugs will be used in the example. Letrozole and Tamoxifen will be used for anti-endocrine purposes with Trastuzamab to address the Her2 over expression and Doxorubicin to eliminate existing cancer cells. It is assumed that the intrinsic apoptotic pathways are viable. RAD001 is not used here since the AKT pathway appears to be intact.

The amendments to the network architecture to accommodate the drug interactions (and most particularly, letrozole whose plasma levels may be reduced by 35% - 40% by tamoxifen) are shown in Figure 10. In the present example, the indicator proteins are Her2, E2α and PR. The total amount of the Her2 is considered in this example. This includes free Her2 as well as the complexes Her2:NRG:Her3, Her2:NRG:Her3:PI-3K, Her2:EGF:Her1 and Her2:EGF:Her1 :GRB2:SOS.

Referring to Figure 1 1 a, the simulated behaviour of the indicator proteins in a healthy non-cancerous cell matches the behaviour that might be expected therefore. More particularly, the level of the pro-proliferative nuclear receptor ERa, both in the phosphorylated and unphosphorylated form is inhibited both by active ERβ and active PR. For comparison, a cancerous cell is simulated by the revised mathematical model, by

- setting the survival and anti-growth signals were set to an ON state;

- setting the mitogenic growth factor stimulation at 10%; - increasing Her2 expression by 400%;

- holding estrogen and gestagen signals at 1 (as in a normal cell);

- increasing ERa signal to a level necessary to generate a 300% increase in Est:ERα and Est:ERα-P;

- increasing the PR signal by an amount necessary to generate a 100% increase in the level of Gest:PR-a; and

- reducing ERβ, the other antagonist of ERa to zero.

Referring to Figure 1 1 b, although the concentration of active PR is increased, it is insufficient to restrict ERa and the concentrations of E2F (not shown) and Cdk2 (not shown) are significantly enhanced, which in a real situation, might activate the intrinsic apoptopic pathways.

As before, the object is to find a combination of drug concentrations which will alter the protein expression of the cancer cell (shown in Figure 11 b) to make it more similar to that of the healthy cell (shown in Figure

1 1 a). The GA is run with a population of 128 chromosomes over a single generation. A metric fitness function is used with a highest fitness selection. Using the above conditions, an optimal drug regime is determined to be:

Tamoxifen - 19,395 units

Letrozole - 43,306 units

Trastuzumab - 1 1 ,728 units Doxorubicin could also be added to this to deal with the elimination of cancer cells. The resulting concentrated profiles of the indicator proteins on treatment with this regime is shown in Figure 1 1 c. In particular, Est:ERα and Est:ERα-P is found to be very satisfactory: Her2 acceptable and, although Gest:PR-a is still high, this is probably unavoidable since it is an antagonist of ERa. The results (not shown) of a simulation of the effect of the drug regime on the concentration of cyclins in the cancer cell indicates that the cancer cell is now unable to proliferate.

The results of a simulation of the effect of the above drug regime on a healthy cell is shown in Figure 1 1 d. In this case, ERa and Her2 are clearly suppressed but are still available. It is likely in this instance that proliferation through the endocrine pathways would not be possible if the healthy cell were required to grow and divide. There may be a problem with threshold levels in this situation in that the drug level in the cancer cell would need to be minimised and still be able to restrain ERa whilst maintaining a high enough level to allow proliferation through the endocrine pathway in healthy cells.

Example 3: ERΗER2^" Breast Cancer

In the ERΗER2^" breast cancer, ERa is over expressed with a score of 6, PR has a score of 5 and the cancer cells are found to be P53 null. Since Her2 is not over expressed, Trastuzumab is not necessary here. However, both anti-endocrine drugs, Tamoxifen and Letrozole, will be required. Since the cell is P53 null it is likely that the intrinsic apoptotic pathway is disabled. Thus, drugs relying on the intrinsic apoptotic pathway (e.g. Doxorubicin) should probably be avoided. The measure for eliminating cancer cells will be by using RAD001 to reduce translation and rely on the immune system to eliminate the cancer cells through fast ligands and the extrinsic apoptotic pathways. Thus, in this example only three of the available drugs will be used, namely Tamoxifen, Letrozole and RAD001.

The amendments to the network architecture (of Figures 7 A to 7E) to accommodate the interactions between these drugs is shown in Figure 12. ERβ is assumed to be zero and thus certain equations need to be changed. For example, a rate equation needs to be changed to construct a link between mTOR-P, reduced translation and subsequent cell death. The indicator proteins are Est:ERα, Est:ERα-P, Gest:PR-a and mTOR-P. In a similar fashion to the previous examples, the algorithm is applied. The simulated behaviour of the indicator proteins in a healthy cell and a cancerous cell are shown in Figures 13a and 13b respectively.

Before running the search algorithms the problem of how much RAD001 to use is considered. Too high a dose will eliminate translation completely and cause high toxicity levels. This is a common problem in cancer drug therapy and here the problem is addressed by using the model to investigate the effect of mTOR-P on the pro-translational proteins S6-P and elF-4G/A/B/E. The situation in a healthy cell is shown in Figure 13c. Ensuring that translation is not totally eliminated but rather reduced to a lower level will result in a minimum quantity of drug being used and the model indicates that mTOR-P should be reduced to the level shown in Figure 13d to achieve this. As a result, the level of RAD001 will be a minimum amount consistent with restricting translation to a minimum but still positive level.

In setting up the model for the cancer cell shown in Figure 13b, the following procedure was adopted:

- the Survival and Anti-growth signals were set at ON; - the mitogenic growth factor stimulation was set at 10%; - the P53 transcription factor was removed;

- the Estrogen and Gestagen signals were held at 1 , as in the case of the normal cell;

- the ERa signal was increased to the level necessary to generate a 500% increase in Est:ERα and Est:ERα-P;

- the PR signal was increased by an amount necessary to generate a 400% increase in the level of Gest:PR-a; and

- ERβ was removed.

The search algorithm was then run using the values in Figure 13a for

Gest:PR-a, Est:ERα and Est:ERα, and Figure 13d for mTOR-P. The results of the search using a population of 128 over 1 generation are shown in Figure 13e where it is seen that the search was successful in all cases except for Gest:PR-a which is unchanged by the drug input. This is due to the fact that, although PR is positively regulated by ERa, it is also the case that PR inhibits ERa. However, since PR does not cause proliferation, in the model at least, then the high level of Gest:PR-a can be tolerated.

In the examples demonstrated here only those parameters related to the amount of drug in a system are allowed to vary while running an algorithm, e.g. algorithm 1. Since most of the parameters in the model are therefore fixed the rate of convergence to an acceptable solution is usually rapid.

If the result shown in Figure 13e is acceptable then the required drug regime as determined by Algorithm No. 1 is as follows:

Tamoxifen - 23,663 units Letrozole - 29,413 units RAD001 - 22,262 units The results can be tested by using the model to generate simulations of the cancer cell attempting to execute the cell cycle, while subject to the above drug input, and by looking at the effect of the drug regime on healthy cells. However, an important element of this example is in demonstrating that the model can be used to calculate the minimum drug dosage providing that the molecular targets of the respective drug are known.

Apparatus/System for Developing a Cancer Model

Referring to Figure 14, an apparatus or system 40 for developing a cancer model may comprise a model generating module 42 adapted to define a generic model of a cell having a conserved dynamic state. The model may comprise: - one or more molecular species, which undergo molecular interactions within the cell with one or more other molecular species;

- a set of parameters associated with the cell.

The model may also comprise a plurality of mathematical equations which substantially represent the dynamic state of the cell. The mathematical equations may, in turn, comprise:

- a one or more phase variables at least partially corresponding to the or each molecular species; and

- parameters to allow variations thereto to be investigated. Furthermore, the model may also comprise one or more input functions which can be used to generate the variations to the mathematical equations.

The apparatus or system 40 may also comprise a modification module 44 adapted to modify the mathematical equations of the generic model to represent a one or more changes caused by the breast cancer to the molecular interactions. Similarly, the apparatus or system 40 may also comprise a determination module 46 adapted to determine an initial condition of the generic model by an iterative process until the initial conditions meet a predetermined performance requirement.

Similarly, the apparatus or system 40 may also comprise a calculating module 48 adapted to solve the mathematical equations of the generic model to produce a plurality of solutions which represent the time- dependent concentrations of one or more selected molecular species. The calculating module 48 solves the mathematical equations based on one or more input functions and, in doing so, determines the effect of the input function on the dynamic state of the cell.

Furthermore, the apparatus or system 40 may also comprise a testing module 50 adapted to test the solutions (from the calculating module 48 to determine the best solution which fits a predetermined requirement at least partially established from qualitative clinical data of a patient. In addition, the apparatus or system 40 may also comprise a selection module 52 adapted to select the best solution (from the testing module 50) to demonstrate the effect of the input functions on the dynamic state of the cell with respect to the breast cancer.

Conclusion

Thus, the methodology of the preferred embodiment enables the construction of models which enable further sub-classification of breast cancer beyond the conventional ER:Her2 system. In particular, the underlying mutations along with the levels of ER and Her2 are used to establish the molecular concentration dynamics of the cancer cells so opening up a range of new sub-classifications which can now be addressed directly by targeting the specific molecular aberrations with targeted drug therapy as identified using the model. It has also been shown that if the molecular targets of a drug are known, it is possible to adjust the dosage to obtain a desired effect at the molecular level using the model. With addition calibration this would add significant optimisation of therapeutic dosage and for minimisation of dosage with regard to toxicity and collateral damage to healthy cells. This in turn would lead to more effective management of a particular disease state.

While particular embodiments of the present invention have been described herein for purposes of illustration, many modifications and changes will become apparent to those skilled in the art. Accordingly, the appended claims are intended to encompass all such modifications and changes as fall within the true spirit and scope of this invention.

Claims

1. A method of developing an interrogatable model of a breast cancer in order to determine characteristics of the breast cancer which may aid in the diagnosis and treatment thereof; the method comprising the steps of:

2. The method of claim 1 , wherein the step of defining the generic model comprises: determining a network for the cell which is a representation of all the molecular species and the molecular interactions there between to enable analysis of all major cellular life cycle processes.

3. The method of claim 1 or claim 2, further comprising defining the conserved dynamic state in terms of molecular concentrations of the one or more molecular species.

4. The method of claim 3, further comprising solving the equations such that the mass of the modelled one or more molecular species is conserved in the dynamic state.

5. The method as claimed in any one of the preceding claims wherein the mathematical equations comprise kinetic equations based on generic expressions for linear synthesis, exponential synthesis, sigmoidal synthesis and exponential degradation.

6. The method as claimed in any one of the preceding claims wherein the step of solving the one or more equations comprises the step of employing a genetic algorithm to solve the equations.

7. The method as claimed in any one of the preceding claims wherein the step of testing the solutions comprises the step of assessing the solutions in terms of a square of a Euclidean norm of a difference between a one or more phase variables of a given solution and a target matrix comprising desired time profiles of the phase variables.

8. The method as claimed in any one of claims 1 to 6 wherein the step of testing the solutions comprises the step of comparing a number of oscillations in the phase variables of a given solution with a desired number of oscillations thereof.

9. The method as claimed in any one of claims 1 to 6 wherein the step of testing the solutions comprises the step of comparing a one or more maxima and non-zero minima in the phase variables of a given solution against a pre-defined qualitative description of an order of events in the cell cycle.

10. The method as claimed in any one of claims 1 to 6 wherein the step of testing the solutions comprises the step of assessing trends and relationships between proteins as expressed in the size and intensity of experimental blot data and comparing the trends and relationships with those expressed in the phase variables of a given solution.

1 1. The method as claimed in claim 10 wherein the step of comparing the trends and relationships comprises the step of comparing the timing of the respective maxima and relevant gradients.

12. The method as claimed in any one of claims 1 to 6 wherein the step of testing the solutions comprises a combination of the steps of claims 8, 9, 10 and 1 1.

13. The method of any preceding claim, further comprising applying one or more calibration protocols to enable the generic model to be adapted to represent a specific cell type.

14. The method of any preceding claim, wherein the generic model is based on an architecture of influencing factors including: mitogenic signalling; cell regulatory control pathways; transcription factors; gene expressions; intervention procedures; or growth and division cycle information.

15. The method as claimed in claim 14 wherein the mitogenic signalling influencing factors comprise mitogenic signalling factors for survival, anti- growth, extrinsic apoptosis and growth and division.

16. The method as claimed in claim 14 or claim 15 wherein the intervention procedures comprise variables representing external damage, internal damage and drug inputs to the cell.

17. The method as claimed in claim 16 wherein the step of modifying a one or more mathematical equations of the generic model comprises the step of modifying the mathematical equations to represent the effect of the drug inputs on the molecular interactions within the cell.

18. The method of any preceding claim wherein the step of testing the solution comprises applying a fitness test to the solutions to find the solution which meets a predetermined level of fitness.

19. The method of any preceding claim, further comprising iterating variations and repeating the steps for solving the one or more equations; testing the solutions, and selecting the best solution in order to test the effects of different inputs or parameters on the disease.

20. A computer program comprising instruction for carrying out the method of any one of the preceding claims, when said computer program is executed on a programmable apparatus.

21. A system of developing an interrogatable model of a breast cancer in order to determine characteristics of the breast cancer which may aid in the diagnosis and treatment thereof; the system comprising:

- a model generating module adapted to define a generic model of a cell having a conserved dynamic state, which model includes: o one or more molecular species, the or each of which undergo molecular interactions within the cell with one or more other molecular species; o a set of parameters associated with the cell; o one or more mathematical equations which substantially represent the dynamic state of the cell, the or each equations comprising: - a one or more phase variables at least partially corresponding to the or each molecular species; ■ parameters to allow variations thereto to be investigated; and o one or more input functions which can be used to generated the variations;

22. A method of developing an interrogatable model of Burkitt's Lymphoma in order to determine characteristics of the Burkitt's Lymphoma which may aid in the diagnosis and treatment thereof; the method comprising the steps of: - defining a generic model of a cell having a conserved dynamic state, which model includes: o one or more molecular species, the or each of which undergo molecular interactions within the cell with one or more other molecular species; o a set of parameters associated with the cell; o one or more mathematical equations which substantially represent the dynamic state of the cell, the or each equations comprising:

23. The method of claim 22, wherein the step of defining the generic model comprises: determining a network for the cell which is a representation of all the molecular species and the molecular interactions there between to enable analysis of all major cellular life cycle processes.

24. The method of claim 22 or claim 23, further comprising defining the conserved dynamic state in terms of molecular concentrations of the one or more molecular species.

25. The method of claim 24, further comprising solving the equations such that the mass of the modelled one or more molecular species is conserved in the dynamic state.

26. The method as claimed in any one of claims 22 to 25 wherein the mathematical equations comprise kinetic equations based on generic expressions for linear synthesis, exponential synthesis, sigmoidal synthesis and exponential degradation.

27. The method as claimed in any one of claims 22 to 26 wherein the step of solving the one or more equations comprises the step of employing a genetic algorithm to solve the equations.

28. The method as claimed in any one of claims 22 to 27 wherein the step of testing the solutions comprises the step of assessing the solutions in terms of a square of a Euclidean norm of a difference between a one or more phase variables of a given solution and a target matrix comprising desired time profiles of the phase variables.

29. The method as claimed in any one of claims 22 to 27 wherein the step of testing the solutions comprises the step of comparing a number of oscillations in the phase variables of a given solution with a desired number of oscillations thereof.

30. The method as claimed in any one of claims 22 to 27 wherein the step of testing the solutions comprises the step of comparing a one or more maxima and non-zero minima in the phase variables of a given solution against a pre-defined qualitative description of an order of events in the cell cycle.

31. The method as claimed in any one of claims 22 to 27 wherein the step of testing the solutions comprises the step of assessing trends and relationships between proteins as expressed in the size and intensity of experimental blot data and comparing the trends and relationships with those expressed in the phase variables of a given solution.

32. The method as claimed in claim 31 wherein the step of comparing the trends and relationships comprises the step of comparing the timing of the respective maxima and relevant gradients.

33. The method as claimed in any one of claims 22 to 27 wherein the step of testing the solutions comprises the steps of claims 29 to 32.

34. The method of any one of claims 22 to 33, further comprising applying one or more calibration protocols to enable the generic model to be adapted to represent a specific cell type.

35. The method of any one of claims 22 to 34, wherein the generic model is based on an architecture of influencing factors including: mitogenic signalling; cell regulatory control pathways; transcription factors; gene expressions; intervention procedures; or growth and division cycle information.

36. The method as claimed in claim 35 wherein the mitogenic signalling influencing factors comprise mitogenic signalling factors for survival, anti- growth, extrinsic apoptosis and growth and division.

37. The method as claimed in claim 35 or claim 36 wherein the intervention procedures comprise variables representing external damage, internal damage and drug inputs to the cell.

38. The method of any one of claims 22 to 37 wherein the step of testing the solution comprises applying a fitness test to the solutions to find the solution which meets a predetermined level of fitness.

39. The method of any one of claims 22 to 38, further comprising iterating variations and repeating the steps for solving the one or more equations; testing the solutions, and selecting the best solution in order to test the effects of different inputs or parameters on the disease.

40. A computer program comprising instruction for carrying out the method of any one of the preceding claims, when said computer program is executed on a programmable apparatus.

41. A system of developing an interrogatable model of a Burkitt's

Lymphoma in order to determine characteristics of the Burkitt's Lymphoma which may aid in the diagnosis and treatment thereof; the system comprising:

^■ a one or more phase variables at least partially corresponding to the or each molecular species; ^■ parameters to allow variations thereto to be investigated; and o one or more input functions which can be used to generated the variations; - a modification module adapted to modify the or each mathematical equation of the generic model to represent a one or more changes caused by the Burkitt's Lymphoma to the molecular interactions;

- a calculating module adapted to solve the one or more equations to produce a plurality of solutions which represent the time-dependent concentrations of one or more selected molecular species based on one or more input functions to determine the effect of the input function on the dynamic state of the cell ;

- a testing module adapted to test the solutions to determine the best solution which fits a predetermined requirement; and

- a selection module adapted to select the best solution to demonstrate the effect of the input functions on the dynamic state of the cell with respect to the Burkitt's Lymphoma.