Morphologies of compressed active epithelial monolayers
EEPJ manuscript No. (will be inserted by the editor)
Morphologies of compressed active epithelial monolayers
Jan Rozman , Matej Krajnc , and Primoˇz Ziherl Joˇzef Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, SloveniaReceived: date / Revised version: date
Abstract.
Using a three-dimensional active vertex model, we numerically study the shapes of strained un-supported epithelial monolayers subject to active junctional noise due to stochastic binding and unbindingof myosin. We find that while uniaxial, biaxial, and isotropic in-plane compressive strains do lead to theformation of longitudinal, herringbone-pattern, and labyrinthine folds, respectively, the villus morphologycharacteristic of, e.g., the small intestine appears only if junctional tension fluctuations are strong enoughto fluidize the tissue. Moreover, the fluidized epithelium features villi even in absence of compressive strainprovided that the apico-basal differential tension is large enough. We analyze several details of the differentepithelial forms including the role of strain rate and the modulation of tissue thickness across folds. Ourresults show that nontrivial morphologies can form even in unsupported, non-patterned epithelia.
The function of many animal organs strongly relies on theshape of the constituent epithelia. The foundations of thevarious organ forms are set in the early embryonic devel-opment, where the initially simple structures gain on com-plexity through different types of deformations includingfolding, wrinkling, in-plane stretching, and shear. Drivenby a complex biochemical machinery, these processes areindividually often discussed within the domain of mate-rials science, but they are generally poorly understood inthe biological context where they usually act together, andtheir synergistic and antagonistic effects remain theoreti-cally largely unexplored. This is quite understandable asthe parameter space increases with each additional pro-cess involved, which makes the analysis rather difficult.Among the best-known epithelial morphologies are theintestinal villi. Villi form in embryos, e.g., upon a sequen-tial differentiation of muscle tissues of the gut, which ex-ert compressive stresses on the growing epithelium. Thesestresses trigger buckling and villus formation through afew intermediate steps involving longitudinal and herring-bone (zigzag) fold patterns [1]. So far, the mechanisticaspects of villus formation were interpreted in terms ofseveral distinct processes such as tissue growth in confinedspace [1,2,3] and mechanical apico-basal polarity [4]. Eachof these processes may give rise to villus morphogenesis byitself although it is also possible that they act jointly soas to ensure robustness. However, the relative role andimportance of the different possible mechanisms of villusformation remains an open question because most theo-retical studies focus on a single one.Many processes that contribute to epithelial morpho-genesis can be theoretically investigated by describing the tissue as an elastic or viscoelastic continuum [1,2,3,5] butit is often advantageous to resort to cell-based compu-tational schemes such as agent-based models [6,7]. Theagent-based models build on cells whose state and re-sponse to signals from the environment are controlled bya set of rules and are thus well suited for the descriptionof growth regulation, cell differentiation, spatial arrange-ment of the different cell populations, etc. At the sametime, these schemes assume a very simple physical formof the cells, which are typically approximated by spheres,and they rely on a rather stripped-down representationof the mechanical cell-cell interaction. As a result, theycan accommodate only a few non-specific morphogeneticmodes such as folding due to rapid growth [8,9].A geometrically and mechanically more accurate com-putational framework is provided by vertex models whereeach cell is represented by a polyhedron carrying a certainmechanical energy. Although still fairly coarse-grained,these models allow one to explore the contributions of thedifferent structural elements of cells and tissues as well asto readily introduce various kinds of cell-level dynamics.While it is plausible that cell-resolution description maynot be needed for many epithelial formations, the abovefeatures render the vertex models quite versatile and oftenadvantageous over continuum theories because the ener-gies and processes involved have a reasonably clear mi-croscopic interpretation. In addition, the vertex modelsbypass the very choice of a continuum theory. The choiceof such a theory poses a challenge since the deformationscharacteristic of epithelia are large, and hence even withinthe domain of elastic materials alone there exist a num-ber of possibilities rather than a single one like in thesmall-deformation harmonic limit. At the same time, acontinuum theory derived from a given vertex model in a a r X i v : . [ c ond - m a t . s o f t ] F e b Jan Rozman et al.: Morphologies of compressed active epithelial monolayers bottom-up manner can be quite illuminating at a qualita-tive level as it may provide a connection between the vari-ous geometric quantities that describe the tissue shape [10,11,12].The main topological dynamic processes of interest invertex models are neighbor-exchange T1 transition, celldivision, and cell expulsion, which all result in in-planerearrangements of cells and involve a certain kind of ac-tivity. Cell division generally takes place after the dividingcell has grown, often so as to double in volume, which evi-dently relies on metabolism, synthesis of biomolecules, etc.Irrespective of how the daughters arrange after the divi-sion, the local topology of the tissue is changed; the sameapplies to cell expulsion. As such, these two events con-tribute to tissue fluidization [13]. A third mechanism offluidization is the T1 transition, where four cells arounda given lateral side locally rearrange such that the twothat initially shared the side no longer do so [14]. Thisrearrangement generally requires climbing across an en-ergy barrier [15], which is evidently an active process andcan be implemented by a fluctuating tension representingbinding and unbinding of junctional myosin [16,17], by amore coarse-grained active T1 transition model [14], or insome other manner.In a recent paper [12], we used an active vertex modelof epithelia to study the form of organoid-like epithelialshells, finding that the most elaborated, branched mor-phologies appear only if cell-level activity leading to in-plane cell rearrangements through T1 transitions is strongenough. The well-developed finger-like protrusions in thebranched shapes are physically very reminiscent of the in-testinal villi and crypts. As it appears possible that theseprotrusions also arise in a planar patch of an active ep-ithelium provided that effect of the fixed-enclosed-volumeconstraint can be substituted by in-plane stress, we nowexplore the vertex model of epithelia under in-plane com-pressive strain. The main elements of the model as wellas the representative morphologies obtained including aclose-up of a villus are shown in Fig. 1. We find thatwhile this strain does produce folded morphologies, villior crypts appear only if two additional requirements arefulfilled: (i) the tissue is completely fluidized and (ii) theapico-basal differential tension exceeds a certain thresh-old. We also find that strain is not essential for the forma-tion of the protrusions, which appear even in its absenceas long as the above two requirements are met. Overall,while most previous studies relate the intestinal villi andcrypts to the differential growth of the epithelium rela-tive to supporting tissues [1,2,3], our work demonstratesthat they can form in an unsupported and non-patternedepithelial tissue.The disposition of the paper is as follows: In Sec. 2we present the computational model used, elaborating allthe details of the vertex model as well as our implementa-tion of the tension-fluctuation scheme as the mechanismof active in-plane cell rearrangements. Section 3 describesthe key results whereas Sec. 4 discusses the state diagramand the relevance of our main findings in the context of experiments and existing theories. Section 5 concludes thepaper.
We use a 3D vertex model where the tissue is representedby a monolayer sheet-like packing of N c polyhedral cellsof identical fixed volumes V c [12,14,18,19]. The apical andbasal sides of the cells are polygons with the same num-ber of vertices whereas the lateral sides are rectangular(Fig. 1). The shape of each cell is parametrized by thepositions of vertices r j = ( x j , y j , z j ), which obey the over-damped equation of motion η ˙ r j = −∇ j W + F ( A ) j . (1)Here η is an effective viscosity associated with the move-ment of vertices and W = N c X i =1 (cid:20) Γ a A ( i ) a + Γ b A ( i ) b + Γ l A ( i ) l (cid:21) (2)is the energy of the epithelial sheet; A ( i ) a , A ( i ) b , and A ( i ) l are the areas of the apical, basal, and lateral sides of the i -th cell and Γ a , Γ b , and Γ l are the corresponding surfacetensions, each of them assumed to be identical in all cells;the factor of 1/2 in the third term arises because eachlateral side is shared between two cells.The second term on the right-hand side of Eq. (1), F ( A ) j , is the active force describing the noise at cell-celljunctions which arises from the stochastic turnover dy-namics in the junctional myosin. These dynamics takeplace at the lateral cell-cell contacts, close to the apicalcell surface [20]. Since the lateral sides of our model cellsare not discretized along the apico-basal axis and T1 tran-sitions always occur simultaneously on the apical and thebasal cell side, we associate these fluctuations with theapical and the basal edge of cell-cell junctions such that F ( A ) j = − N e X k =1 γ k ( t ) ∇ j h L ( k ) a + L ( k ) b i . (3)Here L ( k ) a and L ( k ) b are apical and basal edge lengths ofthe k -th cell-cell junction, respectively, and the sum goesover all N e junctions. Junctional line tensions γ k ( t ) aredescribed by the Ornstein-Uhlenbeck processd γ k ( t )d t = − τ m γ k ( t ) + ξ k ( t ) , (4)where 1 /τ m is the turnover rate of molecular motors and ξ k ( t ) is the Gaussian white noise with properties h ξ k ( t ) i =0 and h ξ k ( t ) ξ k ( t ) i = (2 σ /τ m ) δ k,k δ ( t − t ), σ being thelong-time variance of tension fluctuations [12,16,17].As a mechanism of tissue fluidization, junctional ten-sion fluctuations are one of the key elements of our model:If the magnitude of fluctuations exceeds a critical value, an Rozman et al.: Morphologies of compressed active epithelial monolayers 3 longitudinal folds labyrinthine foldsprevilli villi ( uniaxial, β - α = σ = ) ( isotropic, β - α = σ = )( uniaxial, β - α = σ = ) ( uniaxial, β - α = σ = ) compression1 1.5 2 2.5 300.10.20.3 tissue tension α + β t en s i on f l u c t ua t i on s σ D eff ( - ) fluidsolid a ) b ) c ) d ) e ) f ) g ) h ) Fig. 1.
Schematic of a six-coordinated cell in the model epithelial monolayer (a). The fluctuating-tension edges at the apical(top) and the basal (bottom) side of the lateral surface (cyan) are highlighted by thick green lines. In panel b we plot the diffusioncoefficient of the cell as a function of the tissue tension α + β and the magnitude of tension fluctuations σ . Superimposed onthe color-coded plot is the numerically obtained location of the solid-fluid transition (white points) and the white line showsEq. (10). Panel c shows an oblique view of the initial flat configuration, indicating the direction of uniaxial compressive strain.Also included are apical-surface views of representative α + β = 2 examples of the longitudinal (d) and labyrinthine folds (e) aswell as the previlli (f) and villi morphologies (g); the strain mode, the differential tension β − α , and the magnitude of tensionfluctuations σ are listed in each panel. Cells are colored according to the number of neighbors: Blue for 4, yellow for 5, whitefor 6, brown for 7, and green for 8. Lateral cell sides are shown in cyan. Panel h shows a representative σ = 0 . α + β = 2, β − α = 0 . the monolayer undergoes a transition from solid-like tofluid-like state [17]. In the latter, cells can move withinthe monolayer and this motion has a dramatic and oftenfar-reaching effect on the mechanical features as well ason the shape of the tissue. For example, it can triggerlarge-scale deformations and cell rearrangements duringembryonic development [16,21] but it can also be an es-sential factor in the formation of the more elaborated 3Dtissue morphologies in highly constrained systems such asorganoids [12]. We first non-dimensionalize the model. We use the lateraltension Γ l as the unit of tension whereas the natural choiceof the unit of surface area is based on cell volume V c andreads V / c . The dimensionless energy is given by w = N c X i =1 (cid:20) αa ( i ) a + βa ( i ) b + 12 a ( i ) l (cid:21) , (5)where α = Γ a Γ l (6)and β = Γ b Γ l (7) are the reduced apical and basal surface tensions, respec-tively, whereas a ( i ) a = A ( i ) a /V / c is the reduced area of theapical side (and analogously for a ( i ) b and a ( i ) l ). The reducedapical and basal tensions jointly determine the preferredshape of individual cells. In particular, the tissue tension α + β controls the cell width-to-height ratio: For α + β considerably smaller than 1 this ratio is large which corre-sponds to squamous tissues whereas for α + β much largerthan 1 it is small which represents columnar tissues. Onthe other hand, the apico-basal differential tension β − α favors cells of conical shape which may lead to bucklingat the tissue scale. Unless stated otherwise, we report re-sults for columnar model tissues (Fig. 1a) with α + β = 2,giving the cell width-to-height ratio of about 0.5, and wevary the differential tension β − α so as to explore thecontribution of the apico-basal polarity to the formationof various morphologies. Without loss of generality, werestrict the discussion to β > α so that the differentialtension β − α is positive. The obtained non-trivial mor-phologies that feature basal constrictions and protrusionspenetrating into the lumen can be translated into thosecharacterized by apical constrictions and recesses simplyby swapping the sign of the differential tension. For exam-ple, the villus morphology at β − α = 0 . β − α = − .
8, withFig. 1g showing the basal rather than the apical view ofthe latter.
Jan Rozman et al.: Morphologies of compressed active epithelial monolayers
We investigate tissues consisting of N c cells with periodicboundary conditions along the two in-plane directions de-noted by x and y . The in-plane compression is simulatedby decreasing the linear dimensions of the simulation box L x and L y either along a single in-plane axis so as to ge-nerate uniaxial strain or along both of them by the sameamount, which gives rise to isotropic strain. In our model,these dimensions are decreased linearly in time so that L µ ( t ) = L (0) µ (cid:18) − ε µ tT (cid:19) , (8)where L (0) µ is the initial dimension of the box in question, µ is either x or y , T is the total simulation time, and ε µ is thefinal in-plane strain along the corresponding direction. Ofcourse, ε µ corresponds to the physical definition of strainonly if no out-of-plane buckling occurs, which is not thecase in our simulated morphologies. The competition between cell surface tension and fric-tion with the environment determines the typical timescale of vertex movements given by τ = η/Γ l , which ishere used as the unit of time. After introducing the di-mensionless time t/τ → t , the equation of motion forthe dimensionless vertex positions r j /V / c → r j reads˙ r j = −∇ j w + f ( A ) j , where f ( A ) j is the dimensionless activeforce on the j -th vertex measured in units of Γ l V / c . Inturn, the dynamical equation for dimensionless junctionaltensions γ i / (cid:16) Γ l V / c (cid:17) → γ i reads ˙ γ i ( t ) = − γ i ( t ) /τ m + ξ i ( t ), where the long-time variance of dimensionless linetension fluctuations ξ i is given by σ / (cid:16) Γ l V / c τ − (cid:17) → σ and the dimensionless inverse myosin turnover rate is τ /τ m → /τ m . Here we study the effect of the magni-tude of junctional tension fluctuations σ but we keep themyosin turnover rate fixed at 1 /τ m = 1.To determine the critical magnitude of junctional ten-sion fluctuations σ ∗ where the solid-fluid transition takesplace as previously discussed and shown to hold in a 2Dmodel [17], we simulate fluctuations in a flat non-stressedepithelial patch in absence of differential tension, i.e., at α = β . During these simulations, we affix the basal sidesof the cells to a flat plane so as to prevent buckling outof the plane. Starting from a honeycomb arrangement, weallow the simulation to run at a given σ from t = − to t = 0. Following Ref. [17], we then continue to run thesimulation to t = 10 and we measure the mean squaredisplacement of cell centers given byMSD( t ) = (1 /N c ) N c X i =1 (cid:12)(cid:12)(cid:12) r ( i ) b ( t ) − r ( i ) b (0) (cid:12)(cid:12)(cid:12) , (9) where r ( i ) b ( t ) is the in-plane position of the center of basalside of the i -th cell. Using MSD as a measure of cell move-ment, we can then determine the effective diffusion co-efficient D eff = MSD( t ) / (4 t ) | t →∞ . Similar to what wasobserved in the 2D version of the vertex model [17], thediffusion coefficient only takes finite values above a certainthreshold of tension fluctuations. At the tissue tension of α + β = 2 which pertains to most results presented in thispaper, the threshold is σ = σ ∗ ≈ . α + β between1 and 3, the critical σ ∗ increases with α + β as σ ∗ ( α + β ) ≈ .
123 ( α + β ) . . (10) Within the computational scheme outlined above, we ex-plore the parameter space of our model by performing aseries of simulation runs; rather than embarking on a sys-tematic scan of this space, we focus on the most interest-ing effects at work. In each run, we start with a flat mono-layer consisting of N c identical cells that form a hexagonalprismatic honeycomb, and we simulate its evolution uponcompression until the final strain is reached at a time T .Here we examine a honeycomb that consists of 64 rows of64 cells arranged in a rectangular patch with an aspectratio of 1 : √ / N c = 4096. The compression time measured in units of τ is T = 1000 unless stated otherwise, and α + β = 2,again unless stated otherwise. In the isotropic compressionmode ε x = ε y = 15 %, whereas in the uniaxial compres-sion mode ε x = 27 .
75 % and ε y = 0 so that the projectedareas of the isotropically and the uniaxially strained tis-sues given by L x ( T ) L y ( T ) are identical. To facilitate theevolution of the tissue from the very regular initial hon-eycomb state, we break its apico-basal symmetry by dis-placing each vertex at t = 0 in a random direction suchthat the displacements are normally distributed aroundzero with a standard deviation of 0 . V / c .An essential element of our scheme is the implementa-tion of the T1 transition, which is executed in two steps asdescribed in the following. If the length of a junction (de-fined as the average of the corresponding apical and basaledge length) falls under a threshold value l = 0 .
01 andhas decreased in length since the previous time-step, thejunction is first merged into a fourfold vertex, giving riseto a local 4-cell rosette arrangement. The rosette is thenmaintained for a short time interval of 0.02; this is the sec-ond step. After this time, it is resolved into two three-wayvertices, displacing each from the pre-resolution positionby 0.0005. After the T1 transition, a new line tension isassigned to the junction from a normal distribution withmean 0 and standard deviation σ .Lastly, we note that our scheme does not include thesteric repulsion between the cells; it is therefore possiblefor non-neigboring cells to overlap. In practice, this issuemay arise only on the basal (inner) side of the villus mor-phologies, and this may only happen after the villi are an Rozman et al.: Morphologies of compressed active epithelial monolayers 5 already established. As such, self-overlap may slightly af-fect the detailed shape of the villi but not their existenceor the qualitative conclusions of the paper. Our scheme produces four types of morphologies of com-pressed epithelial monolayers: The longitudinal and thelabytrinthine folds depend on the type of compression(uniaxial and isotropic, respectively) whereas previlli andvilli are independent of the details of compression but re-quire that the tissue is fluidized by junctional fluctuations.
The simplest nontrivial epithelial shape occurs in a uniax-ially compressed tissue with a level of junctional tensionfluctuations below the solid-fluid transition. Much like arod compressed in the lengthwise direction, the solid-likeepithelium undergoes an Euler instability and a bucklingtransition, assuming a corrugated shape that features aseries of parallel folds running perpendicular to strain di-rection (Fig. 1d, Fig. 2a, and Movie S1). In our scheme,buckling does not develop immediately after the strain isapplied but only after a certain finite time. The develop-ment of the deformation can be appreciated by monitor-ing the deviation of the apical vertices from their averageheight δz a , which is close to 0 in a flat tissue before buck-ling and finite in a modulated tissue. Figure 2b shows thatat the magnitude of the initial random displacement of thevertices used here, the deformation of a σ = 0 solid-liketissue rapidly forms around time t ≈
500 as witnessed bythe step-like profile of δz a ( t ); this buckling time is slightlyshorter at larger differential tensions β − α . After that,the deformation continues to grow which shows that theresponse of the tissue to strain is dynamic. As may beanticipated, the dynamics involves fold coarsening; herewe collect and report tissue shape right after compressiontime T = 1000. Notably, the time at which buckling oc-curs does significantly depend on the scale of the initialperturbation of the vertices at time t = 0 (Fig. S1).The σ = 0 fold morphologies obtained at a compres-sion time of T = 1000 shown in the top row of Fig. 2a donot seem to be very sensitive to the apico-basal differentialtension β − α . However, the number of folds decreases withthe compression time T as shown in Fig. 2c. The effect ofthe strain rate is qualitatively similar to folding in a toymodel of a fixed-end growing tissue consisting of circularcells where the number of folds decreases with the dura-tion of the cell cycle; that is, in a slowly dividing tissuethe number of folds is smaller than in a rapidly divid-ing one [8]. The data shown in Fig. 2c also illustrates thesubdominant dependence of the number of folds on thedifferential tension. Despite the large error bars, our datasuggest that the number of folds slightly increases with β − α . This effect is not surprising since the differentialtension promotes the constriction of the basal cell sides,inducing a preferred truncated-pyramid cell shape which - - β - α = β - α = β - α = σ = σ = . σ = . t s t anda r dde v i a t i on δ z a β - α T nu m be r o ff o l d s β - α a ) b ) c ) Fig. 2.
Height profiles of the apical surface of epithelial mor-phologies formed by uniaxial compression along the x axis by27.75 % at σ = 0 , . , and 0 . β − α = 0 , . , and 0 . β − α = 0 and 0.2) and the villi ( β − α = 0 .
8) morphologies.The profiles show the positions of the centers of apical sides,with 0 corresponding to the average height. In panel b we plotthe temporal profile of the standard deviation of the height ofapical vertices in three σ = 0 solid-like tissues, averaged overten simulation runs, which indicates that buckling appears ata time of t ≈ σ = 0solid-like tissue vs. compression time T at three different β − α ,averaged over ten simulation runs; error bars indicate standarddeviation. Note that the precise count of folds is somewhat sub-jective, as some individual folds are branched or do not stretchacross all of the patch. is more easily accommodated in a folded than in a flat tis-sue [4,22]. Lastly, we note that at a finite yet subcriticalmagnitude of junctional tension fluctuations the numberof folds in a uniaxially compressed epithelial monolayer is Jan Rozman et al.: Morphologies of compressed active epithelial monolayers more or less the same as at σ = 0 as illustrated by the σ = 0 .
15 sequence shown in the middle row of Fig. 2a.The longitudinal folds are also interesting because theyallow us to explore the thickness-curvature coupling [4,12]as transparently as possible. The effect of the coupling canbe quantified by thickness modulation ρ = h − h h , (11)where h is the height of a cell defined as the distancebetween the centroids of its apical and basal sides whereas h = (cid:0) √ / (cid:1) − / ( α + β ) / is the equilibrium cell heightin a flat regular honeycomb tissue (Fig. 3b) [14].The mechanism behind thickness modulation goes asfollows. In a tissue with a given apico-basal differentialtension β − α >
0, the preferred minimal-energy shapeof cells is a truncated pyramid with the basal side smallerthan the apical side because of a larger basal tension; sucha shape has a certain curvature. A fold includes cells ofpreferred shape but it must also consist of those with theopposite curvature because the integrated curvature of itscontour must vanish. In the latter, the basal side is largerthan the apical side, which is energetically unfavorable.The local energy increase can be minimized by moving theapical and the basal side closer to each other (accompa-nied by the expansion of the in-plane diameter due to thefixed-volume constraint). In the elastic theory, this phe-nomenon shows as the thickness-curvature coupling termwith a modulus proportional to β − α [10,12].In Fig. 3a we compare three examples of fold mor-phologies at β − α = 0 , . , and 0.8. One can spot the con-strictions on the basal (bottom) sides of crest cells and themore broad and open shape of the grooves present in the β − α = 0 . β − α = 0folds. In the β − α = 0 . ≈
10 % of theequilibrium height in a flat honeycomb at β − α = 0 to ≈
25 % at β − α = 0 .
8, which emphasizes that the effectis proportional to β − α . In epithelia that are completely fluidized by junctionaltension fluctuations, buckling caused by strain takes on adifferent from than in solid-like tissues even at small apico-basal differential tension (Fig. 1f, Fig. 2a, and Movie S2).Instead of the relatively evenly spaced and regular folds,the uniaxially strained fluidized epithelium forms a struc-ture consisting primarily of bumps and pits, i.e., protru-sions pointing up and down. In addition to these spatiallyisolated features, the nearby bumps and pits fuse into un-even short crests and grooves, respectively, essentially dou-blets or triplets of individual protrusions as illustrated bythe σ = 0 . , β − α = 0 morphology in the bottom row of β - α = β - α = β - α = - - ρβ - α = β - α = β - α = σ = a ) b ) Fig. 3.
Side views of the σ = 0 longitudinal folds in a uniaxiallystrained model tissue at apico-basal differential tensions β − α = 0 , . , and 0.8 (a). Thickness modulation in these threefold morphologies (b). Fig. 2a. Interestingly, there seems to be no preferred ave-rage in-plane orientation of these features, which suggeststhat this morphology does not depends on the directionof the applied strain.The bumps and pits are similar to villi and crypts,respectively, in that they are local rather than extendedspatial features of the epithelium, but they are much lessdeveloped. At the same time, as the apico-basal differen-tial tension is increased, the bumps become more sharplydefined to the point where they can be viewed as villi,and the space between them is flattened out. As such thebumps-and-pits morphology may be referred to as the pre-villi (or, equivalently, precrypts if β − α is negative). The third epithelial morphology are the villi observed incompletely fluidized tissues at a large enough apico-basaldifferential tension (Fig. 1g and Movie S3). As shown bythe bottom row of Fig. 2a, at β − α = 0 . β − α = 0 . an Rozman et al.: Morphologies of compressed active epithelial monolayers 7 coordinated cells at the tip. This characteristic distribu-tion has been first reported in the branches of organoid-like epithelial shells [12] but is apparently an intrinsic fea-ture of protrusion-forming vertex models.The results presented here do not depend very muchon the tissue tension α + β in a significant range around α + β = 2. The three discussed morphologies— longitu-dinal folds, previlli, and villi—also arise at significantlylower ( α + β = 1) as well as higher ( α + β = 3) tissuetensions under an appropriate rescaling of the differentialtension β − α and tension fluctuation σ (Fig. S2). Thesame is also true of the labyrinthine morphology intro-duced in Sec. 3.4. Furthermore, our choice to compressthe tissue along the x axis does not significantly affectthe observed morphologies; the results are essentially thesame for model tissues under compression along the y axis(Fig. S3).Our model produces an inconsequential artifact in thevillus morphologies where the basal sides of many cells atthe tips are very small or even vanish. In these pyramidalcells, one or more basal vertices can move towards theapical side so as to increase the length of the basal edgesat unchanged small or zero area of the basal side. Thismay be energetically favorable because of the fluctuating,occasionally negative basal edge tension (Fig. S4a). As aresult, several spikes may appear on the basal side of thevilli tips (Fig. S4b).The basal spikes may be suppressed by including anenergy penalty on the very small basal areas. This auxi-liary term reads P i k b h A ( i ) b − A i H (cid:0) A − A ( i ) b (cid:1) , wherethe sum runs over all cells, k b is the modulus, A ( i ) b is thebasal area of the i -th cell, A is the threshold area underwhich the penalty is activated, and H is the Heavisidefunction [23]; here k b = 200 and A is set to 30 % of theequilibrium basal area in a flat regular honeycomb tissue.With this modification, the villi are almost completely freeof the spikes on the basal side (Fig. S4c). After analyzing uniaxially compressed tissues, we also con-sider two other strain modes. Under isotropic compressionwith ε x = ε y = 15 %, a solid-like σ < σ ∗ monolayer formsa disordered labyrinthine arrangement of folds (Fig. 4a,Fig. S5, and Movie S4). Here the folds do not span acrossthe whole patch but are of finite length, usually bentrather than straight, and lack a preferred orientation. Thedetailed shape of the folds depends only slightly on theapico-basal differential tension. An even more interestingresult of the isotropic-strain case is the σ = 0 . , β − α = 0 . σ = 0. For clarity, we use the ε x = 27 .
75 %, T = 1000 fold state with β − α = 0 . ε y = 22 % in time T = 2200. We find that the se-quential compression as described here indeed leads to aherringbone pattern of zigzagging folds shown in Figs. 4band c. The different morphologies can be arranged in a state di-agram shown in Fig. 5a, where we combine the resultspertaining to the different types of strains so as to bet-ter emphasize the individual roles of strain, apico-basaldifferential tension, and activity. This allows us to groupthe two low-activity fold morphologies—longitudinal andlabyrinthine—in a single domain extending to σ . σ ∗ ≈ .
17 (the threshold σ is slightly higher at lower β − α , butthis effect is not very pronounced). On the other hand,the σ > σ ∗ domain is occupied by previlli and by villiat small and large differential tensions, respectively. Inthe explored region of the phase space, villi appear at β − α ≈ σ is increased, the boundary shifts rapidly to significantly β - α = σ = σ = . a ) b ) c ) Fig. 4.
Height profiles of the apical surface of an isotropicallycompressed epithelial monolayer with a strain of 15 % and amagnitude of junctional tension fluctuations σ = 0 and 0.2 (a),displaying the labyrinthine and the villus morphologies, respec-tively. The color legend is the same as in Fig. 2a. Panel b showsthe σ = 0 , β − α = 0 . y axis by 22 % at a com-pression time of T = 2200; panel c is an oblique 3D view ofthe herringbone pattern. Jan Rozman et al.: Morphologies of compressed active epithelial monolayers β - α t en s i on f l u c t ua t i on s σ longitudinal / labyrinthinefolds p r e v illi villi -
10 0 10apical height ε x = ε y = α + β = β - α = σ = bar2 ε x = ε y = α + β = β - α = σ = a ) b ) c ) d ) e ) Fig. 5.
State diagram of a compressed epithelial monolayerwith α + β = 2 at isotropic strain with ε x = ε y = 15 %combined with that at equivalent uniaxial strain with ε x =27 .
75 % and ε y = 0 (a). Panels b and c show an array ofprevillous bumps in a solid-like tissue with a large differentialtension at ε x = ε y = 0 with σ = 0 at time t = 10000, whereasin panels d and e we plot an example of the villus morphologyinduced by tension fluctuations in an unstrained ε x = ε y = 0flat tissue at t = 3000. smaller differential tensions down to β − α ≈ .
4; fromthis point on, the boundary is relatively vertical, reaching β − α ≈ . σ = 0 .
3, the largest magnitude of junctionaltension fluctuations considered. We note that the distinc-tion between the previlli and the fully developed villi isquite arbitrary, hence the blurred boundary in the statediagram. More specifically, as the differential tension in-creases, the symmetry between the previllous bumps andpits breaks so that the former exceed the latter. After-wards, the bumps become more clearly defined, produc-ing the villi (Fig. S6). Similarly, the boundary betweenthe solid-like folds and the previlli can also be somewhat unclear, because, e.g., some of the longitudinal orientationis retained at σ slightly above 0.165.The state diagram in Fig. 5a emphasizes that villus for-mation relies on a strong enough junctional tension fluc-tuations as well as on a large enough apico-basal differen-tial tension. The diagram corresponds to isotropic in-planestrain with ε x = ε y = 15 %, and to the equivalent uni-axial strain with ε x = 27 .
75 % and ε y = 0, but we notethat it remains qualitatively the same at any moderatecompression.The natural question to ask is to what extent are theseresults affected by compression. To examine this effect, wescanned the parameter space so as to find that an initiallyflat σ = 0 model tissue placed in a perfectly fitting boxwith ε x = ε y = 0 did undergo a deformation after a largeenough differential tension was turned on, forming a pre-villous array; an example with α + β = 3 and β − α = 2 . ε x = ε y = 0 demonstrates that this morphology maybe induced by differential tension alone; note that whenstudying the effects of compression, we only focused onmoderate differential tensions up to β − α = 1, which areinsufficient to produce previllous bumps at α + β = 2 inabsence of strain. Likewise, compression is not essential forthe existence of villi. This is illustrated by Figs. 5d and eshowing villi that developed at ε x = ε y = 0 from a fluid-like flat tissue with σ = 0 . α + β = 2 , and β − α = 0 . an Rozman et al.: Morphologies of compressed active epithelial monolayers 9 shape of the mesenchyme resulting from the compressivestrain due to the muscle around it.Our results suggest that a very similar progression ofthe morphologies can conceivably be obtained in a sub-strate-free epithelium. We show that the first two stepsof the process—the formation of longitudinal folds andzigzags—can be obtained by a suitable sequential biaxialstrain in the solid-like regime. The third step, the forma-tion of villi from zigzags, was not explicitly tested in outmodel, but based on the above findings it seems possi-ble that a suitable combination of tissue fluidization andan increase in differential tension would produce such anoutcome. In a similar manner, our results also relate to thebuckled morphologies generated by differential isotropic [3]and anisotropic [5] growth of an epithelium on a support-ing tissue as such growth produces strain.More importantly, our model also provides an interpre-tation for the development of the mouse-gut villi whichform directly from the smooth epithelium. In Ref. [1],this particular morphology was attributed to a softer en-doderm compared to that in the chick embryo, but ourresults suggest that the mouse-gut-like villi may also re-sult from a strong enough activity, a suitable compressivestrain, and a large enough differential tension. We alsonote that the mouse-gut villi do not form an ordered ar-ray but are randomly distributed across the epithelium.This agrees with our villus morphologies as illustrated bybottom-right profile in Fig. 2a but not with the regularpatterns obtained in models that build exclusively on elas-tic instability [1,3,5].The predictions of our model depart from those of theexisting theories in several respects; here we mention twokey differences. The first one is the detailed shape of thevilli as the most elaborated nontrivial epithelial morpho-logy. Our model villi are considerably taller and moretubular than those seen in the buckled elastic sheets [1,3,5], which are shorter and more similar to our previlli.As far as their height-to-width ratio is concerned, our villiare quite close to many real ones. Secondly, our modelepithelium is devoid of an explicit support such as anelastic foundation. As such, it provides a basis for theunderstanding of villus-like formations characteristic oforganoids as well as an insight into certain embryonicstructures such as the archenteron in the sea urchin. Geo-metrically, these formations are closely reminiscent of theintestinal villi, yet they develop in absence of surroundingtissues, which implies that their morphogenesis can onlybe related to the intrinsic mechanics of the epitheliumand to global constraints such as those fixing the lumenor blastocoel volume [12].Of course, a generalized version of our model where theepithelium would rest on an elastic foundation could wellbe very interesting. We expect that such a substrate maydefine the wavelength of the folds [2] as well as the lengthand the detailed shape of the villi, but it may also giverise to novel epithelial forms. In addition, the substratewould distinguish between the villi and the crypts. Ourpresent model applies equivalently to both of these mor-phologies; the only difference between the apical and the basal surface is due to the differential tension, and thusa villus at a given β − α > β − α <
0. However, in anepithelium supported by an elastic substrate this equiva-lence no longer exists, with the villus growing out of thesubstrate and the crypt protruding into it.Lastly, we emphasize that in our model, the existenceof villi relies on the topological activity, the villi them-selves having a nontrivial topological structure (Fig. 1h).This element of mechanics is absent in the continuum the-ories of epithelial form and it is plausible that the interplayof substrate elasticity and cell-level activity could resultin novel morphologies if the former would interfere withthe topology of the tissue.
The most important aspect of this study may well be theexplicit comparison of the results of selected external andinternal factors responsible for epithelial morphogenesiswithin a single-cell-population tissue. The strain arisingfrom the competition between the tissue and any adjacentstructure, be that it results from the growth of the epithe-lium attached to a substrate or from stress generated bythe substrate, is a very plausible external source of buck-ling, yet we find that the buckled shape induced by it in-cludes folds and bumps, but not true, well-developed villi.On the other hand, a high enough level of activity withinthe tissue and its ensuing fluidization is enough to promotevillus and crypt formation provided that the apico-basaldifferential tension is large enough. Naturally, in a differ-entiated tissue consisting of biologically and mechanicallydistinct cell populations there may exist additional inter-nal processes that contribute to the development of bothfold and villus morphologies.The potential internal origin of villi and crypts alsopoints to an intriguing hypothesis that certain other mor-phogenetic process in epithelial tissues may too be in-trinsic rather than induced or driven by the environment.Branching, for example, relies on the formation of a Y junc-tion at the tip of the epithelial tube, and the saddle-shapedbases of the two growing buds are quite similar to thoseseen in the villi. From this perspective, budding is just avillus growing from the tip of an existing villus. Shouldthis be possible, our vertex-model framework would pro-vide the micromechanical interpretation of the elementaryprocess of branching morphogenesis [24].We thank P. Mrak and M. Rauzi for helpful discus-sions, and we acknowledge the financial support from theSlovenian Research Agency (research core funding No. P1-0055 and projects No. J2-9223 and Z1-1851).
Author contribution statement
All authors were involved in the preparation of the manu-script. All authors have read and approved the final manu-script.
Data availability statement
The datasets generated during and/or analyzed during thecurrent study are available from the corresponding authoron reasonable request.
Conflict of interest
The authors declare that they have no conflict of interest.
Publisher’s Note
The EPJ Publishers remain neutral withregard to jurisdictional claims in published maps and in-stitutional affiliations.
References
1. A.E. Shyer, T. Tallinen, N.L. Nerurkar, Z. Wei, E.S. Gil,D.L. Kaplan, C.J. Tabin, L. Mahadevan, Science , 212(2020).2. E. Hannezo, J. Prost, J.-F. Joanny, Phys. Rev. Lett. ,078104 (2011).3. P. Ciarletta, V. Balbi, E. Kuhl, Phys. Rev. Lett. ,248101 (2014).4. N. ˇStorgel, M. Krajnc, P. Mrak, J. ˇStrus, P. Ziherl, Bio-phys. J. , 269 (2016).5. M. Ben Amar, F. Jia, Proc. Natl. Acad. Sci. USA ,10525 (2013).6. J. Galle, M. Loeffler, D. Drasdo, Biophys. J. , 62 (2005).7. P. Buske, J. Galle, N. Barker, G. Aust, H. Clevers, M. Lo-effler, PLoS Comput. Biol. , e1001045 (2011).8. D. Drasdo, Phys. Rev. Lett. , 4244 (2000).9. D. Drasdo, S. H¨ohme, Phys. Biol. , 133 (2005).10. M. Krajnc, P. Ziherl, Phys. Rev. E , 052713 (2015).11. P.A. Haas, R.E. Goldstein, Phys. Rev. E , 022411(2019).12. J. Rozman, M. Krajnc, P. Ziherl, Nat. Commun. , 3805(2020).13. J. Ranft, M. Basan, J. Elgeti, J.-F. Joanny, J. Prost,F. J¨ulicher, Proc. Natl. Acad. Sci. USA , 20863 (2010).14. M. Krajnc, S. Dasgupta, P. Ziherl, J. Prost, Phys. Rev. E. , 022409 (2018).15. D. Bi, J.H. Lopez, J.M. Schwarz, M.L. Manning, Soft Mat-ter , 1885 (2014).16. S. Curran, C. Sandkvist, J. Bathmann, M. de Gennes,A. Kabla, G. Salbreux, B. Baum, Dev. Cell , 480 (2017).17. M. Krajnc, Soft Matter , 3209 (2020).18. S. Okuda, Y. Inoue, M. Eiraku, T. Adachi, Y. Sasai,Biomech. Model. Mechanobiol. , 413 (2015).19. M. Misra, B. Audoly, I. G. Kevrekidis, S.Y. Shvartsman,Biophys. J. , 1670 (2016).20. T. Lecuit, P.-F. Lenne, Nat. Rev. Mol. Cell Biol. , 633(2007).21. A. Mongera, P. Rowghanian, H.J. Gustafson, E. Shelton,D.A. Kealhofer, E.K. Carn, F. Serwane, A.A. Lucio, J. Gi-ammona, O. Campas, Nature , 401 (2018).22. M. Krajnc, N. ˇStorgel, A. Hoˇcevar-Brezavˇsˇcek, P. Ziherl,Soft Matter , 8368 (2013). 23. L. Sui, S. Alt, M. Weigert, N. Dye, S. Eaton, F. Jug,E.W. Myers, F. J¨ulicher, G. Salbreux, C. Dahmann, Nat.Commun. , 4620 (2018).24. E. Hannezo, C.L.G.J. Scheele, M. Moad, N. Drogo,R. Heer, R.V. Sampogna, J. van Rheenen, B.D. Simons,Cell , 242 (2017). upplementary information Morphologies of compressed active epithelial monolayers
J. Rozman , M. Krajnc , and P. Ziherl Jožef Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia a r X i v : . [ c ond - m a t . s o f t ] F e b upplementary figures t s t anda r dde v i a t i on δ z a δ r = - β - α t s t anda r dde v i a t i on δ z a δ r = - β - α t s t anda r dde v i a t i on δ z a δ r = - β - α t s t anda r dde v i a t i on δ z a δ r = - β - α t s t anda r dde v i a t i on δ z a δ r = - β - α t s t anda r dde v i a t i on δ z a δ r = - β - α t s t anda r dde v i a t i on δ z a δ r = - β - α t s t anda r dde v i a t i on δ z a δ r = β - α Figure S1:
Standard deviation of the height of apical vertices δz a as a function of time at eight different values ofthe standard deviation of the initial perturbation of the vertices δr .1 - β - α = β - α = β - α = σ = σ = . ( / ) / σ = . ( / ) / - - β - α = β - α = β - α = σ = σ = . ( / ) / σ = . ( / ) / α + β = α + β = a ) b ) Figure S2:
Height profiles of the apical surface of model epithelia formed by uniaxial compression in the horizontaldirection by 27.75 % for tissue tensions α + β = 1 (a) and 3 (b). To facilitate comparison, we show the morphologiesat differential tensions β − α such that the ratio ( β − α ) / ( α + β ) in the left, middle, and right column in panelsa and b as well as in Fig. 2a is 0, 0.1, and 0.4, respectively. We find that in the studied range of tensions, areasonably good way to select equivalent tension fluctuations for comparison at different α + β is to rescale thevalue of σ by the characteristic in-plane dimension and the effective mean line tension. The characteristic in-planedimension is given by the linear size of the equilibrium apical or basal side of a hexagonal-prismatic cell equal to a = (cid:0) √ / (cid:1) / ( α + β ) − / . To determine the effective mean line tension, we follow Ref. [1]. In a flat tissue witha constant cell height h and consequently a constant cell apical and basal area a , we recast the non-active partof the energy [Eq. 5 of the main text] as w = P N c i =1 (cid:2) αa + βa + h p ( i ) / (cid:3) = const . + P N e j =1 h l ( j ) , where p ( i ) is theperimeter of the i -th cell and l ( j ) is the length of the j -th junction. From this expression we can see that the cellheight h plays the role of the mean effective line tension; in turn, h can be estimated by the equilibrium value ina hexagonal-prismatic cell h = (cid:0) √ / (cid:1) − / ( α + β ) / . Therefore, we chose σ values for comparison by the relation σ ( α + β ) / ( α + β ) / ≈ σ ( α + β ) / ( α + β ) / . The morphologies obtained at the appropriately rescaled differentialtensions and tension fluctuations for α + β = 1, 2, and 3 are very similar. Note that while the profiles shown hereare plotted in same-size rectangles so as to facilitate comparison, the actual area of the more columnar tissue with α + β = 3 is smaller than in the tissue in Fig. 2a, whereas the area of the more cuboidal tissue with α + β = 1 islarger than in that in Fig. 2a. 2 - β - α = β - α = β - α = σ = σ = . σ = . Figure S3:
Height profiles of the apical surface of epithelial morphologies at α + β = 2 formed by uniaxial compressionalong the y axis by 27.75 %. The morphologies obtained are equivalent to those that appear under compression alongthe x axis. a ) b ) c ) Figure S4:
Example of a pyramidal cell at the tip of the villus, with vanishingly small basal side and several basalvertices shifted towards the apical side (a). The apical side is shown in white, the lateral sides are cyan, and the basalvertices are indicated by blue points. Panel b shows an example of a villus with spikes at the basal side of cells at thetip, whereas that in panel c is free of such spikes as it was computed using a generalized version of our vertex modelwhich includes an auxiliary term that prevents the appearance of the spikes by imposing an energy penalty on smallbasal areas. For simplicity, the villi presented in panels b and c were computed in a patch of 32 ×
32 = 1024 cellswith α + β = 2, β − α = 0 .
8, and σ = 0 . ε x = 27 .
75 % and ε y = 0.3 - β - α = β - α = β - α = σ = σ = . σ = . Figure S5:
Height profiles of the apical surface of epithelial morphologies at α + β = 2 formed by isotropiccompression by ε x = ε y = 15 %. β - α = β - α = β - α = β - α = β - α = a ) b ) c ) d ) e ) Figure S6:
Apical surface view of computed epithelial morphologies following uniaxial compression with α + β = 2, σ = 0 . ε x = 27 .
75 %, and ε y = 0, at differential tension β − α = 0 . β − α = 0 . β − α = 0 . β − α = 0 . β − α = 0 . upplementary movies Movie S1.
Temporal evolution of longitudinal folds with α + β = 2, β − α = 0, and σ = 0 under uniaxial compressionwith ε x = 27 .
75 % and ε y = 0. Movie S2.
Temporal evolution of previlli with α + β = 2, β − α = 0, and σ = 0 . ε x = 27 .
75 % and ε y = 0. Movie S3.
Temporal evolution of villi with α + β = 2, β − α = 0 .
8, and σ = 0 . ε x = 27 .
75 % and ε y = 0. Movie S4.
Temporal evolution of labyrinthine folds with α + β = 2, β − α = 0, and σ = 0 under isotropic compres-sion with ε x = ε y = 15 %. Supplementary references [1] M. Krajnc, S. Dasgupta, P. Ziherl, J. Prost, Phys. Rev. E.98