Deformable Cell Model of Tissue Growth

This paper is devoted to modelling tissue growth with a deformable cell model. Each cell represents a polygon with particles located at its vertices. Stretching, bending and pressure forces act on particles and determine their displacement. Pressure-dependent cell proliferation is considered. Various patterns of growing tissue are observed. An application of the model to tissue regeneration is illustrated. Approximate analytical models of tissue growth are developed.


Introduction
Mathematical and computer modelling of tissue growth is used in various biological problems such as wound healing and regeneration, morphogenesis, tumor growth, etc. Tissue growth can be described with continuous models: partial differentiation equations for cell concentrations (reaction-diffusion equations), Navier-Stokes or Darcy equations for the velocity of the medium, and elasticity equations for the distribution of mechanical stresses.There is a vast literature devoted to these models (see, e.g., [1][2][3][4][5] and references therein).Another approach deals with individual-based models where cells are considered as individual objects.It allows a more detailed description at the level of individual cells though analytical investigation of such models becomes impossible and their numerical simulations are often more involved than for the continuous models.There are various lattice and off-lattice cell-based models which can be regrouped as spherical particle models, cellular automata and deformable cell models.
Spherical particle models originate from molecular dynamics.Each cell is considered as a sphere which can interact with the surrounding cells.The force acting between them depends on the distance between their centers.Additional random and dissipative forces can be introduced, as is the case in dissipative particle dynamics [6].The motion of particles is determined by Newton's second law.Moreover, cells can divide, die and change their type.Different cell types can have different behavior from the point of view of their motion, division, and death.Already in the 1980s, a similar model, with the equation of motion replaced by some algorithm of cell displacement, was used to study embryogenesis [7].The spherical particle method to model morphogenesis was developed in a recent work [8].In [9], one-cell layered tissues were studied under the assumption that there were additional adhesion forces between cells.Hematopoiesis and leukemia development were investigated in [10][11][12][13].Consecutive cell proliferation and differentiation allowed the description of normal and leukemic hematopoiesis in the bone marrow.Dynamics of cell populations and the conditions of leukemia development were studied.In more complex models, cells can be represented as two connected spherical particles [14].
In the cell-center models, cells are considered as polygons with their centers connected by springs, and cell shape determined by Voronoi tessellation.This approach was used for crypt modelling in [15,16] and brain cancer in [17,18].Particle models can be combined with continuous models, ordinary and partial differential equations for various intracellular and extracellular concentrations.This hybrid method was used for modelling of cancer [19,20], hematopoiesis and blood diseases [21][22][23][24] and blood clotting in flow [25,26].In some cases, it is possible to justify transition from particle models to continuous models [27][28][29].
Another approach to tissue growth is based on cellular automata (see [5,30,31] and references therein).This method can be easier to implement and it is better adapted for the transition from discrete to continuous models.On the other hand, it imposes a square lattice geometry which can influence the results.This method is used to study tumor growth [32,33] and cell invasion [34,35], angiogenesis [36,37], tissue engineering [38], and development of avian gut tissue [39].This approach was further developed in the cellular Potts model [30,40] (see also [41,42] and references therein).
Deformable cell models take into account cell geometry and possible deformation due to mechanical forces.The cell membrane can be considered as an ensemble of particles with various forces acting on them (Section 2).This approach is developed for modeling blood flows where blood cells are considered as individual objects (see [43] and references therein).It is close to the spherical particle models though the particles here do not correspond to individual cells.This method gives more detailed information about cell geometry, adhesion and deformation but it is more difficult to realize, especially in the case of cell division (not considered in the works on blood flows).The questions about the direction of cell division and the mechanical properties of the cell wall become of primary importance here.If they are not sufficiently well described, then the model will fail because of geometrical constraints or because of an accumulation of mechanical stresses.
There are various realizations of the deformable cell model, including the immersed boundary method [44], the subcellular element method [45], and the level set method [46] (see [30] and references therein).The deformable cell model coupled with Navier-Stokes equations was used to study the development of epithelial acini [47,48] and ductal tumor [49].The interaction of cell geometry and migration for the development of cancer was studied in [50].A phase field model for collective cell migration was developed in [51].Deformable cell models for plant growth were developed in [52,53].Growing shoots and roots have a particular cell organization where dividing cells are located in outer layers of the meristem.Differentiated cells do not divide.The direction of cell division depends on cell location.A precise description of the location of dividing cells and of the direction of their division allowed us to reproduce a self-similar growing structure of the plant root [52].Directions of cell division and visco-elasto-plastic properties of cell walls are important in modelling root growth with a deformable cell model.
In this work, we continue to develop the deformable cell model of tissue growth.The main assumptions of the model are similar to those in [52] except for the location of dividing cells and the direction of their division.Instead of a particular cell structure and division specific for root and shoot meristem, we consider uniform tissue growth where dividing cells can be located everywhere.The direction of their division is determined by a special algorithm described in the next section.We will describe the deformable cell model in Section 2. Some properties of growing tissues and applications are discussed in Section 3. Section 4 is devoted to an approximate analytical model of tissue growth.

Forces
We begin the description of the model with mechanical forces acting in cellular structures.An individual cell at its mechanical equilibrium (without deformation) represents a regular polygon with n vertices.Initially, it has an area S 0 , the length of the sides l 0 and the angle between any neighboring sides α 0 .These parameters can be different for different cells.We place a particle of mass m at each vertex of the polygon.
If the polygon changes its shape, then there are three forces which tend to return it to its original configuration.They depend on the change of side length, angles and volume.Similar to other deformable cell models, we define them as follows.If the length l of a side changes, then there is a force acting on each of the two vertices located at this side: (Figure 1b).Here, − → τ is the unit vector along the side.If the angle α between two sides changes (Figure 1c), then there is momentum which creates a pair of forces applied to the vertices and acting on each side in the direction perpendicular to it, α 0 = 180 • .Finally, if the cell area S changes, then there is a force (pressure) applied to each side and acting in the direction perpendicular to it (Figure 1d).Let us note that the first terms in ( 1)-( 3) describe elastic force and the second term viscous dissipation.They correspond to visco-elastic properties of biological cells.The problem is considered in dimensionless variables.Coefficients k i are material parameters which determine the resistance to stretching, bending and compression.Parameters µ i represent artificial viscosity introduced in order to avoid numerical instability.The total force − → F i acting on the i-th vertex is a sum of these three forces.The position − → r i of this vertex is determined by the equations where − → v i is the vertex velocity, µ is the coefficient which determines viscous friction from the surrounding medium.
For further description of the model, we denote by l 0 the side length at mechanical equilibrium and by l its length in the process of cell growth and deformation.Due to plastic deformations, l 0 can change.So, we will distinguish the initial value l i 0 which corresponds to the moment of cell appearance, and the current value l 0 , which can be different from l i 0 .We introduce two additional features in the model: If the length of a side becomes twice as large as its initial length l i 0 , then an additional vertex is introduced in the middle of the side.At the moment of side division, the stretching force is preserved.Hence the number of vertices depends on cell deformation.This allows us to better describe cell shape in the case of large deformations.

2.
If the length l of a side becomes sufficiently large, that is l i 0 /l ≤ γ, with some γ, 0 < γ < 1, then the initial length l 0 increases irreversibly in such a way that it satisfies the equality l 0 = γl.In this case, if the stretching force is removed, then the spring will not return to its initial length but to some greater length.This corresponds to the irreversible deformation of the cell wall when it grows.Here, γ is a material parameter characterizing plastic deformation observed for biological tissues.In the simulations presented below, we set γ = 0.8.

Cell Growth
The cell area S 0 (without deformation) of growing cells changes in time with the rate: where σ(t) is some given function.We set σ = 0, if a cell does not grow and σ(t) = const in the case of linear growth.Other growth rates can also be considered [52].When S 0 increases two-fold, the cell divides.The area of the daughter cell (without deformation) returns to the two-fold lower value.When studying root growth, we also introduced cell growth without division [52].

Cell Division
One of the important characteristics of growing tissue is the directions of cell division.If these directions are not coordinated, then cells can become strongly deformed.This contradicts biological observations and can result in failure of the numerical method.We will choose the direction of cell division according to the algorithm explained below.This algorithm is specific for unpolarized cells [47].
Let us recall that cells are represented as polygons.They grow, increasing their volume, and divide into two cells with approximately equal areas.Cell division occurs along a section connected to vertices (diagonal) of the polygon.Only such sections are considered.For each cell, the algorithm of choice of the section consists of two steps.In the first step, we choose all sections which divide the polygon into two parts with approximately equal area.Since these two areas are determined by the vertices of the polygon, in general they are not exactly equal to each other.Therefore, we introduce a maximal possible difference .Hence we choose all sections such that the areas S 1 and S 2 of two resulting post-division polygons satisfy the condition |S 1 − S 2 | < .Here, S 1 + S 2 = S 0 .In the second step of the algorithm, we choose the shortest section among all of those chosen in the first step.Let us note that this algorithm depends on the location of cell nodes.
Thus, we can summarize this algorithm as a choice of the shortest section by dividing the polygon into two subpolygons with approximately equal areas.Due to this algorithm of cell division, cells cannot be strongly deformed (elongated).We will see that it provides stable tissue growth in a wide range of parameters.Other approaches to the direction of cell division were used in [52] in order to model plant root growth and other polarized cells in [47].

All Cells Grow and Divide
Modelling of tissue growth with a deformable cell model implies certain restrictions on the location of dividing cells and on the direction of their division.If these condition are not satisfied, then it will lead to a rapid accumulation of mechanical stresses and to failure of the cellular structure.One possible cell organization was suggested in [52] for the description of root growth.According to the biological observations, dividing cells are located in this case close to the outer surface of the growing root.Directions of cell division are also precisely determined.In this work, we consider a uniform tissue growth where all cells of the tissue can divide.
Let us recall that each cell is characterized by its equilibrium area S 0 and current area S. The equilibrium area grows here as a linear function of time.An example of a growing tissue is shown in Figure 2.Each cell grows and divides.The direction of division is determined by the shortest diagonal algorithm described above.Initially, there are 20 nodes at the surface of the cell in this simulation.Other nodes can be added if the distance between two neighboring nodes exceeds some given value.We will use the values of parameters the algorithm, we choose the shortest section among all of those chosen in the first step.Let us note that this algorithm depends on the location of cell nodes.
Thus, we can summarize this algorithm as a choice of the shortest section by dividing the polygon into two subpolygons with approximately equal areas.Due to this algorithm of cell division, cells cannot be strongly deformed (elongated).We will see that it provides stable tissue growth in a wide range of parameters.Other approaches to the direction of cell division were used in [52] in order to model plant root growth and other polarized cells in [47].

All Cells Grow and Divide
Modelling of tissue growth with a deformable cell model implies certain restrictions on the location of dividing cells and on the direction of their division.If these condition are not satisfied, then it will lead to a rapid accumulation of mechanical stresses and to failure of the cellular structure.One possible cell organization was suggested in [52] for the description of root growth.According to the biological observations, dividing cells are located in this case close to the outer surface of the growing root.Directions of cell division are also precisely determined.In this work, we consider a uniform tissue growth where all cells of the tissue can divide.
Let us recall that each cell is characterized by its equilibrium area S 0 and current area S. The equilibrium area grows here as a linear function of time.An example of a growing tissue is shown in Figure 2.Each cell grows and divides.The direction of division is determined by the shortest diagonal algorithm described above.Initially, there are 20 nodes at the surface of the cell in this simulation.Other nodes can be added if the distance between two neighboring nodes exceeds some given value.We will use the values of parameters  the algorithm, we choose the shortest section among all of those chosen in the first step.Let us note that this algorithm depends on the location of cell nodes.Thus, we can summarize this algorithm as a choice of the shortest section by dividing the polygon into two subpolygons with approximately equal areas.Due to this algorithm of cell division, cells cannot be strongly deformed (elongated).We will see that it provides stable tissue growth in a wide range of parameters.Other approaches to the direction of cell division were used in [52] in order to model plant root growth and other polarized cells in [47].

All Cells Grow and Divide
Modelling of tissue growth with a deformable cell model implies certain restrictions on the location of dividing cells and on the direction of their division.If these condition are not satisfied, then it will lead to a rapid accumulation of mechanical stresses and to failure of the cellular structure.One possible cell organization was suggested in [52] for the description of root growth.According to the biological observations, dividing cells are located in this case close to the outer surface of the growing root.Directions of cell division are also precisely determined.In this work, we consider a uniform tissue growth where all cells of the tissue can divide.
Let us recall that each cell is characterized by its equilibrium area S 0 and current area S. The equilibrium area grows here as a linear function of time.An example of a growing tissue is shown in Figure 2.Each cell grows and divides.The direction of division is determined by the shortest diagonal algorithm described above.Initially, there are 20 nodes at the surface of the cell in this simulation.Other nodes can be added if the distance between two neighboring nodes exceeds some given value.We will use the values of parameters   The first division occurs along a diameter of the circular cell.The second division takes place along the perpendicular diameter.These first two divisions form the axis of the growing cell structure.They determine the next several divisions.However, during further growth of the tissue, the axes are not preserved and radial structure begins to emerge.The outer cells divide either in the direction parallel to the outer surface or perpendicular to it.This cell organization becomes more obvious during further growth (Figures 3 and 4a).Let us also note that the shapes of the cells represent mostly polygons with four or five vertices.The cell structure symmetry is lost after several divisions because cells have approximately equal area but not exactly equal.Therefore, even after the first division, the two daughter cells are slightly different.

Pressure-Dependent Proliferation
Cell division leads to cells' compression and to internal pressure growth.Initially, internal pressure in a separate cell equals zero because the cell has its equilibrium shape.When it grows, its equilibrium area S 0 increases linearly in time.Its actual area also grows but slower because the cell membrane resists elongation.This difference creates internal cell pressure.When the cell divides, its equilibrium area S 0 is divided by 2 and its actual area is divided approximately by 2. Therefore, the value of pressure remains approximately the same.Moreover, since cells are attached to each other, they cannot return to the circular shape for which the area is maximal.Cell deformation increases the pressure.
Maximal pressure (among all cells) is reached at the center of the domain filled by cells.It grows exponentially in time (Figure 5a) because cells at the center need to push surrounding cells more and more in order to expand.This is related to the exponential growth of cell number (Section 4).For a fixed time, the maximal pressure depends on parameters.It grows approximately as a linear function of viscosity (Figure 5b).It can be compared with fluid flow in a porous medium where pressure difference is proportional to the viscosity (Section 4).
They determine the next several divisions.However, during further growth of th are not preserved and radial structure begins to emerge.The outer cells divide eith parallel to the outer surface or perpendicular to it.This cell organization becomes mo further growth (Figure 4 Left).Let us also note that the shapes of the cells represen with four or five vertices.The cell structure symmetry is lost after several divisi have approximately equal area but not exactly equal.Therefore, even after the first daughter cells are slightly different.

Pressure-Dependent Proliferation
Cell division leads to cells' compression and to internal pressure growth.Initially in a separate cell equals zero because the cell has its equilibrium shape.When it grow area S 0 increases linearly in time.Its actual area also grows but slower because th resists elongation.This difference creates internal cell pressure.When the cell divid area S 0 is divided by 2 and its actual area is divided approximately by 2. There pressure remains approximately the same.Moreover, since cells are attached to each return to the circular shape for which the area is maximal.Cell deformation increas Maximal pressure (among all cells) is reached at the center of the domain filled exponentially in time (Figure 5 Left) because cells at the center need to push surrou and more in order to expand.This is related to the exponential growth of cell numbe fixed time, the maximal pressure depends on parameters.It grows approximately a of viscosity (Figure 5 Right).It can be compared with fluid flow in a porous mediu difference is proportional to the viscosity (Section 4).Biological cells cannot divide when they are too compressed (see [54] and re In order to describe this effect, we introduce the maximal pressure beyond which c Let us recall that pressure P inside cells is determined by Equation (3).We imp condition that if P is greater than some critical value p max , then the cell does not div shows a snapshot of the growing cells' structure for the maximal division pre We obtain a regular structure close to a circle with radial cell organization near Cells there divide either in the direction parallel to the boundary or perpendicula They determine the next several divisions.However, during further growth of are not preserved and radial structure begins to emerge.The outer cells divide ei parallel to the outer surface or perpendicular to it.This cell organization becomes m further growth (Figure 4 Left).Let us also note that the shapes of the cells represe with four or five vertices.The cell structure symmetry is lost after several div have approximately equal area but not exactly equal.Therefore, even after the fir daughter cells are slightly different.

Pressure-Dependent Proliferation
Cell division leads to cells' compression and to internal pressure growth.Initia in a separate cell equals zero because the cell has its equilibrium shape.When it gro area S 0 increases linearly in time.Its actual area also grows but slower because resists elongation.This difference creates internal cell pressure.When the cell divi area S 0 is divided by 2 and its actual area is divided approximately by 2. The pressure remains approximately the same.Moreover, since cells are attached to eac return to the circular shape for which the area is maximal.Cell deformation incre Maximal pressure (among all cells) is reached at the center of the domain fille exponentially in time (Figure 5 Left) because cells at the center need to push surr and more in order to expand.This is related to the exponential growth of cell numb fixed time, the maximal pressure depends on parameters.It grows approximately of viscosity (Figure 5 Right).It can be compared with fluid flow in a porous medi difference is proportional to the viscosity (Section 4).Biological cells cannot divide when they are too compressed (see [54] and In order to describe this effect, we introduce the maximal pressure beyond which Let us recall that pressure P inside cells is determined by Equation (3).We im condition that if P is greater than some critical value p max , then the cell does not d shows a snapshot of the growing cells' structure for the maximal division p We obtain a regular structure close to a circle with radial cell organization ne Cells there divide either in the direction parallel to the boundary or perpendicu hen they are too compressed (see [54] and references therein).troduce the maximal pressure beyond which cells do not divide.cells is determined by Equation ( 3).We impose an additional me critical value p max , then the cell does not divide.Figure 4 Left cells' structure for the maximal division pressure p max = 10. to a circle with radial cell organization near its outer border.ction parallel to the boundary or perpendicular to it.The latter outer layer and provides growth of the perimeter.Similar cell ss-section of growing plant branches.Biological cells cannot divide when they are too compressed (see [54] and references therein).In order to describe this effect, we introduce the maximal pressure beyond which cells do not divide.Let us recall that pressure P inside cells is determined by Equation (3).We impose an additional condition that if P is greater than some critical value p max , then the cell does not divide.Figure 4a shows a snapshot of the growing cells' structure for the maximal division pressure p max = 10.We obtain a regular structure close to a circle with radial cell organization near its outer border.Cells there divide either in the direction parallel to the boundary or perpendicular to it.The latter increases the number of cells in the outer layer and provides growth of the perimeter.Similar cell structures can be observed in the cross-section of growing plant branches.
Computation 2016, 5, x   ivision pressure can influence the properties of the growing cell wn in Figure 4 Middle.Here, the maximal division pressure is set at d when p max = 5 (otherwise, the first divisions do not occur).In this divide even at the outer surface.This results in the appearance of the th a curved surface and "leaves" separated by an interface.Cells at he pressure there is greater than the critical value p max .structures can be considered as instability of symmetric circular ts from cell competition for space.More compressed cells do not row and take even more space.As happens in many other examples, tinction limit where the cell structure stops growing because of the (c) The value of the maximal division pressure can influence the properties of the growing cell structure.Another example is shown in Figure 4b.Here, the maximal division pressure is set at p max = 4 after a short initial period when p max = 5 (otherwise, the first divisions do not occur).In this case, there are cells which do not divide even at the outer surface.This results in the appearance of the specific asymmetric structure with a curved surface and "leaves" separated by an interface.Cells at the interface do not divide since the pressure there is greater than the critical value p max .
The appearance of curved structures can be considered as instability of symmetric circular structures.This instability results from cell competition for space.More compressed cells do not grow, and less compressed cells grow and take even more space.As happens in many other examples, the instability occurs near the extinction limit where the cell structure stops growing because of the excessively high pressure.
Computation 2016, 5, x  aximal division pressure can influence the properties of the growing cell ple is shown in Figure 4 Middle.Here, the maximal division pressure is set at tial period when p max = 5 (otherwise, the first divisions do not occur).In this h do not divide even at the outer surface.This results in the appearance of the cture with a curved surface and "leaves" separated by an interface.Cells at de since the pressure there is greater than the critical value p max .curved structures can be considered as instability of symmetric circular ity results from cell competition for space.More compressed cells do not d cells grow and take even more space.As happens in many other examples, ar the extinction limit where the cell structure stops growing because of the e. ore example where cells die by apoptosis if they do not divide during some e check points during the cell cycle.In this case, they are removed from the In this case, there is a curved front of cells slowly propagating outwards and nside (Figure 4 Right).This is a non-stationary structure where cells divide enough space; they stop dividing and disappear if there is not enough space.We consider one more example where cells die by apoptosis if they do not divide during some time.This is related to the check points during the cell cycle.In this case, they are removed from the computational domain.In this case, there is a curved front of cells slowly propagating outwards and numerous cell clusters inside (Figure 4c).This is a non-stationary structure where cells divide and expand if they have enough space; they stop dividing and disappear if there is not enough space.

External Supply of Nutrients
If nutrients come from outside of the growing tissue, as is the case for tumor growth, then exterior cells have more nutrients than the cells inside the tissue.They will grow and divide faster than cells inside the tissue.They will also compete for nutrients.This competition can lead to nonuniform tissue growth and to the emergence of asymmetric structures.This effect was also observed in hybrid models of tumor growth [55].
In this section, we illustrate it with the deformable cell model.We assume that the growth rate of exterior cells is proportional to the length of their exterior boundary.This assumption implies that the concentration of nutrients in the medium surrounding the tissue is constant and their consumption by cells is proportional to the length of the cell boundary.
If we assume that nutrients can be transported inside the tissue, then interior cells will also have some nutrients.In order to simplify the model, we will suppose that the concentration of nutrients in the extracellular matrix at some point inside the tissue is the same as in the cell located at this point.The flux of nutrients q between cells is proportional to the difference of concentrations c 1 and c 2 in the neighboring cells and to the length L of the boundary between them, q = λ(c 1 − c 2 )L.Here, λ is a parameter which determines the flux intensity.
Figure 6 shows examples of numerical simulations of tissue growth with an external supply of nutrients only to the outer cells and to the interior cells in the case of nutrient diffusion in the tissue.We observe asymmetric growth which results from cell competition for nutrients.A cell which has a longer exterior boundary will grow and divide faster and will compress surrounding cells, reducing their exterior boundary and growth rate.Moreover, the pressure-dependent proliferation increases this competition because compressed cells can stop their growth completely.Emerging structures are more pronounced here than in the model without nutrients (Section 3.2).nutrients only to the outer cells and to the interior cells in the case of nutrient diffusion in the tissue.We observe asymmetric growth which results from cell competition for nutrients.A cell which has a longer exterior boundary will grow and divide faster and will compress surrounding cells, reducing their exterior boundary and growth rate.Moreover, the pressure-dependent proliferation increases this competition because compressed cells can stop their growth completely.Emerging structures are more pronounced here than in the model without nutrients (Section 3.2).

Wound Healing
Together with tumor growth, wound healing is one of the main applications of tissue growth models.An extensive literature review on this subject can be found in [56].A conventional approach to study wound healing consists in the analysis of travelling wave solutions of the reaction-diffusion system of equations for the cell density and for the concentration of some substances which influence cell proliferation (e.g., oxygen).The wave of cell proliferation fills the nutrients only to the outer cells and to the interior cells in the case of nutrient diffusion in the tissue.We observe asymmetric growth which results from cell competition for nutrients.A cell which has a longer exterior boundary will grow and divide faster and will compress surrounding cells, reducing their exterior boundary and growth rate.Moreover, the pressure-dependent proliferation increases this competition because compressed cells can stop their growth completely.Emerging structures are more pronounced here than in the model without nutrients (Section 3.2).

Wound Healing
Together with tumor growth, wound healing is one of the main applications of tissue growth models.An extensive literature review on this subject can be found in [56].A conventional approach to study wound healing consists in the analysis of travelling wave solutions of the reaction-diffusion system of equations for the cell density and for the concentration of some substances which influence cell proliferation (e.g., oxygen).The wave of cell proliferation fills the

Wound Healing
Together with tumor growth, wound healing is one of the main applications of tissue growth models.An extensive literature review on this subject can be found in [56].A conventional approach to study wound healing consists in the analysis of travelling wave solutions of the reaction-diffusion system of equations for the cell density and for the concentration of some substances which influence cell proliferation (e.g., oxygen).The wave of cell proliferation fills the spatial domain, which corresponds to the wound and which is surrounded by the remaining tissue.In this section, we illustrate modelling of wound healing with the deformable cell model.We will see that pressure-dependent proliferation is an appropriate mechanism to control wound closure.Figure 7 shows an example of numerical simulations of tissue regeneration.The original tissue (left image) is obtained from a single cell as described above.The cells of this tissue are fixed-they do not grow or divide any more.After that, we remove the internal part of the tissue.The cells at the boundary of the removed tissue can grow and divide.They fill the available place.They stop growing and divide when the tissue is regenerated due to pressure-dependent proliferation.Cells push each other, internal pressure grows and they no longer grow and divide.The new tissue fits the form of the wound well.However, some damage remains at the place where two parts of the tissue meet each other.Its size depends on the mechanical properties of cells.It decreases if cells are more deformable and better fit each other.
Another example of tissue regeneration is shown in Figure 8.The main difference here in comparison with the previous example is that the region where a part of the tissue is removed is not bounded by the remaining tissue (Figure 8b).If we use the same approach as in the previous case, we obtain a different form of the regenerated tissue compared with the initial tissue (Figure 8c).In order to reproduce the form of the original tissue, we need to introduce different growth rates for different cells.Consider cells located at the wound surface in Figure 8b.We impose a particular growth rate for each of these cells.Their descendants will preserve these growth rates.The cells located closer to the outer surface grow and divide slowly.The growth rate and division increase as much as the cells go deeper into the damaged tissue.The maximal growth rate is reached at the center.Such distribution of the growth rate allows us to approximate the initial tissue form by the regenerated tissue (Figure 8d).It should also be noted that this tissue continues to grow even when the initial form is regenerated.Therefore, there are additional mechanisms which control the size of the tissue [57][58][59].
not grow or divide any more.After that, we remove the internal part of the tissue.The cells at the boundary of the removed tissue can grow and divide.They fill the available place.They stop growing and divide when the tissue is regenerated due to pressure-dependent proliferation.Cells push each other, internal pressure grows and they no longer grow and divide.The new tissue fits the form of the wound well.However, some damage remains at the place where two parts of the tissue meet each other.Its size depends on the mechanical properties of cells.It decreases if cells are more deformable and better fit each other.Another example of tissue regeneration is shown in Figure 8.The main difference here in comparison with the previous example is that the region where a part of the tissue is removed is not bounded by the remaining tissue (Figure 8b).If we use the same approach as in the previous case, we obtain a different form of the regenerated tissue compared with the initial tissue (Figure 8c).In order to reproduce the form of the original tissue, we need to introduce different growth rates for different cells.Consider cells located at the wound surface in Figure 8b.We impose a particular growth rate for each of these cells.Their descendants will preserve these growth rates.The cells located closer to the outer surface grow and divide slowly.The growth rate and division increase as much as the cells go deeper into the damaged tissue.The maximal growth rate is reached at the center.Such distribution of the growth rate allows us to approximate the initial tissue form by the regenerated tissue (Figure 8d).It should also be noted that this tissue continues to grow even when the initial form is regenerated.Therefore, there are additional mechanisms which control the size of the tissue [57][58][59].
Let us note that if the critical value p max is high enough, it will not influence cell division.On the other hand, there is a biochemical regulation of the cell division rate, implying that some cells divide faster than others.This regulation is not explicitly introduced in the model but it acts implicitly through the cell location.not grow or divide any more.After that, we remove the internal part of the tissue.The cells at the boundary of the removed tissue can grow and divide.They fill the available place.They stop growing and divide when the tissue is regenerated due to pressure-dependent proliferation.Cells push each other, internal pressure grows and they no longer grow and divide.The new tissue fits the form of the wound well.However, some damage remains at the place where two parts of the tissue meet each other.Its size depends on the mechanical properties of cells.It decreases if cells are more deformable and better fit each other.Another example of tissue regeneration is shown in Figure 8.The main difference here in comparison with the previous example is that the region where a part of the tissue is removed is not bounded by the remaining tissue (Figure 8b).If we use the same approach as in the previous case, we obtain a different form of the regenerated tissue compared with the initial tissue (Figure 8c).In order to reproduce the form of the original tissue, we need to introduce different growth rates for different cells.Consider cells located at the wound surface in Figure 8b.We impose a particular growth rate for each of these cells.Their descendants will preserve these growth rates.The cells located closer to the outer surface grow and divide slowly.The growth rate and division increase as much as the cells go deeper into the damaged tissue.The maximal growth rate is reached at the center.Such distribution of the growth rate allows us to approximate the initial tissue form by the regenerated tissue (Figure 8d).It should also be noted that this tissue continues to grow even when the initial form is regenerated.Therefore, there are additional mechanisms which control the size of the tissue [57][58][59].
Let us note that if the critical value p max is high enough, it will not influence cell division.On the other hand, there is a biochemical regulation of the cell division rate, implying that some cells divide faster than others.This regulation is not explicitly introduced in the model but it acts implicitly through the cell location.not grow or divide any more.After that, we remove the internal part of the tissue.The cells at the boundary of the removed tissue can grow and divide.They fill the available place.They stop growing and divide when the tissue is regenerated due to pressure-dependent proliferation.Cells push each other, internal pressure grows and they no longer grow and divide.The new tissue fits the form of the wound well.However, some damage remains at the place where two parts of the tissue meet each other.Its size depends on the mechanical properties of cells.It decreases if cells are more deformable and better fit each other.Another example of tissue regeneration is shown in Figure 8.The main difference here in comparison with the previous example is that the region where a part of the tissue is removed is not bounded by the remaining tissue (Figure 8b).If we use the same approach as in the previous case, we obtain a different form of the regenerated tissue compared with the initial tissue (Figure 8c).In order to reproduce the form of the original tissue, we need to introduce different growth rates for different cells.Consider cells located at the wound surface in Figure 8b.We impose a particular growth rate for each of these cells.Their descendants will preserve these growth rates.The cells located closer to the outer surface grow and divide slowly.The growth rate and division increase as much as the cells go deeper into the damaged tissue.The maximal growth rate is reached at the center.Such distribution of the growth rate allows us to approximate the initial tissue form by the regenerated tissue (Figure 8d).It should also be noted that this tissue continues to grow even when the initial form is regenerated.Therefore, there are additional mechanisms which control the size of the tissue [57][58][59].
Let us note that if the critical value p max is high enough, it will not influence cell division.On the other hand, there is a biochemical regulation of the cell division rate, implying that some cells divide faster than others.This regulation is not explicitly introduced in the model but it acts implicitly through the cell location.Let us note that if the critical value p max is high enough, it will not influence cell division.On the other hand, there is a biochemical regulation of the cell division rate, implying that some cells divide faster than others.This regulation is not explicitly introduced in the model but it acts implicitly through the cell location.

Approximate Model
Passing from the equation of motion of individual particles to an approximate averaged equation of motion of the medium, we get ∂v ∂t where v is the velocity of the medium and p is the pressure (see [27] for more detail).We consider this equation in the one-dimensional case in order to obtain an approximate analytical solution.
It corresponds to the two-dimensional case with plane geometry or to the circular tissue with a sufficiently large radius.Equation ( 5) is a macroscopic analogue of Equation ( 4) with an average force proportional to the pressure gradient.It should be noted, however, that force in (4) also contains other components.So, the analytical continuous model is an approximation.Parameter µ in both equations corresponds to viscous friction.We will consider this equation in the quasi-stationary approximation, which takes the form of the Darcy equation for the porous medium.Next, we will assume that the medium is incompressible.This assumption is justified for sufficiently large values of the pressure force coefficient k 3 in the cell model.In this case, cells preserve their volume and, in the limit of large coefficients, the medium becomes incompressible.In this case, the mass conservation (or continuity) equation can be written as follows: where W is the rate of production of mass.In the one-dimensional case, divergence is the space derivative ∂v/∂x.Production of mass is determined by the rate of cell division W = f (p).We will suppose that it depends on the pressure.When biological cells are strongly compressed, they do not divide.This assumption is consistent with pressure-dependent proliferation in the cell model (Section 3), which can be considered as discretization of the continuous model.From Equations ( 6) and ( 7), we get where prime denotes the space derivative.

Constant Proliferation
If the proliferation rate does not depend on pressure, then we set f (p) = b, where b is some positive constant.We solve Equation (8) in the interval −ξ < x < ξ with the boundary condition p(±ξ) = p 0 .We obtain Since dξ/dt = v(ξ), then taking into account (6), we have Hence the length of the interval grows exponentially.If we consider a circular tissue with the radius ξ, then neglecting the influence of the curvature on the growth rate, we get the formula for the tissue area A = πξ 2 : where A 0 is the initial area.Figure 9a shows an example of numerical simulations of tissue growth with a constant proliferation rate.The linear dependence of the logarithm of tissue area on time, log A(t), shows that the growth rate is exponential.Three different proliferation rates are considered, dS 0 (t)/dt = 0.1, 0.05, 0.025.The slopes of the straight lines are in the same proportions.For large amounts of time, the assumptions of the model may not be satisfied, and the growth rate slows down in comparison with the exponential one.
Let us suppose that the medium fills the half-axis x ≤ ξ.We should complete the formulation of the problem by the boundary condition for the pressure where p 0 is some non-negative constant, p 0 < p m .
We can now solve problem ( 8) and (9).Suppose that the value p = p m is reached at x = 0. Taking into account the continuity of pressure and of its first derivative at x = 0, we obtain where α = bµ/(2k).From the boundary condition ( 9), we find the value of ξ: In the case of circular tissue growth, neglecting curvature, we can approximate the growth rate of irs radius R by the formula dR dt = v.
Hence we obtain the formula for the tissue area A = πR 2 : Figure 9b shows the results of numerical simulations of tissue growth and comparison with the analytical approximation.At the beginning of growth, when the pressure-dependent proliferation is not yet established, all cells divide, and the growth rate is exponential.Therefore, instead of Formula ( 10), we approximate it by the expression where the coefficients a and b are chosen to fit the simulations presented in Section 3.1 with two different values of maximal pressure: p m = 5 and p m = 6.The curves are well approximated by Formula (11) with a = 0.0033, b = −1142 and a = 0.0045, b = −2227.
In order to compare more precisely the values of a obtained in numerical simulations with the coefficient a = π 2bk(p m − p 0 ) µ (12) in (10), we need to specify the value of pressure p 0 at the boundary.This is the pressure in the cells at the outer layer.It varies during cell growth.Its average value depends on the values of parameters.
For the set of p 0 ≈ 4. Hence the right-hand side in (12) increases twice when we compare p m = 5 and p m = 6 while the coefficient a on the left-hand side obtained from numerical simulations in 0.0045/0.0033≈ 1.5 times.This difference is related to the fact that intracellular pressure in the numerical model (Section 3) is only a part of the total pressure in the analytical model.The latter also takes into account the contribution from elastic forces.If we decrease the value of the stretching coefficient, k 1 = 30 (instead of k 1 = 40 above), then pressure in cells at the outer layer decreases.Therefore, growth accelerates and we get a = 0.004.This example shows how mechanical properties of cells influence the growth rate of the tissue and how they influence macroscopic phenomenological parameters.

Discussion
In this work, we developed a deformable cell model of growing tissue where cells are attached to each other.Two neighboring cells have the same boundary and cannot be separated.Another possible approach is to consider separated cells which can adhere due to additional forces acting between them [47,48].The direction of cell division is chosen along the shortest diagonal which divides the cell into two cells with approximately equal areas.This algorithm appears to be very robust and allows us to simulate tissue growth with more than 10 4 cells.Pressure-dependent proliferation, which corresponds to biological observations, can result in the appearance of asymmetric tissue structures.
The spring model of cells (Section 2) can be considered as an approximation of an elastic medium.Together with viscosity in the equation for particle displacement, we obtain a viscoelastic medium.Using the porous medium equation to describe it, we simplify the model but we obtain a good approximation of the tissue growth rate.
Modelling of tissue growth with the deformable cell model can be applied to study tumor growth.In the case of pressure-dependent proliferation, it describes the linear growth rate of the tissue radius specific for tumors after a short initial stage of exponential growth [54].The growing tumor can lose its radial symmetry, resulting in the appearance of spatial structures.This instability can be related to cell competition for nutrients, and it can be modeled with partial differential equations [60][61][62][63][64] or with cell-based models [32,55] (see also references in [54]).However, as is indicated in [54], in vitro experiments show a smoother shape of the growing tumor.It is suggested that another possible mechanism of pattern formation can be related to mechanical stresses: cells diffuse along the surface searching for a free space.Results presented in this work show that indeed, competition for nutrients is not the only possible mechanism of instability.Mechanical constraints lead to the emergence of spatial structures.The introduction of nutrients increases their roughness.Hence pressure-dependent proliferation, a linear growth rate (exponential rate at the beginning) and various mechanisms of pattern formation of growing tissues observed in numerical simulations are in agreement with the experiments in vitro.
Another classic example of tissue growth is wound healing.It is conventionally modelled with reaction-diffusion equations for the concentrations of cells and some substances which determine their proliferation.A wave of cell proliferation propagates in the wound region and fills it.The deformable cell model of tissue growth is an appropriate tool to study wound healing.Even without additional concentrations, pressure-dependent proliferation allows us to control wound closure.When two borders of the wound meet and cells touch each other, their proliferation stops.However, if we consider the wound cross-section perpendicular to the epidermis, the wound area is not bounded by the remaining tissue from one side.Therefore, the problem here is different: how to control the form and the size of the growing tissue in an open domain.In this case, we need to introduce the rate of cell proliferation which depends on its position.This question requires further investigation.
The model presented in this work has certain limitations.First of all, we consider a two-dimensional tissue (cell sheet) while most biological tissues are essentially three-dimensional.Numerical simulations in this case become more involved but we expect that the main properties discussed above remain the same.Cell division and death are determined by the intracellular

Figure 1 .
Figure 1.Schematic representation of the model: a structure which consists of four deformable cells (a); stretching force acts between two neighboring vertices along the line connecting them (b); bending force acts on neighboring sides as the angle between them differs from the angle α 0 (c); pressure acts on all sides in the direction perpendicular to them (d).
8, unless other values are specified.
8, unless other values are specified.

Figure 2 .
Figure 2. Example of simulations where all cells grow and divide; S 0 is a linear function of time; direction of cell division is chosen along the shortest diagonal.First divisions are shown on the left, and a snapshot of the growing cell structure is shown on the right.

Figure 2 .
Figure 2. Example of simulations where all cells grow and divide; S 0 is a linear function of time; direction of cell division is chosen along the shortest diagonal.First divisions are shown on the left, and a snapshot of the growing cell structure is shown on the right.

Figure 2 .
Figure 2. Example of simulations where all cells grow and divide; S 0 is a linear function of time; direction of cell division is chosen along the shortest diagonal.First divisions are shown in (a); and a snapshot of the growing cell structure is shown in (b).

Figure 3 .
Figure 3. Snapshots of growing tissue with the maximal division pressure p max = 15 in t moments of time.Non-dividing cells with pressure greater than p max are shown in dar

Figure 3 .
Figure 3. Snapshots of growing tissue with the maximal division pressure p max = 15 in moments of time.Non-dividing cells with pressure greater than p max are shown in d

Figure 3 .
Figure 3. Snapshots of growing tissue with the maximal division pressure p max = 15 in the consecutive moments of time.Non-dividing cells with pressure greater than p max are shown in dark grey.

Figure 4 .Figure 4 .
Figure 4.A snapshot of the growing cell structure with the maximal division pressure p max = 10 p max = 4 (Middle) and with the disappearance of non-dividing cells (Right).

7 of 17 ing
cell structure with the maximal division pressure p max = 10 (Left), disappearance of non-dividing cells (Right).the growing cell tissue.(Left) maximal pressure as a function of time; e t = 700 for different values of viscosity coefficient µ.Dimensionless d.

Figure 4 .
Figure 4.A snapshot of the growing cell structure with the maximal division pressure p max = 10 (a); p max = 4 (b) and with the disappearance of non-dividing cells (c).

Figure 4 .
Figure 4.A snapshot of the growing cell structure with the maximal division pr p max = 4 (Middle) and with the disappearance of non-dividing cells (Right).

Figure 5 .Figure 5 .
Figure 5. Maximal pressure in the growing cell tissue.(Left) maximal pressu (Right) maximal pressure at time t = 700 for different values of viscosity coeffi time and pressure units are used.

Figure 5 .
Figure 5. Maximal pressure in the growing cell tissue.(a) maximal pressure as a function of time; (b) maximal pressure at time t = 700 for different values of viscosity coefficient µ.Dimensionless time and pressure units are used.

Figure 6 .
Figure 6.A snapshot of growing tissue with an external supply of nutrients only to the outer cells (Left) and with nutrient diffusing inside the tissue (Right).

Figure 6 .
Figure 6.A snapshot of growing tissue with an external supply of nutrients only to the outer cells (Left) and with nutrient diffusing inside the tissue (Right).

Figure 6 .
Figure 6.A snapshot of tissue with an external supply of nutrients only to the outer cells (a) and with nutrient diffusing inside the tissue (b).

Figure 7 .
Figure 7. Numerical simulations of wound healing.Original tissue (Left), damaged tissue (Middle), tissue regeneration (Right).Grey cells are fixed, while white cells divide and fill the missing part of the tissue.

Figure 7 .
Figure 7. Numerical simulations of wound healing.Original tissue (Left), damaged tissue (Middle), tissue regeneration (Right).Grey cells are fixed, while white cells divide and fill the missing part of the tissue.

Figure 7 .
Figure 7. Numerical simulations of wound healing.Original tissue (Left), damaged tissue (Middle), tissue regeneration (Right).Grey cells are fixed, while white cells divide and fill the missing part of the tissue.

Figure 7 .
Figure 7. Numerical simulations of wound healing.Original tissue (a); damaged tissue (b); tissue regeneration (c).Grey cells are fixed, while white cells divide and fill the missing part of the tissue.

Figure 8 .
Figure 8. Numerical simulations of wound healing.Original tissue (a); damaged tissue (b); tissue regeneration with equal rates of cell growth (c); cell growth rate depends on the initial position of the first dividing cell in its lineage (d).