A level set-based approach for modeling cellular rearrangements in tissue morphogenesis

Mathematical models and numerical simulations can provide an essential insight into the mechanisms through which local cell-cell interactions aﬀect tissue-level cell morphology. Among such morphological phenomena, cellular patterns observed in developing sensory epithelia have gained keen attention of researchers in recent years, because they are thought to be of utmost importance for accurate sensory functions. However, most of current computational approaches to cellular rearrangements lack solid mathematical background and involve experimentally unreachable parameters, whereby only weak and ambiguous conclusions can be made based on simulation results. Here we present a simple mathematical model for tissue morphogenesis together with a level set-based numerical scheme for its solution as a tool to rigorously investigate evolving cellular patterns. This combined framework of a model and a numerical method features minimum possible number of physical parameters and guarantees reliability of simulation results, including correct handling of topology changes, such as cell intercalations. In this framework, we adopt the viewpoint of free energy minimization principle, and take cellular rearrangement as a gradient ﬂow of a weighted surface energy associated with cell membrane, where the weights are related to physical parameters of the cells, for example, cell-cell adhesion and cell contractility. We present the applicability of this model to a wide range of tissue morphological phenomena, such as cell sorting, engulfment or internalization. In particular, we stress that this method is the ﬁrst one to be successful in computationally reproducing the experimentally observed development of cellular mosaic patterns in sensory epithelia. Thanks to its simplicity and reliability, the model is able to capture the essence of biological phenomena, and may give a strong helping hand in deciphering unsolved questions of morphology.


Introduction
of elements to solve mechanical equations characterizing tissue dynamics. In the two-dimensional setting, each cell-cell junction is represented as a piecewise linear curve and interfacial tensions are modeled using 64 constant force on point masses along each cell-cell junction. In this sense, this formulation acts like the 65 vertex dynamics model. 66 Although each specific morphogenetic phenomenon involves a large number of biological and physical 67 factors, it is a well accepted understanding in the modeling community that it is not reasonable to con-  without relying on artificial procedures that inevitably involve unphysical parameters. Hence, the presented 86 model is well suited for testing various hypotheses related to developmental morphology. We support this 87 fact in Section 6 using the example of formation of cellular patterns in sensory epithelia. As mentioned 88 above, such pattern formation entails complex curved shapes of cell junctions, largely different cell sizes and 89 frequent topology changes, which make its numerical simulation challenging. Let us now present the mathematical model and the numerical algorithm for its solution. 92 We start with the formulation of the mathematical model. We represent an aggregate of cells as a 93 bounded domain Ω ⊂ R d (d = 2 or 3) partitioned into N closed sets C 1 , . . . , C N , representing cells (see 94 Figure 2 for the basic notation). It naturally follows that cell-cell junction γ ij := C i ∩ C j = ∂C i ∩ ∂C j is the 95 boundary between cells C i and C j . We consider cellular rearrangement as the L 2 -gradient flow (see (Laux 96 constrained by each cell's prescribed volume V 0 ( = 1, 2, . . . , N ). Here, the weights σ ij = σ ji > 0 for 98 i = j may be related to cell-cell adhesion and/or cell contractility (see Section 3 for more). When i = j, we 99 formally set σ ij = σ ii = 0.

100
In materials science, this problem (considered without volume constraint to begin with) is widely known 101 as the Mullins (1999) model for normal grain growth where C i denotes a crystallite/grain in polycrystalline dimensions. To alleviate this drawback, Merriman et al. (1994) introduced the MBO thresholding scheme for diffusion-generated curvature-dependent motion of multiple junctions, which is based on the level set for-120 mulation of Osher and Sethian (1988) for propagating fronts with curvature-dependent speed. This scheme 121 tracks interfaces implicitly by following level sets -allowing it to naturally handle topological changes. In the configuration u k+1 at the next time step t k + δt by minimizing a linear functional: and G δt is the Gaussian kernel, over the convex constraint set  To numerically implement (2), a straightforward scheme is to incorporate volume constraints with La-135 grange multipliers, however, a more efficient approach employs auction algorithms (Jacobs et al., 2018).

136
The idea is to assign cell membership to each point of the discretized domain Ω m = {x 1 , . . . , x m } ⊂ Ω by 137 simulating an auction, so that each cell C j contains v j points, where the integers v j are chosen in such a way 138 that the mass ratios v j / k v k are as close to V 0 j /|Ω| as possible.

139
The Esedoḡlu-Otto scheme is simple and efficient but there are some issues that need to be tackled, with the lowest bid by x.

171
(e) Update the price p i * = min

172
We now briefly comment on the parameters related to the numerical implementation. The discretization This condition is often satisfied in materials science but there is no guarantee that it will hold in biological i.e., the minimization problem (2). Stability can be obtained by a simple modification of the original proof 201 (see (Xu et al., 2017) for the basic idea) and convergence has been established in (Laux and Swartz, 2017).

202
Finally, the convergence of the auction algorithm is well-known (see the references in (Jacobs et al., 2018)).

203
It follows that the localized auction scheme is justified locally, but analysis of its global behavior is still 204 missing.

205
In this section, we have presented an algorithm that can express arbitrary curved shapes of cellular 206 junctions thanks to its being rooted in the level set approach. In the following two sections, we address the 207 remaining two benefits of our algorithm announced at the end of the introductory section: small number of 208 parameters (Section 3), and natural treatment of topology changes (Section 4).

Parameters of the Algorithm
Besides the parameters m, δt, and ε, used in the numerical implementation of the model, the only 211 parameters of the model itself are the weights σ ij , showing a major advantage over other similar models.

212
For example, the vertex dynamics model includes, besides these weights, several further parameters, such as 213 the elastic stiffness coefficients for the volume and perimeter, minimal distance of two vertices to initiate T1 214 transformation, a parameter for the resulting distance of the vertices, etc. Consequently, our model allows 215 us to focus on the role of the target force represented by the parameters σ ij , without having to analyze the 216 impact of other effects. This is essential as such analysis often presents an infeasible task. energy in a suitable manner that has to be determined as a part of the particular model), or a combination 239 of both. In Section 6, we give a specific example of the design of the weights σ ij that express the adhesion 240 energy and are quantified through the experimentally measurable quantity of β-catenin intensity. 241 We would like to bring attention to the fact that the weights σ ij in the algorithm can be time-dependent, 242 which is expressed by the index k in equation (3). This is essential, as almost all morphogenetic phenomena 243 are driven by temporally changing forces. Moreover, this time-dependence turns our model from a mere 244 energy minimizing gradient descent system into an out-of-equilibrium one, as expected for a model of a 245 living system.

246
The minimal set of parameters, i.e., the weights σ ij , can be augmented by new parameters according to 247 necessity to express various additional aspects of the target biological phenomenon. One example of such 248 additional parameters are the mobilities µ ij mentioned in the beginning of Section 2. In fact, the algorithm 249 as presented in Section 2 advances each junction γ ij with the mobility µ ij = 1/σ ij . If one wishes to prescribe 250 mobilities in a different way, it is possible to modify the algorithm by introducing so-called retardation 251 terms, as explained in Section 5.1 of (Esedoḡlu and Otto, 2015). This method was improved in (Salvador 252 and Esedoḡlu, 2019) by replacing the computation of retardation functions by a more efficient convolution 253 step.

254
Another optional set of parameters that the discrete cell volumes v k i (i.e., the number of grid points in cell C i at a given time t k ) satisfy It is possible to modify the auction algorithm to extend it to this type of constraint. The implementation   Consider a wetting case for a 4-cell aggregate of three types: 2 blue, 1 red, and 1 gray cell (see Figure   269 3), with σ BB = 1.5, σ BG = 0.5, and σ RR = σ GG = σ BR = σ RG = 1.0. Here, B, R, G denote the cells of 270 type blue, red and gray, respectively, so that, for example, σ BR means the interfacial energy weight for the 271 junction between blue and red cell. Note that this violates the triangle inequality, since σ B1G + σ B2G = 272 1 < 1.5 = σ B1B2 . We evolve the initial configuration using three different algorithms: the original Esedoḡlu Observe that for the Esedoḡlu-Otto scheme, a new gray cell grows at the BB-junctions -wetting occurs.

279
When implemented with usual auction dynamics, we see that the gray cell splits and some of its parts 280 appear in the BB-junction. In these cases, optimization (2) is taken over all possible cells, and thus σ- and σ BR = 0.530, we evolve the initial configuration using the same numerical schemes as above. Note that 288 since the σ ij 's satisfy the triangle inequality condition, wetting cannot take place. However, Figure 4 and 289 Movie S2 show that in the evolution computed by Esedoḡlu-Otto algorithm alone, red cells grow in the 290 vicinity of the blue tricellular junction -an unnatural cell dynamics. This phenomenon persists even when 291 auction algorithm is incorporated, but is completely eliminated upon introducing the localization. complex whose representative marker, β-catenin indicates the differential adhesions required to drive self-297 organized cell movements in OE. Following the differential adhesion hypothesis (Steinberg, 1963; Foty and  Using the vertex dynamics model with σ ij = α −1 ij , they were able to confirm that the first case leads to 301 cellular intercalation, while the second does not; thereby, supporting their idea that differential adhesion in 302 heterotypic cell-cell junctions drives cell intercalations.

303
In order to compare the performance of the standard vertex dynamics algorithm and our scheme, we  To confirm whether this aspect of vertex dynamics is pertinent to topology changes such as intercalation,  This, on the contrary, does not occur in the level set-based approach, as it is theoretically proved to satisfy 342 the tricellular angle condition. Lastly, the vertex dynamics model can only take a restricted set of paths 343 to minimize its energy. Indeed, since the level set-based method allows for arbitrary shapes of junctions, 344 contrary to the vertex model where only polygonal shapes are allowed, it is able to follow the gradient descent by having to take only polygonal deformations.
We elaborate on the first difference mentioned in the previous paragraph. In our model, the energy 348 consists only of surface energy, and the volume preservation is realized through the constraint condition. 349 Thus, the cell sizes are preserved precisely. On the other hand, in the vertex model, the energy is defined 350 as the sum of surface energy and volume penalty term, while it is necessary to use a large value of penalty 351 coefficient ρ in order to preserve the cell sizes to some extent. The effect of the surface energy term is then 352 relatively weakened when ρ is large. In this way, the level set-based approach allows us to focus on the 353 effect of surface or adhesion energy, unlike the vertex model where it is difficult to separate and correctly 354 understand this effect (see Figure 6). This provides a typical example of the importance of reducing model  To conclude this section, we remark on some computational aspects of the level set-based method. First, 362 we discuss the jump-like behavior of the junction length in the level-set based method, which the reader might 363 have noticed in Figures 5 and 6. This is due to fact that the level set method is performed on a fixed grid  only when interfacial tensions satisfy σ BR > 1 2 (σ BB + σ RR ); meanwhile a sufficient condition for mixing is 401 σ BR < 1 2 σ BB or σ BR < 1 2 σ RR .

402
We recreate these simulations using our level set-based approach. Consider an initial aggregate of 50  pattern as in embryonic E18 stage (see Figure 11C and Movie S9).

490
We conclude the list of applications with two morphogenetic phenomena including a medium that were 491 already treated by other methods, and show that our method is able to obtain at least comparable outputs.  2002)) and generate its evolution using our level set-based scheme with time step size δt = 0.001. Observe 501 that the final configuration results in total engulfment of the blue cell mass by the red cell mass (see Figure   502 12 and Movie S10). This is consistent with theoretical considerations that in order to reduce the energy of 503 BR-junction, the engulfed blue mass becomes perfectly circular; as well as, the outside engulfing red mass.

504
Note that this is not easily achieved by finite element-based simulations, cf., Figure 8   The area (corresponds to "length" in a two-dimensional model) of a cell-cell junction γ ij can be estimated using the so-called "heat content approximation", which says that the area of the junction is proportional 669 to the heat that flows from cell C j to cell C i in a short time δt: Here G δt (x) = (4πδt) The approximate energy E δt is still nonlinear, so to devise a simple minimization scheme, we adopt 679 the approach of iterative minimization, which is justified by Lemma 5.2 of (Esedoḡlu and Otto, 2015).

680
Namely, we assume that we have an approximation of the minimizer u k of E δt in K and compute a updated Here L E δt is the linearized energy given by The functions ϕ k i correspond to the functions ψ k i in the main algorithm. It was proved in (Esedoḡlu and Otto, approximate energy E δt in every step and correctly approximates the L 2 -gradient flow of the original energy 687 E in the limit δt → 0.

688
In the case when there is no volume constraint, the minimization problem (A.6) becomes a problem of 689 minimizing a linear function over a simplex set K and thus the solution is obtained immediately as