Pinning-depinning transition in a stochastic growth model for the evolution of cell colony fronts in a disordered media

We study a stochastic lattice model for cell colony growth, which takes into account proliferation, diffusion, and rotation of cells, in a culture medium with quenched disorder. The medium is composed both by sites that inhibit any possible change in the internal state of the cells, representing the disorder, as well as by active medium sites, that do not interfere with the cell dynamics. By means of Monte Carlo simulations we find that the velocity of the growing interface, which is taken as the order parameter of the model, strongly depends on the density of active medium sites (ρA). In fact, the model presents a (continuous) second-order pinning-depinning transition at a certain critical value of ρcrit A , such as for ρA > ρ crit A the interface moves freely across the disordered medium, but for ρA < ρ crit A the interface becomes irreversible pinned by the disorder. By determining the relevant critical exponents, our study reveals that within the depinned phase the interface can be rationalized in terms of the Kardar-Parisi-Zhang universality class, but when approaching the critical threshold, the non-linear term of the Kardar-Parisi-Zhang equation tends to vanish and then the pinned interface belongs to the quenched Edwards-Wilkinson universality class. PACS numbers: 87.17.Aa, 02.70.-c, 64.60.De, 64.60.-i, 68.35.Ct


I. INTRODUCTION
The study and understanding of growing interfaces is of great interest in many areas of science (physics, chemistry, biology, etc.), as well as for a wide variety of technological applications (micro-and nano-electronics, development of new materials, storage devices, etc.). From the theoretical point of view, the development of both the dynamic scaling framework for the rationalization of growth interfaces [1], and the formalism of continuous growth equations [2], such as that proposed by Kardar-Parisi-Zhang (KPZ) [3], provide powerful tools that have largely stimulated a broad spectrum of multidisciplinary activity [2,4,5].
In particular, the experimental and theoretical characterization of growth interfaces has been carried out for various biological systems such as tumor growth and cell cultures [6][7][8][9][10][11], as well as bacterial and fungi colonies [12][13][14][15][16][17]. The understanding of the collective behavior of interacting living cells and organisms, which is rather complex and somewhat intriguing, can be considerably enhanced by the characterization of its interfacial behavior, despite of that different points of views have also been used, such as active matter [18][19][20] and collective motion [21][22][23][24]. In fact, the dynamic scaling theory allow us to classify many types of growth interfaces in terms of a reduced number of universality classes that can be identified by means of few critical exponents, as e.g. in the study of critical phenomena in condensed matter and statistical physics [2]. Within this context, Huergo et al. [7][8][9] have performed in vitro experiments of growth interfaces of Vero and Hela cell cultures, which are two mammalian cell lineages that can be replicated through many cycles of division without becoming senescent. These interfaces have been characterized by means of the dynamic scaling theory with roughness and growth exponents that agree with those predicted by the well known KPZ universality class, i.e., α = 1/2 and β = 1/3, respectively [3]. Subsequently, we have proposed a minimal model for cell colony growth that gives growth interfaces in the same universality class [10]. More recently, Huergo et al. [25] studied the 2D front spreading dynamics of Vero cell colonies in a standard culture medium that has been gelled by the addition of methylcellulose, i.e. a hydrophilic polymer that increases the viscosity of the solution without being toxic to cells. It has been found that some cells become considerably enlarged due to the presence of methylcellulose. Furthermore, the colony population exhibits a quite small duplication rate as compared to experiments performed in a standard culture medium (without methylcellulose) [7][8][9]. The presence of rather large cells that behave as slow-moving obstacles tends to quench the displacement of smaller neighboring cells.
Consequently, it is found that the colony front displaces at a constant velocity that is about one half that the value already reported for the growth in absence of methylcellulose [25].
By applying the dynamic scaling theory to the digitalized images, Huergo et al. [25] have concluded that the obtained exponents, within a certain range of colony age and size, are consistent with those expected from the KPZ equation with quenched noise and a positive non-linear term, i.e. the PQKPZ equation [2,26,27]. While the authors reported the existence of pinning effects due to the presence of methylcellulose in the culture medium, an actual pinning-depinning transition has not been identified.
Motivated by those experiments and based on our experience gained by modeling the growth of a cell colony in an homogeneous media [10], the aim of this paper is to propose and study the generalization of the above-mentioned model for the case of a growing medium with quenched disorder. It is well known that the study of the displacement of interfaces in the presence of disorder introduces a very interesting analogy with the study of critical phenomena. In fact, the stationary motion of the interface requires the application of a driving force F , such as below a certain critical threshold the interface becomes pinned. So, one has a pinning-depinning transition where the velocity of the interface plays the role of the order parameter [2]. Our study reveals that within the depinned phase the interface can be rationalized in terms of the KPZ universality class. However, when approaching the critical threshold, the non-linear term of the KPZ equation tends to vanish and the pinned interface belongs to the quenched Edwards-Wilkinson (QEW) universality class [28].
The manuscript is organized as follows: in Section II we provide the definition of the model and the description of the simulation method, while in Section III we give an overview of the relevant theoretical background. The results are presented and discussed in Section IV, and our conclusions are stated in Section V.

II. STOCHASTIC CELL GROWING MODEL IN A DISORDERED MEDIUM
We use a lattice model recently proposed by us [10] that considers proliferation, diffusion, and rotation of cells. The model is formulated in a discrete two-dimensional square lattice, with four nearest-neighbors, and each site can be occupied by a cell, or by the culture medium. In the model, cells can be in three different states (I 1 , I 2 , and M ) since the proliferation cycle [29] is modeled in three main steps. The first part of the interphase (I 1 −cells) represents the G 1 phase, where the cellular contents, excluding the chromosomes, are duplicated. The second part (I 2 −cells) corresponds to the S and G 2 phases, where the chromosomes are duplicated and the cell checks the duplicated chromosomes, respectively.
Finally, we considered the mitotic phase (M −cells) where the cellular energy is focused on the orderly division into two daughter cells. The proliferation cycle as considered in the model is shown in Figure 1-a. The cell dynamics, which depends on the state of each particular cell, can be summarized as follows: (i) A cell of type I 1 spontaneously becomes an I 2 -cell, with probability P I 1 = 1, remaining at the same site. (ii) A cell in the state I 2 can either grow becoming a cell in state M with probability P I 2 = 0.99, or diffuse into a nearest-neighbor site with probability P d I 2 = 0.01, so that P I 2 + P d I 2 = 1. Since an M −cell occupies two neighboring sites, the growth process can only proceed if a randomly selected neighboring site of the original I 2 −cell is actually occupied by the culture medium and not by another cell. Note that the M -cells can be placed either horizontally or vertically, since we consider four nearest-neighbors. In the same fashion, diffusion of an I 2 −cell is only allowed to proceed to a neighboring medium site. (iii) Finally, a double cell of type M splits out into two I 1 −cells, occupying the same sites as the original ones, with probability P M = 0.98. Also, this type of cell can diffuse and rotate (in steps of ±π/2 degrees in the square lattice) with probabilities P d M = 0.01 and P r M = 0.01 (independently of the direction of rotation), respectively; such as P M + P d M + P r M = 1. Of course, these two latter processes involve randomly selected neighboring sites and they are constrained by the excluded volume of the cells. Previously, by studding the stochastic cell model described above in an homogeneous media, we have shown that the universal behavior of the growing interface does not depend on the relative values of the grow and diffusion probabilities [10]. Nevertheless, for higher values of the diffusion probabilities, the model dynamics is noisier, and greater lattice sizes should be used in order to obtain reliable results. This aspect is particularly sensitive in the present work, since large simulation times are required in order to determine the pinning transition, as we will see later. The rates of cell diffusion and proliferation in real biological systems can be very variable. In case of normal cells it depends on several aspects, like the cell type and function, density of cells, injury, aging, etc. In case of tumor cells, due to the lack of control mechanisms, usually both replication and movement are modified in an unexpected way. In several cases, diffusion is faster than proliferation, as in the case of highly motile cells as glioma [22,30], or for Vera and Hela cells [7][8][9]25]. Despite of that, in order to avoid prohibitive computing times, we choose to work with very low diffusion probabilities (and therefore, high values of grow probabilities), since we expect that this aspect does not interfere in the universal behavior of the growing interfaces.
As already mentioned, our aim is to propose a numerical approach to actual experiments of cell cultures where the culture medium is gelled by the addition of methylcellulose [25].
In order to formulate our model we considered that methylcellulose contributes to decrease the cell duplication rate. Also, the additive favors the occurrence of large cells that hinder cell motility in the colony, actually acting as pinning sites [25]. So, in our model, we chose at random a certain density of blocking sites (ρ B ) belonging to the medium, such as when the growing, diffusion, or rotation of a cell occurs in one of those sites, the cell also becomes blocked for future changes in its internal state or movement. Under these conditions we say that the cell becomes pinned (note that for a cell of type M it is enough that one of its sites becoming blocked in order to pinned the cell). On the other hand, in active medium sites (whose density we call ρ A ) the cell dynamics evolves according to the previously described rules for the homogeneous system, see [10]. Note that ρ A + ρ B = 1, so that for ρ A = 1 one has an homogeneous medium, as considered by us in [10,31], and the growth interface belongs to the KPZ universality class.
Simulations are performed in a square lattice with sides D × L and four nearest-neighbor sites. Along the D-direction (say the vertical direction) we define the lines of our system and identify the spatial location of a cell in a given line by the index i (0 ≤ i ≤ D − 1). Similarly, the L-direction defines the columns and the positions along them are identified by the index j (0 ≤ j ≤ L−1). The initial condition is selected by locating I 1 −cells filling the first column (j = 0) of the substratum. In this way the culture can grow along the L-direction and a growth interface that essentially runs parallel to the D-direction is naturally established.
The two-dimensional cell growing model setup represents a sector of the glass substratum of Petri dishes as used in actual growth experiments [7][8][9]25], where empty (parallel) strips surrounded by cell-cultured ones are obtained by using microfabricated stencils, as it has also been used for the study of scratch-wound-healing experiments [22,32], collective cell migration patterns [23,24], and interface characterization [7,9]. Also, we take periodic boundary conditions along the D-direction, while for the L-direction, which corresponds to the cell growth direction, the boundaries are closed on the left-hand side (because of the first column of the substratum is filled by cells), and at the right-hand side the boundary is left open, being this later assumption irrelevant since the interface of the culture is never allowed to reach the boundary. Figure 1-b shows a typical snapshot configuration obtained for ρ A = 0.9, that gives insight into the geometry used and the shape of the growing culture.
Therefore, in order to implement the simulations we first set the concentration of blocking sites (ρ B ) that are distributed at random on the substratum with probability ρ B . In this way, the presence of blocking sites represents the disorder of our growing medium and can be identified as a quenched noise. Subsequently, we implement the initial condition, as discussed above. Then, the actual dynamics process starts, namely a cell is selected at random and depending on its state (I 1 , I 2 or M ), we proceed according to the above-defined cell dynamic rules ((i), (ii) or (iii), respectively). After each random selection the simulation time t  Along this paper we will focus our attention on the properties of the growth interface of the cell culture. Despite that there are many possible definitions of such an interface, a multi-valuated interface (MVI) is the adequate one for our system, due to the presence of overhangs, especially near the pinning transition. Under these circumstances a singlevaluated interface does not provide a realistic description of the system [10]. In order to locate the MVI interface one first determines the biggest cluster of cells (connected by means of nearest neighbors) in contact with the initially filled first column of the culture, i.e. the so-called percolation cluster [33][34][35], that is denoted as land. On the other hand, empty sites in contact with the last column of the sample, which is never reached by the culture, are linked through both nearest and next-nearest neighbor sites and form a large cluster that is termed the sea. The sites that are not connected neither to the largest cluster of land nor sea are identified as islands and lakes, respectively, but they are irrelevant for the definition of the interface. In fact, the interface is given by the seashore where land and sea are in contact [33][34][35] (see Figures 1-c to 1-h). We call the set of sites belonging to the interface, {k}, k = 1, ..., N int , as the MVI, where the way of order these sites is irrelevant and N int ≥ D is the number of sites of the MVI. For further technical details on the method used for the location of the interfaces, see e.g. [33][34][35]. It is worth mentioning that the theoretical framework for the rationalization of interfaces (see next Section) has been developed for the case of single-valued interfaces, so we will discuss the use of a MVI below.

III. THEORETICAL OVERVIEW
The KPZ equation [2,3,5] is a nonlinear continuous growth equation for the time and spatial dependence of the interface height (h(x, t)) that evolves due to random deposition, diffusion, and lateral growth, namely The first term is the simplest possible form of accounting for the smoothing of the interface, so that ν is the diffusivity or viscosity. Also, the second (nonlinear) term accounts for the growth in the direction perpendicular to the interface, while F is the driving force. Finally, η is the noise that is used in order to introduce the disorder affecting the motion of the interface. Two main types of disorder are usually considered: (a) the thermal or annealed one that depends on time, i.e. η(x, t); and (b) the quenched noise that is frozen in the medium and depends on the local position on interface, i.e. η(x, h). Notice that for λ = 0, equation (1) yields the well known Edwards-Wilkinson growth equation [2,5].
On the other hand, for the study of discrete models (e.g. by means of numerical simulations) it is customary to measure the interface width given by the rms of the interface position, namely where D is the vertical size of the sample (see Figure 1-b), the brackets denote a spatial average, while the overbar denotes average over different realizations of the disorder. In order to determine the average horizontal position of the interface, h , we compute the horizontal coordinate (L-direction, see Figure 1-b) of each site belonging to the MVI. Therefore, h is defined as the mean value of being N int the number of sites of the interface, such as for a MVI one has N int ≥ D.
According to the standard dynamic scaling approach of Family-Vicsek [1], at early times the interface width grows with time as a power-law, namely where β is the growing exponent, and t x is the crossover time given by being z the dynamic exponent. Subsequently, one observes the saturation of the interface width that scales with D according to where α is the roughness exponent. The dynamic scaling theory [1] is based in the following scaling Ansatz: Here, f (u) is an universal scaling function such as where K is a constant. It is easy to show that only two exponents are independent and the relationship z = α/β, holds [2]. By considering annealed noise in (1+1)-dimensions, the KPZ universality class is represented by the following exponents: β KP Z = 1/3, and α KP Z = 1/2, which yields z KP Z = 3/2 [2]. Furthermore, due to Galilean invariance, a single exponent remains independent since the relationship α + z = 2 also holds [2,5].
As already mentioned in the Introduction, by considering quenched noise the interface undergoes a pinning-depinning transition that takes place at a certain critical value F c of the driven force (see equation (1)). The proper order parameter is the velocity of the interface, which for second-order transitions is expected to vanish when approaching criticality as where f = (F − F c )/F c is the reduced driving force, and θ is the order parameter critical exponent [2]. Note that one has a critical transition if the nonlinear term of equation (1) is positive (λ > 0), i.e. the positive quenched KPZ (PQKPZ) case, while for λ < 0, i.e. the negative quenched KPZ (NQKPZ) scenery, the pinning-depinning transition is of first-order [26]. Furthermore, for λ = 0 one has the quenched Edwards-Wilkinson (QEW) critical behavior. It is worth mentioning that equation (11) holds independently of the value of λ.
For the case of a critical transition (that is for λ > 0 and λ = 0) one has a diverging correlation length (ξ ) given by the characteristic length of the pinned domains (in the direction parallel to the interface) for F → F + c , which as in the case of standard critical phenomena, behaves as where ν is the parallel correlation length exponent.
Of course, not all the introduced exponents are independent. By using equations (5) and (6); and the fact that the growth correlation length becomes saturated by the sample width at criticality (ξ ∼ D), near the the pinning transition one expects that [28] v ∼ which also yields (c.f. equations (9) and (11)) where ν ⊥ is the correlation length exponent in the direction perpendicular to the interface, see also equation (33) below.
In order to determine the universality class of the pinned interface it is essential to investigate the behavior of the nonlinear term (λ) of equation (1) close to the critical driving force. For this purpose it is useful to impose an average tilt to the interface m = ∇y, and perform the transformation [28] h( Then it can be shown that the average velocity of the interface given by v(m) ≡ ∂y ∂t depends on the tilt according to [28,36] v In this way the coefficient of the nonlinear term of growing equations, such as e.g. equation (1), can be evaluated by considering specific models growing on tilted substrata. In particular, if λ diverges as [37] λ in the vicinity of F c , where φ ′ > 0 is an exponent, one has that the nonlinear term dominates yielding a QKPZ behavior, whose most studied case is the directed percolation depinning (DPD) universality class [38]. However, as already anticipated in the Introduction, we are interested in the case that λ vanishes when approaching criticality, namely with φ > 0. In this case, the nonlinear term of the KPZ equation (equation (1)) has to be dropped and one has the Edwards-Wilkinson equation with quenched noise, namely which leads to the QEW universality class. Pointing again our attention to the tilt dependent velocity of the interface given by equation (16) and by using equations (11) and (18) one It is worth mentioning that most exponents of the pinned phase corresponding to the QEW universality class can be derived from renormalization group calculations [39][40][41], yielding to first order and where ǫ = 4 − d, being d the dimension. Furthermore, according to reference [42] equations (21) and (22) may be exact for all orders in ǫ. Additionally, functional renormalization group calculations yields where α = (2σ − 1)/3. So, since our simulations are performed in (d + 1) dimensions with d = 1 (see Section II) the renormalization group predictions for the QEW universality class are α QEW = 1, z QEW ≈ 4/3 and β QEW ≈ 3/4 (see equation (9)). Also, one has ν QEW = 1 Further insight on the behavior of the interface close to criticality can be obtained by performing time dependent measurements of the interface velocity (v(t)), which resemble the well known epidemic analysis frequently used in order to study irreversible phase transitions [43]. In fact, for F > F c (f > 0) one expects that v(t) will approach a stationary value according to equation (11). On the other hand, just at criticality (f = 0) a power-law should hold, where q is the exponent of the critical velocity. Furthermore, within the subcritical regime (f < 0) the interface becomes pinned, so that the propagation velocity decays where the characteristic pinning time ξ t diverges at criticality according to where ν t is the corresponding critical exponent. Now it is useful to derive relationships among the above defined exponents. In fact, since some part of the interface is pinned at or below F c , and due to the random nature of the quenched noise, one has that the average height of the interface H(t) = h(x, t) (as measured by taken the pinned sites as a reference frame) is proportional to the width of the interface w [44,45], namely It is worth mentioning that equation (28) holds at the critical pinning transition, so that the exponent β corresponds to that of the universality class of the pinned interface, which may differ from that of the free interface. Then, v(t) = dH(t) dt ∼ t β−1 , so that at criticality one recovers equation (25) with and by considering that β QEW ≈ 3/4 (see equations (9), (21) and (23)), one has q QEW ≈ 1/4.
Further relationships among exponents are obtained by using finite-size scaling arguments.
In fact, as in the case of their equilibrium counterpart, the nonequilibrium pinning-depinning transition is affected by finite-size effects causing rounding of the transition and shifting of the critical points. In order to analyze those effects we recall that close to the critical point the correlation length in the direction parallel to the interface is constrained by the lateral size, so that the transition is observed at a certain "effective" size dependent driven force f c (D) such as (see also equation (12)), which yields a relationship for the shift of the critical point, namely where f c (D = ∞) ≡ f c is the true critical point and A is a constant.
Furthermore, the correlation length in the direction perpendicular to the interface is of the order of its width, which is also constrained by the system size. So, by using a similar line of reasoning one can write where we have used equation (6). So, by comparing equations (31) and (32) one obtains Since α QEW = 1 and ν QEW = 1 (see equations (21) and (22), respectively) one determines ν QEW ⊥ = 1. Let us now analyze the relationship among the characteristics lengths (ξ and ξ ⊥ ) and time (ξ t ) already defined. Close the pinning transition ξ gives the typical length of the interface attached to two pinning points. As time evolves, some region of the interface may start to move suddenly during a certain characteristic time of order ξ t , so that it moves a typical distance ξ ⊥ in the direction perpendicular to the interface. In this scheme, the nearly vanishing velocity of the interface is given by where we have used equation (11), so one has which by comparing to equation (14) and using equation (33) gives By remembering that z QEW ≈ 4/3 and ν QEW = 1 (see equations (23) and (22), respectively) one has ν QEW t ≈ 4/3. Finally, we recall that by using z = α/β, as well as equations (33) and (36), one obtains As already mentioned in the previous Section, the above described theoretical framework has been developed for Single-Valued Interfaces (SVIs), while in the present work we will evaluate MVIs. One expects that the interface could be described by means of a SVI when only the site located at the rightmost side of the interface is relevant for the dynamics.
However, in the present work the overhangs are important, specially when approaching the pinning transition (see Figure 1-d,f,h), and the use of a MVI seems to be the adequate tool to describe the system. The reason for the successful description of our data in terms of a MVI (see below), despite the fact that KPZ or EW equations are single-valued functions, is not totally clear. In the literature one can find a discussion about the relevance of the overhangs and its implication on the universality behavior of different models described by KPZ and EW equations [27,28,46,47], but scarce information is found about the use of MVIs or SVIs. Recently, we have reported a healthy agreement between simulation results of a cell colony growth model in an homogeneous media by using a MVI and the analytical results obtained from the KPZ equation [10].

IV. RESULTS AND DISCUSSION
We simulate the growth of cell colonies in a disordered medium, as described in Section II, on samples of different sizes, e.g. 32 ≤ D ≤ 512, and for different values of the driving force that is identified as the density of active medium sites (F ≡ ρ A ). The dimension of the horizontal size of the square lattice L is not relevant since we always ensure that the growth interface never reaches the right-hand border of the lattice. Earlier, we have shown that the cell colony model simulated in an homogeneous medium (ρ A = 1) gives growth interfaces that belong to the KPZ universality class [10].
In order to verify the effect of the disordered medium on the cell growth interface, we At a certain critical concentration (ρ crit A ) the interfaces do not move anymore (v = 0) by the excessive presence of inactive sites, so that one observes that all the interfaces become pinned below the critical threshold. The pinning-depinning transition of the interface takes place at a critical value of density of active medium sites, which can roughly be estimated from the Figure 2-(b) as ρ crit A ∼ 0.74. In fact, as already discussed in Section III, at criticality one expects that the interface velocity would decay according to equation (11). Therefore, a more accurate estimation of ρ crit A can be obtained by plotting in a log-log scale v versus , for different trial values of the critical density of active medium sites ρ * A , as shown in Figure 2-(c). Here, one has that the best straight line is obtained for ρ * A = ρ crit A = 0.740 (2), in full accordance with the result of Figure 2-(b). Also, from the slope of this line the exponent θ = 0.21 (3) is also determined (see equation (11)). The study of the growth front velocity as a function of time can give us further insight on the universal behavior of the system. In fact, for ρ A > ρ crit A the interfaces grow at constant speed (after a short transient), whereas for ρ A < ρ crit A the interfaces become pinned, reaching zero velocity for a finite time (see Figure 3-(a)). According to equation (25), we expect a power-law decay of v just at the critical point, which allow us to determine the critical density as ρ crit A = 0.74(1) (see Figure 3-(a)), in full accordance with the results already discussed in the context of Figures 2-(b) and (c). From the short-time slope of the data presented in Figure 3-(a) for ρ crit A = 0.74 we obtain q = 0.17(1) for the exponent of the critical velocity (see equation (25)). Also, within the subcritical regimen we can define a cutoff time of the velocity (t * ) as the time such as v(t * ) ≡ v * = 0.01, where the value of v * is arbitrary, with v * ≪ v(t = 0). The cutoff time t * should diverges as ρ A → ρ crit A , since just at criticality the velocity exhibits a power-law behavior. In fact, t * is proportional to the characteristic pinning time ξ t , that is t * ∼ ξ t ∼ |f | −νt (see equation 27), in agreement with the fact that pinning causes the drop of the interface propagation velocity. Therefore, from the density dependence of t * (inset of Figure 3-(a)) we determine the critical exponent ν t = 1.32 (5). On the other hand, also within the subcritical phase we expect an exponential decay of the velocity of the interface (c.f. equation (26)), as can be seen from the log-linear plot of v versus t shown in Figure 3- From the fit of the data we obtain the characteristic pinning time (ξ t ) that scales according to equation (27), as verified in the inset of Figure 3-(b), which yields ν t = 1.36(2) for the temporal correlation length exponent, in full accordance with the result obtained from the cutoff time of the velocity (inset of Figure   3-(a)). In order to test both the critical velocity exponent q and the temporal correlation length exponent ν t , we rescale the data presented in Figure 3-(a) for the subcritical regime according to equations (25) and (27), by using q = 0.17 and ν t = 1.32. We can observe a good data collapse in Figure 3-(c).
. From the best fit of the data we obtain ν t = 1.32 (5). (b) Log-lineal plot of v as a function of time (t) as measured within the subcritical phase (ρ A < ρ crit A ), as indicated. Inset: Log-log plot of the characteristic pinning time (ξ t ) as a function of . From the best fit of the data we obtain ν t = 1.36 (2). (c) Scaling plots of the data already shown in (a) according to equations (25) and (27) Also, the power-law behavior of the interface velocity at the critical point (see equation (25)) is observed for larger times when the size of the sample is increased, as shown in Figure   4-(a) for samples of different sizes, while the drop of the curves is due to finite-size effects. In the inset of Figure 4-(a) we show the cutoff time of the velocity (t * ) determined for v * = 0.01 as a function of 1/D. As expected, t * diverges as D → ∞, since a full power-law behavior is expected to be found only at the thermodynamic limit. We propose that the cutoff time is proportional to the crossover time of the interface width t x , that is t * ∼ t x ∼ D z (equation 5), since the drop of the interface velocity within the subcritical phase and the subsequent pinning causes the saturation of the interface width. Therefore, from the best fit of the data presented in the inset of Figure 4-(a), we obtain z = 1.2 (2). Also, we rescale the data presented in Figure 4-(a) according to equations (25) and (5), by using q = 0.17 and z = 1.2 (see Figure 4-(b)). An excellent data collapse is found, in particular note that the scaled velocity has slope approximately zero for a broad time window that corresponds to the power-law decay of the raw data with exponent q = 0.17, and subsequently it drops due to the finite-size effects, as expected. In order to better characterize the pinning transition we define the pinning probability, p pin (ρ A ), as the occurrence frequency of pinned interfaces at a certain density of active medium sites ρ A : where n pinning (ρ A ) is the number of pinned interfaces observed and n total is the number of realizations performed. Figure 5 shows the pinning probability as a function of the density of active medium sites for different lattice sizes. For high values of ρ A the interfaces move freely so that p pin = 0, whereas for small values of ρ A all the interfaces are pinned and one has p pin = 1. The size-dependent transitions between these two regions occurs for effective values of the critical density of active medium sites. Note that the transition between p pin = 1 and p pin = 0 becomes rather stepped for the larger lattice sizes. Also, we define the "effective" lattice size-dependent critical density ρ 50% A (D) as the density of active medium sites such as p pin = 0.5. As in the standard percolation problem [48], we expect that ρ 50% A (D) should scale as (c.f. equation (31)) where ρ crit A (∞) = ρ crit A is the critical density in the thermodynamic limit, B is a constant, and ν is the correlation length exponent in the direction parallel to the interface. From the extrapolation of ρ 50% A (D) to infinite lattice sizes (inset of Figure 5) we obtain an independent measurement of the critical active medium density ρ crit A = 0.741(3) and also ν = 0.98(2). After locating the pinning-depinning transition by means of both dynamic and stationary measurements, we study the cell growth interface by using the standard dynamic scalingapproach of Family-Vicsek [1] in the region where it moves freely, in order to verify the effect of the disordered medium sites on its properties. Therefore, we performed simulations studies well above the critical point (ρ crit A ≃ 0.74), that is for ρ A = 0.90. Figure 6 shows loglog plots of the interface width (w) versus time for different lattice sizes. One observes that after some transient period of the order of ≃ 50 MCS, the interface width scales according to equation (4), e.g. w ∼ t β , with β ≃ 0.30 (3), being this behavior more evident for the largest samples, as expected due to finite-size effects. Subsequently, w saturates as predicted by equation (6), so that w sat ∼ D α , with α = 0.47(4) (c.f. left-hand side inset of Figure 6). This figure of the roughness exponent is also consistent with the determination of the structure factor [10] (see right-hand side inset of Figure 6), which yields α = 0.45 (5). Based on these results we conclude that the free interface with ρ A = 0.9 belongs to the KPZ universality class (α = 1/2, β = 1/3 [3]). This finding agrees with our previous determination for the growth interface in absence of disorder (ρ A = 1), that also belongs to the KPZ universality class [10]. Results obtained for the free interface but for ρ A < 0.9 are affected by large crossover effects so that a dedicated study would require a computational effort beyond our capability. For this reason and based on the fact that one expects that the universality class of the free interface would not depend on ρ A above the critical point (that is for ρ A > ρ crit A ), we conclude that the free interface belongs to the KPZ universality class. Now, we focused our effort in the characterization of the interfacial properties of the pinned interface just at the critical point. We will show that these results are in full agreement with the epidemic measurements presented and discussed previously. Therefore, we have studied the interface width at the critical point, that is ρ A = ρ crit A = 0.74. In that way, Figure 7-(a) shows that the time behavior of w for t ≪ t x allows the determination of the growing exponent β (according to equation (4)) as 0.83 (3). From the lattice size dependence of the interface width during the saturation regime of the interface w sat (as given by equation (6)) we obtain the roughness exponent α = 1.02(3), as can be seen from the inset of Figure 7-(a). Also, from Figure 7-(b) we observe the expected scaling behavior (see equation (7)) when using the exponents obtained from Figure 7-(a), with z = α/β = 1.23 (6).
Of course, departure from scaling is observed for the short-time regime due to the influence of the initial condition selected given by a flat interface. It is worth mentioning that the obtained exponents β = 0.83(3), α = 1.02(3), and z = 1.23(6) of the growing interface at the critical point are far-from the corresponding ones for the KPZ universality class (namely β KP Z = 1/3, α KP Z = 1/2, and z KP Z = 3/2 [2]). They present, in fact, a better accordance with the values of the QEW universality class predicted by renormalization group calculations, e.g. β QEW ≈ 3/4, α QEW = 1, and z QEW ≈ 4/3 [39][40][41]. Also, note that only α QEW is expected to be exact (see equation (21)), while β QEW and z QEW are approximations of order ǫ 2 , with ǫ = 4 − d and d = 1 (see equations (9) and (23)).  (7), and applied to the data already shown in Figure 7-(a). The data collapse is achieved by using the exponents α = 1.02 and z = α/β = 1.23. More details in the text.
Since the interface at the critical point presents a clear departure of the KPZ behavior, it is worthwhile to investigate the dependence of the non-linear term of the KPZ equation (λ in equation (1)), in order to determine if it remains relevant or it vanishes when approaching the critical pinning-depinning transition. In this way, we impose different average tilts m to the starting interfaces and study the velocity v as a function of the squared tilt m 2 for different values of the density of active medium sites (with ρ A > ρ crit A ), as can be seem in Figure 8-(a). Note that the results for m = 0 are the same as those shown in Figure   2-(b). By using equation (16) we are able to determine λ, which clearly decreases as the critical point ρ crit A = 0.74 is approached (see Figure 8-(b)), indicating that it vanishes at the critical point. Therefore, when the density of active medium sites is greater than the critical value (ρ A > ρ crit A ) equation (1) holds, and the growing interface belongs to the KPZ universality class. However, at criticality the behavior of the system is described by the Edwards-Wilkinson equation (i.e. with λ = 0) with quenched noise (equation (19)).
Summing up, in this Section we have already presented and discussed epidemic like studies of the interface velocity, static scaling results of the pinning probability, as well as a dynamic scaling analysis of the interface width close to the critical pinning-depinning transition, that takes place at ρ crit A = 0.74(1). All those obtained results allow us to determine independently eight different critical exponents that are summarized in Table I. Also, the relationships between the determined exponents allow us the derivation of several exponents, as presented in Table I, summing up a total of nine different exponents obtained in the present work. Exponents reported by other authors and renormalization group results for pinning-depinning transitions in the QEW universality class are also listed in Table I for the sake of comparison.

V. CONCLUSIONS
In the present work we have studied a cell colony growth model that considers proliferation, diffusion, and rotation of cells [10] in a culture medium with quenched disorder.
Actually, we are inspired in experiments of cell culture where the culture medium is gelled by the addition of methylcellulose, whose presence contributes to decrease the cell duplication rate, as well as the occurrence of large cells that hinder cell motility in the colony,   (1)) as obtained in the present work both from direct calculations (2nd column) and from relationships between the obtained exponents (4th column). In both cases the equations used in order to determine them are indicated, in the 1st and 3rd columns, respectively. For the sake of comparison we also show previous results from others authors (5th column) and results obtained with the renormalization group (RG) for pinning-depinning transitions in the QEW universality class (6th column). The symbol * indicates the references to RG calculations: [39][40][41]. Exponents identified with the symbol † were obtained from relationships, as indicated in the text.
acting as pinning sites [25]. We have modeled the disordered culture medium by the presence of blocking sites, i.e., sites that inhibit any possible change in the internal state of the cells, as well as active medium sites, that do not interfere with the cell dynamics. Our aim was to study and characterize the behavior of the interface of growing cells in this gelled medium for different concentrations of active medium sites, verifying the existence of a pinning-depinning transition.
As expected, we find that the velocity of the growing interface, that constitutes the order parameter of the model, depends strongly on the density of active medium sites ρ A . In fact, the model presents a second-order pinning-depinning transition for ρ crit A = 0.74 (1). That is, if ρ A > ρ crit A the interface can move freely, whereas for ρ A < ρ crit A the interface becomes pinned and its velocity vanishes for a finite time. It is worth mentioning that we determined the critical point by means of different dynamical and stationary methods, such as from the measurement of the velocity of the growing interface, from the power-law behavior of the velocity at criticality both as a function of time and of the reduced force ( ρ A −ρ A crit ρ A crit ), as well as from the definition of a pinning probability. We observed a full accordance between the obtained values.
Also, we shown that when the interface can move freely (e.g. for ρ A > ρ crit A ) its scaling behavior can be rationalized in terms of exponents belonging to the KPZ universality class (we obtain β = 0.30(3) and α = 0.47(4) for ρ A = 0.9, in accordance with the expected values α KP Z = 1/2, β KP Z = 1/3 [3]). This result is in agreement with our previous analysis for the cell model in a medium without disorder (that is ρ A = 1) [10] and with in vitro studies of Vero and Hela cells in an homogeneous culture medium [7,9,31], that also belong to the KPZ universality class. However, at criticality, the interface is very rough, the overhangs are clearly not negligible (see Figure 1-f) and the growing interface is described by exponents that are far-from the corresponding ones for the KPZ universality class. An analysis of the dependence of the non-linear term of the KPZ equation (λ in equation (1)), performed by studying different tilts applied to the starting interfaces, allows us to determine that λ → 0 when the critical point is approached. This result indicates that at the critical point the EW equation with quenched noise (see equation (19)) should provide the correct description of the system, which is then expected to belong to the QEW universality class. In fact, when comparing the exponents obtained by us at the critical point and those reported for the QEW universality class, we can find an excellent agreement for α = 1.02(3) and ν = 0.98(2), namely α QEW = 1 and ν QEW = 1 [39][40][41], whose values, obtained by means of renormalization group calculations, are expected to be exact. Also, the determination of z = 1.2(2) and ν t = 1.32(5) − 1.36 (2) are in very good agreement with the theoretical calculations, z QEW ≈ 4/3 and ν QEW t ≈ 4/3 [39][40][41]. Our determinations of β = 0.83(3), θ = 0.21(3) and q = 0.17 (1) show some differences with the theoretical values that are considered as approximate up to order ǫ 2 , with ǫ = 4 − d and d = 1, namely β QEW ≈ 3/4, θ QEW ≈ 1/3 and q QEW ≈ 1/4 [39][40][41]. As expected, our determination of the critical points is affected by finite-size effects. In fact, the closeness of the critical density rises the probability that the interfaces become pinned (even in the region where the interface moves freely), particularly for smaller lattice sizes. Therefore, we expect that the agreement between simulation and exact results for the critical points would be enhanced by the use of bigger lattice sizes and larger simulation times. It is worth mentioning that different models considered in the QEW universality class present scarce agreement with theoretical results, particularly for β and θ, as in our case. The critical exponents obtained at the critical point are summarized in Table I. On the other hand, Amaral et al. [28] reported a pinning-depinning transition for a random field Ising model, in which the non-linear term λ tends to zero just at the critical point, as in our case.
Interestingly, we have also proposed an interplay between the drop of the velocity of the interface within the subcritical phase, and the crossover time of the interface width observed in the context of the standard dynamic scaling approach of Family-Vicsek [1]. In that way, we determine the exponent z from measurements of the size dependence of the cutoff time of the velocity of the interface (see equation (5)), in excellent agreement with the expected theoretical value (see Table I).
Summarizing, by means of a numerical simulations we found a second-order pinningdepinning transition in a cell growing model, which is driven by the density of active medium sites ρ A . The universality class of the transition changes from KPZ, in the region where the interfaces can move freely (ρ A > ρ crit A ), to QEW universality class at the critical point ρ crit A . Recent experimental results of Vero colonies cells in a medium gelled by the addition of methylcellulose suggest that the moving interfaces may belong to the QKPZ universality class [25], since they found α = 0.63(4), β = 0.75(5) and z = 0.84(5) for a concentration of 2.5% of methylcellulose in the medium culture. However, no attempts to perform a study of the dependence of the interface properties on the concentration of the additive are reported.
Consequently, an actual pinning-depinning transition could not be identified. According to our results, the precision in the determination of the critical exponents depends on the closeness to the critical point, due to the presence of both finite-size and crossover effects that may hinder the actual physical behavior. Consequently, it would be very interesting perform further experiments for different concentrations of methylcellulose, in order to determine the possible existence of a true pinning-depinning transition as well as the actual universality class of such a transition.