Clonal dynamics of surface-driven growing tissues

The self-organization of cells into complex tissues relies on a tight coordination of cell behavior. Identifying the cellular processes driving tissue growth is key to understanding the emergence of tissue forms and devising targeted therapies for aberrant growth, such as in cancer. Inferring the mode of tissue growth, whether it is driven by cells on the surface or cells in the bulk, is possible in cell culture experiments, but difficult in most tissues in living organisms (in vivo). Genetic tracing experiments, where a subset of cells is labeled with inheritable markers have become important experimental tools to study cell fate in vivo. Here, we show that the mode of tissue growth is reflected in the size distribution of the progeny of marked cells. To this end, we derive the clone-size distributions using analytical calculations in the limit of negligible cell migration and cell death, and we test our predictions with an agent-based stochastic sampling technique. We show that for surface-driven growth the clone-size distribution takes a characteristic power-law form with an exponent determined by fluctuations of the tissue surface. Our results show how the mode of tissue growth can be inferred from genetic tracing experiments.


I. INTRODUCTION
The self-organization of cells into tissue relies on the coordination of cell proliferation and differentiation in space and time.Broadly, tissue growth can be driven by the spatially homogeneous proliferation of cells (bulk growth).This mode of growth is characteristic of softer tissues like tendroins, arteries, or brain [1].Alternatively, tissues may grow by the preferential proliferation of cells on the surface, for example, because these cells have access to signaling molecules or vasculature.Surface-driven growth is often found in shells, horns, some bones [1], or tumors, where cells on the tumor surface have better access to nutrients [2].As a specific example of surfacedriven growth, in some fish and amphibians the eyecup forms by cell division in the outer part of the eye, the ciliary marginal zone [3].Understanding whether a given tissue grows by cell proliferation on its surface or in its bulk is important for targeting treatments during aberrant growth, such as cancer, it can form a template for developing synthetic tissues, and for understanding pathological development scenarios.In the example of the eyecup, cell divisions outside of the ciliary marginal zone, in the bulk, lead to the formation of additional blood vessels and scar tissue, and eventually to a disorder called proliferative retinopathy and to a complete loss of the eye's functionality [4].
The regulation of cell proliferation and the ensuing spatial distribution of proliferation patterns is governed on the one hand by complex biochemical signaling networks and cell-to-cell communication [5].On the other hand, it relies on how microscopic mechanical parameters, such as stresses, translate to a macroscopic scale.The connection between both is not well understood [6], such that an a priori inference of the mode of tissue growth is infeasible from a tissue mechanics perspective [5,7].Live imaging gives access to spatio-temporally resolved cell kinetics and allows for the estimation of tissue mechanical parameters [8][9][10][11].However, live imaging is highly challenging in vivo and it is usually limited to specific cases of embryonic development [12] or to studies of epithelium or other surface tissues [13,14].
The recent advent of genetic tracing experiments allows studying cell-fate behavior in vivo.In these experiments, a subset of cells is genetically labeled with fluorescent markers or genetic barcodes [15].As cells divide, these labels are passed on to all progeny of a labeled cell, termed a clone.The probability density of the sizes of such clones provides indirect information about the history of cell division, differentiation, and cell death events between labeling and the time point of analysis [16,17].For example, the first moment of the clone size distribution, the average clone size, reflects the rate of proliferation and whether both daughter cells remain proliferative or not.The functional form of the clone size distribution reflects how the fate of individual cells is decided [9,10,18], and the presence of mechanical forces [19] leading to clone fragmentation and merging [8,11,20].Therefore, the combination of genetic tracing experiments and tools from statistical physics has become a standard method for unveiling cell-fate behavior in vivo [21][22][23][24].
Here, we derive a theory that can help identify the mode of tissue growth from genetic tracing experiments.By studying the stochastic dynamics of clone boundaries and employing ideas from the range expansion process [25,26], we show that, for surface-driven tissue growth, the clone size distribution follows a characteristic powerlaw decay.The decay exponent depends on the rough-ness of the surface, which in turn is determined by the mechanisms regulating the tissue interface.We confirm our theoretical predictions with Monte Carlo stochastic lattice simulations with forward labeling.
This paper is organized as follows: In the following section, we introduce the model for surface-driven tissue growth and our mean-field scaling argument for the clone size distribution.We also show that clonal dynamics of surface-growing tissues can be mapped to the first-passage problem of a random walk in the presence of the absorbing boundary.In Sec.III, we present results for the agent-based lattice simulations of modified Eden cluster growth with label forwarding.We conclude in Sec.IV with a summary and discussion of our main results.

A. Model description
We consider a tissue in d spatial dimensions which are separated from other tissues by a d − 1 dimensional boundary.At a time t = 0 a random subset of cells is labeled, and when cells divide, this label is passed on to all progeny of a labeled cell.We are interested in how the size distribution of the number of cells that carry a given label at time t relates to the growth mode of the tissue.To begin our analysis of the clonal dynamics in surface-driven growing tissues, we make a simplifying assumption that cells labeled with the same marker remain spatially segregated.For this to hold true, the rate of cell death needs to be small compared to the rate of cell proliferation, which is generally the case for growing tissues.Moreover, the typical length scale associated with cell migration also needs to be small compared to the spatial extension of clones that result during the time course of the experiment.Under these conditions, cells that share the same marker form spatially-segregated clones such that boundaries that separate clones with different markers are well-defined.Such spatially segregated domains have indeed been observed in experiments of the growing retina of medaka fish [27].Under these assumptions, we show that the clone size distribution can be obtained from stochastic and geometric arguments alone without making further assumptions about tissue mechanics.
As the tissue grows, the boundaries of clones are subject to stochastic fluctuations, which originate from the random nature of cell divisions at the tissue surface [25,26].In the following, we will first derive expressions for the clone-size distribution in two spatial dimensions and then extend these results to three spatial dimensions.To this end, we will consider the stochastic wandering dynamics of clone boundaries -an approach that was first applied in the context of the random genetic drift of the range expansion process [25,26].
In two dimensions, the clone boundary dynamics stemming from the stochasticity of cell divisions can be de- fined by a stochastic process, X(t).The difference in distance between two adjacent clone boundaries, ∆X, has, as the label does not influence proliferation, a vanishing mean, ⟨∆X⟩ = 0, and the time evolution of its variance is described by a wandering exponent ζ, As an example, a wandering exponent ζ = 1/2 describes the Brownian motion of the distance between clone boundaries.The presence of tissue surface fluctuations may alter the wandering exponent, e.g., for tissue interfaces that are described by the Kardar-Parisi-Zhang equation, the wandering exponent takes a value ζ = 2/3 [25,28].
As the tissue grows, adjacent clone boundaries are subject to stochastic coalescence events.Thereby, a clone that is enclosed by two merging boundaries loses its access to the growing tissue surface (see Fig. 1).As a result of this merging event of the domain boundaries, the number of persisting clones, N p , i.e., the number of clones that have access to the front and that continue to forward their label and grow in size, decreases with time as

B. Mean-field scaling argument
In order to derive the size of persisting clones, we note that, in contrast to non-growing or bulk-driven tissues, the clonal dynamics in surface-driven tissue growth depend on the proximity of the clone to the surface of the growing tissue.Only clones containing cells at the tissue's surface can continue growing and contributing to the asymptotic shape of the clone-size distribution.If the linear extension of the tissue surface stays constant at a value L, we get an approximate expression for the average size of persistent clones, ⟨n p ⟩, by dividing the total tissue area by the number of persistent clones at time t, Here, v is the growth rate of the tissue that we assume to stay constant for a given cell division rate [29,30].Asymptotically, the fraction of clones with access to the surface vanishes.Therefore, the clone size distribution, P (n), is well approximated by collecting the sizes of clones that have lost their access to the front, i.e., by collecting the clones that have reached their terminal size.To calculate P (n), we therefore first calculate the number of clones that have lost their access to the moving front in a time interval dt, Then, using the mean-field argument that all persisting clones grow with the same average rate n(t) = ⟨n p (t)⟩ in Eq. ( 3), we obtain the clone-size distribution for surfacedriven growing tissues, The clonal size distribution has a unique, previously unobserved, power-law form, with an exponent that only depends on the wandering exponent ζ that describes the stochastic motion of clone boundaries.This result is in contrast to log-normal distribution and exponential distributions observed for bulk-driven growing tissues and in homeostatic tissues, respectively [9][10][11].
In three-dimensional tissues, clone boundaries are defined by stochastic surfaces.For a given clone, consider a slice along the direction of the growth.Within this slice, we expect the distance between the clone boundaries, ⟨(∆X) 2 ⟩ to scale like t 2ζ .In the absence of anisotropies, this scaling holds for all d − 1 directions perpendicular to the growth direction.We now consider the number of cells that share the same marker in a given slice perpendicular to the growth direction, A. Its deviation from the average, ∆A, fluctuates as If the number of cells in a given slice remains constant, the number of clones that have access to the surface decreases as (cf.Eq. ( 2)) For growth with a constant growth rate v, the average size of persistent clones increases as Finally, utilizing the same line of argument as for twodimensional tissues, the clone size distribution in d = 3 reads Taken together, these scaling arguments predict that the clone size distributions follow characteristic power laws.The associated exponents depend on the spatial dimension and the wandering exponent of clone boundaries, which is again influenced by the roughness of the tissue surface.For flat surfaces, where ζ = 1/2, the clone size distribution decays with an exponent of 4/3 for planar tissues and 3/2 for volumnar tissues.We derived these results in the limit of negligible curvature.For curved tissue surfaces, clone boundary coalescence halts asymptotically if the mean squared displacement of clone boundaries increases slower than the expansion of the tissue surface [26,27,31].Therefore, the results are strictly valid if 2 ζ > d − 1 for surfaces with constant curvature.Even if this is not the case, our results are applicable if the linear extension of clones, ∆X, is much smaller than the length scale associated with the curvature.This is generally the case not too long after labeling.Since genetic tracing experiments typically trace clones over several rounds of cell divisions, we expect our results to be broadly applicable.

C. Analogy with the first-passage problem
The power-law form of clonal size distribution can be obtained by associating coalescence events of clone domain boundaries that follow ⟨(∆X) 2 ⟩ ∼ t with the first passage events of a Brownian walker that hit the origin (see Fig. 2 (a)).Specifically, if the distance between two clone boundaries performs a random walk as the tissue front advances, then, since the position of the tissue interface h depends explicitly on time (e.g., h ∼ vt), the final clone size can be associated with the area X(t)dt that the random walker would cover before it hits the absorbing boundary at the origin.As such, the size distribution of clones that reached their final size is equal to the size distribution of the areas subtended by a randomwalk trajectory where t f is the first-passage time of the stochastic process X(t) hitting the absorbing boundary.
For diffusive motion, in the continuum limit, the firstpassage Brownian functional can be written as It has been evaluated analytically in [32], where, for large areas, the authors demonstrated that the probability density function P (A) follows a power law This result is in agreement with what we obtained with the mean-field argument for a Brownian case ζ = 1/2.We also confirmed this result by running Monte Carlo simulations for the one-dimensional random walk in the presence of an absorbing boundary at X = 0, for which we computed the areas A = t f t=0 X t covered by the Brownian walker before it hits the absorbing boundary (Fig. 2 (b)).For anomalous diffusion, ζ ̸ = 1/2, exact solutions for the first-passage functional A are not known apart from specific stochastic processes [33].

III. AGENT-BASED LATTICE SIMULATIONS
To test the validity of our analytical predictions, we performed numerical simulations of surface-driven growth in d = 2 and d = 3.For these simulations to verify the predicted power-law exponents they need to generate clones spanning orders of magnitude in size.Simulations of such extent are impossible when considering tissue mechanics and the details of biochemical processes underlying cell fate regulation.However, for surface-driven growth, if the rate of cell motility and loss are significantly smaller than the rate of cell division, the clone-size distribution is expected to depend only on the wandering and coalescence statistics of clone domain boundaries, and not on the underlying tissue mechanics or regulatory biochemical signaling network.Therefore, we use computationally efficient lattice simulations that capture the stochastic dynamics of clone domain boundaries and their relation to the clone-size distribution without necessarily being accurate microscopic representations of the tissue mechanics and regulatory processes.Specifically, we employ a modified version of the Eden cluster growth model, which is a minimal agentbased model that produces surface-driven cluster growth [29,34].In addition to the diffusion-limited branching process A λ − → A + A that increases the size of clusters by 1 with a rate λ, we randomly label agents A at the beginning of the simulation and allow them to pass their label upon replication.To produce and sample clone statistics, we employ Monte Carlo simulations of the described diffusion-limited birth process with label forwarding on two-and three-dimensional lattices.All of our simulations are initialized with a fully-occupied line (d = 2) or plane (d = 3) at x = 0, while the rest of the lattice is empty.Initially, each agent is endowed with a unique label.We update the system state using the Monte Carlo random sequential updating scheme.For a randomly chosen agent, we select at random an empty nearest neighbor lattice site and generate a new agent with the same label with a rate λ.
First, we simulated systems where clones have boundaries that follow Brownian dynamics, i.e., ζ = 1/2 in Eq. ( 1) and Eq. ( 6), respectively.This situation is realized for tissues with sharp and smooth interfaces, where surface fluctuations can be neglected.We achieve that in our simulations by choosing a space-dependent growth  rate, λ = (1 − tanh [α(x − x 0 )])/2, where the coefficient α sets the surface sharpness, x 0 = vt determines the surface position, and v sets the growth velocity.Choosing λ to have a sigmoidal functional form that varies only along the growth direction prohibits the development of surface fluctuations in directions perpendicular to the growth.As shown in Fig. 3 (a)-(b) for d = 2 and d = 3 cluster growth, respectively, the cluster interface stays flat at all times for this choice of λ, and clone boundaries perform a random walk.We then collect the number of different labels N p that have access to the front and measure the size of these clones to compute ⟨n p ⟩.To obtain the clonal size distribution, we collected the sizes of the clones that have lost their access to the advancing surface.
As shown in Fig. 3 (c)-(d), our simulations reflect the predictions made by the scaling arguments given above with some slight deviations for d = 3 cluster growth.Specifically, all our results for two-dimensional growth follow the predicted values after the early-time transient.We attribute this early-time transient to relaxation from the initial flat, sharp interface with zero width to the steady-state interface with a finite 1/α width.In three spatial dimensions, as shown in Fig. 3 (c), the number of persistent clones N p and the average size of persistent    clones ⟨n p ⟩ deviate slightly from our mean-field analysis.
It is plausible that the slower decay in the number of persistent clones N p for d = 3 stems from the fragmentation and merging clones, which can occur in d = 3 and is not considered in the mean-field theory.Nevertheless, if we substitute measured power laws for N sim p ∼ t −0.85 and ⟨n sim p ⟩ ∼ t 1.9 into Eq.( 4) and Eq. ( 5), the value for the clonal size distribution exponent that is predicted by the mean-field argument agrees well with our simulation data.This fact supports the connections that we have established when deriving the mean-field expression for the clone size distribution.
If the tissue interface is rough in the sense that the interface gradients are pronounced, tissue surface height fluctuations can no longer be neglected.As has been shown in [25,35], the boundaries between clones receive an additional drift contribution that comes from a surface tilt, rendering their dynamics superdiffusive, i.e., ζ > 1/2.In situations when the tissue interface evolution is described by the KPZ equation [28], which is often the case for rough interfaces, the expression for the wandering exponent becomes ζ = 1 + (χ − 1)/z [25].In this expression, the dynamical exponent z describes how the characteristic linear extension L x of the surface height fluctuation grows with time L x ∼ t 1/z , while the roughness exponent χ determines the scaling ratio of the fluctuation height to fluctuation's linear extension L z ∼ L χ x .For a one-dimensional tissue interface (d = 2, d surf = 1) whose dynamics is described by the KPZ equation, the exponents that characterize the scaling of surface fluctuations would belong to 1 + 1 KPZ universality class and take χ = 1/2, z = 3/2 values [28].Consequently, the wandering exponent will be equal to ζ = 2/3 and, according to our mean-field argument, P ∼ n −7/5 in d = 2.For three-dimensional tissues with two-dimensional KPZ interface, χ ≈ 0.382 and z = 1.618 [36], which results in ζ = 0.618 and P ∼ n −1.55 .We first confirm that the Eden cluster growth with constant agent division rate λ produces a rough traveling front with KPZ fluctuations [37][38][39].We do that in our d = 2 and d = 3 simulations by computing the width of the interface height fluctuations W (L, t) at the front of the growing Eden cluster: and then, employing the Family-Viscek scaling W (L, t) = L −χ W (t/L z ), we confirm that the surface height fluctuations are characterized by the KPZ scaling exponents.
After that, we use the same simulation algorithm to generate clones with different sizes, i.e., we begin with only the line or plane at X = 0 being occupied with agents that have unique labels, and we allow agents to pass their labels as they reproduce.Once the tissue interface reaches the other end of the simulation box X = L x , we stop the simulation and collect the sizes of the clusters that share the same label.As shown in Fig. 4 (a)-(b), in comparison to simulations with a flat front, pronounced surface gradients make clones lose their access to the surface at a faster rate.Measuring the mean-squared displacement of the clone domain boundaries from d = 2 simulations, we confirm that ⟨(∆X) 2 ⟩ ∼ t 4/3 and that the number of the persisting clones drops as N p (t) ∼ t −2/3 (Fig. 4 (c)).Similarly to the flat interface scenario, the clonal size distribution for simulation with KPZ surface dynamics follows a power-law decay with exponents that come very close to our mean-field predictions, as shown in Fig. 4 (d).For d = 3 simulations, we again observe slight deviations for N p and ⟨n p (t)⟩ quantities from our meanfield predictions, which we attribute to possible branch-ing and merging of the clones.

IV. CONCLUSION
In summary, we have studied the dynamics of clones for both d = 2 and d = 3 surface-driven growing tissues.We found that the clone-size distribution takes a power-law form with exponents depending on the tissue dimension and the nature of fluctuations in the surface.The power laws in the clone size distribution translate to associated power laws in the time evolution of the average clone size and the number of clones with access to the surface.These results suggest how the mode of tissue growth can be identified using genetic tracing experiments.While genetic tracing experiments using fluorescent markers typically do not produce a sufficiently high number of clones to confidently identify such power laws, recently developed technologies using genetic barcodes produce millions of unique clones in tissue [40] and can, therefore, be used to infer the mode of tissue growth as well as distinguish different kinds of surface fluctuations.
Throughout this work, we have been assuming that the tissue growth by cell divisions at its surface would lead to the formation of clonal sectors, i.e., the cells that share the same label would stay grouped up, and the domain walls that separate different clones would be clearly defined.While we have excluded cell migration and cell turnover from our analysis in order to be able to use this clone domain-wall dynamics approach, both of these processes are essential for tissue growth and remodeling and cannot be excluded from consideration entirely [41][42][43].As a potential direction for future research, it would be interesting to build a model that does not rely on the clonal domain wall description and would allow label mixing by including cell death and migration processes.Finally, the model of surface-driven tissue growth that we considered in this work leaves out the tissue mechanics and the effects that may come from it.For example, it would be interesting to consider a situation where cells that have lost access to the tissue interface and are still not far from the surface may continue dividing and regain access to the tissue interface once again.

FIG. 1 .
FIG. 1. Schematic illustration of expected clonal dynamics in surface-driven growing tissues -tissues where the cell divisions occur predominantly at tissue's surface.Cells that share the same color belong to the same clone.The arrow indicates the growth direction, and ∆X indicates the distance between two clone boundaries.

FIG. 2 .
FIG. 2. (a)Sample trajectories of one-dimensional Brownian walker in the presence of an absorbing boundary at X = 0.The enclosed areas can be obtained by evaluating the integral of first-passage Brownian functional A = t f 0 X(t)dt, where t f is the first time the process X(t) crosses the origin.(b) Probability distribution P (n) that area A takes a specified value n, when the Brownian motion begins at X = 0: comparison of the analytical prediction by[32] with our Monte-Carlo simulations.

3 FIG. 3 .
FIG.3.Snapshots of Monte Carlo lattice simulation of the birth process with label forwarding for (a) two-dimensional 1000 × 500 lattice and (b) three-dimensional 400 × 200 × 200 lattice.In both cases, the simulation begins at X = 0 with only the first perpendicular layer being filled with agents, each having a unique label that they can when they reproduce.The snapshots are taken right before the front reaches the other end of the lattice X = Lx.For both (a) and (b) we have kept periodic boundary conditions in directions perpendicular to the growth.(c) Time evolution of the number of persisting clones Np(t) divided by its initial value Np(0) and the average size of persisting clones ⟨np(t)⟩ (inlay).Both quantities are obtained from two-and three-dimensional Monte Carlo lattice simulations of a simple birth process with label forwarding, and time is measured in Monte Carlo steps.(d) The clonal size distribution and its local decay exponent.The inlay shows the local exponent.Error bars depict ± standard deviation.The data in (c) was obtained from simulations on 1000 × 500 and 1000 × 100 × 100 latices and were averaged over 10 4 independent realizations.The data for the clonal size distribution in (d) was obtained from simulations on 500 × 200 and 100 × 50 × 50 lattices and were averaged over 10 6 independent simulation runs.

3 FIG. 4 .
FIG.4.Snapshots of Monte Carlo lattice simulation of the birth process with constant rate λ = const.with label forwarding for (a) a two-dimensional 1000 × 500 lattice and (b) a three-dimensional 400 × 200 × 200 lattice.In both cases, the simulation begins at X = 0 with only the first perpendicular layer being filled with agents, each with a unique label they can forward when they reproduce.The snapshots are taken right before the front reaches the other end of the lattice X = Lx.For both (a) and (b), we have kept periodic boundary conditions in directions perpendicular to the growth.(c) Time evolution of the number of persisting clones Np(t) divided by its initial value Np(0) and the average size of persisting clones ⟨np(t)⟩ (inlay).Both quantities are obtained from two-and three-dimensional Monte Carlo lattice simulations of a simple birth process with label forwarding with λ = const., and time is measured in Monte Carlo steps.(d) The clonal size distribution and its local decay exponent, obtained from simulations with a rough KPZ front.The inlay shows the local exponent.Error bars depict ± standard deviation.The data in (c) was obtained from simulations on 1000 × 500 and 1000 × 100 × 100 latices and were averaged over 10 4 independent realizations.The data for the clonal size distribution in (d) was obtained from simulations on 500 × 200 and 100 × 50 × 50 lattices and were averaged over 10 6 independent simulation runs.