Two handy geometric prediction methods of cancer growth

In present day societies, cancer is a widely spread disease that affects a large proportion of the human population, many research teams are developing algorithms to help medics to understand this disease. In particular, tumor growth has been studied from different viewpoints and different mathematical models have been proposed. Our aim is to make predictions about shape growth, where shapes are given as domains bounded by a closed curve in R2. These predictions are based on geometric properties of plane curves and vectors. We propose two methods of prediction and a comparison between them is shared. Both methods can be used to study the evolution in time of any 2D and 3D geometrical forms such as cancer skin and other types of cancer boundary. The first method is based on observations in the normal direction to the plane curve (boundary) at each point (normal method). The second method is based on observations at the growing boundaries in radial directions from the ”center” of the shape (radius method). The real data consist of at least two input curves that bind a plane domain.


Introduction
The evolution in time of some objects is the subject of study of many researchers worldwide. Special attention has been given to cancer, and a way to understand this disease is to know how it evolves over time.
We can find, in reading, many studies about mathematical modeling of tumors, which attempt to predict the growth of tumors from a mathematical point of view.
Williams and Bjerknes [14] introduced a stochastic model for the spread of cancer cells. Cells, both healthy and diseased, are situated on a planar lattice. Tumor extension is through cell division, one daughter keeps his position, while the other usurps the position of a neighbor; abnormal cells reproduce at a faster rate than normal cells.
One model that has been used to describe tumor growth is the exponential growth model [15] and another model used to describe tumor dynamics is a Gompertz curve or Gompertz function [11]. This is a type of mathematical model for a time series, where growth is slowest at the end of a time period [15]. A proposed mathematical model based on energy conservation (Universal Law model) was derived to model tumor growth [13]. This model was tested against empirical data and the results fit a variety of in vitro and in vivo data [5].
The development of tumor models is important as they offer a way to better understand the kinetic growth of malignant tumors which may lead to the development of successful treatment strategies.
Our aim is to propose two simple prediction methods of the shape growth, where shapes are given as bounded domains by a closed curve in R 2 .
Let us assume that the cells, both normal and abnormal, are situated on a planar lattice R 2 and let D t (D t ∈ S , where S = sample space or Region of Interest) denote a bounded domain occupied by cancer cells at time t. Given the domain D t ∈ R 2 we consider that it is bounded by a simple closed curve parameterized by arc length α t : [0, L] −→ R 2 . So, the original cancerous population occupies the planar domain D 0 . For simplicity, Williams and Bjerknes [14] restricted their attention to the process α 0 starting with a single abnormal cell at the origin. In this context, the following question arises: how fast does α 0 grow, and what is the geometric nature of α t for time t > 0 (but not too large).
To make the prediction we use two geometrical methods and start from the hypothesis that the speed of variation in time (growing) is constant in each direction (but not equal). So the evolution in time of the tumor can be expressed by the variation of each vector. Starting from this, the problem is resumed to determine the values of each vector for a future time (t + ∆t).
In the first method (normal method) we construct the vectors in the normal direction to the curve and in the second one (radial method) we construct the vectors as an extension of the radius of a circle in which the tumor can be enrolled.
The usual data consist of in at least two curves that bound two growing planar domains. The first curve can be the contour of one tumor when it was discovered (at time t), and the second curve is the same tumor after a while ∆t. For real data analysis both curves can be provided by two analysis with computed tomography (CT), at a time interval.
Although these curves are continuous (parameterized) curves, all our calculations, which are based on the comparison of these curves, are done on digital curves, that is, a discretized version of the parameterized curves using contour points of each curve.
The values of the vectors from a specific time can be calculated if we know the parametric function or the discretized step. From both methods of calculus the input data represent the coordinates of the contour points.
In section 2 we present a brief theory of growth shapes and curve deformations, and we propose our two prediction methods of tumor growing, depending on the direction of growing chosen from each point of the curves.
In Section 3 we implement all the mathematical calculations in Matlab software, and we build a library of functions to run both methods. We present the results of stimulation study based on random and parametric curves. We also present the analysis of a real data set. The paper ends with some conclusions and open ideas.

Shape and growth description
At present, there exists a variety of growth models for objects in discrete space; see for instance [2], [11], [9], and [8]. In the Richardson model, introduced in [12], the growth is described by a Markov process. For a growing object in the plane, the state at time t is a random subset Y t of Z 2 consisting of the "infected sites", and Y 0 (initial tumor) consists of a single site.
So, it is deduced from the preceding results that the tumor shape Y t at present time t, depends on the structure of the initial tumor shape Y 0 . Then the tumor shape in a future Y t+∆t , is a function which depends on the edge and structure of the cancer in the present time Y t , and also on some external factors like mitosis, nature of cancer (benign or malign), density, etc all these factors can be included in a function g(t) ; then, In our study we consider the boundary of the tumor at different times (not the entire tissue). This boundary is represented by a closed curve which may be star-shaped or of any other random form. In real studies, a tumor is discovered time after starting its development and, for this reason, the initial domain D 0 does not consist in just one cell but in a closed domain bounded by a curve α 0 . So, each domain D t ∈ R 2 is bounded by a closed curve α t and we propose two simple methods to predict the tumor growth.
To predict the tumor growth means to find the curve α t+∆t based on observations at times t i < t + ∆t. Both methods are based on geometric properties of curves and vectors and the main differences between them consist in the directions of the vectors chosen at points P i ∈ α t (s i ).
The first method (normal method) can be applied to general curves (not necessarily star-shaped) but the observed curves must be close enough to avoid self-intersections between normal lines. The second method (radial method) can be applied only to star-shaped domains with respect to a point in D 0 .

Normal method
The first method is based on observations in the normal direction to the plane curve (boundary) at each point.
We consider negatively oriented planar closed curves α : I −→ R 2 (that is, when traveling on α one always has the curve interior (tumor) to the right). Moreover, when the curve α is parameterized by arc length s; then, the signed curvature function of α(s) is defined as [4] dT where T (s) is the unit tangent vector and N (s) is the unit normal vector oriented to the exterior of D (like the sunlight). The sign of the signed curvature κ indicates the direction in which the unit tangent vector rotates as a function of the parameter along the curve. If the unit tangent rotates counterclockwise, then κ > 0. If it rotates clockwise, then κ < 0 (see Figure 1).
For a plane curve given by a parametrization α(s) = x(s), y(s) , the signed curvature is expressed as (3) Figure 1: Signed curvature κ for a negatively oriented planar closed curve α To describe growth shape, our data consists of a discrete number of shapes D 0 , D 1 , . . . D n obtained at times t i = t 0 + i∆t, respectively, and bounded by closed simple curves α i , i = 1, 2, . . . n, where each curve α i is defined from the preceding one as where N i−1 (s) is the unit normal to α i−1 (s) and f i−1 (s) is a differentiable function.
In order to make a prediction, that is, in order to construct a curve α n+1 (s) from the information provided by α 0 , α 1 , . . . , α n and therefore a function f n (s) from f i−1 (s), i = 1, 2, . . . n, we consider the curve: where α n (s) and N n (s) are defined in (4) and we suppose that f n (s) = f n−1 (s), ∀s ∈ I. Then, we suppose that f i (s) = f i−1 (s), ∀s ∈ I and for i = 1, 2, . . . n which means that the function Note that if f i (s) = k, constant, then α i+1 (s) is the parallel curve to α i (s) at a distance k [6]. This is the simplest case when the evolution in time in each direction is constant and equal to k. Note that to ensure that the curve α n+1 is well defined it is necessary that, for all s ∈ I, the distance from the curve α n to the point α n+1 (s) will be f n (s).
Therefore, in our prediction, we take into consideration that the growing rate is different for each normal direction but constant for equal periods of time. In real studies the tumor is discovered after some time t and instead of a parameterized curve the study is based on digital curves.
In Figure 2a we show the digitization of a brain tumor at time t 1 which corresponds to curve α 1 (s) plotted at points P t 1 i (x i , y i ), with coordinates x i and y i , i = 1, . . . , n.
After a time ∆t, for t 2 = t 1 + ∆t, the same growing tumor has a different contour (Figure 2b), α 2 (s), represented at points P t 2 j (x j ), j = 1, . . . , m. The number of points n and m are factors which represent the resolution or the step of digitization for each brain tumor contour.
Note that since the curves are digital we only know a discrete number of points P t 1 i , i = 1, 2, . . . n and P t 2 j , j = 1, 2, . . . m for each curve: α 1 (s i ) and α 2 (s j ). Then, in general, it is possible that the estimated point P t 2 j ∈ α 2 (s j ) with j ∈ {1, . . . , m} does not correspond to a point P j ∈ α 2 (s j ). For this reason, we join the points P t 2 j (x j , y j ) by straight segments and we obtain an approximated polygonal curve to α 2 (s j ) that, as there is no confusion, we will denote it also by α 2 (s). Now, from the n points P t 1 i (x i , y i ) in curve α 1 and an approximation of the values of N (s i ), we will compute the points in the polygonal curve α 2 . We will suppose that we have two curves α 1 and α 2 and we will predict, under our assumptions, the curve α 3 .
We proceed as follow. We consider three consecutive points P t 1 i−1 , P t 1 i and P t 1 i+1 and the segment P t 1 i−1 P t 1 i+1 ; then we draw the perpendicular line to this segment which passes through the point P t 1 i . The intersection between this line and the polygonal curve α 2 corresponds to the calculated point P t 2 j and the distance between P t 1 i and P t 2 j is f 1 (P t 1 i ). We repeat this process for i = 1, 2, . . . , n and we obtain all the points of P t 2 j , for i = 1, 2, . . . , n. Other methods to approximate the normal direction to a point P t 1 i , like the bisector direction or the median direction, can be found in [1]. Figure 3: Calculus of the predicted points P t 3 k with the normal method At each point P t 2 j we determine the normal vector N 2 (P t 2 j ) to the polygonal curve α 2 and, following that direction, at a distance f 1 (P t 1 i ), we plot the predicted points P t 3 (x j , y j ), j = 1, 2, . . . , m, of the predicted curve α 3 (s), at time t 0 + 2∆t (see Figure 3).
The condition that must be required to apply this method is that vectors − −−− → P t 2 i P t 3 j do not intersect between them for i = 1, 2, . . . , n. This can be accomplished when the distance f 1 (P t 1 i ) for each point P t 1 i does not exceed the value of the κ(P t 1 i ) [6]. There exist several different points approximation methods for the discrete approximation of the curvature; for instance, circle approximation or angle approximation (see, for instance [10]). However, we will apply the Archimede's theorem of the area of a parabolic segment to approximate the curvature from a parabola [7].
Given P t 1 i and the segment P t 1 i−1 P t 1 i+1 whose length is c, we have that [7], where A( i ) is the area of the triangle of vertices P t 1 i−1 , P t 1 i and P t 1 i+1 .
Since the distance d from the point P t 1 i to the line defined by P t 1 i−1 P t 1 i+1 (see the triangle 2 = P t 1 1 P t 1 2 P t 1 3 in Figure 3) is: where: The approximation of the curvature when we know the coordinates of the points P t 1 i−1 , P t 1 i and P t 1 i+1 is, from (6), Moreover, if we have the approximation of κ 1 (P t 1 i ) via (8) and that of f 1 (P t 1 i ) we may know if there exists a relation between local curvature and the growth shape. A discussion in this line that has been considered in interesting applications [3].

Curve evolutions
This method of predicting growth is a particular case of curve evolution, where each point of a curve α moves in the normal direction with speed equal to the function f (s) at that point. Consider a family of smooth closed curves α(s, t) where t means time and s is the parameter of the curve ( s and t are independent) and suppose that α(s, t i ) = α i (s), for i = 1, . . . , n.
The mathematical formulation in this case is where N (s, t) is the normal vector to the curve α(s, t) and F (s, t) is the speed function. In principle the function F may depend on many factors like local and global properties of the growing curve and time t. However, in our method F does not depend on t.
For numerical implementation we approximate Then, since ∆t is constat for each i, we have: Several examples of functions F (s, t) can be found in [1], for instance, when the function is the curvature of the initial curve, we have a curvaturedriven evolution of the initial curve α 0 .

Radius method
In this case we consider star-shaped domains from an origin O(x 0 , y 0 ) and vectors which are represented by the radius line from this origin to the each point P t 1 i of first curve α 1 in that direction. This method is quite simple to understand from the theoretical viewpoint and also from the point of view of calculus. Using this method the choice of the point of origin of the tumor is important O = O(x 0 , y 0 ) (see Figure 4a).
Starting from this point we construct a line from each contour point P t 1 i of the curve α 1 and continue to the intersection of the second contour α 2 which give us the points of the second contour P t 2 j , for each i = 1, 2, . . . , m. These lines with the direction from the center point O to the points P t 1 i will be called radius vectors and denoted by − → R t 1 i . The estimated points P t 2 j of the second curve α t 1 +∆t , as it is a second time, are: P t 2 1 , P t 2 2 , P t 2 j , here j = 1, 2, . . . , n. The third curve (the simulated curve) is α t 1 +2∆t because the time intervals between t 1 , t 2 and t 3 are supposed to be equal. These three curves are denoted by α 1 , α 2 and α 3 respectively.
For each j = 1, 2, . . . , n the distances between P t 1 i and P t 2 j are the values of the function f 1 (P t 1 i ) and we can calculate the spatial coordinate of the contour points P t i k with j = 1, 2, . . . , n at time t 1 + 2∆t.

Simulated data: random curves
The simulation consist in creating three discrete closed curves with different areas, keeping the same center. Curves are created by joining a predetermined number of points arranged in a circle of radius defined by the user plus a uniformly distributed random variable. Each curve is simulated with a different number of contour points. After generating the curves we apply both prediction methods, and calculate the predicted contour for each method. Figure 4: Simulated tumors at time t, t + ∆t and t + 2∆t using random curves: a) prediction with radius method; b) prediction with normal method Radius method: we start to construct each vector − −−− → R i (P i ) from the center O(x 0 , y 0 ) of the first tumor α t to each contour point P i (x i , y i ) at time t, and we continue to the intersection of the second contour P j (x j , y j ) at time t + ∆t and so on by the predicted point P k (x k , y k ).
So as not to overload the picture, we plot just the first two iterations and the corresponding points of calculus P t 1 1 , P t 2 1 , and the predicted point P t 3 1 (see Figure 4a). Normal method: in order to begin the simulation, the first vector is the normal vector − −−−− → N 1 (P t 1 1 ) to surface α t in point P 1 (x 1 , y 1 ), and in that direction we find the point P j (x j , y j ) i.e P t 2 1 at the intersection with the second contour α t+∆t . Starting from this point in the direction of the normal vector to second contour, we calculate the spatial coordinates of the predicted point ( Figure  4b).
To show the precision of these methods we repeated simulation 100 times and calculated the values of absolute error (ε a ) and relative error (ε r ) between the area of the contour curve at time t + 2∆t names as t 3 i.e α t+2∆t (s) or α 3 (P t 3 i ) (the area of the polygon determined by the points P t 3 i ) and the area of the predicted curve (the polygon determined by the estimated points P t 3 k ). The results are given in Table 1 did not require too many computing resources. In a machine with i3 processor and 4 GB of RAM the time for simulating the growing tumor based on three curves, each one with resolution n = m = 50 points, based on normal method was 1.897642 sec, and it is very accurate.
The running time to simulate the growing of the tumor with the radius method, in the same conditions as in the preceding method, was 0.909440 sec, which means a decrease of the computational effort.
When star-shaped tumors are considered, both methods are accurate in prediction and are computationally fast, and can be used with success in lower computational machines.
Comparing the results we obviously note that the radius method is more effective, and the user requires less recourses from the computing machine, so it is very clear that it is a better method.

Simulated data: parametric curves
To check for the accuracy of each algorithm, we now consider a parametric irregular form and its growth at regular intervals of time. The curves come from a Fourier series expansion of a sine and cosine function.
In Figure 5 we plot in red the curve α t 1 , in green the curve α t 2 and in blue the curve α t 3 . Taking into account the first two boundaries of the shape we predict the third curve (drawn in magenta in Figure 5) and we compare with the real evolution of the function.
We also calculated the area of the surface drawn by the predicted contour and the errors absolute and relative reported to the area of surface bounded by curve α t 3 : a) relative error : radius method: ε r = 0.005278 normal method: ε r = 0.005389 b) absolute error : radius method: ε a = 0.002709 normal method: ε a = 0.002766 As can be noted, the absolute and relative errors for the normal method are larger than these for the radius method.

Real data
We now show the performance of the two methods over real data. It was necessary to implement, in one library, a series of functions to process images from a magnetic resonance tomography, computed tomography or any other analysis or image of a tumor. The input data is represented by a complete set of images of a brain tumor taken at intervals of one month, to which we apply the methods to predict its growth.
This particular tumor is found in the Central Nervous System and is called glioblastoma multiform. In conformity with the World Health Organization, this tumor is the most aggressive tumor and the type and grade is according to IV-th classification. In the scans that we are studying it is clear that there is a presence of multiple tumors in the body which is known as metastasis. These metastatic tumors are children of primary tumors from the breast, lung, colon, stomach and skin (melanoma), but in our case, the first was the brain tumor.
From this complete set of analysis we selected three images taken in the same plane: one from November 9 th , 2009 (see Figure 6a), one from December 8 th , 2009 (Figure 6b) and the last from January 10 th , 2010 (Figure 6c).
The input data for this analysis was the first two curves.The boundary tumor of November was digitalized in 64 points: α t 1 = {P t 1 1 , P t 1 2 , . . . , P t 1 64 }. The curve is closed, so P t 1 1 ≡ P t 1 64 . The second one (corresponding to December) was approximated with 61 points, then: α t 2 = {P t 2 1 , P t 2 2 , . . . , P t 2 61 }. We applied both prediction methods starting from the first two images and we calculated the prediction of growth of this tumor, which we scheduled for January 10 th , 2010. We can directly compare the results of the prediction methods, with the third image from the set of analysis.
The reference curve (plotted in blue in Figure 7) is the approximation in 60 points of the real contour of brain tumor. Because we start from the first curve and all the calculus are related to the number of points P t 1 i , (in this case i = 64), the number of points of the prediction is identical, i.e k =64 points.
The area of the surface domain D 3 representing a polygon with 59 sides is Area(D 3 ) = 7852.13479 and the predicted areas are: a) radius method: Area(D 3 ) = 7774.37985; b) normal method: Area(D 3 ) = 7776.43387.
The absolute and relative errors are: a) radius method :(see Figure 7a) -absolute error: ε a = 77.7549 -relative error: ε r = 0.00990 b) normal method :(see Figure 7b) -absolute error: ε a = 75.7009 -relative error: ε r = 0. 00964 We note that the results are directly affected by the digitization of the curves by the physicians when they choose the contour points of the tumor. Also we must take into account that the resolution of these kinds of images is low (512x512 pixels). We would obtain better results if the input images

Conclusions and future research
Mathematical modeling always tries to find a compromise between simplicity of analysis and requirements of realism. On the one hand, we have extremely complex natural and biological systems; on the other hand, we need to formally address some quantitative issues about these systems which can often be done only through the use of mathematical models. For most of the realistic problems, the solution of the corresponding exact equation is, in practice, impossible, so we need to make approximations. Making approximations to solve very difficult problems is not a new idea. Appropriate models enable accurate prediction of future behavior which can be used to control and optimize various aspects of the system in question.
Both prediction methods proposed here are simple computing, fast, and provide good results. They can help medics cure and better understand the propagation of cancer. For a better adjustment of the prediction models in the very near future, we are trying to develop new models without restricting the growing velocity. In these cases, the function f (x) which transforms the curve α t into the curve α t+∆t is a second-order degree function and from each point of the curve we must solve an equation system in order to find the parameters of each function f i (x). We hope that this will lead us to remove the predictions error in comparison with real data.