GAUSSIAN RADIAL GROWTH

The growth of planar and spatial objects is often modelled using one-dimensional size parameters, e.g., volume, area or average width. We take a more detailed approach and model how the boundary of a growing object expands in time. We mainly consider star-shaped planar objects. The model can be regarded as a dynamic deformable template model. The limiting shape of the object may be circular but this is only one possibility among a range of limiting shapes. An application to tumour growth is presented. A 3D version of the model is presented and an extension of the model, involving time series, is briefly touched upon.


INTRODUCTION
Modelling of biological growth patterns is a rapidly developing field of mathematical biology.Its state-of-the-art was explored at the successful conference On Growth and Form, held in 1998 in honour of D'Arcy Thompson (1860Thompson ( -1948) ) and his famous book, cf.Thompson (1917).Out of the conference grew a monograph which contains substantial biological material and an overview of mathematical modelling of spatio-temporal systems, cf.Chaplain et al. (1999).Examples of growth mechanisms studied are growth of capillary networks, skeletal growth and tumour growth.
Modelling of tumour growth has attracted particular interest in recent years.Tumour growth was one of the high priority topics of the recent multidisciplinary conference arranged by the European Society for Mathematical and Theoretical Biology in July 2002.More than 500 scientists from a wide range of disciplines participated.One of the subjects discussed was pattern formation problems, relating to tumour formation and progression, in particular the question of tumour shape.
The models suggested for tumour growth are either continuous or discrete.In Murray (2003), the continuous approach is explained in relation to brain tumours.The simplest models involve only total number of cells in the tumour, with growth of the tumour usually assumed to be exponential, Gompertzian or logistic (Swan, 1987;Marusic et al., 1994).More powerful deterministic models describe the change of the spatial arrangement of the cells under tumour growth.The discrete models are most often cellular automaton models, cf.Qi et al. (1993) and Kansal et al. (2000).
The growth literature contains very few examples of statistical modelling and analysis of growth patterns.An exception is the paper by Cressie and Hulting (1992).Growth of a planar star-shaped object is here modelled, using a sequence of Boolean models.The object Y t+1 at time t + 1 is the union of independent random compact sets placed at uniform random positions inside the object Y t at time t.More formally, where {x i } is a homogeneous Poisson point process in the plane and Z(x i ) is a random compact set with position x i .Note that this model is Markov since Y t+1 only depends on the previous objects via Y t .The model is applied to describe the growth pattern of human breast cancer cell islands.Practical methods of estimating the model parameters, using the information of the complete growth pattern, are devised.A related continuous model has recently been discussed in Deijfen (2003).The object Y t is here a connected union of randomly sized Euclidean balls, emerging at exponentially distributed times.It is shown that the asymptotic shape is spherical.
In the present paper, we propose a Gaussian radial growth model for star-shaped planar objects.The model is a dynamic version of the p−order shape model introduced in Hobolth et al. (2003).The object at time t + 1 is a stochastic transformation of the object at time t such that the radius vector function of the object fulfils where Z t is a cyclic Gaussian process.The coefficients of the Fourier series of Z t θ ∈ [0, 2π), have important geometric interpretations relating to the growth process.The overall growth from time t to t + 1 is determined by the parameter µ t .The coefficients A t,1 and B t,1 determine the asymmetry of growth from time t to t + 1, while A t,k and B t,k affect how the growth appears globally for small k ≥ 2 and locally for large k ≥ 2. Under the proposed p-order growth model , where the variances satisfy the following regression model The organization of the paper is as follows: in the first section, we introduce the Gaussian radial growth model.Then, we study the induced distributions of object size and shape under the radial growth model.After that, an application to tumour growth is discussed and a statistical analysis of a growth pattern is presented.An extension of the model, involving time series, is then briefly described.Finally, a 3D version of the model is presented.

THE GAUSSIAN RADIAL GROWTH MODEL
Consider a planar bounded and topologically closed object with size and shape changing over time.The object at time t is denoted by Y t ⊂ R 2 , t = 0, 1, 2, . . . .We suppose that Y t is star-shaped with respect to a point z ∈ R 2 for all t.Then, the boundary of Y t can be determined by its radius vector function Hobolth et al. (2003), a deformable template model is introduced, describing a random planar object as a stochastic deformation of a known star-shaped template, see also the closely related models described in Hobolth and Jensen (2000), Kent et al. (2000) and Hobolth et al. (2002).We use this approach here and describe the object at time t + 1 as a stochastic transformation of the object at time t, such that Here, {Z t } is a series of independent stationary cyclic Gaussian processes with Z t short for {Z t (θ ) : θ ∈ [0, 2π)}.The process Z t is stationary if the distribution of Z t (θ + θ 0 ) − Z t (θ ) does not depend on θ , while the process is said to be cyclic if Z t (θ + 2πk) = Z t (θ ), for all k ∈ Z.The initial value R 0 of the radius vector function is assumed to be known.
Note that Y t is used as a template in the stochastic transformation, resulting in Y t+1 .The increment process Z t can be written as where µ t ∈ R represents a constant radial addition at time t and U t a stochastic deformation with mean zero of the expanded object with radius vector function R t + µ t , cf.Fig. 1. (The object with radius vector function R t + µ t is in geometric tomography known as the radial sum of Y t and a circular disc of radius µ t , cf.Gardner (1995).)Fig. 1.The object Y t+1 is a stochastic transformation of the object Y t (grey), using a constant radial addition (shown stippled) followed by a deformation.
Because of the independence of the Z t s, the model is Markov in time, in the sense that it uses information about the object at the immediate past to describe the object at the present time.More specifically, under Eq. 2 the conditional distribution of R t+1 given R t , . . ., R 0 depends only on R t .The model suggested in Cressie and Hulting (1992) possesses a similar Markov property.
If Y t is non-circular, it can be natural to extend the model (Eq.2), using an increasing time change function where L t (θ ) is the distance travelled along the boundary of Y t between the points indexed by 0 and θ .Note, however, that if R 0 ≡ 0, then the boundary of Y t is expected to be approximately circular, since E(Z t (θ )) = µ t does not depend on θ ∈ [0, 2π).
The following result is important for the construction of parametric models in the framework of model (Eq.2).The result also implies a simple simulation procedure for a stationary cyclic Gaussian process on [0, 2π).
Proposition 2.1 The process Z t is a stationary cyclic Gaussian process on [0, 2π) with mean µ t ∈ R if and only if there exist λ The proof of Proposition 2.1 is not complicated.Given a stationary cyclic Gaussian process Z t on [0, 2π) with mean µ t , consider the stochastic Fourier expansion of Z t .Then, simple calculations involving the stochastic Fourier coefficients give the result.The other assertion is trivial.
Note that in Proposition 2.1, the λ t,k s are allowed to be zero, meaning that A t,k ≡ 0, almost surely.

The Fourier coefficients
k = 1, 2, . .., have interesting geometric interpretations relating to the growth process.It is clear that the coefficient A t,0 determines the overall growth from Y t to Y t+1 .The Fourier coefficients A t,1 and B t,1 play also a special role.Numerically large values of the coefficients will imply an asymmetric growth from Y t to Y t+1 .In order to interpret geometrically the remaining Fourier coefficients A t,k and B t,k , k = 2, 3, . .., let us consider an increment process for which all Fourier coefficients except those of order 0 and k are zero, θ ∈ [0, 2π), does not depend on i.Therefore, A t,k and B t,k affect how the growth appears globally for small k and locally for large k.The variances λ t,k control the magnitude of the spread of the Fourier coefficients.
Since the zero-and first-order Fourier coefficients play a special role in relation to the growth process and may in applications well depend on explanatory variables, we shall desist from specific modelling of these coefficients.In the following we will assume that A t,0 = µ t is deterministic.Furthermore, we suppose that A t,1 = B t,1 = 0 or, equivalently, we concentrate on modelling A special case of the Gaussian radial growth model is the p-order growth model.This model is inspired by the p-order model described in Hobolth et al. (2003), where the stochastic deformation process is a stationary Gaussian process with an attractive covariance structure, described below.The model is called p-order because it can be derived as a limit of discrete p-order Markov models defined on a finite, systematic set of angles θ , cf.Hobolth et al. (2002).
The parameters α and β determine the variance of lower order and higher order Fourier coefficients, respectively.Furthermore, p determines the smoothness of the curve Y .In fact, the curve Y is k − 1 times continuously differentiable where k is the unique integer satisfying p ∈ (k et al., 2003).Note that the first Fourier coefficients of Y are set to zero.
We can now give the definition of the p−order growth model.

Definition 2.3
The series Z = {Z t } follows a p−order growth model if the Z t s are independent and Z t ∼ G p (µ t , α t , β t ) for all t.
The parameters α t and β t determine, respectively, the global and local appearance of growth from Y t to Y t+1 .As before, p determines the smoothness of the curves Z t .The overall growth pattern is specified by the µ t s.Their actual form depends on the specific application.Tumour growth has often been described by a Gompertz growth pattern where κ t is the average radius at time t and η and γ are positive parameters determining the growth, implying that For more details, see e.g., Steel (1977).
Note that the p-order growth model allows for negative values of R t (θ ).However, the parameters µ t , α t and β t will be chosen such that this will practically never occur.Fig. 2 shows simulations of the increment process Z t from time t to t + 1 for different values of α t and β t under the second-order growth model, i.e., p = 2.A large value of α t gives increments that are fairly constant while a small value of α t provides a more irregular growth on a global scale.The parameter β t controls the local appearance of the increment process, the smaller β t the more pronounced irregularity on a local scale.Fig. 2. Simulated objects under the second-order growth model.The object at time t is fixed while the object at time t + 1 is simulated under the indicated values of α t and β t .

DISTRIBUTIONAL RESULTS
In this section, we study the induced distribution of object size and shape under the p-order growth model.The limiting shape may be circular but, as we shall see, there is a whole range of possibilities.
Unless otherwise explicitly stated, we assume that R 0 ≡ 0. We then have for θ ∈ [0, 2π) where and The shape of the object at time T will be represented by its normalized radius vector function which can be regarded as a continuous analogue of the standardized vertex transformation vector in shape theory, cf.Hobolth et al. (2002).
Under the assumption of independent increments, the distribution of the area of the object at time T , A(Y T ), is known, provided that the radius-vector function R T is positive.
Proposition 3.1 Assume that the radius vector function R T of the object Y T is positive and that it satisfies Eqs.6-8.Then, where V k , k = 2, 3, . . ., are mutually independent exponentially distributed random variables with mean 1. Proof.The area of the object at time T is we have that A(Y T ) < ∞, almost surely.Using Eq. 6 and Parseval's equation, we get that where V k , k = 2, 3, . . ., are mutually independent exponentially distributed random variables with mean 1.
The distribution of the area of Y T is thus a sum of independent Gamma distributed random variables.The saddlepoint approximation of such a distribution is easily derived, cf.Jensen (1992).
It does not seem possible to get a correspondingly simple result for the distribution of the boundary length of Y T .This seems apparent from the expression for the boundary length of which is valid in the case where R T is differentiable.
As we shall see now, the class of p-order growth models is quite rich in the sense that the shape of the limiting object, represented by its normalized radius vector function, may be distributed according to any porder model G p (1, α, β ) with mean 1.For large values of α and β , the shape is close to circular.
Let us consider the p-order growth model with proportional parameters, i.e., α t = γβ t .Equivalently, we assume that there exists a sequence {τ t } of positive real numbers such that and {X t } are independent and identically G p (0, α, β ) distributed.If σ 2 = Var(X t (θ )), then Z t (θ ) ∼ N(µ t , τ 2 t σ 2 ) under Eq. 9. Examples of choices of τ t are τ t = 1, √ µ t or ρ t+1 , cf.Eq. 7. If τ t = 1, the variance of the increment Z t (θ ) is constant in time.If τ t = √ µ t , we obviously need that µ t ≥ 0 for all t and we have that Var(Z t (θ )) ∝ E(Z t (θ )) such that the variance of the increment Z t (θ ) is proportional to the average increase in the radius at time t.If τ t = ρ t+1 , then the distribution of the shape of the object defined by the radius vector function is constant in time, i.e., the distribution of does not depend on t.
In the proposition below, we show that under Eq. 9 the shape of Y t is distributed according to a p-order model.
Below, we study examples of different limiting shapes under the model (Eq.9).

Example 3.3 (Constant increment growth)
Let the situation be as in Proposition 3.2 with µ t = µ and τ t = 1 in Eq. 9.The increment processes Z t are thereby independent and identically distributed.It follows from Proposition 3.2 that where ᾱT = T µ 2 α and βT = T µ 2 β .Since ᾱT → ∞ and βT → ∞ for T → ∞, the boundary of the object becomes more circular and smooth as T increases.An example is shown in Fig. 4. The limiting object has circular shape.
Example 3.4 (Wiener growth) Let the situation be as in Proposition 3.2 with µ t arbitrary and τ t = √ µ t in Eq. 9.This special case is called a Wiener growth model since Var(R T ) ∝ E(R T ).If µ t = µ such that ρ T = T µ, the process is called a Wiener process with linear drift.If ρ T = δ T ψ for some δ , ψ > 0, then with parameter H = ψ 2 , which is a discrete analogue of self-similarity, cf.Sato (1999).Notice that If ρ T → ρ < ∞, the limiting object can have any stochastic shape determined by G p (1, αρ, β ρ).
Example 3.5 Let the situation be as in Proposition 3.2 with µ t arbitrary and τ t = ρ t+1 in Eq. 9.The normalized radius vector function is distributed as If ρ 2 T / ∑ T t=1 ρ 2 t → 0 as T → ∞, the objects become more irregular both globally and locally as T increases.An example is shown in Fig. 4.

AN APPLICATION
For illustrative purposes, we consider a data set consisting of human breast cancer cell islands, which have been observed in vitro in a nutrient medium on a flat dish.This data set has earlier been analysed in Cressie and Hulting (1992).Three profiles of cancer cell islands are available.The data set is presented in the upper left corner of Fig. 5.  Fig. 5.The tumour growth data (upper left corner) and simulations under the second-order growth model with µ t , α t and β t replaced by the maximum likelihood estimates.
The centre of mass of Y 0 is used as reference point.The data consist of increments in n t directions, equidistant in angle, t = 0, 1.For convenience, z t is normalized with the average radius of Y 0 .Only digitized images are available.As n t , we have used approximately 25% of the number of pixels on the boundary of the digitized image of Y t , t = 0, 1.
Under the p-order growth model, the mean value parameters µ t can be estimated by the average observed increment at time t.The variance parameters can be estimated using the likelihood function where L t (α t , β t ) is the likelihood function based on the Fourier coefficients the likelihood becomes where c t,k = [a 2 t,k + b 2 t,k ]/2 are the observed phase amplitudes.In applications, a t,k and b t,k are replaced by discrete versions of the integrals in Eq. 4.
The choice of the cut-off value K t is very important.Clearly, K t must not be too large in order to avoid that the estimates are influenced by the digitization effects.On the other hand, if the cutoff value K t is too small information about the growth pattern is lost.The choice of K t should be an intermediate value for which the estimate of the local parameter β t is stable.Whether a specific choice of K t is appropriate can also be judged from visual inspection of simulated growth patterns under the estimated model.

The estimated regression curves
are shown in Fig. 6, together with 95% confidence limits for the logarithm of the phase amplitudes.The model fits the data well which can also be seen from the fractile diagrams (QQ plots) for the normalized Fourier coefficients, also shown in Fig. 6.
Simulations under the second-order growth model with µ t , α t and β t replaced by the maximum likelihood estimates are shown in Fig. 5.
Since the data set only contains two increments, it is not meaningful to try to evaluate the Markov assumption.Note also that the Z t s are assumed independent but not necessarily identically distributed.
Fig. 6.The two upper figures show the observed phase amplitudes (full-drawn lines) together with the estimated regression curves (stippled) and 95% confidence limits, t = 0, 1.The two lower figures show fractile diagrams (QQ plots) for the normalized Fourier coefficients a t,k / λt,k , b t,k / λt,k together with 95% confidence limits, t = 0, 1.
we have that β t (Z t − µ t ) ∼ G p (0, γ, 1) are independent and identically distributed.Thus, under the assumption (Eq.12) of proportionality and with sufficient number T of time points, we can examine the independence of for selected values of θ ∈ [0, 2π), using a runs test, for instance.

A TIME SERIES EXTENSION
Let us suppose that where X = {X t } is a stationary time series of cyclic Gaussian processes satisfying the ARMA model equation We assume that W = {W t } is a sequence of i.i.d.stationary cyclic Gaussian processes on [0, 2π) with If φ i = 0, i = 1, . . ., r, and ψ j = 0, j = 1, . . ., s, Z follows the p-order growth model with independent increments, treated in the previous sections.
Under the general ARMA model (Eq.13), the Fourier coefficients of X and W of a given order follow a one-dimensional ARMA model.Furthermore, for fixed θ ∈ [0, 2π), X t (θ ) follows a one-dimensional ARMA model.Aspects of this time series approach has earlier been discussed in Alt (1999).An early example concerning year ring widths is discussed in Kronborg (1981).
Note that in the special case of a MA model (φ 1 = • • • = φ r = 0), the marginal distribution of Z t belongs to the class of p-order models .
Note also that in this case Z t and Z t are independent if |t − t | > s.

EXTENSION TO THREE DIMENSIONS
The p-order growth model for planar objects can easily be extended to three dimensions.Consider a spatial bounded and topologically closed object Y t ⊂ R 3 which is star-shaped for all t with respect to z ∈ R 3 .Clearly the boundary of the object can be determined by In the same way as in the planar case we let the object Y t+1 be a stochastic transformation of the object Y t , such that where {Z t } is a time series of Gaussian procesess on [0, 2π) × [0, π].Writing the stochastic process Z t in terms of its Fourier-Legendre series expansion we get, cf.Hobolth (2003), where φ n,m are the spherical harmonics and A t,n,m are random coefficients.Using a similar reasoning as in Hobolth (2003) it can be seen that A t,0,0 determines the overall growth from Y t to Y t+1 .The coefficients A t,1,m , m = −1, 0, 1, control the asymmetry of growth, and the remaining coefficients A t,n,m for n ≥ 2, m = −n, . . ., n, affect how the growth appears globally for small n and locally for large n.A p-order growth model can be defined by assuming that A t,0,0 = µ t , A t,1,m = 0 for m = −1, 0, 1 and A t,n,m ∼ N(0, λ t,n ) , n = 2, 3, . . ., m = −n, . . ., n, independent, where λ −1 t,n = α t + β t (n 2p − 2 2p ) .
As in the planar case, the increment processes may be chosen to be normal after a transformation.A simulation from such a model, where {Z t } is a series of log-Gaussian processes, is shown in Fig. 6.

DISCUSSION
The p−order growth model has mainly been suggested as a general tool for analyzing observed radial growth patterns.The model may, however, also be of interest as a building block in other modelling situations, for instance in models for tessellations where cells are created by radial growth from each point of a point process.
The p−order growth model can be extended in various ways.It is obviously easy to modify the model such that the increments are Gaussian after a transformation.An example is log-Gaussian increments.If the number of increments observed is not too small it is also of interest to try to model the dependency in the series Z = {Z t }.We have discussed a time series approach.Another alternative is to look at Lévy based models, Z t (θ ) = A t (θ ) h t (a; θ )Z(da) , A t (θ ) ∈ B, where B is the Borel field of [0, 2π) × R and Z is a Lévy basis on [0, 2π) × R. A detailed study of the Lévy based growth models is ongoing research in our group, cf.Schmiegel et al. (in preparation).These models can also be formulated in continuous time.
The likelihood used in the application is correct if the increments are independent.If the marginal distributions of the Z t s belong to the class of p-order models but the Z t s are dependent, the likelihood may still be used as a pseudo-likelihood.
In relation to tumour growth in particular, it will also be of interest in the future to try to embed specific mathematical models in a stochastic framework.A starting point could here be a study of dynamic point process models with a specified time-dependent intensity function.

Fig. 3 .
Fig. 3. Top: Simulated growth pattern under the constant increment second-order growth model.Bottom: The corresponding normalized profiles, representing the shape of the object.

Fig. 4 .
Fig. 4. Top: Simulated growth pattern under the model described in Example 3.5.Bottom: The corresponding normalized profiles, representing the shape of the object.