VISUALISATION OF THE NUMERICAL SOLUTION OF PARTIAL DIFFERENTIAL EQUATION SYSTEMS IN THREE SPACE DIMENSIONS AND ITS IMPORTANCE FOR MATHEMATICAL MODELS IN BIOLOGY

Numerical analysis and computational simulation of partial differential equation models in mathematical biology are now an integral part of the research in this field. Increasingly we are seeing the development of partial differential equation models in more than one space dimension, and it is therefore necessary to generate a clear and effective visualisation platform between the mathematicians and biologists to communicate the results. The mathematical extension of models to three spatial dimensions from one or two is often a trivial task, whereas the visualisation of the results is more complicated. The scope of this paper is to apply the established marching cubes volume rendering technique to the study of solid tumour growth and invasion, and present an adaptation of the algorithm to speed up the surface rendering from numerical simulation data. As a specific example, in this paper we examine the computational solutions arising from numerical simulation results of a mathematical model of malignant solid tumour growth and invasion in an irregular heterogeneous three-dimensional domain, i.e., the female breast. Due to the different variables that interact with each other, more than one data set may have to be displayed simultaneously, which can be realized through transparency blending. The usefulness of the proposed method for visualisation in a more general context will also be discussed.


(Communicated by Qing Nie)
Abstract.Numerical analysis and computational simulation of partial differential equation models in mathematical biology are now an integral part of the research in this field.Increasingly we are seeing the development of partial differential equation models in more than one space dimension, and it is therefore necessary to generate a clear and effective visualisation platform between the mathematicians and biologists to communicate the results.The mathematical extension of models to three spatial dimensions from one or two is often a trivial task, whereas the visualisation of the results is more complicated.The scope of this paper is to apply the established marching cubes volume rendering technique to the study of solid tumour growth and invasion, and present an adaptation of the algorithm to speed up the surface rendering from numerical simulation data.As a specific example, in this paper we examine the computational solutions arising from numerical simulation results of a mathematical model of malignant solid tumour growth and invasion in an irregular heterogeneous three-dimensional domain, i.e., the female breast.Due to the different variables that interact with each other, more than one data set may have to be displayed simultaneously, which can be realized through transparency blending.The usefulness of the proposed method for visualisation in a more general context will also be discussed.

1.
Introduction.Mathematical modelling of spatio-temporal biological problems often results in the generation of coupled multidimensional data sets.Solutions of one-dimensional models of the form u(x) can simply be visualised using a twodimensional ux-plot (i.e., a graph of u against x).Two-dimensional functions of the form u(x, y) can be represented using either a three-dimensional uxy-plot or a coloured two-dimensional image, in which different colours represent different function values (Figure 1).The mapping of a function which is the solution to a partial differential equation in three spatial dimensions u(x, y, z) is difficult, as a "fourth dimension" has to be introduced for the function values at each point in time, and to imagine and visualise structures in four and higher dimensions is nontrivial.One way to visualise the structure of those four-dimensional data sets at different points in time is to present two-dimensional slices through the volume at discrete spatial intervals.Figure 2 shows an example of 2D slices through a volume data set in MATLAB.This technique has the major disadvantage that if the data set is highly spatially heterogeneous, interesting structures might be lost between these discrete slices or hidden by other slices.The visualisation of the interaction of coupled equations with each other is even more difficult but often necessary, especially in mathematical biology.The extension of mathematical models to a third spatial dimension is important, to provide an ideal platform on which both biologists and mathematicians can communicate with each other.For cancer modelling, in particular, three-and four-dimensional visualisation can be useful for surgeons to localize the estimated tumour position within the domain for surgery and treatment planning.
In this paper we discuss the marching cubes visualisation algorithm combined with transparency blending, and present an adaptation and its usefulness in displaying coupled four-dimensional data sets generated from a mathematical model in three spatial dimensions.The mathematical model discussed here is an established model that describes solid tumour growth and invasion into surrounding healthy tissue, but until now it has only been solved in one and two spatial dimensions, as the visualisation of three-dimensional solutions remained challenging.Here we present first simulation results in three spatial dimensions and then change the computational domain to represent the female breast.Despite this specific example it should be emphasised that this non-commercial technique is easy to implement and adapt and therefore useful for displaying any four-dimensional data set produced from the numerical/analytical solutions of PDE models in three spatial dimensions.
2. The marching cubes visualisation technique.When solving spatial models numerically using finite difference schemes, we subdivide space into equal line segments (1D), squares on a lattice (2D), and cubes (3D).In three-dimensional space, to generate such equal cubes (referred to as voxels), we use orthogonal sets of parallel planes that have equal spacing between them [4].The resulting vertices P x,y,z at the intersections of the planes are assigned different values as determined by the model solution at position (x, y, z).Thus the solution of a model in three space dimensions results in four-dimensional data sets.
Most three-dimensional volume rendering and surface construction algorithms have their origin in medical imaging [10] to find structures with specific properties within the datasets.These properties are usually grey values in computer tomography or magnetic resonance images.
In most mathematical models the information we want to visualise is either density values for different variables or Boolean information, i.e., whether an object is present at a certain point in space.From this information we have to construct a surface which shows us the exact location and structure of the objects of interest.Similar objects are classified by similar grey values in medical images or similar density values resulting from our mathematical model.To obtain a graphical representation of objects of the same class, so-called isosurfaces are computed [7].Isosurfaces in medical images are defined by the same grey value, and respectively in mathematical models areas defined by the same density value v: f (x, y, z) = v may form a continuous object.
The marching cubes algorithm [10] uses a divide-and-conquer approach to locate surfaces in a data set.The divide-and-conquer design paradigm is well established in computing.Large problems are divided into subproblems that are conquered by solving them recursively.The solutions of all subproblems are combined to obtain the solution for the original problem.
As the numerical grid is subdivided in cubes, the same subdivision can be used to find surface structures.Instead of analyzing the entire data set at once, we analyse each cube independently.All surface patches that are found in this way are combined to obtain the overall surface.Each cube has 12 edges that connect 8 vertices which hold information on the calculated equations.If vertices on one edge have a value inside and a value outside the structure, i.e., below or above the isovalue threshold, we know that the surface intersects this edge.The intersection of the surface with the edge is interpolated using the data values at both vertices.The interpolation is important to obtain a smooth surface approximation with only a little error, especially when there is a big difference between the values on both vertices of the edge, as shown in Figure 3. From the intersections of the isosurface with the cube's edges we can than approximate the surface that intersects the cube.With eight vertices in a cube and two possible stages for each vertex, there are 256 different cases of how the surface intersects the cube.Rotation and inverting vertex values (the two possible stages for each vertex) reduce all possibilities to 15 basic ways of intersection.These basic cases as well as the triangulation of the obtained surface patch is shown are Figure 4.The marching cubes algorithm comprises the following steps: • define threshold for isosurface • consider each voxel • reduce problem to one of the basic cases • approximate intersection with edge • draw triangular surface patches.
To visualise surfaces in a data set with different structures we have to know the value ranges of the structure that we want to visualise.For example, to cite a problem arising in clinical practice, to display and visualise bones in a computer tomography data set, we have to know the range of grey values that represent bones in this imaging technique.This information is used as the threshold to filter relevant information.
The mathematical model that we will consider emphasises that there are different densities of tumour cells within the solid tumour as well as in the tissue.The density value may also be interpreted as the probability there are cancer cells at a certain point in space.It is very unlikely that surgeons find all cancerous cells in the tissue, given that single tumour cells are not detectable.The computational results from the model can therefore help to predict how far into the tissue the tumour has invaded and where to make the necessary surgical cuts.The display of density variations within the tumour is also useful for researchers to analyze the internal structure and changes due to a heterogeneous environment Figure 4. Fifteen basic cases for the surface intersecting a cube and display of the triangular surface parts.In the simplest case the cube is either completely outside the object or completely inside, and thus the surface does not intersect the cube (case 1).The other fourteen cases (cases 2 − 15) describe the possibilities of intersection.Adapted from [10].and proliferation.Different density levels are therefore used as thresholds for the isosurface calculation.

Transparency blending.
To study the internal structure of different objects and to visualise the interaction of structures with each other, different isosurfaces or different data sets have to displayed at the same time.This is only possible in three spatial dimensions with an implementation of transparency blending.Each structure has a degree of opacity which enables the observation of objects hidden by other structures.Figure 5 shows the importance of transparency blending to observe different densities and structures in the interior of a displayed object.The cubes in the left and right panel look similar from outside, but through transparency blending we can visualise not only the outer surface, but also structures that are hidden inside the object.While the left cube is homogeneous, the right cube features density variations inside.With the opportunity to look inside different objects, we can again use a colour map to describe different density or concentration values inside the structures that are shown.

4.
Marching cubes adaptation to reduce calculation time.The marching cubes algorithm is a powerful technique to visualise three-dimensional spatial data sets.However, there are various adaptations that focus on speeding up the computation time [3,5,8].In this paper we present a straightforward extension of the marching cubes algorithm to speed up the calculation and visualisation time of mathematical models, as the computation time is especially important for the interaction with large data sets.For the visualisation of isosurfaces, we assume that only subregions of the data set have to be rendered.Hence, the extraction of Through transparency blending we can observe not only the outer surface of one isocontour (left) but also the interior of a body (right).The visualisation of only the outer surface of the cube hides internal structures, whereas with the introduction of transparency we can observe difference isocontours inside the cube.Different pink colours represent different isovalues, with less intense pink meaning lower isovalues than intense pink isocontours.an isocontour does not require searching all the voxels if we pre-analyse the data and store it in an efficient data structure [3].We introduce a "search list" in which we store all cubes that contain relevant information on certain structures with all other cubes being neglected.To ensure fast access to voxels with the same isovalue, we sort the list with respect to the minimum vertex value.With established search algorithms we can ensure efficient access to voxels containing relevant information which need to be displayed.
Figure 6 shows a scene in which only a small spherical object is located in the center of the domain.In this case we consider a 25 × 25 × 25 domain (number of voxels n, n = 25 3 = 15625), with the object intersecting only nine voxels.The object is assumed to be heterogeneous with a gradual density fall-off towards the boundary (right).When reading the data set, every voxel has to be touched once.
With the standard marching cubes algorithm, the calculation of 10 different isocontours as seen in Figure 6 (right) would then cost another 10 × n voxel considerations.In contrast with the adaptive marching cubes, if we store the 9 voxels in a separate list, we reduce the costs to 10 × 9 considerations, and therefore reduce the number of considerations by 99.94%.If the object has a heterogeneous structure, not all cubes in the "search list" have to be considered for each isovalue.If we sort the cubes in the list with respect to their minimum vertex value, we can then using search algorithms to faster find the cubes of interest for each isosurface.

Mathematical model of cancer growth and invasion.
In this section we apply the adaptive marching cubes rendering technique with transparency blending to computational results of mathematical models for cancer growth and invasion.The specific model we consider in this paper is a development of the invasion model by Anderson et al. [1,2] and Enderling et al. [6].Consisting of a system of three coupled PDEs, it describes the invasion of host tissues extra-cellular matrix (ECM) by a solid tumour.Cancer cells are known to secrete matrix-degrading enzymes (MDE) which destroy the tissue upon contact to make space for the tumour to grow into.The three equations model the interactions between tumour cells (n), host tissue (f ) and matrix-degrading enzymes (m) produced by the tumour cells.We recently extended the original model by a proliferation term to model breast cancer specific tumour dynamics [6].Previous analysis of the models has only been done in one and two spatial dimensions as the 3D interpretation of the results remained challenging.In this paper we present the extension to a third spatial dimension and an irregular domain to model the female breast shape.
Tumour cell migration is assumed to consist of two components: unbiased migration through random motility and biased migration through haptotaxis.Haptotaxis is migration up a fixed or bound gradient of ECM molecules, such as fibronectin.In our model, we include a term to simulate the effect of proliferation of tumour cells.For simplicity we neglect cell death and consider proliferation to be modulated only by available space.To create space, cancerous cells are known to produce matrixdegradative enzymes that diffuse into the surrounding ECM and degrade it upon contact [11].Hence the system of equations in non-dimensionalised form is where d n , d m are diffusion coefficients of the cells and MDEs, respectively, and λ, γ, η, α, β are positive constants.The estimation of these parameters is discussed in detail in [1] and [6].The proposed model ( 1) is now solved in three spatial dimensions.Initially a spherical tumour is placed in the center of the domain.For the visualisation of the domain boundaries, we use a box.(All simulation videos are available online: http://www.maths.dundee.ac.uk/%7Eheikman/visualisation) Figure 7 shows the growth of a three-dimensional tumour in homogeneous and in heterogeneous tissue (top and bottom row, respectively).To emphasize the effect of haptotaxis, the ECM is modelled as a heterogeneous tissue.The heterogeneity is obtained by a combination of sine functions with frequency variations in x, y and z directions.Figure 8 shows the initial setup of the spherical tumour in the center of the domain as well as different densities in the surrounding tissue.More intense blue means areas of high tissue density, whereas less intense blue represents areas in the tissue of lower density.1)) grows smoothly and uniformly in homogeneous tissue (f in (1)), but with spatial variations in ECM density, the tumour shapes) becomes irregular.In these simulation results only the outer tumour surface is visualised; i.e., tumour density threshold n is n > 0.
We used the aforementioned rendering and transparency blending techniques to visualise the tissue and tumour data sets at the same time.This is important when we want to study the interaction between different variables in mathematical models.Due to haptotaxis, as time evolves the tumour cells migrate directed towards areas with high tissue density (Figure 8).
The effects of proliferation upon the internal structure of the tumour is perhaps one of the more interesting questions we can examine with our model.With the aid of transparency blending, we can observe the density variations within the tumour.Figure 9 shows the results of the tumour evolution both with and without proliferation.The observable shape of the tumour is in both cases similar.The size of the tumour when proliferation is included is slightly larger.However, differences become clearer when examining the different isosurfaces.Without proliferation, the highest tumour cell densities are found close to the boundaries of the tumour close to the tissue.The location of the tissue is shown in Figure 8.In the current model, proliferation is modulated by available space.When the ECM is degraded, there is more space for tumour cells to proliferate.Thus the internal structure of the tumour changes rapidly.High densities of tumour cells are found not only nearby the tumour boundary but inside as well.6. Irregularly shaped domains.The behaviour of the solution of mathematical models and biological structures in different shaped domains is an important topic of research.Natural domains are unlikely to be cubic, so we applied our tumour growth model to breast cancer.The shape of the breast dictates that modelling should be undertaken using an irregularly shaped domain.In particular we examine how tumour invasion is modulated by the domain shape especially near the domain boundaries.
The breast shape was modelled using a modified Gaussian distribution function in three dimensions.The heterogeneous ECM density is modelled using sine functions with frequency variations in x, y and z directions.The tumour is placed close to the upper outer boundary where most early breast cancer tumours are located [9].
With aid of transparency blending, we can observe the interaction of the tumour with the surrounding tissue as well as the impact of the domain boundaries. of the ECM results in asymmetric growth of the tumour with areas of high density near areas of high tissue density.Figure 11 shows different time steps of the simulation and different view points to observe the scene.The exact location of the tumour is important for both before surgery and after surgery during therapy planning.
7. Discussion.The visualisation of the solution of system of couples PDEs in three spatial dimensions is a challenging task for mathematical modelling.The scope of this paper was to apply the established marching cubes volume rendering technique combined with transparency blending in OpenGL to mathematical models of solid tumour growth and invasion.We have introduced a third spatial dimension into a previously discussed mathematical model and changed the computational domain to model a breast-like shape.In this model, tissue heterogeneity and domain shape have a big influence on the growing tumour, especially on its geometry.Numerical analysis and computational simulation of equations describing tumour invasion in three space dimensions can be carried out using well-known numerical techniques, whereas the presentation of the results, especially the visualisation, offers considerable scope for creativity.We have presented an extension to the marching cubes technique to decrease the algorithm computing time for large  coupled data sets as produced by mathematical models.The resulting algorithm implemented in Java and OpenGL is a useful and flexible method of visualisation which will work just as well on a PC as it does on a Mac or other Unix/Linux based systems.The shaded outer surface of the tumour can be visualised within the domain and offers the possibility of finding its exact location in the domain.The introduction of transparency blending also allows an observation of density variations within both the tissue and the tumour simultaneously.The proposed method of simulating and visualising three-dimensional mathematical models of biological problems may create an ideal communication platform in the dialogue between surgeons and mathematicians.An exact study of different objects is possible given by various options in interacting with the output.Masking of different structures, zooming options and transparency blending are the major features amongst standard options like rotation.In addition, the proposed technique can be used for the visualisation of any four-dimensional data set, including the numerical solution of three-dimensional mathematical models.

Figure 1 .
Figure 1.Different visualisation options for a two-dimensional function u(x, y).The function values can either be plotted against a third axis (left) or be represented by a colour map (right).The figures are adapted from MATLAB.

Figure 2 .
Figure 2. Visualisation of 3D data with MATLAB using 2D slices through the volume.Because parts of each slide are covered by other slides, information in spatially heterogeneous data sets may be hidden (adapted from MATLAB).

Figure 3 .
Figure 3. Interpolation of the intersection of an isosurface of value 0.2 with an edge of a cube.One of the cube's vertices has the value 0.1 below the threshold value of 0.2, and the other one has a value of 0.9 above the threshold (left panel).Instead of assuming the intersection of the isovlaue 0.2 to be in the centre of the edge (centre panel), linear interpolation improves the approximation of the point of intersection (right panel).

Figure 5 .
Figure 5. Visualisation of an object in three spatial dimensions.Through transparency blending we can observe not only the outer surface of one isocontour (left) but also the interior of a body (right).The visualisation of only the outer surface of the cube hides internal structures, whereas with the introduction of transparency we can observe difference isocontours inside the cube.Different pink colours represent different isovalues, with less intense pink meaning lower isovalues than intense pink isocontours.

Figure 6 .
Figure 6.The visualisation of a small object emphasizes the benefit of an adaptive algorithm.Instead of going through the whole domain every time to look for certain isocontours (right), only a restricted search area would speed up the visualisation greatly.Different isocontours are shown by different colours with bright pink representing a high isovalue and a darker pink a lower isovalue.

Figure 7 .
Figure 7. Solid tumour growth and invasion into a homogeneous (top) and a heterogeneous (bottom) tissue.The tumour (n in model (1)) grows smoothly and uniformly in homogeneous tissue (f in (1)), but with spatial variations in ECM density, the tumour shapes) becomes irregular.In these simulation results only the outer tumour surface is visualised; i.e., tumour density threshold n is n > 0.

Figure 8 .
Figure 8. Solid tumour growth and invasion into a heterogeneous tissue.As time evolves (left to right, top to bottom), the initially spherical tumour (red, n in (1)) evolves to a heterogeneous shape due to haptotaxis towards areas of high tissue density (intense blue, f in (1)).

Figure 9 .
Figure9.Three-dimensional tumour (n in(1)) growing in heterogeneous tissue (f in (1)), both without proliferation (first row) and with proliferation (second row) included in the model.In both simulations the outer surface is similar, but the visualisation of isosurfaces with transparency blending shows huge variations in the internal structure.Shaded isosurfaces enable the exact study of areas of equal density.In the right column, tumour densities above a 0.8 threshold are visualised.

Figure 10 .
Figure 10.Tumour growth in a heterogeneous, irregularly shaped domain.The tumour shown in red in the top row is initially located close to the nearby domain surface, which is shown in dark yellow.Different tissue densities are shown in less intense blue, which denotes lower density, and intense blue, which denotes higher density, respectively.Different tumour densities are shown in the bottom row.High densities are represented by an intense pink colour and lower densities in a less intense pink, respectively.

Figure 11 .
Figure 11.Different time steps visualised from different viewpoints to enable an exact location of the tumour in the domain.The tumour is again shown in red, ECM densities with different blue colours, and the domain boundaries in yellow.