Growth dynamics of cancer cell colonies and their comparison with noncancerous cells.

The two-dimensional (2D) growth dynamics of HeLa (cervix cancer) cell colonies was studied following both their growth front and the pattern morphology evolutions utilizing large population colonies exhibiting linearly and radially spreading fronts. In both cases, the colony profile fractal dimension was d(f)=1.20±0.05 and the growth fronts displaced at the constant velocity 0.90±0.05 μm min(-1). Colonies showed changes in both cell morphology and average size. As time increased, the formation of large cells at the colony front was observed. Accordingly, the heterogeneity of the colony increased and local driving forces that set in began to influence the dynamics of the colony front. The dynamic scaling analysis of rough colony fronts resulted in a roughness exponent α = 0.50±0.05, a growth exponent β = 0.32±0.04, and a dynamic exponent z=1.5±0.2. The validity of this set of scaling exponents extended from a lower cutoff l(c)≈60 μm upward, and the exponents agreed with those predicted by the standard Kardar-Parisi-Zhang continuous equation. HeLa data were compared with those previously reported for Vero cell colonies. The value of d(f) and the Kardar-Parisi-Zhang-type 2D front growth dynamics were similar for colonies of both cell lines. This indicates that the cell colony growth dynamics is independent of the genetic background and the tumorigenic nature of the cells. However, one can distinguish some differences between both cell lines during the growth of colonies that may result from specific cooperative effects and the nature of each biosystem.


I. INTRODUCTION
The study of cell colony dynamics has received increasing attention due to its role in comprehending complex biological processes such as tumor growth. To advance the knowledge of these processes, complemented research work from both biochemical and biomechanical standpoints became particularly attractive [1][2][3]. This multidisciplinary problem has been approached in different ways, for instance by extending the conceptual framework of statistical and fractal analyses to study the evolution of rough interfaces in living systems, as was done to investigate the dynamics of rough interfacial profiles of different biosystems including bacteria [4] and cell colonies [5][6][7]. Concepts such as fractality and dynamic scaling analysis then became of increasing interest in biology and medicine [8,9] to establish, for instance, a possible relation between the fractal dimension of tumor contours and their degree of malignancy [10].
The roughness dynamics of condensed phases can be characterized by a set of scaling exponents (α, β, and z = α/β, which denote the roughness, the kinetic, and the dynamic exponent, respectively), derived from the scaling analysis of profile data. In the past decade, a number of publications have shown that the dynamic behavior of systems of a very different nature obeys common scaling relations that can be linked to a certain universality class [11], at least over a certain spatiotemporal range. For a set of coherent scaling exponents, whose validity extends over a relatively wide spatiotemporal range, it is possible to interpret the growth of condensed phases by either a single or various complex mechanisms [12]. Accordingly, the scaling relations of multiconstituent systems have been mimicked either by discrete models or continuous differential equations involving a relatively small number of essential relationships [12]. In addition, this knowledge of the overall growth process makes it feasible to draw out information about multiple interactions among phase constituents as well as the participation of cooperative effects [13]. As recently reported [6,7], for two-dimensional (2D) growth of Vero (African Green Monkey kidney) cell colonies, the influence of those effects, which depend on the colony population, is reflected in the age-dependent heterogeneity of biosystems produced by changes in both the cell shape and size distribution in the colony. Accordingly, the evolution of biosystems becomes much more complicated than that of inanimate material systems because of their larger number of specific variables. These variables are, in general, more difficult to control and to study separately than those involved in inanimate materials.
The scaling analysis of tumor growth fronts, as well as those from different cell line colony 2D profiles, has been used to determine the universality class involved in their growth processes [5], and scaling data were interpreted by the molecular beam epitaxy (MBE) model. This conclusion was revisited to assure whether the validity of extending the mathematical formalism for fixed size colony growth experiments was directly applicable to radially expanding systems [14][15][16]. However, in contrast with the MBE model, scaling exponents derived from a cellular automaton model for 2D colony growth dynamics appeared to be much likely to be represented by the Kardar, Parisi, and Zhang (KPZ) -type continuous equation [17].
A first attempt to clarify the above situation was made utilizing experimental designs developed for the growth of cell colonies starting from either quasilinear [6] or quasicircular [7] spreading 2D fronts using Vero cells, a nontumorigenic cell line that exhibits a null or almost null contact inhibition and continues growing and dividing indefinitely, which is a typical property of cancer cells [18]. Data from those experiments yielded the scaling exponents α = 0.50 and β = 0.32 for width fronts above 100 μm [6,7]. These exponents fulfilled the expectation of both the KPZ continuous equation and the Family-Vicsek relation [19]. Those experiments also showed that the morphology of growth patterns became more complex than the assumptions involved in deducing both the standard KPZ equation [20] and the automaton model [17]. However, despite the specific characteristics of the biosystem, the above set of scaling exponents essentially captured the biosystem scaling behavior. Seemingly, a KPZ-like dynamics was accomplished irrespective of the fact that linearly 2D fronts were constrained to advance perpendicular to the initial front, whereas for radial fronts a progressive expansion of the colony front was involved. The comparison between both growth geometries also showed that radially growing experiments became useful to distinguish the exponential to linearly front displacement transition along the colony growth [7].
As the preceding conclusions were derived from a noncancerous cell line, it appeared important to extend the approach to cancer cells to attempt to determine whether the growth dynamics of cell colonies could be specifically related to the tumorigenic nature of the cells or to their genetic background. For instance, it was recently reported that nontumorigenic and cancer cells exhibited differential morphology and motility responses to changes in substrate rigidity and microtopography [21]. For this purpose, the HeLa cell line, one of the oldest human cell lines that is derived from cervix cancer [22], was selected. These cells have been used for research into cancer, AIDS, the effect of radiation and toxic substances, gene mapping, as well as many other scientific pursuits [23][24][25].
In this work, we report on the evolution of large population linearly and radially growing HeLa cell colonies. In general, the morphology evolution of colony patterns shows small changes in the average size of irregular cells as the colony age increases, and their distribution becomes rather homogeneous without the formation of well-ordered domains. However, at advanced growth stages, some large, multinucleated cells located at the outermost colony regions and 3D cell domains at the innermost ones are formed. The dynamic scaling analysis shows that the growth dynamics can be interpreted in terms of a KPZ-like continuous equation. The dynamic scaling for both colony geometries renders the scaling exponents in the width range from a lower cutoff of about 60 μm upward. This cutoff is related to the average cell size at the colony front and determines a characteristic width at which the influence of overhangs on the colony front becomes negligible, although the influence of the latter became slightly larger for quasicircular colonies. Comparing dynamic data obtained for HeLa cells with those reported for Vero cells, one can conclude that the 2D growth of cell line colonies in a standard culture medium can be explained by the KPZ dynamic universality class.

A. Cell colony growth procedures
HeLa cells (passage 44) were cultured using a Roswell Park Memorial Institute (RPMI 1 640) medium containing 10% fetal bovine serum (FBS) maintained in a 5% carbon dioxide controlled atmosphere at 37 • C, changing one-half of the culture medium every 2 days.
Colonies with linearly growth fronts were obtained by first covering the central region of the Petri dish bottom with a 2.2cm-wide and 100-μm-thick sterilized Teflon tape [6]. Then, disaggregated HeLa cells (30 000-40 000 cell mL −1 ) were shed and left to grow for about 2 days until confluence in the Teflon-free region was reached. Eventually, the Teflon tape was removed, leaving a cell-free central region with the formation of two facing linear colony fronts of width L that started to shift in opposite directions perpendicular to L. From this stage on, the morphology evolution of the colony pattern and the 2D growth front displacement were followed.
Colonies exhibiting radially spreading growth fronts were prepared by shedding disaggregated cells in the Petri dish and were left growing until a 3D cell cluster with an average radius of about 250-300 μm at the colony center was formed. Subsequently, the 3D cluster was carefully removed with a micropipette and transferred to a second Petri dish containing fresh culture medium. Eventually, the colony continued spreading from the 3D seed rim as a 2D domain.
Each colony growth was followed until it contacted the border of a neighbor colony, a fact that commonly occurred after 5-8 days. Most of the neighbor colonies were presumably formed by occasional cell detachment from 3D colony domains, as cell-cell adherence energy became weaker than cell-substrate interaction energy [26][27][28][29][30].
Colony fixation and staining with May-Grünwald Giemsa was occasionally performed to improve the detection of large multinuclear cells and cell filopodia, and to evaluate the cell size statistical distribution at 2D growth patterns. The viability of the cells was routinely checked using the exclusion Tripan-Blue test. Cell duplication was determined by labeling with proliferating cell nuclear antigen (PCNA).

B. Colony imaging and data processing
Colonies were imaged using a Canon digital camera coupled to a Nikon TS100 phase-contrast inverted microscope with a CFI flat field ADL 10× objective at a resolution of 0.88 μm/pixel, where the field covered by the microscope was about 1000 μm. As the colony grew, an increasing number of partial images were taken to cover the whole colony . Subsequently, these images were stitched to compose the complete image of the colony. Colony fronts were manually traced using a Wacom graphic tablet with a tracing error on the order of the pixel. The typical evolution of colony fronts from 2D radially and linearly spreading colonies is depicted in Fig. 1.
The analysis of colony fronts was performed employing an in-lab developed program to obtain their fractal dimension, the front roughness, and other parameters related to the front displacement, with each front consisting of N points (pixels). For the sake of simplicity, both the radial and the linear front distances are denoted as h. For radially spreading colonies, the program allowed us to determine the center of mass (c.m.), the instantaneous distance h i (t) from the c.m. for the ith point at the front (i = 1,2, . . . ,N), and the corresponding arc (s i ), whereas for linearly spreading colonies, one obtained the instantaneous colony height h i (t) from the starting linear front of width L for each ith front site. The mean colony distance at time t and the mean front displacement velocity were calculated as hh(t)i = P h i /N and hvi = dhh(t)i/dt, respectively.
The instantaneous colony front roughness was determined from the standard deviation of the front height fluctuations, where L is the width of the growing front, which is constant for linear fronts and increases with time as L = 2π hh(t)i for radially expanding fronts [7]. The dependence of the local roughness w(l,t) on the front width l (l 6 L) was evaluated from the standard deviation of h i for different front lengths l at time t. While this procedure is straightforwardly applied for linearly growth fronts, radially spreading colonies require the evaluation of the roughness at different arc lengths s (s = l < L) [5,7,31]. For this purpose, the coordinate system was transformed from angle radius into arc radius. Then, the location of each point at the front was determined by coordinates h i ,s i , with the arc s i being measured along a circle of radius hh(t)i [7].
Global and local roughnesses were evaluated for both experimental and overhang-corrected data. Overhang-corrected colony profiles were obtained by taking the maximum value of h i (t) at each ith position at the front. As shown below, the influence of overhang correction, required for a proper scaling of the front roughness, becomes practically null when the length of the colony front exceeds about 2-3 average cell diameters (20 μm). rather compact and exhibit a fairly homogeneous distribution of cells with an average 20 μm cell diameter at the front, although there are not distinguishable ordered domains. The 2D colony front can be described as a tessellation of cells with different irregular shapes (Fig. 2). As time increases, one can observe a relatively small number of large cells, some of them multinucleated, principally located at the outer region of the colony, and the appearance of 3D domains at the innermost colony region (Figs. 3 and 4). According to PCNA labeled growing colonies images, cell proliferation takes place at both the colony front and the bulk (Fig. 5).
On the other hand, 2D radially spreading HeLa colonies started from a 3D seed and continued growing as a 2D phase around the border of the seed cluster (Fig. 4). For both linearly and radially spreading colonies, the cell size at the colony front appears rather homogeneously distributed around a maximum at 1000 μm 2 (Fig. 6), i.e., a mean cell radius of about 20 μm. On the other hand, as the colony size increases, the average cell  area decreases from the colony front inward (Fig. 7). At the longest time, the area of cells located at the innermost part of the colony becomes about 25% lower than those of the initial colony front. In this case, the cell area at the colony border can exceed 2000 μm 2 .
The colony front mean displacement, hhi, increases linearly with time (Fig. 8), irrespective of the colony geometry. For radially expanding fronts, the normalized value of hhi, i.e., hhi − R 0 , where R 0 is the mean radius of the initial 3D seed cluster, converges to the same initial point [ Fig. 8(b)]. The average colony front displacement velocity derived from these plots results in (0.09 ± 0.02) × 10 −4 μm min −1 , a figure that is about one-half the value reported earlier for Vero cell colonies under comparable conditions [6,7].

B. Fractality and roughness exponents from 2D front colony patterns
The fractal dimension of HeLa colony fronts, evaluated by the box-counting method, resulted in d f = 1.20 ± 0.05 (Fig. 9), a figure that agrees with that obtained for Vero cell colonies under both quasilinear [6] and radial [5,7] growth conditions. These values of d f indicate complex colony contours with a low degree of ramification. In addition, plots resulting from radially spreading colonies shift upward parallel with time due to the increase of the colony front width [ Fig. 9(b)].
According to the dynamic scaling analysis [12], the interface roughness of a 2D condensed system of width L is expected to increase with time for t ¿ t s , where t s denotes the roughness saturation time, as follows: and when t > t s , the roughness saturation w s is attained. The latter should increase with L according to where the value of t s depends on the system width (t s ∝ L z , z = α/β) [12]. Therefore, from the slopes of the log w(L,t) versus log t and log w(l,t) versus log l plots, one can obtain the growth and roughness exponents β and α, respectively, and from their ratio, the dynamic exponent z. Roughness data of both experimental and overhangcorrected colony fronts, evaluated from Eq. (1), fulfill linear log-log plots with time with the slope β = 0.32 ± 0.04, irrespective of the colony geometry (Fig. 10). This value of β agrees with that previously reported for Vero cell colonies approaching either quasilinear [6] or quasicircular [7] growth fronts.
The evaluation of the local roughness w(l,t) was performed for both the experimental and the overhang-corrected linearly and radially spreading colony fronts. Thus, the w(l,t) versus l (l 6 L) log-log plots collapse into a single curve for l above a certain characteristic length l c , irrespective of the colony geometry (Fig. 11). These plots exhibit a linear portion with the slope approximately 0.5 over about one decade in l, and they tend to attain a limiting value as l → L. Consequently, for  l above l c ≈ 60 μm, the application of scaling relationships to overhang-corrected data becomes feasible.
The scaling exponents derived from the dynamic scaling analysis (Figs. 10 and 11) were utilized to test the fulfillment of the Family-Vicsek relation [32], where f (x = tL −z ) is the scaling function that depends on x as follows: Equation (4) brings about the collapse of scaling data for an appropriate consistent set of exponents α and β. Log-log plots of Eq. (4) for both linearly and radially spreading colonies show collapses approaching a reasonable single universal curve (Fig. 12).
The scaling analysis in the real space was complemented by the Fourier analysis of front profiles [33], where only longwavelength modes contribute to scaling. In contrast, real-space scaling Eq. (4) involves all wavelength modes, including short ones, and consequently stronger finite-size effects are expected [12]. Thus, the structure factor S(k,t) of overhang-corrected growth fronts is [12,17,33,34] S(k,t) = hĥ(k,t)ĥ(−k,t)i, whereĥ is the kth Fourier mode of the profile height fluctuation around its average value at t. Then, one can obtain the corresponding Family-Vicsek relation in terms of S(k,t), with α s being the spectral roughness exponent and f (x = kt 1/z ) the scaling function that depends on x as follows: f (x) = ½ const for x À 1,  Therefore, from the slope of the logS(k,t) versus logk plots, one can obtain the value 2α s + 1. Structure factor data obtained from both linearly and radially spreading colonies approach a straight line with a slope close to −2.0 ± 0.1 (Fig. 13). Then, in the range 0.02 6 k 6 1, it results in α s = 0.50 ± 0.05, irrespective of both t and the colony-spreading geometry.
From the value of α s derived from the slope in Fig. 13, one can then plot Eq. (8) for both linearly and radially spreading colony data (Fig. 14). Both log-log plots exhibit a reasonable collapse into a single curve with a straight line portion of slope 2, as expected for α s = 0.5. Considering the proper noise of biological systems, the above results from the dynamic scaling analysis are fairly good.

IV. DISCUSSION
The growth dynamic data of HeLa cell colonies with a relatively large population can be compared with those reported for similar Vero cell colonies [6,7] to find out whether there are any specific characteristics that could be related to the different tumorigenic nature and genetic background of the cell lines. Accordingly, from such a comparison, we get the following: (i) The colony growth patterns from HeLa cells depict relatively less ordered domains with respect to Vero colonies [6,7].
(ii) The density of relatively large cells in the HeLa colonies is smaller than that reported from Vero colonies.
(iii) HeLa colonies show the formation of more compact colony patterns (Fig. 2) as compared with Vero ones. This fact presumably produces a lower contribution of overhangs at the colony front.
(iv) The constant displacement velocity of HeLa colony growth fronts (Fig. 8) results in almost one-half the value that has been found for Vero colonies [6]. This could be attributed to changes in cell-cell and cell-substrate interactions.
(v) The fractal dimension of colony fronts from both cell lines is about 1.20 ± 0.05, indicating not far extended front ramifications.
(vi) The scaling analysis of colony fronts gives the same set of exponents, i.e., α = 0.5, β = 0.32, and z = 1.5, for both cell lines, covering a front width from a lower cutoff l c upward.
(vii) The growth dynamics for both cell lines becomes practically independent of the colony growth geometry.
The scaling exponents α = 0.50 ± 0.05 and β = 0.32 ± 0.04, derived from the dynamic scaling analysis in the l c to L width range, agree with what one would expect from both the standard KPZ continuous equation [20] and the cellular automaton model [17] in 2D space. This fairly good coincidence is accomplished despite the fact that some specific features of the biosystems-such as the colony age dependence of both the average cell size and morphology, the  formation of overhangs at the front, and the cell distribution in growth patterns-are variables that were not considered in the standard KPZ model. These findings are interesting for further improving theoretical models. Furthermore, the 2D front displacement in relatively large populated colonies involves the inward formation of 3D phase domains. Under constant cell nutrition, this contribution, combined with others such as contact inhibition, apoptosis, cell death, cell detachments from the colony itself, and cooperative phenomena involving biological and physicochemical interactions [35], produces a retardation in the 2D front displacement and, as referred to elsewhere [7,17], is responsible for the transition from an exponential growth regime when colonies consist of a small number of cells to a linearly front displacement regime when a large population in the colony is reached [7,36].
The dynamics of cell colony fronts is assisted by a number of cooperative interactions that result from cell compression, i.e., a progressively shorter average cell-cell distance in going from the front toward the colony bulk, which should affect the displacement of cells in the colony. Then, in contrast to a single cell movement that has been outlined as consisting of four steps [37], namely pseudopodial protrusions, the formation of focal adhesion, the development of contractile forces (translocation), and the detachment of all adhesions (retraction), in a colony these steps are temporally and spatially integrated through regulated intercellular responses, such as a signaling pathway and cycloskeletal reorganization [38]. Therefore, the colony growth global process becomes extremely complex as it involves the concerted actions of biological constituents of the colony and the physicochemical properties of the whole biosystem. However, despite the large number of variables that influence the biological behavior of cell colonies, the dynamic scaling interpretation apparently mimics the essential macroscopic and physical aspects of the growth process.
The KPZ behavior of both cell lines appears to be mainly related to those contributions that affect the generation of protrusion forces, a complex process in which chemical signaling results in the production of mechanical energy [1]. In general, the standard KPZ dynamics of spreading 2D fronts of condensed phases involves the participation of processes whose kinetics can be accounted for by a general equation that contains linear, nonlinear, and stochastic terms [12]. Each one of these terms is related to a certain characteristic contribution to the overall process. For instance, the linear term is related to a process that tends to smooth the surface, such as those related to bulk diffusion, surface diffusion, surface tension, or convective transport. On the other hand, biased displacement of the colony front is accounted for by the nonlinear term, which explains the roughening of the front due to the appearance of a lateral growth, related to the cell's own duplication at the front, local cell deformations, and cell motility [36]. Finally, the stochastic term accounts for those stochastic processes inherent to the evolution of a cell colony, such as cell duplications and random walk (cell motility) at the front, and cell membrane fluctuations.
Therefore, although a colony growth process involves concerted biological actions and physicochemical properties of the whole biosystem, which in turn would depend on the age of the colony and the size and shape distribution of the cells therein, the global growth process can be explained by the KPZ model. By comparing the scaling analysis of HeLa and Vero cell colonies, it also results that the growth dynamics is independent of the genetic backgrounds and the tumorigenic nature of cells. The universal dynamics is compatible with a linearly growth regime, and the lateral growth at the colony front is the determining factor for the colony-spreading process. However, it should be noted that, as stated elsewhere [39], the velocity of the front displacement cannot be explained by the simple formation of new cells at the front. The motility of cells apparently plays a key role in the advance of the colony front, and experiments tracking the paths of marked cells will help to advance the development of cell colony growth models [36].

V. CONCLUSIONS
(i) 2D HeLa cell colonies with a relatively large population result in patterns that display a rather disordered distribution of irregular cells of about the same average size. The morphology evolution of HeLa cells involves a smaller average size and a narrower cell size distribution at the colony front as compared with Vero cell colonies. As the colony size increases, a few large cells are formed, particularly at the outer regions. In addition, the progressive formation of 3D cell domains at inner colony regions takes place. These observations are independent of the colony growth geometry.
(ii) The HeLa colony front displaces at a constant velocity of about (0.09 ± 0.02) × 10 −4 μm min −1 , irrespective of the spreading geometry. This figure is about one-half the value that has been reported for Vero colonies under comparable conditions.
(iii) The rough colony profiles reveal an appreciable contribution of overhangs for front widths shorter than 100 μm, a figure that slowly increases with the colony age. The influence of overhangs at HeLa colony fronts becomes less remarkable than for Vero cells.
This set of scaling exponents is predicted by the KPZ equation for a condensed phase 2D growth.
(v) The Family-Vicsek relation is reasonably fulfilled by the above set of dynamic scaling exponents.
(vi) The comparison between HeLa and Vero cell colony dynamic data, under similar conditions, shows no clear-cut dependence of the scaling exponents on either the tumorigenic nature or the genetic background of the cells.