## Abstract

The biology of cancer is a complex interplay of many underlying processes, taking place at different scales both in space and time. A variety of theoretical models have been developed, which enable one to study certain components of the cancerous growth process. However, most previous approaches only focus on specific aspects of tumour development, largely ignoring the influence of the evolving tumour environment. In this paper, we present an integrative framework to simulate tumour growth, including those model components that are considered to be of major importance. We start by addressing issues at the tissue level, where the phenomena are modelled as continuum partial differential equations. We extend this model with relevant components at the cellular or even sub-cellular level in a vertical fashion. We present an implementation of this framework, covering the major processes and treat the mechanical deformation due to growth, the biochemical response to hypoxia, blood flow, oxygenation and the explicit development of a vascular system in a coupled way. The results demonstrate the feasibility of the approach and its applicability to *in silico* studies of the influence of different treatment strategies (like the usage of novel anti-cancer drugs) for more effective therapy design.

## 1. Introduction

Cancer is one of the leading causes of death in most Western societies, responsible for nearly 30% of deaths in Switzerland (BFS 2007). The progression of cancer involves several stages, defined by changes in the behaviour of the tumour cells. Most tumours are caused by abnormalities in the genetic material of the transformed cells. Cancer progression seems to involve the accumulation of such abnormalities, making the tumour cells more invasive (Hanahan & Weinberg 2000). Initially, the transformed cells proliferate and form a cluster that gradually increases in size. In this early phase, tumours growing *in vivo* tend to develop without a dedicated vascular network, relying on diffusion from the neighbouring healthy vascularized tissue for the supply of oxygen and nutrients and the removal of wastes (e.g. CO_{2}). Diffusion from the boundary alone, however, does not suffice to support larger tumours, which often develop hypoxic and eventually necrotic regions. In response to acute hypoxia, the developing neoplasm produces a cascade of specialized agents, collectively called tumour angiogenesis growth factors (TAFs). These biochemicals promote the proliferation and migration of endothelial cells, thereby creating a specialized vascular network that is able to supply the tumour. This formation of a vascular system is called tumour-induced angiogenesis (Carmeliet & Jain 2000), enabling extensive growth that leads to increasing pressure and traction on the tumour's micro-environment (Gordon *et al*. 2003). There is growing evidence of the influence of pressure on the transport of nutrients (and drugs), tissue composition and proliferation. In this paper, we restrict our attention to solid tumours, although some concepts could also be applied to invasive ones. A solid tumour is an abnormal mass of tissue that usually does not contain cysts or liquid areas. They are named for the type of cells that form them (e.g. sarcomas, carcinomas and lymphomas) and may be benign or malignant.

The complexity and quantity of biological knowledge, combined with rapidly increasing computational power and the vision of more efficient therapy and drug design, has generated an increasing interest in the mathematical and computational models of tumour growth, angiogenesis, tissue mechanics, blood flow and other relevant fields in the past three decades. At the different spatial scales, various types of models have been developed, including continuous, discrete, deterministic, stochastic, analytical and numerical methods. A review of all the approaches is beyond the scope of this paper and exhaustive surveys can be found in Araujo & McElwain (2004), Byrne *et al*. (2006) and the very recent paper by Bellomo *et al*. (2008). Additional literature on specific aspects relevant to our model was presented in our previous works (Szczerba & Székely 2002, 2005; Lloyd *et al*. 2007) and will not be repeated here.

One of the main limitations of many current modelling approaches is that they neglect the important interactions between the individual underlying processes, covering a broad spectrum of spatial and temporal scales. As documented in Bellomo *et al*. (2008), there is a long list of different models, operating at the characteristic (natural) scale of the phenomenon they represent. If we are to build models that capture the interdependence between different phenomena, it is inevitable that we develop methods to integrate the corresponding scales into one consistent model. One of the first tumour models incorporating multiple scales in a systematic way was presented by Alarcón *et al*. (2004, 2005). It included three modules: oxygen and TAF diffusion; cell division dynamics; and blood flow. The cell cycle dynamics is represented by a cellular automaton, which shares the same domain as a pre-existing hexagonal vascular network. They do not include the process of sprouting angiogenesis, instead replacing it by remodelling due to (flow induced) shear stress and TAF, based on the work by Pries *et al*. (1998). More recently, they have extended their approach to include cell crowding (Betteridge *et al*. 2006). The model described by Frieboes *et al*. (2006) treats the cell dynamics and the transport of oxygen and growth factors using convection–diffusion–reaction equations solved using sophisticated numerical methods (Macklin & Lowengrub 2007). Cell-to-cell adhesive forces are modelled by a surface tension, treating the tissue as an incompressible fluid. Their recent work (Sanga *et al*. 2007) includes angiogenesis using a model based on the approach by Plank & Sleeman (2004). The recently published paper by Kim *et al*. (2007) presents a hybrid model of avascular tumour growth, which treats the actively proliferating rim using a discrete approach, while using a continuum viscoelastic material response for the remaining regions. This reduces the computational burden compared with a purely discrete approach.

The framework that we present here allows one to consider multiple phenomena in a coupled way by interfacing components both horizontally (component wise) and vertically (scale wise). A first attempt for computational implementation is presented in Lloyd *et al*. (2007). Here, we continue by including a detailed cellular-level simulation of oxygen transport to the tissue. This is an important model component also with respect to the transport of drugs (Jang *et al*. 2003). Additionally, we have added several mechanisms of vessel remodelling due to shear stress, including blood flow resistance remodelling and vessel regression.

Horizontally, any finite-element-based module can be added to the framework by giving access to the global underlying information like the domain representation (geometry, topology and constraints) and expecting its contribution to a global coefficient matrix, or performing a requested number of iterations in case no global assembly is required or the problem is nonlinear. The dialogue between the different components does not need to be restricted to the same resolution. For example, the oxygen source in a tissue-level reaction–diffusion model can be estimated from a cellular scale blood flow simulation. The interaction is realized via constitutive laws, e.g. effective responses to loads, tractions or constraints modulated by other components of the model. We have currently implemented and tested structural mechanics, chemical transport/reaction and blood flow modules.

Technological innovation has provided biological and clinical research with new insights and large amounts of data to study tumour growth. It has, however, very quickly become clear that cancer is a complex multifactorial process, extremely variable between different types of cancer cells, host tissues or progression stages. In order to be predictive, we need to understand not only individual mechanisms but also the interaction between physiological processes at many different spatial and temporal scales. This calls for an integrative approach, analogous to those proposed by the Europhysiome initiative (STEP Consortium 2007), which intends to establish a better coordination between various European Physiome projects and develop the so-called virtual physiological human. This methodological and technological framework aims to enable collaborative investigation of the human body as a single complex system.

There are a few core areas of interest in this initiative, for example, addressing cardiac physiology and mechanics, in order to understand the fundamental aspects of treating various heart malfunctions (Hunter & Nielsen 2005). Another example is the investigation and multifunctional modelling of pathological tissue growth for enabling the development of more effective tumour therapy strategies. The Physiome project addressing tumour growth requires the use of a broad range of theoretical and numerical tools, including stochastic processes, chemical kinetics and transport, finite-element modelling of fluid flow, structural mechanics and many more. In addition, these tools have to be applied over a range of spatial and temporal scales allowing one to capture the structure–function interactions on all relevant levels (sub-cellular, cellular, tissue and macroscopic). The benefits of such a model-based understanding of physiological mechanisms involved in tumour development can hardly be overestimated for both clinical care and basic biological research. In the context of the above-outlined Physiome project paradigm, our objective is to investigate various aspects of tumour growth, neovascularization and remodelling in one coherent modelling framework. A comprehensive package will allow the simulation of the underlying physical and biochemical processes, covering mutual relationships between mechanical and biochemical stimuli, transport phenomena like flow, diffusive transport or active cell migration, or possibly even gene expression, resulting in tissue formation, deformation or removal. Such an *in silico* model can greatly contribute to the development of novel treatment strategies such as optimizing dose delivery during radiation therapy or systematic selection of anti-angiogenic drugs depending on tumour phenotype.

## 2. Description of framework

Our model combines several previously identified components in a consistent framework, including neoplastic tissue growth, blood and oxygen transport and angiogenic sprouting. We use a hybrid approach combining continuum scale reaction–diffusion equations and microscale transport for addressing the oxygen delivery problem. Sprouting angiogenesis is treated by explicitly modelling sprouting vessel tips, proliferating and migrating up the gradient of angiogenic growth factors that are produced in hypoxic regions of the tumour. Finally, a continuum model of tissue behaviour and growth is added to the multi-physics, multi-scale simulation. The equations describing individual model components are solved in each time step, using finite-element methods on body-fitting meshes.

The time scales at which the diffusion of molecules takes place (nanoseconds to seconds) and at which the cells proliferate (days) allow us to treat the processes in a quasi-static manner, i.e. at a given time step, we compute a concentration map of oxygen according to the sources and consumption by the tissue based on the actual tumour and vascular geometry. In the following, we give a more detailed description of the individual model components.

### (a) Tissue growth

The growth model was presented in our previous paper (Lloyd *et al*. 2007) and will only briefly be described here in order to keep this paper self-contained. Tumour development is assumed to be relatively slow; therefore, it is reasonable to surmise that over time the neighbouring tissue will accommodate for the stretching by increased growth or cell migration. The volumetric growth due to cell proliferation is prescribed by an initial strain condition *ϵ*_{0} (Zienkiewicz & Taylor 2000). We compute the new deformation caused by a volume expansion using the finite-element method. In our implementation, we neglect the accumulation of stress during growth, instead assuming that the residual stress caused by inhomogeneous volume increase dissipates completely, constituting an elastoplastic growth law. A discussion of alternative biomechanical tissue material laws is given in §4*a*. From biological research, it is known that cells need certain conditions in order to undergo mitosis. One mechanism, which regulates proliferation through changes in cell cycle dynamics, is the dependency on oxygen and glucose (Jiang *et al*. 2005). At low concentrations, oxygen consumption and cell proliferation are reduced according to Freyer & Sutherland (1985) who explained this observation by the increased number of cells in a quiescent state. Under acute hypoxia, cells will eventually die, leading to necrotic regions in the tumour. Based on these observations, we have defined a functional relationship between the prescribed strain (i.e. tissue growth) and the availability of oxygen(2.1)where *f* is selected as a piecewise linear function. For high oxygen concentrations (>*h*_{1}), the cells reach the maximum proliferative potential, which decreases at lower levels. In necrotic regions with oxygen concentrations below a threshold *h*_{2}, we allow for the removal of tissue debris by prescribing a negative strain leading to volume loss.

### (b) Oxygen transport

Oxygen is transported by the blood and diffuses through the vessel walls and into the tissue. At a macroscopic scale, this process can be modelled by a reaction–diffusion equation(2.2)where *c* is the oxygen concentration and is the diffusion constant. The oxygen consumption rate depends on the tissue type and can vary depending on the metabolic activity. The source term in healthy tissue usually matches the consumption term. Its magnitude is governed by the blood flow inside the vascular system and the local tissue oxygen content, which influences how much oxygen diffuses through the vessel wall.

For computational reasons, and since we are mainly interested in the development of the tumour and not the healthy tissue, we choose to use a dual representation of the oxygen source. Inside the tumour, an explicitly modelled vasculature supplies the tissue with oxygen, while in the healthy tissue, we postulate a homogeneous source of oxygen from an implicit pre-existing regular vascular system. In the implicit representation, oxygen diffuses into the tumour from the boundary. Under the assumption that the oxygen levels in the healthy tissue are not affected by the presence of the tumour, we can model this as a reaction–diffusion problem (2.2) with constant boundary conditions at the tumour surface. The computation of oxygen transported by the blood vessels and diffusing into the tissue is quite involved and calls for a special treatment discussed in §2*c*.

### (c) Haemodynamic conditions in the blood vessel system

In order to estimate the oxygen source term , several issues need to be addressed. For computing the oxygen transport and partial pressure distribution in the vessel system, we first need to solve the blood flow problem. The partial pressure difference is a driving force behind oxygen delivery through the vessel wall; therefore, the local tissue oxygen concentration and the corresponding partial pressure should be known.

The vessel network generated by our angiogenesis simulation consists of straight vessel segments (pipes) connected in a network. Following our previous work on tumour-induced angiogenesis, we compute the pressure distribution and flow in the network (Szczerba & Székely 2002, 2005). For laminar steady flow in stiff and straight tubes, the flow is computed according to Poiseuille's law(2.3)where *Q*_{ij}, *R*_{ij}, *L*_{ij} and *μ*_{ij} are the flow, the radius, the length and the viscosity, respectively, between nodes *i* and *j*. The variables *P*_{i} and *P*_{j} are the pressure values at the nodes and *G*_{ij} is the conductance. The flow fulfils mass balance at the junctions between pipe segments, i.e. for all nodes *i*. Mass balance and Poiseuille's law lead to a system of linear equations, with pressure as an unknown quantity. For large vessels, the blood viscosity *μ*_{ij} is approximately constant, while for small diameters *D*, the Fahraeus–Lindqvist effect causes non-Newtonian behaviour and has to be taken into account. Therefore, in order to correct for small vessels, we use an apparent viscosity, empirically estimated under *in vivo* conditions by Pries *et al*. (1996).

Usually, the vascular structures we generate consist of many hundreds of thousands of vessel segments. For such large systems, it is necessary to take advantage of the sparsity of the system in order to cope with the excessive memory requirements.

### (d) Oxygen supply by the vascular system

As blood flows through the vessels, oxygen diffuses through the vessel walls into the tissue. Using a mass balance principle for the blood oxygen content and Fick's first law, we can compute the amount of oxygen which is delivered to the tissue(2.4)where *Q* is the blood flow; is the oxygen concentration (ml O_{2} per ml blood); d*S* is the surface of the vessel between *x* and *x*+d*x*; *w* is the wall thickness; and is the diffusion constant of oxygen in the vessel wall. The position *x* is upstream of *x*+d*x*.

Oxygen concentration in the tissue is proportional to the oxygen partial pressure according to Henri's law(2.5)This linear relationship is only valid for dilute solutions as in the case of oxygen dissolved in the tissue interstitial fluid. The parameters and are called Henri's constant and solubility, respectively. Assuming that Henri's law is valid inside the vessel wall, the gradient term in equation (2.4) can be replaced and approximated by a difference quotient(2.6)where is the pressure on the interface between the vessel and the tissue. The partial pressure on the surface of the vessel lumen is given by . The product is also known as Krogh's oxygen diffusion coefficient (Ursino *et al*. 1989). Basically, oxygen is removed downstream, eventually leading to a small difference in the oxygen partial pressure between the tissue and the blood, i.e. at low blood , very little oxygen will be delivered to the tissue. Additionally, we can see that the blood gradient along the vessel will be smaller if the flow is higher.

Because only a small proportion of blood oxygen content is dissolved and the major part is bound to haemoglobin, the shape of the – curve is dominated by the oxygen–haemoglobin dissociation curve. Its sigmoidal shape results from the influence of bound oxygen molecules on the conformation of the haemoglobin, which facilitates the binding of additional oxygen molecules. In general, the oxygen content in the blood can be given by(2.7)where is the solubility of oxygen in blood and is the oxygen saturation of the haemoglobin, i.e. the percentage of haem sites that are combined with oxygen. We use Hill's equation for the saturation function (Ji *et al*. 2006). The gradient of the function depends on the Hill exponent *n*_{Hill}. The value at which haemoglobin is saturated by 50% is given by *P*_{50}. In tissue regions with high carbon dioxide concentrations or high pH, *P*_{50} increases and the saturation decreases, i.e. more oxygen is passed to the tissue. Currently, however, we ignore this effect. Each gram of haemoglobin can combine with 1.36 ml of oxygen. With a normal haemoglobin concentration of 150 g l^{−1} (Ward 2006), and neglecting the dissolved oxygen, fully saturated blood has an oxygen content of *β*=1.36 ml O_{2} g^{−1}×0.150 g ml^{−1}=0.204 ml O_{2} ml^{−1}.

Using the same nomenclature as for the blood flow calculation, we can rewrite (2.6) as(2.8)Because, according to equation (2.7), the blood partial pressure depends on the oxygen concentration in a nonlinear fashion, the above relation is a nonlinear function in *c*_{i} and *c*_{j}. If *c*_{i} is known, the unknown downstream value *c*_{j} can be calculated using a Newton–Raphson scheme. In the case where more than one input vessel segment joins in a point, we make the assumption that the blood is perfectly mixed in the junction and all downstream pipes have equal blood oxygen concentrations.

Finally, if the oxygen concentration in the arterial inflow pipes is known, we can compute the distribution in the entire vascular system by propagating through the network in the direction of the blood flow. By integrating over the surface of the vessels, we can estimate the local amount of oxygen delivered by the explicit vessel network. This result is used in the macroscopic-level reaction–diffusion equation (2.2), which we solve using a Bubnov–Galerkin finite-element method (Huebner *et al*. 2001). Note that the actual amount of oxygen flux across the vessel walls depends on the oxygen concentrations in the tissue via . Therefore, the macroscopic and microscopic models are coupled and must be solved iteratively.

### (e) Sprouting angiogenesis

The formation of blood vessels by sprouting angiogenesis is a process where capillary sprouts depart from the pre-existing parent vessels in response to externally supplied chemical stimuli. By means of endothelial cell proliferation, migration and remodelling, the sprouts then organize themselves into a branched connected network structure. The preliminary capillary plexus does not necessarily require blood flow to form.

It is widely agreed on that a major cause of elevated angiogenesis growth factors (TAFs) is hypoxia. Hypoxia upregulates many genes, but the induction of vascular endothelial growth factors (VEGFs) is perhaps the most remarkable, up to 30-fold within minutes (Carmeliet 2003). As the growth factors diffuse into the surrounding tissue, there is some uptake and binding by the cells (Ausprunk & Folkman 1977). Therefore, the process can be modelled by a reaction–diffusion equation with a natural decay term(2.9)where *c*_{TAF} is the TAF concentration; *D*_{TAF} is the diffusion coefficient; and *Θ*_{0} is the decay rate. The source term depends on the local tissue oxygen concentration. If the oxygen tension is below a threshold, growth factors are produced by the affected cells. We model this by a function , where *H*(.) is the Heaviside function; *S*_{0} is a constant production rate; and *h*_{2} is the oxygen concentration threshold, below which the tissue becomes hypoxic. Currently, we treat the multitude of angiogenesis growth factors as a single generic growth stimulus. For the TAF diffusion constant, we use the estimated value for VEGF (*D*_{VEGF}=1.0×10^{−6} cm^{2} s^{−1}; Gabhann & Popel 2007).

Once the initial distribution of the growth agents due to secretion by the hypoxic cells has been established, endothelial cells respond to the stimulus by sprouting. In our previous work, we developed a model of sprouting angiogenesis, which generates functional flow networks, consisting of interconnected straight pipes (Szczerba & Székely 2002, 2005).

Similar to our previous method, we model the dynamics of the sprouting vessel tips using a Lagrangian approach. The major innovations include the coupling of the vascularization to the tissue growth via the time-dependent changes in the growth factor distribution, which typically decreases with increasing vascularization, since the initially hypoxic region receives a new source of oxygen. Additionally, we allow for the capability of the vascular system to adapt to current haemodynamic conditions by remodelling the vessels' diameters in order to decrease the flow resistance. The details of this process are presented in §2*f*.

### (f) Vascular remodelling

Once vessels have sprouted and penetrated the tissue, they undergo various changes, usually described as remodelling. Regulatory processes, which are believed to take place, include modification of the radii of vessels under high shear stress in order to reduce the flow resistance, adaptation of the vessel wall thickness due to pressure and general adaptations under the influence of vasoactive factors or metabolic stimuli. We propose a set of rules that take various biophysical and haemodynamic quantities such as flow, shear stress, growth factors and vessel maturation into account. Within this context, it is reasonable to apply a hybrid modelling approach, where physical quantities are computed and used as input to a rule-based simulation of the investigated processes. A rule is applied according to a probability, which depends on the physical quantity. This is fundamentally similar to the approach by Gödde & Kurz (2001) who proposed two probability functions for adding and deleting a vessel segment in a hexagonal grid. The following rules have been implemented.

Endothelial cell recruitment (migration and proliferation) depends on the shear stress

*τ*at the vessel wall and a local concentration of TAF. Shear stress in pipe flow is given as(2.10)with the same variable definitions as in §2*e*. Shear stress varies remarkably little within large ranges of vessel sizes and on average is approximately 1–3 Pa (Kamiya & Togawa 1980; Miyashiro*et al*. 1997).High shear stress

*τ*(above the optimal shear stress*τ*_{opt}) generally will call for the recruitment of endothelial cells in order to increase the vessel diameter and reduce the flow resistance.When the flow drops below a threshold and

*c*_{TAF}is low, vessel retraction occurs. Otherwise, if*c*_{TAF}is high it will allow the vessel to persist.When the tissue receives a sufficient amount of oxygen, no more TAF is produced and the remodelling eventually ceases.

Maturation of vessels will gradually reduce the responsiveness of a vessel to undergo remodelling.

These remodelling guidelines based on different physical and chemical stimuli and their interaction are depicted in table 1. Note that table 1 has been implemented in a probabilistic way using heuristic functions as described previously. Currently, we ignore the effect of pressure on the vessel wall thickness and instead use an empirical functional relationship between radius and thickness like in Szczerba & Székely (2002).

## 3. Results

In this section, we present results corresponding to a number of different stages in the development of our tumour model. To initialize the pathology growth inside the hosting tissue, we place a small avascular ball inside a larger domain representing the healthy tissue. The initial tumour cluster does not contain any necrotic tissue and is small enough to receive enough oxygen for rapid growth. At a certain distance from the tumour, we define several pre-existing vessels, which will be the origin of the new sprout formation. The growth process is iteratively solved and updated, in each time step solving quasi-static problems using the most recent solution from the other simulation components as described in §2 (figure 1).

Figure 2 presents the simulation results at two different stages of the development. The figure depicts a volume rendering of the angiogenesis growth factor concentration and the developing vascular system. The first stage represents the avascular phase, with diffusion from the neighbouring healthy tissue as the main source of oxygen. In figure 2*b*, the vasculature has penetrated hypoxic area and is already dense enough to support a dramatic increase in growth rate. The formation of a vascular shell around the tumour characterized by the rather chaotic geometric arrangement of the individual vessel segments can be well explained by the diminishing gradient of TAF concentration resulting in decreasing influence of the directed motility during tip migration. The resulting geometry resembles qualitative experimental observations made on casts of tumour neovascularization (Walocha *et al*. 2003). The spatial distributions of some governing factors and characteristic system descriptors are shown in figure 3, corresponding to the same time step as in figure 2*b*. The oxygen concentration, the local growth expansion in per cent and the TAF concentration are shown in figure 3*a–c*, respectively. In figure 4, we have computed various entities that quantitatively describe the development process and allow the identification of characteristic events during the process. Initially, the growth seems to be more or less linear until the angiogenic switch is turned on. The last third of the process shown exhibits a nearly exponential volumetric growth. This correlates well with the onset of the availability of additional oxygen supply driven by the flow through the newly formed vessels, as shown in figure 2*b*. In figure 4*d*, we see that the volume of the hypoxic region of the tumour increases at the beginning, but as the new oxygen source becomes available it gradually decreases again until it reaches zero.

We have verified the numerical accuracy of our results in several complementary ways. The numerical properties of the used Galerkin finite-element methods on Lagrangian linear/quadratic tetrahedral meshes are well documented to be *O*(*h*^{p}) convergence of the energy norm for a polynomial order *p* finite-element method (Zienkiewicz & Taylor 2000). In order to find a suitable mesh resolution, we have performed the same simulation with higher polynomial order finite-element methods. This has the advantage that it does not change the performance of the search algorithms in the hybrid oxygen delivery module, where we need to locate a vessel segment within a finite element. The results suggest that the chosen mesh resolution is sufficient to accurately solve the involved equations (reaction–diffusion and elasticity).

The verification of the independence of time steps within the quasi-static algorithm is complicated by the fact that angiogenic growth is constituted by discrete events. We suggest that time steps should be chosen at the same order of magnitude as eventually prescribed by the angiogenic growth time scale.

## 4. Discussion

### (a) Alternative tissue models

The treatment of the tissue response to volumetric growth using a linear constitutive law is certainly a major simplification and relying on a nonlinear hyperelastic growth model as discussed by Ambrosi & Mollica (2004) could be considered. Nevertheless, tissue is a highly complex material, which depending on the time scale and applied loads, can be hyperelastic, viscoelastic or plastic (Gladilin *et al*. 2003). For instance, larger deformations lead to irreversible energy dissipation, for example, through relaxation or adaptive remodelling of the tissue. Inhomogeneities should also be taken into account, as parts of the solid tumour are known to be stiff and can well be approximated as elastic, while parts, such as the necrotic regions, are nearly fluid, possibly leading to convective cell mixing. This could possibly explain often highly penetrative patterns formed by the tumour–tissue interface. Such considerations are currently absent in existing theories of tissue growth, largely due to a lack of accurate experimental observations on the physiological reaction of the cancerous and the surrounding healthy tissue to stress exerted by the growing neoplasm. Under these circumstances, one cannot formulate reliable nonlinear tissue models that could improve the realism of our simulations. Once the necessary findings are available, we can easily incorporate them into our model, which is a major advantage of our framework.

### (b) Parameter estimation and validation

Our ultimate goal is to develop a comprehensive framework allowing creating, testing and validating models of tumour development, which can be used as a research tool to test the hypotheses and perform *in silico* experiments for drug and therapy design. We are currently working on setting up a corresponding experimental environment based on animal models. In order to be able to collect comprehensive spatio-temporal *in vivo* data, mouse tumour models will be investigated using multimodal non-invasive imaging approaches addressing various aspects involved in the formation of neovasculature. Two types of models will be used: subcutaneously implanted tumours for development of imaging tools and *in situ* tumour models (e.g. gliomas), which also consider the tumour–host tissue interaction. Tissue hypoxia will be assessed using fluorine-18-labelled fluoromisonidazole (Wyss *et al*. 2006), which under hypoxic conditions will be reduced to form reactive metabolites, which bind to intracellular macromolecules. The PET label is therefore effectively trapped in hypoxic cells, local tracer activity being a measure for the degree of hypoxia. Reduced cellular oxygen levels will lead to stabilization of hypoxia-inducible factor-1α (HIF-1α), which after dimerization with the β-analogue will form an active transcription factor that prompts the expression of a large number of downstream genes, several of them related to angiogenesis. The level of HIF-1α and selected downstream products will be assessed using reporter gene assays (Lehmann *et al*. 2008). A result of the activation of the HIF signalling cascade is the development of neovasculature, which can be characterized by a number of structural and physiological parameters derived from magnetic resonance imaging (MRI) studies. Immature vessels display a high vascular permeability, which can be assessed by analysing the extravasation of extracellular MRI contrast agents such as Gd chelates. This method is used both clinically and preclinically as biomarker for assessing the angiogenic potential of neoplastic tissue (as potential indicator of malignancy) and the therapeutic efficacy of anti-angiogenic drug candidates (Rudin 2007). The rate of contrast agent uptake by tissue depends on the blood flow rate, the vascular permeability and the vascular surface, important parameters characterizing the vasculature. The use of intravascular contrast agent (e.g. iron oxide nanoparticles) allows the estimation of the local tissue blood volume (Rudin *et al*. 2005) and of the local average vessel diameter (Tropres *et al*. 2001). These longitudinal *in vivo* studies will be complemented by high-resolution X-ray tomographic analysis of the tumour vascular system. The obtained data will not only provide the experimental environment of model building and simulation but also serve as a basis for the quantitative validation of the developed methods.

The resulting models will depend on a number of parameters, which have to be estimated as precisely as possible. While some of them will be directly measurable by the *in vivo* and *in vitro* measurements described above, in many cases, these may not be directly observable with available technology and must be estimated using model identification techniques. In this context, we are working on a parameter estimation framework. Biological hypotheses will be translated into the form of a mathematical model, usually depending on several unknown parameters. It will be implemented relying on the simulation framework described in this paper showing some similarity with its iterative structure. The first step is the generation of the anatomical model, which is fed as input to the main optimization loop. Based on the morphological information, which is eventually augmented with physiological properties, the simulation of the investigated process is carried out. The results of the computation are then compared with actual measurements, whereas the goal function of the minimization is calculated based on the detected deviations. The unknown parameters are then adjusted iteratively, until maximal similarity between the measured and calculated values can be reached. The framework is shown in figure 5. The high computational cost per iteration—in each step, a complete simulation must be performed—will force us to use as few variables as possible. Simultaneously, minimizing the number of free parameters also reduces the potential of overfitting to the measured data.

An important aspect of this framework is the possibility to identify outliers and systematic deviations between the model and the reality (e.g. by residual analysis; Ljung 1999). These can be analysed and interpreted by the user in order to modify and refine the overall model, then restarting the parameter estimation, and thus closing the highest level of a nested iterative analysis loop. Once carefully validated, these methods will serve as a powerful and efficient framework for formulating and testing hypotheses about various aspects of tumour growth. In addition, this approach offers the most generic way to bridge over arbitrary large gaps in scale by naturally and seamlessly integrating input and output variables as well as model parameters connected to all relevant components of the system and processes investigated, from single molecules to whole organisms.

### (c) Sensitivity analysis

As part of the process of parameter identification and model validation, it is good practice to test the sensitivity of the model behaviour with respect to its parameters. This is helpful in providing guidance as to how an experiment should be performed, by identifying measurable variables that are sensitive to certain parameters. Additionally, it provides information on how confident we can be about identified parameters or how accurate they need to be in order to make the model sufficiently useful and valid. Such analysis of the model described in this paper should test the sensitivity of entities like local volumetric growth with respect to parameters of the lumped proliferation model (e.g. thresholds) or the importance of vessel wall permeability, network architecture or connectivity on oxygen concentrations, drug delivery and other transport mechanisms.

## 5. Conclusion and outlook

We have presented an integrated framework to model solid tumour development, handling some of the major processes, i.e. the mechanical deformation due to growth, the biochemical response to hypoxia, blood flow, oxygenation and the explicit development of a vascular system in a coupled way. The results demonstrate the feasibility of the approach. Unlike previous methods to simulate angiogenesis, we do not treat the vascular sprouting remodelling as a process that takes place within a static domain, or as a static distribution of growth factors. Instead, our angiogenesis and remodelling approach is coupled to a dynamically growing domain, with evolving TAF concentrations. Whereas in a fixed domain vascular growth and remodelling eventually reach a static equilibrium, this is not the case for vascularized tumours. With our modelling framework, we are now able to perform *in silico* studies to identify factors that could possibly stabilize this process.

We are, however, aware of the limitations of the current model. The description of the tissue as elastoplastic is a simplification, as discussed in §4. Within our framework, it is, however, straightforward to exchange the current material law with any other behaviour. This could be done in a horizontal fashion, by using a different tissue-level growth model, or vertically by replacing the functional relationship between incremental volume expansion and oxygen concentration (lumped model) by a cellular-level simulation of the cell population dynamics. This could, for instance, be implemented by running a separate simulation at each Gauss integration point (Nickerson *et al*. 2005). Several other extensions and improvements are planned, including a more detailed representation of the metabolic activity, treatment of individual growth factors and the role of different cell types, including endothelial cells and macrophages.

Finally, it has to be noted that the level of detail that can be handled is limited by the computer resources that are available. For instance, the simulation of the vessel system as an explicit network structure allows us to study important aspects such as blood transport, remodelling and leakiness of vessels. However, this does not come without a significant computational burden. Overall, the number of unknowns necessary to capture local variations in our computations are in the range of 10^{6}, which is close to current constraints imposed by computer memory and the simulations need computation times from hours to several days. Therefore, we are convinced that for any Physiome-type approach to tumour modelling it will be necessary to rely on current advances in high-performance computing. It is inevitable that the envisioned *in silico* modelling environment is implemented in a grid-enabled, massively parallel, distributed fashion, similar to what is used today for weather forecasts.

## Acknowledgments

This work has been performed within the frame of the Swiss National Centre of Competence in Research on Computer-Aided and Image-Guided Medical Interventions (NCCR Co-Me) supported by the Swiss National Science Foundation.

## Footnotes

One contribution of 11 to a Theme Issue ‘The virtual physiological human: building a framework for computational biomedicine II’.

- © 2008 The Royal Society