PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Math Biosci. Author manuscript; available in PMC 2010 May 28.
Published in final edited form as:
PMCID: PMC2878366
NIHMSID: NIHMS201417

Biological Growth on a Surface

Abstract

A new biological growth model is introduced. After random selection of possible growth sites, cells are added in an allowable random direction. Unlike the Eden model, the current model is grid independent and can be adapted to any curved surface. Growths of colonies of as many as 105 cells are simulated on planar, cylindrical, and spherical surfaces. It is found that the interior density is constant, whereas the boundary is fractal.

1. INTRODUCTION

In 1961, Eden [1] introduced a simple but important stochastic growth model that mimics biological growth such as that of a tumor. The algorithm is essentially as follows. A cell of a planar square grid is labeled “infected.” Then any one of the four possible adjacent cells is randomly chosen to be infected. The pair now has six possible growth surfaces. The process continues until a cluster is formed. Many papers have been written on the Eden model [2]. It was found that the interior core is solid [3] and that the boundary is fractal [4].

Tumors or diseases also may be found on curved surfaces such as those on the skins of organs and fruits. In these cases, the Eden model is difficult to apply because a suitable grid cannot be found. On the other hand, the Eden model is not isotropic; that is, the growth tends to lengthen along the axes of the grid [5]. The purpose of the present paper is to introduce a new random branching growth model that is equivalent to the Eden model but without the aforementioned disadvantages.

1.1 TILING A CURVED SURFACE

The square grid of the Eden model constitutes the boundaries of the square cells, which “tiles” the surface. There are many other ways to tile a planar surface with similar shapes, such as triangles or rectangles [6]. Considering minimizing the boundary length, which is true for most biological units, the best tile for a plane is probably the regular hexagon, a shape closely related to the equilateral triangle.

A cylindrical surface is developable; that is, it can be cut along the axis and spread fiat. One must keep in mind, of course, that the cut edges are actually connected. Thus the planar Eden model can be adapted to a cylindrical surface if the circumference contains an integer number of cells [4].

Another basic geometry is the cone. Although the conical surface is developable, there is almost no way to tile the cone successfully. One can construct a rhombic grid by cutting along a ray, developing the surface, and using lines parallel to the edges as grid lines. There are, however, distinct changes in orientation along the ray when the edges are rejoined.

The problem is much more complicated for the basic shape, the sphere. The spherical surface is not developable and any tiling must be done in three dimensions. Considering our requirements of similar shapes of equal area and minimal boundary lengths, the tiling would be a projection of a regular polyhedron onto a circumscribing sphere. The best tiling is the icosahedron, the platonic polyhedron with most surfaces [7-9] (Figure 1a). The icosahedron (composed of 20 equilateral triangles) still has too few cells for modeling a tumor spread. The n-frequency icosahedron almost satisfies all our requirements. Each basic triangle of the icosahedron is partitioned into smaller, identical triangles. There can be as fine a subdivision as needed by increasing the frequency (number of cells = 20n2). The small triangles are then projected onto the circumsphere of the icosahedron (Figure 1b). Such a tiling is still not perfect, because owing to projection, the small triangles are not the same size, nor they are equilateral. Those near the middle of the basic icosahedron triangles are slightly larger. Also the relations between neighboring cells are not the same everywhere. Note that all nodes are at the confluence of six cells except at the 12 icosahedron vertices where there are only five cells.

Fig. 1
(a) The icosahedron. (b) The three-frequency icosahedron projected onto a sphere.

We conclude that, except for the cylinder, a homogeneous grid on a curved surface does not exist.

1.2 DIRECTIONAL PREFERENCES OF A GRID

Another problem of the Eden model is its preference for growth in directions along the grid axes. This anisotropy is inherent in any algorithm based on a grid or lattice. Although the Eden model looks somewhat circular at first, when the cluster is large, the anisotropy becomes obvious because the shape is distorted into a diamond and, for much larger clusters, the shape becomes a cross. The reason is due to the fact that the distances to the layer of cells diagonal to the axes are larger by a factor of 2 than distances to the layers along the axes [2, 5]. If triangular or hexagonal tilings were used, the asymptotic shape would be like an asterisk. Thus the Eden model could not have described a natural biological growth that is isotropic or grid-independent. To alleviate the grid dependency in one direction, the Eden model has been applied to a strip equivalent to a developed cylindrical surface, starting from a line of infected cells across the strip [2, 4]. Even with this adjustment, the growth is still biased axially, and it is uncertain whether the results would reflect the properties of even a small boundary segment of a large isotropic natural cluster.

This grid dependency is also found in diffusion limited–aggregation (DLA) models, which are based on a process entirely different from that of the Eden model. Instead of growth by infecting neighboring cells as in the Eden model, the DLA cluster enlarges by capturing successive free cells that are in Brownian motion. The cluster shape is lichenlike without a discernable core. Meakin [10] used both a square grid and a gridless (off-lattice) algorithm and found that the density-density correlations are indeed different. A similar study for the Eden model does not exist.

Because the Eden model is based on a grid that is anisotropic and because any grid on a curved surface (except cylindrical) is furthermore nonhomogeneous, we need to construct an algorithm that is gridless but retains the essence of the Eden growth. Such a model is described in the next section.

2. THE RANDOM GROWTH ALGORITHM

The lattice-free growth model is described as follows.

  1. Each cell is represented by a node (point) on the surface. Each node excludes the growth of any new nodes within a normalized distance of 1 from that node.
  2. The nodes can be either “living” or “dead.” Only living nodes are capable of growing a new node at a distance of 1.
  3. A node is selected randomly from the existing living ones for possible growth. Its neighbors within a distance of 2 are identified. Each neighbor excludes a certain range of growth directions of the chosen living node.
  4. A random direction from the chosen living node is tried. If this direction is not excluded, a new (living) node is grown. Go to (3). If this direction is excluded, another random direction is tried. The angles of these directions may be related by multiples of 2π/K, where K is the number of directions tested. The process is like randomly choosing each direction of the spokes of a K-spoked stopped wheel.
  5. If all K directions are excluded, the node is declared “dead.” Go to (3).
  6. The growth process stops until a given number of nodes are grown or a given domain is filled.

The result is a distribution of nodes on the surface. To visualize the cells, one can surround each node by a Voronoi polygon (enclosing a Dirichlet domain). The polygon contains the set of all points that is closest to the node inside the polygon. An example is shown in Figure 2, where the boundary nodes are enclosed by partial circles. We shall not be concerned with the properties and the computer construction of the Voronoi polygons, which are described elsewhere [11].

Fig. 2
Nodes surrounded by Voronoi polygons and partial circles.

3. PROPERTIES OF PLANAR GROWTH

To determine the growth properties of the model, large-scale simulations (resulting in clusters of as many as 105 nodes) are performed. A small colony of 1000 nodes (with K = 12 directional trials) is shown in Figure 3. As in the Eden model, the interior density is relatively constant. This can be shown by filling a circular region of large diameter and counting the nodes enclosed within different-sized concentric circles. A typical result is shown in Figure 4. We see that local inhomogeneities become insignificant for regions of dimensions larger than about 40. The variation of density ρ with the number of trials K is shown in Figure 5. The density increases slightly as K is increased, reaching an asymptote of ρ = 0.8276 for large K, where each selected living node must grow a new node unless all directions are entirely excluded by its neighbors.

Fig. 3
A colony of 1000 nodes, K = 12.
Fig. 4
Density versus diameter of domain, K = 12.
Fig. 5
Density as function of K.

Now let us consider the boundary of the colony consisting of the living nodes capable of further growth. (The boundary of Figure 3 is shown in Figure 6.) Let the locations of the boundary nodes be at (xn, yn), n = 1, N. The center of mass is then

(x,y)=1N1N(xn,yn).
(1)

The mean radius is

R=1N1N[(xnx)2+(yny)2]12.
(2)

We define the thickness σ of the boundary by

σ2=1N1N{[(xnx)2+(yny)2]12R}2.
(3)

Using samples of clusters of as many as 105 nodes, we find that the number of living nodes N is proportional to the mean radius R, which in turn is proportional to the square root of the total number of nodes. The proportionality constant is slightly dependent on the number of trials K (Figure 7). Because the density is constant, the number of living nodes is proportional to the mean radius R. On the other hand, the thickness of the boundary is fractal:

σ~Rβ.
(4)

Using 90 samples and a confidence level of 95%, we find for K = 6, β = 0.167±0.047; for K = 12, β = 0.410±0.054; for K = 24, β = 0.396±0.051, which is not significantly different from that of K = 12. Because β < 1, larger colonies have smaller relative boundary thicknesses.

Fig. 6
The boundary (living) nodes of Figure 3.
Fig. 7
The number of living nodes N as a function of total number of nodes M for planar growth.

Because of the fact that our algorithm is grid independent and thus istropic, the shape of the boundary is guaranteed to be near circular.

Let M be the total number of nodes and N be the number of boundary nodes. Define one unit of time as the time for growing a new node. Of course, new nodes are grown concurrently from all living nodes instead of consecutively as in our algorithm. Thus

dMdt=N.
(5)

Because density is constant, the rate of increase in infected area is proportional to the number of living nodes. For planar growth from a single seed,

M=πρR2.
(6)

Thus from Equation (6), the rate of increase of the mean radius is

dRdt=12πρRdMdt=N2πρR.
(7)

Because N ~ R, we find that the rate of increase of the mean radius is constant, or the size (diameter) of the growth increases with constant speed. Figure 8 shows a typical growth of the boundary of a colony at equal time intervals.

Fig. 8
Boundaries of planar growth at equal intervals (M = 2500. 10,000, 22,500).

4. GROWTH ON A CYLINDRICAL SURFACE

As we mentioned before, any cylinder (not necessarily circular) can be opened and spread flat after a longitudinal cut. Let the boundary of the spread-out strip be at x = ± xm. The cut edges are actually connected. The algorithm remains the same as in the planar case except for the following modifications.

  1. If a new growth is at (x1, y1) and |x1| < xm -2, the growth will be oblivious to the cut boundaries and (x1, y1) is added to the living pool.
  2. If xm -2 < |x1| < xm, the point (x1, y1) becomes a living node and a dead node is added just outside the other boundary at (x1 - 2xmsign(x1), y1).
  3. If |x1| > xm, the point (x1, y1) becomes a dead node with a living node added at (x1 - 2xmsign(x1), y1).

Figure 9 shows the growth on a cylinder from an infected node. At small times, the boundary is simply connected and almost circular. For large times, the boundary becomes two rings encircling the cylinder. Figure 10 shows the number of boundary (living) nodes N versus the total number of nodes M for a typical case. With fixed cylinder circumference (2xm), N increases linearly with M as in Figure 7. As the boundary is reached, N decreases and then tends to a constant for large M. A more consolidated universal graph, plotting (N/xm) versus (Mxm) is shown in Figure 11. Theoretically, the maximum value of N is (π/2) times its asymptotic value for large M. Because the rate of increase of total area is proportional to the number of living nodes, we find that the fastest spread is when the growth just begins to link itself by wrapping around the cylinder.

Fig. 9
Boundaries of growth on a cylinder with surface developed to a strip (M = 2500, 10,000, 16,000).
Fig. 10
Number of living nodes as a function of total number of nodes for various different cylinder circumference xm (K = 12) [big up triangle, open], ○, □, [diamond with plus]: x = 50, 100, 150, 200.
Fig. 11
Consolidated graph of Figure 10, using normalized parameters. Dashed line is asymptotic value for large M(N/xm = 8.116).

For large M, the statistical properties of the boundary of the growth on a cylinder become independent of its initial development. The situation is similar to the growth on a strip studied by Jullien and Botet [4] for the Eden model. For different xm (10, 20, 50, 100, 200, 500), we simulated the growth from a string of living nodes starting from one end of the strip. The process is stopped when the height of the growth is about 40 times its width (as many as 2 × 105 nodes). Using 60 samples, we find for K = 12, the number of living nodes is proportional to xm; that is, N = 4.058xm. Considering that there are two propagating ring boundaries, N/xm = 8.116, which is the asymptotic line plotted in Figure 11. The mean thickness of the boundaries is fractal:

σ=Axmβ,
(8)

where for K = 12, A = 1.0714 and β = 0.1185±0.0314. For K = 6, the value of A is 1.671 and β = 0.0473±0.022. In contrast, the fractal dimension for the Eden model is higher; depending on the source, β ranges from 0.25 to 0.5 [4].

5. GROWTH ON A SPHERICAL SURFACE

The spherical surface is not developable. In contrast with the cylindrical or planar surfaces, a homogeneous grid does not exist. Thus any growth simulation on a sphere must be gridless. Our random growth algorithm can be adapted to this curved surface as follows.

Consider a sphere of radius 1. Any Cartesian point (x, y, z) can be parametrized by spherical coordinates (ϕ, θ) by using

x=sinϕcosθ,y=sinϕcosθ,z=cosϕϕin[0,π],θin(π,π].
(9)

Because the radius is 1, the arc length from the north pole to any point Pi(ϕi, θi) is ϕi (Figure 12a). The sides of the spherical triangle linking the points P0, P1 and the north pole are great circles (shortest distances). Let d be the (great circle) distance between P0 and P1. From spherical trigonometry [12],

cosd=cosϕ0cosϕ1+sinϕ0sinϕ1cos(θ1θ0).
(10)

Let the distance of each growth be h(1). As before, a living node P0 is selected at random and its neighboring nodes (of distance < 2h) are identified by using Equation (10). Compute the half exclusion angle w from the spherical triangle in Figure 12b:

w=cos1[1cosdsindtanh].
(11)

The angle between P0 P1 and the vertical is α:

α=π+cos1[cosϕ1cosdcosϕ0sindsinϕ0]sign(θ0θ1).
(12)

Thus the range of angles excluded from growth is [αw, α + w]. As in the planar case, the directions of the K spokes of a stopped wheel are tested. New growth commences from P0 if some spike direction is not excluded. Figure 13 shows the boundaries of a growth corresponding to approximately same-time intervals. The longitudes and latitudes are drawn in for reference only.

Fig. 12
(a) The spherical system. (b) Exclusion angle w.
Fig. 13
Growth on a sphere (h = 0.01, K = 12, M = 2500, 10,000, 22,500).

The spherical surface is special because its area is finite and can be filled completely with nodes in finite time. Thus there is no “asymptotic state,” as in the infinite planar and cylindrical surfaces. Therefore the fractal properties of the boundary are not studied in the spherical case.

We simulated the growth on a sphere (K = 12) for different values of h. For h = 0.01, 0.015, and 0.02, there are 99,781, 44,402, and 24,493 nodes, respectively, to completely cover the sphere. For comparison with the planar case, we set h = 1 and R = 1/h. The densities (number per unit square) are 0.7940, 0.7950, and 0.7956—very close to the 0.805 of the planar case. Because h is small, the curvature of the sphere has negligible effect on the density based on h2.

Of interest is the number of living nodes N as a function of the total number of nodes M. Guided by our discussion on density, we found that a plot of Nh versus Mh2 would be universal—that is, independent of h (Figure 14). The number of living nodes reaches a maximum (near the equator of the sphere) and then decreases to zero as the sphere is completely covered. In comparison, for large colonies, the number of living nodes increases as M for planar growth and tends to a constant for the growth on a cylinder. The symmetry of the data about (Mh2)max/2 suggests reversibility of growing outward and growing inward. Because the rate of increase of total nodes is proportional to the number of living nodes [Eq. (5)], the fastest area increase is when the growth covers about half the sphere.

Fig. 14
Number of living nodes Nh versus Mh2. □, [big up triangle, open], ○: h = 0.01, 0.015, 0.02.

If we assume that the boundary is circular and smooth, the area of a spherical cap of growth is

2π(1cosϕ)=Mρ,
(13)

because in unit time,

N=dMdt=2πρsinϕdϕdt.
(14)

Again, for a smooth boundary, dϕ/dt can be supplanted by h. From Equations (13) and (14), we obtain a theoretical formula:

(Nh)=2π(ρh2){1[1(Mh2)2π(ρh2)]2}12.
(15)

Now (ph2) is a constant approximately 0.795, and Equation (15) shows the qualitative behavior of Figure 14, except that maximum Nh is underestimated to be 5.0 instead of the actual 6.5. This means that the boundary is not smooth, as assumed, but fractal, as in the planar case.

6. DISCUSSIONS AND CONCLUSIONS

It is clear that the grid-based Eden model could not have succesfully simulated boiogical growth. Not only it is inapplicable to curved surfaces such as a sphere, but more importantly its results are very much in doubt owing to the inherent anisotropy of a grid. In this paper, we have introduced a gridless random growth model that does not have the deficiencies of the Eden model.

Using our model, we can simulate the spread of disease on any curved surface. We find that the interior has constant density, whereas the boundary is fractal. The thickness of the boundary has fractal dimension of 0.4 for planar growth from a point and 0.12 for the growth along a strip. The latter value is much smaller than those predicted by the Eden model.

In our algorithm, there is a need to identify neighbors near a growing node. We grouped all nodes spatially into labeled bins and thus very easily retrieved those containing possible neighbors. This method greatly facilitated the search.

The parameter K in the present model represents K directions tested for possible growth of a living node. We find K = 12 is adequate because the results are within 2% of the ideal of K = ∞.

The number of living nodes at each stage is related to the rate of increase of infected area. Our Figure 7 for planar growth and our universal plots in Figures 11 and and1414 for growth on cylinders and spheres should be highly useful in the prediction and control of biological growth.

Acknowledgments

The research was partially supported by NIH Grant No. RR01243.

Contributor Information

C. Y. WANG, Departments of Mathematics and Physiology, Michigan State University, East Lansing, Michigan.

J. B. BASSINGTHWAIGHTE, Center for Bioengineering, University of Washington, Seattle, Washington.

REFERENCES

1. Eden M. Proc. Fourth Berkeley Symp. Math. Stat. Prob. 1961;4:233.
2. Vicsek T. Fractal Growth Phenomena. World Scientific; Singapore: 1989. p. 182.
3. Richardson D. Proc. Camb. Philos. Soc. 1973;74:515.
4. Jullien R, Botet R. J. Phys. 1985;A18:2279.
5. Dhar D. In: On Growth and Form. Stanley HE, Ostrowsky N, editors. Nijhoff; Dordrecht: 1986. p. 288.
6. Grunbaum B, Shephard GC. Tilings and Patterns. Freeman; New York: 1987.
7. Wenninger MJ. Spherical Models. Cambridge, England: 1979.
8. Pugh A. Polyhedra: A Visual Approach. Berkeley, CA: 1976.
9. Gasson PC. Geometry of Spatial Forms. Wiley; New York: 1983.
10. Meakin P. In: Fractals in Physics. Pietronero L, Tosatti E, editors. Elsevier; Amsterdam: 1986. p. 213.
11. Preparata FP, Shamos MI. Computational Geometry. Springer; New York: 1985.
12. Sperry P. Short Course in Spherical Trigonometry. Johnson; Richmond, VA: 1928.