Dynamics of anisotropic tissue growth

We study the mechanics of tissue growth via cell division and cell death (apoptosis). The rearrangements of cells can on large scales and times be captured by a continuum theory which describes the tissue as an effective viscous material with active stresses generated by cell division. We study the effects of anisotropies of cell division on cell rearrangements and show that average cellular trajectories exhibit anisotropic scaling behaviors. If cell division and apoptosis balance, there is no net growth, but for anisotropic cell division the tissue undergoes spontaneous shear deformations. Our description is relevant for the study of developing tissues such as the imaginal disks of the fruit fly Drosophila melanogaster, which grow anisotropically.


Introduction
Multicellular organisms develop from a fertilized egg by repeated rounds of cell divisions. As a result of this process, spatial cell packings are formed, mostly in two-dimensional (2D) sheets such as epithelia, or in 3D mesenchymal tissues. Besides cell division, cell death (apoptosis) can play an important role during development [1,2]. As the tissue grows by cell divisions, its proper size and patterning are ensured by cell signaling events which confer specific cell identities at certain positions in space. A prominent example of such signaling molecules involved in tissue patterning are morphogens, which are secreted from localized sources and form graded concentration profiles in the target tissues. Thus, positional information is conveyed to target cells based on their position within the morphogen concentration gradient [3]- [5]. Understanding the precise mechanisms which underlie tissue patterning and growth, however, remains an important challenge in the field of developmental biology [6]- [9].
Developing tissues can be considered as soft materials with visco-elastic properties [10]. Since active processes take place in cells, such as the cytoskeleton dynamics and cell division, such tissues can be described as active complex fluids. The generic physical properties of such active fluids have recently been discussed in the hydrodynamic limit by continuum descriptions [11,12]. The most striking feature of a developing tissue that results from active processes is growth. Cells undergo a cell cycle during which they double in size and then divide into two daughter cells. Cell division thus involves forces and mechanical work performed to move the neighboring cells in order to create space for the newly produced daughter cells [13]- [16]. In addition, apoptosis leads to the removal of cells, and the liberated space is then occupied by neighboring cells. As a consequence of cell division and apoptosis, cells move relative to each other so that cellular packings are remodeled and cells change their nearest neighbors [17].
Important model systems for the study of tissue growth and patterning during development are the imaginal disks of the fruit fly Drosophila melanogaster. These essentially 2D larval structures are precursors of adult organs, such as the wings, legs or eyes [18,19]. It has been observed that cell division in the imaginal disks is oriented. This implies that the cell division axis has a preferred orientation and can be characterized by an angular distribution [20,21]. An important open question is the role of oriented cell division for shape changes of a growing tissue.

3
Here, we present a coarse-grained physical description of cell movements in growing tissues which takes into account oriented cell division and the physical properties of cells in the tissue, in particular, tissue viscosity. By first developing a continuum description, we study the general features of the system and discuss the flow profiles of cell movements that result from oriented cell division and apoptosis in 2D epithelia. We furthermore perform numerical simulations of a discrete model in which individual cells are described as elastic objects that can slide relative to each other subject to friction forces. In the discrete model, fluctuations resulting from cell division and apoptosis are taken into account. Finally, we use the numerical results for the discrete cellular system to test the applicability of the continuum description.

Continuum description
Our continuum theory of anisotropic tissue growth is based on balances of cell number and forces, and extends earlier approaches [13,14] by taking anisotropic stresses into account. This description allows us to study the flow profiles of cell rearrangements in growing tissues with oriented cell division.

Balances of cell number and forces
We introduce the cell density ρ(r, t) at position r and time t as the number of cells per area (volume) in d = 2 (d = 3) dimensions within an area (volume) element. These elements are large compared to the size of a single cell, but small compared to the size of the tissue. The velocity v(r, t) is defined as the averaged velocity of cells situated in the corresponding area (volume) element at r. In a growing tissue, cell number balance is given by Here, k g (ρ, r, t) and k a (ρ, r, t) are the growth and apoptosis rates which account for cell division and cell death, respectively. In general, both rates can depend on cell density, space and time. Forces in the tissue are captured by the stress tensor σ ik and pressure P, which depends on cell density according to P(r, t) = χ(ρ(r, t) − ρ p )/ρ p , where χ is the bulk elastic modulus and ρ p a reference cell density. Force balance implies that where f ext i is an external force density. Inertial forces are small as compared to other forces and have been neglected.

Constitutive relation for anisotropic tissue growth
The stress tensor is related to the velocity field by a constitutive material relation. We focus on the long time limit where the tissue behaves as a viscous fluid. In addition, we take into account the anisotropic active stresses which are on average generated by oriented cell division. The stress tensor can be written as Here, η and ζ are the shear and bulk viscosity, respectively. The vector p is a unit vector which describes the preferred axis of cell divisions, see figure 1. Hence, σ ik is invariant with respect to p → −p. The anisotropic stress results from oriented cell division. Therefore, we assume this stress to be proportional to the growth rate k g . The strength of the anisotropic stress is characterized by the coefficient µ 0 which has units of viscosity. The special case µ = 0 describes isotropic cell divisions. Considering a constant external pressure P ext , we furthermore impose the boundary conditions σ nn = P − P ext and σ nt = 0, where the indices n and t denote the components of the stress tensor normal and tangential to the tissue boundary, respectively. The equations (1)-(4) describe the full dynamics of growing tissues.

Growth of incompressible tissues
We now focus on 2D epithelia. In order to provide some general insights into the anisotropic growth, we consider the incompressible limit in which the cell density is constant, ρ(r, t) = ρ 0 . Furthermore, we do not consider external forces acting on the tissue, i.e. f ext i = 0, and we assume for simplicity that the growth rate k g , the apoptosis rate k a , the preferred orientation of cell division p, the magnitude of the anisotropic stress µ and the viscosities η and ζ are independent of position but could be functions of time. We choose a coordinate system such that cell division is oriented preferentially along the x-axis, i.e. p = e x . In that case, the dynamic equations (1) and (2) together with the constitutive relation (3) are given by Cell division and cell death generate a non-vanishing divergence of the velocity field. The pressure P plays the role of a Lagrange multiplier to impose the condition (5) for ∇ · v. Note that the anisotropic stress of equation (3) disappears in the force balance (6) since it is homogeneous in space so that its divergence vanishes. However, the boundary conditions for the stress given by equation (4) do involve the anisotropy of the stress and depend on p. For arbitrary tissue shapes, the resulting flow field is given by where k 0 = (k g − k a )/2 and k 1 = k g µ/(4η), and the pressure P = P ext + ζ (k g − k a ) which is independent of position. Interestingly, we find that the anisotropy of tissue growth characterized by k 1 is not only determined by the anisotropy of cell division, but also depends on the tissue viscosity η.
In general, the shape of the tissue is deformed under this flow. For the simple case of an elliptical tissue boundary, the shape stays elliptical during growth. However, the lengths l x and l y of the two main axes change with time as l x and l (0) y are the initial lengths. For k g , k a and µ constant in time, the area A of the tissue grows as A(t) = A 0 e (k g −k a )t . The growth rate k g − k a is therefore related to the effective cell doubling time For arbitrary tissue shapes and time independent k g , k a and µ, average cell trajectories follow flow lines that are described by power-laws where (x 0 , y 0 ) denotes a reference position on the trajectory. There are two important cases which we want to discuss: (i) k a = 0, where no apoptosis occurs and (ii) k g = k a , where cell division and cell death are balanced. In case (i), equation (8) also holds for time-dependent k g , because in this case the ratio (k 0 − k 1 )/(k 0 + k 1 ) is independent of time. Several cases can be distinguished. For µ = 0, it follows that k 1 = 0, so that growth is isotropic and flow lines are radial. For 0 < µ < 2η, the tissue grows at a higher rate along the x-axis than along the y-axis. In the special situation, where µ = 2η the rates k 1 and k 0 balance, and thus the tissue grows only in one dimension. Finally, for µ > 2η, the tissue shrinks along the y-axis and grows rapidly along the x-axis. In case (ii), the velocity field is v = (k 1 x, −k 1 y) and the flow lines obey y = y 0 (x/x 0 ) −1 . This result also holds for time-dependent k a . In this case, the 2D tissue undergoes so-called convergence-extension rearrangements which imply a spontaneous shear deformation [22].

Discrete model
In order to test our continuum theory for anisotropic tissue growth, we define a discrete representation of tissue growth in two dimensions following related approaches [16]. This allows us to obtain solutions to equations (1)-(4) in more complex situations where parameters are position-dependent. This discrete model generates robustly the large-scale features of cell rearrangements in growing tissues, but it is not intended to capture details of deformations on the cell scale as, e.g., discussed in [23,24].

Dynamic equations
We represent cells as elastic objects with the center of the i-th cell located at position x i . Static forces between the N cells are described by the potential function where V (r ) is a pair potential for the two cells with centers at distance r = |x i − x j |. The pair potential describes adhesive forces as well as elastic forces which keep cell centers at a preferred distance a. The potential force f = −dV /dr is chosen to be piecewise linear, For r d the force vanishes. It is attractive for a < r < d with strongest attractive force f 1 < 0 at r = b. It is repulsive if r < a with maximal force f 0 > 0 at r = 0. The dynamics is described by balancing potential forces with friction forces that account for tissue viscosity. Neglecting for simplicity differences between compressional and shear viscosities, we write the balance of forces acting on cell i as Here,η denotes the tissue viscosity on the scale of a cell, and the sum is over the cells j which are neighbors of cell i. The neighbors of cell i are defined as the n nearest cells and, in addition, all cells for which i is within the n nearest neighbors. Equation (11) is a discrete form of the force balance equation (6). We can therefore relate the parameters of the discrete model to the parameters and phenomenological quantities used in equations (2) and (3). In the vicinity of the preferred cell density ρ p 2/( √ 3a 2 ), the pressure can be estimated as P − √ 3V /a. The bulk elastic modulus is thus given by χ √ 3/2(V − V /a).

Oriented cell division and apoptosis
Cell division is implemented as a stochastic process. If cell i is dividing at time t, a new cell is created and both cells are positioned with their centers on opposite sides of a circle with radius ε = a/4 and center at the original position x i , see figure 1. The axis of cell division is characterized by the angle ϕ ∈ [ − π/2, π/2] with respect to the x-axis. The cell division angle ϕ is a random variable drawn from a distribution Q(ϕ). In the case of isotropic cell division, Q(ϕ) = 1/π. Anisotropies of cell division are captured by a distribution Q with a peak at a preferred orientation ϕ = 0, given by the x-axis, see figure 1. For simplicity, we use a piecewise constant distribution function to describe anisotropies of cell division with Q(ϕ) = 1/( ϕ) for − ϕ/2 < ϕ < ϕ/2 and Q = 0 otherwise, i.e. ϕ describes the spread of the division angles around the x-axis. The anisotropic repositioning of the cell pair during division induces a force dipole via the potential U (see equations (9) and (10)). On average, these force dipoles generate the anisotropic active stress in the continuum limit described by equation (3). The division and death events occur at stochastic times. Each cell has an internal clock which measures its lifetime t L after which the cell divides with probability p and dies with probability 1 − p. If cell i undergoes apoptosis, it is simply removed from the system. In our simulations, the cellular lifetime obeys a Gaussian probability distribution R(t L ) with averaget L and variance σ L t L . The effective cell doubling time ist 2 =t L ln 2/ln(2 p). The growth rate k g and the apoptosis rate k a in the continuum limit are related to the probability p and the average 7 cell lifetimet L by Therefore, the probability p for cell divisions obeys p = k g /(k g + k a ).
After the repositioning or removal of cells corresponding to division and apoptosis, the system relaxes according to the dynamic equation (11). In the absence of cell division and apoptosis, there are no fluctuations in this model. As a consequence, the system relaxes in this case to a stable configuration with elastic properties. As soon as cell division and apoptosis are introduced, fluctuations appear as a consequence of stochastic events. The growing and dividing system then behaves like a fluid where cells can change their neighbors. The resulting flow field can on large scales be described by the continuum equations introduced above.
We solve the force balance equation (11) numerically for all cells i = 1, . . . , N with stochastic cell division and cell death events drawn randomly from the probability p and the distributions R(t) and Q(ϕ). At each time step, the velocities dx i /dt are calculated from equation (11) using a matrix inversion, and the positions x i of all cells are updated. At the boundary of the tissue, no external forces are imposed. All forces in equation (11) are internal forces. Therefore, the motion of the center of massx = (1/N ) i x i needs to be specified, and x is fixed in our simulations. As parameter values, we choose for all simulations the average cell lifetimet L = 8 h, with standard deviation σ L = 0.5 h, length ratios b/a = 1.25, d/a = 1.75, cell diameter a = 2.6 µm, | f 0 / f 1 | = 10 (see equation (10)), and nearest-neighbor number n = 6. The spread ϕ of the division angles around the x-axis is varied in our simulations and specified in the figures below. Note that the viscosityη as well as the forces f 0 and f 1 are not themselves simulation parameters, because only the ratios f 0 /η and f 1 /η occur in the dynamic equation (11). Together with the values given here, these ratios are specified by the values for the dimensionless parameter ξ =ηa/(t L f 0 ), which are given in the figure captions. Simulations in the absence of apoptosis ( p = 1) start with a single cell at time zero. For simulations with p < 1, we start with an initial configuration of 512 cells generated by isotropic growth.

Results of numerical simulations
We discuss shape changes and average cell trajectories obtained in our growth simulations. We first consider the case where no apoptosis occurs ( p = 1). The anisotropic shape of an epithelium can be characterized by the variances of cell distributions a denotes the components a = x and y of the vector x i . Figure 2 shows the increase of the linear dimensions l x = √ I x x and l y = I yy as functions of time for two different angular variations of cell division ϕ. The tissue grows exponentially with different growth rates k x and k y which depend on ϕ, consistent with equation (7). Indeed, k x + k y k g where k g = ln 2/t 2 , independent of ϕ, see figure 3 (inset). The dependence of k 1 = (k x − k y )/2 as a function of ϕ is shown in figure 3 for three different values of the parameter ξ . Since χ f 0 /a, the dimensionless parameter ξ η/(χt 2 ) characterizes the ratio of growth rate and cellular relaxation rate. The anisotropic component of growth is given by k 1 /k g α(π − ϕ), where α is the slope of the curves in figure 3. The coefficient α(ξ ) depends weakly on ξ in a nonmonotonous manner, see figure 3. Since k 1 k g µ/(4η), we can determine µ/η 4α(π − ϕ) x y π/2 ∆ϕ π/6 π/3 2π/3 5π/6 (k +k )/k x y g π ∆ϕ π/6 0 1.0 π Figure 3. Anisotropic growth rate k 1 = (k x − k y )/2 normalized by the total growth rate k x + k y as a function of ϕ for ξ = 0.00165 (red), ξ = 0.0165 (green) and ξ = 0.165 (blue). No apoptosis occurs ( p = 1). The data points and standard deviations are obtained from ten independent simulations for each set of parameters. The inset shows the total growth rate k x + k y normalized by the rate k g = ln 2/t 2 defined by the average cell doubling timet 2 =t L .
of the continuum limit. Figure 4 represents average trajectories of cells and their descendants in the x y-plane. The double logarithmic plot reveals that the average positions exhibit a power law as described by equation (8). We find that for each set of parameters, the slope of the linear fit is indeed given by k y /k x .
We now consider the effects of apoptosis in tissues with anisotropic cell division. Figure 5 shows the shape changes of tissues for p = 0.5 and two different choices of ϕ. In this case, no net growth occurs since proliferation and apoptosis are balanced. Displayed are the relative x,y are plotted as a function of time t relative to the average lifetimet L . For each value of ϕ, ten independent realizations are shown. The growth rates k x and k y are determined from linear fits (black lines), and ξ = 0.00165 for all simulations. changes in length of the principal axes l x /l (0) x and l y /l (0) y , where l (0) x,y are the corresponding values of the initial configuration. The logarithmic plot shows that the lengths l x and l y grow and shrink exponentially with the rates k x and k y as described by equation (7). The total growth rate k x + k y 0, consistent with equal rates for cell division and cell death. As in the case without cell death, k 1 = (k x − k y )/2 depends linearly on the variation of the cell division angle ϕ (not shown).

Discussion
In summary, we have shown that oriented cell division in a developing organ leads to anisotropic tissue growth. The anisotropy of growth rates depends on biophysical properties of cells, in particular, on tissue viscosity. Our continuum theory of an incompressible tissue predicts flow fields and cell trajectories which describe well the average behaviors observed in stochastic simulations of anisotropic growth. In our simulations, small differences between the observed total growth rate k x + k y and the rate k g defined by the average cell doubling time arise as a result of tissue compressibility.
The parameter values used in our growth simulations are motivated by studies of the wing imaginal disk of the fruit fly. Key parameters are the bulk elastic modulus χ and the viscositȳ η of the 2D tissue. We estimate χ χ 3D h 6 10 −3 N m −1 , where χ 3D 200 Pa is the shear modulus of a cell and h 30 µm is the height of the epithelium. The choice ξ = 0.0165 thus corresponds to a local two-dimensional tissue viscosity ofη 3 N s m −1 . Using the tissue height h, this corresponds to a viscosity of η 3D 10 5 Pa s which is a typical value that has been reported in experiments [10].
Our work shows that oriented cell division has interesting consequences for cell rearrangements in a growing tissue. In the special case where apoptosis balances cell division, the tissue does not grow but spontaneously undergoes a shear deformation similar to the so-called convergence-extension transformations. Oriented cell division can therefore control shape changes of tissues during development.