Interface dynamics of competing tissues

Tissues can be characterized by their homeostatic stress, i.e. the value of stress for which cell division and cell death balance. When two different tissues grow in competition, a difference of their homeostatic stresses determines which tissue grows at the expense of the second. This then leads to the propagation of the interface separating the tissues. Here, we study structural and dynamical properties of this interface by combining continuum theory with mesoscopic simulations of a cell-based model. Using a simulation box that moves with the interface, we find that a stationary state exists in which the interface has a finite width and propagates with a constant velocity. The propagation velocity in the simulations depends linearly on the homeostatic stress difference, in excellent agreement with the analytical predictions. This agreement is also seen for the stress and velocity profiles. Finally, we analyzed the interface growth and roughness as a function of time and system size. We estimated growth and roughness exponents, which differ from those previously obtained for simple tissue growth.


Introduction
The mechanics of growing tissues has received much attention recently from physicists and biologists alike. For example, in the case of cancer tissue growth is determined by the competition for space between a tumour and healthy host tissue. A key characteristic determining which tissue can grow as the expense of the other was proposed to be the homeostatic stress [1]. It is defined as the value of the stress for which cell division and cell death balance, corresponding to a homeostatic state with constant cell number and cell density.
For larger (expansive) stresses, the balance is shifted towards cell division, while for lower (compressive) stresses, it is shifted towards cell death. Therefore, if homeostatic stresses are isotropic and if we ignore surface effects, the tissue with the lower homeostatic stress in general grows while the one with larger homeostatic stress disappears by cell death.
Tissue growth is affected by boundary effects and by external mechanical forces [2][3][4][5]. Such effects also play a role at interfaces that separate different tissues [6,7]. A simple mechanical argument shows how the interface between a tissue and surrounding liquid, or a tissue-air interface, could affect tissue growth: a cell has to increase its volume in order to grow and divide. Considering a growing cell for simplicity as a strain dipole, we can discuss the work associated with the insertion of the strain dipole into the tissue. It can be shown that to insert a strain dipole near the surface (where a part of the strain field does not contribute to the work) requires less mechanical work as compared to the case where the dipole is inserted in the bulk. As a consequence, cell division and growth occur at a larger rate near the surface. This effect has been demonstarted by mesoscale simulations of growing tissue spheroids [8].
The population dynamics of two competing tissue layers growing on a substrate has been studied recently by a generalized Fisher-Kolmogorov equation that takes into account the effects of tissue stresses [9]. The Fisher-Kolmogorov equation was originally proposed to describe the kinetics of spreading of advantageous alleles in a population [10]. The resulting traveling-wave solutions are driven by diffusion. Therefore, for vanishing Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. diffusion coefficient, the velocity of these waves approaches zero. In this work, we explore the properties of pressure-driven tissue competition by combining mesoscale simulations of a cell-based model with continuum theory. We find that the tissue with the lower homeostatic stress takes over the system via the propagation of a growth front that moves with constant velocity. The main characteristics of the moving tissue-tissue interface in the simulations, namely the stress and velocity profiles, agree very well with analytical calculation in the continuum theory. Over time, a finite interface width develops. This interface width exhibits a power-law scaling with system size that characterizes interface roughness.

Cell-based mesoscopic model
In order to study the competition between tissues, we used the simulation approach introduced in [11,12]. A growing cell is represented by two point particles i and j that are pushed apart by a growth force Here G is the strength of this repulsion, r ij is the distance andr ij the unit displacement vector between the centers of particles i and j, and r 0 is a length scale at which the force saturates. This growth force is opposed by external stresses, and cell-internal friction g c . This internal friction sets the growth velocity together with the growth strength G. Upon reaching a certain critical size r c , the cell divides. Division implies that new particles are added in very close proximity to i and j which then define two separate cells. Apoptosis is introduced as a constant rate k a at which particles corresponding to a cell disappear. These are the only active processes in the model; the remaining passive interactions resemble soft sticky spheres. The impenetrability of cells is ensured through the volume exclusion force with the cell-cell potential coefficient f 0 and the range R of all pair potentials acting on all pairs of particles i j , that do not belong to the same cell. Adhesion is represented by a simple constant force= - , not belonging to the same cell. Additionally, cells dissipate momentum and energy on the tissue scale by a friction g t with neighbouring cells. All dissipative forces are accompanied by appropriate random forces. Noise intensity is strong enough to break local minima, but small enough to have otherwise no noticeable effect. In particular it is not strong enough to lead to significant cell diffusion by itself. So far the model is identical to the one presented in [11][12][13]. Here however, we focus onthe competition between cell sheets, growing on a substrate. Thus, the simulations are restricted to two-dimensions and a friction force bg i is introduced describing the interaction with a substrate. In the supporting information (SI) we provide a comprehensive list of all parameters. This in-silico tissue grows in number of cells until it fills all available space. When the compartment is filled, the pressure increases and slows down divisions to the point where they are exactly balanced by apoptosis. The pressure at this stationary state is the homeostatic pressure P H , or equivalently the homeostatic stress s = -P H H . The homeostatic stress s H of the tissue depends on the choice of model parameters. In particular, it can be varied by changing the growth-force strength G, as shown in [13]. In the remainder of this work, we measure space in units of the particle interaction range R, time by the apoptosis rate k a (per definition equal to the division rate in the homeostatic state), and force in units of G R 0 2 , where = G G 0 for the weaker tissue in all simulations of competing tissues. Thus we report stresses normalized by R G 4 0 . All rescaled quantities are denoted by an asterisk, e.g. * s s = R G 4 0 . The spatially resolved stress tensor at any point P is obtained through a virial method similar to the one described in [14] (see SI for details). Note that in our case all forces, in particular the growth force, arise from pair interactions and contribute to the stress. This is different from simulations of active micro swimmers, where an additional swim stress contributes to the overall stress of the system [15,16].
For computational efficiency, we designed a set-up with a co-moving simulation box, which we call treadmilling. The simulation box consists of bounce-back boundaries in x-and periodic boundaries in y-direction, and is initially filled with one type of cell (A) which is grown until it reaches its homeostatic state. Then, all cells with < r L 2 x x are changed to another type of cell (B) that has the same interaction parameters, but an increased growth force G, balanced by an increased intra-cell friction coefficient g c to generate approximately the same free growth rate. Thus, B has a smaller homeostatic stress than A. The system size L x is chosen such that both tissues reach their homeostatic state sufficiently far from the interface. Since the tissue with the lower homeostatic stress (larger homeostatic pressure) pushes aside the other, the interface is kept at L 2 x by displacing all cells every 1000 integration steps by h refers to the current averaged interface position. In order to remove the excess number of cells of the stronger tissue, all cells entering a small region near the hard wall are taken out of the simulation. The background friction near this region is continuously increased to ten times its bulk value, such that cell velocities become negligible in accordance with a bulk situation. The weaker tissue simply replenishes itself by growing into the free space that opens up after the cell displacement step (see SI for further details). We verified the validity of this approach by comparing the results to those of fixed boundaries simulations (see SI).
Due to overhangs and island formation (see figure 1(a) and SI movies S1 and S2) a unique interface position is not easily defined. Apart from utilizing the cell number fraction, defined below, to characterize the interface position and width, we also employ Voronoi tesselations, which are determined by the program Triangle [17]. All Voronoi edges that separate particles of different tissue types are therein defined as part of the interface.

Interface dynamics in the mesoscopic simulations
Simulating two tissues that only differ by their homeostatic pressure, we find a stable steady-state interface in a co-moving simulation box. Starting from an initially flat interface, we observe a roughening over time, with interface width reaching a finite plateau on longer time scales (see figure 1). The interface moves with a constant ). The initially flat interface roughens over time, while propagating to the right. Note that the scale is the same as for (c) and velocity v 0 that depends on the homeostatic stress difference s D H . The local average cell velocity v x in the lab frame depends on the distance from the interface. It decays from the interface velocity v 0 to zero on a finite length scale l » R 10 for increasing distance from the interface. The coordinate s refers to cell position in a comoving coordinate system =s x v t 0 , with s=0 describing the interface position. The stress profile ( ) s s xx varies on the same length scale from one homeostatic stress to the other as shown in figure 1(e). Note that both tissues indeed reach their respective homeostatic stresses far away from the interface. In the sections below, we study the competition and interface properties using our simulations and analytic theory. We focus on competition of tissues growing on a substrate in two-dimensions, and analyze the emerging stress and velocity profiles as well as the interface propagation and interface width.

Velocity and stress profiles
The velocity and stress profiles can be understood using a simple continuum model in one-dimension. The competition between two incompressible tissues A and B with constant density is governed by the cell-number balance equation where we have neglected for simplicity terms corresponding to cell diffusion. Here which relates the velocity gradient to the growth rates k c . The growth rates k c , in turn, depend on the local stress σ where we have expanded to linear order and neglected higher order terms, see [1]. H . We find the prediction for the stress profile to be in very good agreement with the simulation results (see figure 1). In order to make such a comparison, we characterized the tissues used in the simulations by the methods outlined in [13], which result in the growth rate coefficients k c , the homeostatic stresses s c H , and the bulk density r b (see SI). The velocity profile, which can be obtained analytically from equations (5) and (6), also agrees very well with the simulations (see figure 1(b)) without adjusting any parameters. Furthermore, the linear dependence of the interface velocity v 0 on the homeostatic pressure diffrence s D H in the simulations can be understood by the analytical calculation (see figure 1(c)).

Interface width
Over the course of the simulations, we observe a finite steady-state width of the interface that mainly depends on the channel width L. Figure 2(a) shows a typical snapshot of the interface for a relatively wide channel. To quantify the interface width, we define ( ) = á ñ -á ñ w t h h 2 2 , with the interface function ( ) h y t , . The first two moments of this function are determined by (see SI) We can then follow the evolution of the interface width w(t) over time, as shown in figure 2(b). For small times t, w displays a power-law increase with time, ( )~b w t t (with growth exponent β), as might be expected by analogy with othergrowing interfaces. For two identical tissues, one would expect an ever growing interface width. However, a difference solely in their respective homeostatic stresses, is enough for interface width to saturate over time (see figure 2(b) and inset). This saturated interface width grows with the system size L, and seemingly obeys a scaling relation~a w L sat [18], where α is called the roughness exponent (see figure 2(c)). We obtain α from a fit of w sat (L) and a fit of the structure factor in the steady state (see [19][20][21][22] for details). Here,˜( ) h k t , denotes the Fourier transformation of the height

Island formation
The simulations revealed another interesting phenomenon: the interface is not a single-valued function, but we furthermore observe islands of weaker tissue left behind in the stronger tissue, (see figure 3(a) and SI movie S1). These islands consist of cells of the weaker tissue that eventually die due to their higher homeostatic stress. We characterize the island formation using Voronoi tesselation (see figure 3(b)), as outlined in the methods section.
The average number of islands á ñ n I per unit length is independent of the system size L but seems to vary nonmonotonically with the homeostatic stress difference s D H (see figures 4(a) and (b)). In the limit of s D = 0, both tissues are the same and thus would mix completely on large time scales, i.e. the number of islands is expected to diverge as s D approaches zero. For larger s D , on the other hand, the increasing interface propagation velocity may explain the higher detachment rate and thereby an increased number of islands.

Stress across interface
So far, we have only considered the principal stress s xx in the direction perpendicular to the interface. The interface, on the other hand, breaks the x-y isotropy, which results in a surface tension (see e.g. [23]). In our simulations, however, the interactions between all cells (even of different type) are exactly the same and, thus, no 'static' interface tension should be observed. On the other hand, the anisotropic active growth of the cells have to be considered for the overall stress. We write the stress equations for an anisotropically growing material in 2d under the assumption of incompressibility (see [24]) . Note the islands that detach from the interface and move away. (b) Interface determined by a Voronoi tessellation for the same data as in (a).
Here, -P represents the isotropic part of the stress tensor, η is the shear viscosity, k g is the growth rate, μ is the strength of the anisotropic growth, and = ab a b Q pp 2 1is the nematic-order tensor that is often used in the context of liquid crystals to characterize nematic order (p is a unit vector describing the preferred axis of growth; in the simulations it is the vector connecting the particles constituting one cell). We consider k g to be a constant since in our simulations all cells always grow until they divide; there is no distinction between the different phases of the cell cycle. Thus, all cells contribute to anisotropic part of the stress, independent of the actual division rate. Hence we shorthand m m ¢ = k g in the remainder. In the integral (see equation (11)) all but the anisotropic part of the active stress drop out. We thus arrive at an interface tension, which is active in nature because it is caused by the anisotropic active growth of the cells.
We measured the difference in the principal stresses s s xx yy as well as the order parameter Q xx as a function of the distance s to the interface in our simulations (see figures 5(a) and (b)). Both seem to show the same behavior. Neglecting shear viscosity and assuming ( ) to be proportional to the concentration profile ( ) j s A , where m¢ c are the directed growth coefficients that differ for both tissues, we write the difference in principal stresses as A least-squares fit of s s yy xx as a function of Q xx (s) and ( ) j s A results in coefficients m¢ c that are roughly constant for different system sizes L, but differ considerably for different growth force strengths G, as would be expected (see SI). This leads to an effective surface tension that arises from the directed growth of cells. We find the surface tension to grow with the homeostatic stress difference, while it is roughly independent of the system size (see figures 5(c) and SI figure S7).
It has been shown previously that velocity gradients can align cell divisions [25] and, indeed, we find Q xx to be proportional to ¶ v x x . The proportionality constant, however, depends on the absolute value of the stress. Similar to the island formation, vanishing homeostatic stress difference leads to vanishing surface tension, but diverging interface width (see figure 5(d)). With increasing homeostatic stress difference the interface tension increases, and the with first decreases, but than increases again for larger differences in homeostatic stress. It is important to note that the saturated interface width w sat and its scaling is not dominated by the aforementioned interfacial tension Γ. We observe no systematic decrease of w sat with Γ as would be expected for surface tension dominated scaling, but instead a slight increase of the interface width for higher tension (see figure 5(d)).

Summary and conclusions
We have investigated the interface dynamics of two competing tissues with different homeostatic stresses. In the context of the continuum theory of growing tissues, this difference leads to a take-over of the tissue with the lower homeostatic stress [9]. We have shown by analytical calculations for the interface dynamics in one dimension that an interface propagates at a constant velocity even in the simple case of vanishing diffusion. The interface velocity depends linearly on the difference in the homeostatic stresses of the tissues.
We employ an established simulation method to model the competition of tissues growing on a substrate in two-dimensions. We find the analytical predictions to capture the stress and velocity profiles of our simulations qualitatively and even quantitatively. For the latter, a few tissue characteristics are obtained from independent single-tissue simulations. Furthermore, the linear take-over and, in particular, the linear dependence of the interface velocity on the homeostatic stress difference are reproduced in our simulations as well.
The initially flat interface profile grows in witdth with a power law for short times, saturating on longer time scales. This saturation width obeys a power-law growth with the system size. Lacking published experimental data of tissue competition on substrates, we compare the scaling with colonies growing freely, i.e. not in competition with second tissue. The obtained exponents clearly differ significantly from those previously determined for free in vitro growth of 15 different cell lines as well as in vitro growth of 16 different tumors [22] (b = 3 8, a = 3 2). In order to calculate the roughness exponent for real tumors, slices through surgically removed tumors were analyzed. Furthermore, the roughness exponent found in freely growing bacterial colonies [26] (a » 0.78) is still much larger than the one we estimated. Instead, our results are more similar to those obtained by a cellular automaton model of a single growing tissue [21] which shows Kardar-Parisi-Zhang (KPZ) scaling [27] (b = 1 3, a = 1 2). KPZ is one of the standard theories for growing interfaces, like in ballistic deposition.More recent experiments on growing microbial colonies [28] also indicate a KPZ scaling. However, it is important to note that the typical length scale on which the stress changes across the interface in our simulations is of the order of l » R 10, such that even the largest feasible simulations with = L R 400 are only one order of magnitude above this length scale.
Although cell-cell interactions in our simulations are identical even between different cell types, a non-zero interface tension is found in our simulations. This interface tension can be understood qualitatively by the anisotropic active cell growth. It is thus fundamentally different from thermodynamic interface tension and results from anisotropic active stresses at the interface. Futhermore, our findings show that this interface tension does not control the scaling of the saturation width of the interface.
These findings suggest that competition experiments provide a powerful tool to study the growth dynamics of tissues. In particular we find that interface velocity, velocity decay length, island formation, and the stress profile are to a large degree independent of the system size. We show that they depend in a well-defined, algebraic maner on tissue properties (homeostatic stress, background friction, growth coefficient). Traction force microscopy of growing and competing tissue colonies [29][30][31] could shed light on the role of homeostatic pressure on growth.
Our work brings the by-now classical field of nonequilibrium interface dynamics [20] to intrinsically growing systems. It results in novel dynamical interface features yet to be studied on a fundamental level.