A diffuse interface framework for modelling the evolution of multi-cell aggregates as a soft packing problem driven by the growth and division of cells
AA DIFFUSE INTERFACE FRAMEWORK FOR MODELLING THEEVOLUTION OF MULTI - CELL AGGREGATES AS A SOFT PACKINGPROBLEM DRIVEN BY THE GROWTH AND DIVISION OF CELLS P REPRINT
J. Jiang ∗ , K. Garikipati † & S. Rudraraju ‡ January 9, 2019 A BSTRACT
We present a model for cell growth, division and packing under soft constraints that arise from thedeformability of the cells as well as of a membrane that encloses them. Our treatment falls withinthe framework of diffuse interface methods, under which each cell is represented by a scalar phasefield and the zero level set of the phase field represents the cell membrane. One crucial element in thetreatment is the definition of a free energy density function that penalizes cell overlap, thus giving riseto a simple model of cell-cell contact. In order to properly represent cell packing and the associatedfree energy, we include a simplified representation of the anisotropic mechanical response of theunderlying cytoskeleton and cell membrane through penalization of the cell shape change. Numericalexamples demonstrate the evolution of multi-cell clusters, and of the total free energy of the clustersas a consequence of growth, division and packing.
Formation of multi-cell aggregates is a foundational process in the evolution of multicellular organisms. Beginningwith a single cell or a small cluster, the growth of aggregates is driven by cell division, differentiation, migrationand cell-cell interactions. Understanding the processes underlying the formation of these aggregates is central tomany phenomena in cellular biology and physiology, including embryogenesis, regeneration, wound healing, tissueengineering, and the growth and metastasis of cancerous tumors. As in most areas of cellular biology, a large body ofwork has focused on understanding the signalling pathways that control the evolution of multi-cell clusters, and it iswidely accepted that these pathways are triggered by the chemical environment and mechanical interactions (intra-cell,cell-cell, aggregate-matrix and other external stimuli). However, the understanding of the spatial and temporal variationsand effects of the chemo-mechanical environment in these cell aggregates remains at a nascent stage. Even in this earlystage, because of its complexity, the field is increasingly leaning on computational models.Early work on modeling growth and interactions in cell clusters included lattice models [10, 11, 23]. These treat cellsas sites on a square or hexagonal lattice and evolve multi-cell configurations through free energy-minimizing cell pairexchanges. These highly reduced order representations have provided significant insights into the effect of cell-cellinteractions on the evolution of cell aggregates. However, the cell-cell exchange processes assumed in these models arenot universally observed in real cell aggregates. Such treatments also limit the incorporation of sub-cellular growthdynamics. Improvement of the lattice models in the form of sub-cellular lattice models using high-Q Potts models havedelivered better geometric representation of cell structure [13, 9]. In sub-lattice models, the cells are represented by acluster of lattice sites rather than a single lattice site, and cell migration is achieved by switching parent cell identity ofthe lattice sites at the boundary. A single such switch allows the cell boundary of one cell to advance by one latticelength into the neighbouring cell. These methods have allowed for a finer representation of the cell geometry, but resultin unrealistic jagged cell boundaries and complex cell shapes that are not simply connected (e.g. cells within cells). ∗ Mechanical Engineering, University of Michigan † Mechanical Engineering, and Mathematics, University of Michigan; corresponding author ‡ Mechanical Engineering, University of Wisconsin-Madison a r X i v : . [ q - b i o . CB ] J a n REPRINT - J
ANUARY
9, 2019The drawbacks of the lattice models consisting of their non-representative cell geometries or a jagged representation ofthe cell boundaries was partially addressed with the development of cell-centric/center-dynamics models [15, 16, 12, 24]and vertex dynamics models [16, 17, 7, 2, 25]. The center dynamics models approximate cell shapes as polygonsgenerated though Voronoi tesselation of a collection of forming points, and the evolution of cell boundaries is achievedthrough a free energy minimizing movement of the forming points. A major drawback of this method is the restrictionof the boundary to a polygonal shape dictated by the underlying tessellation. This restriction was addressed in vertexdynamics models that allowed for the cells to be represented as general polygons defined by the connectivity of thevertices. This connectivity evolves dynamically and is driven by the free energy minimizing pair wise movement ofvertices that conserve the cell volume.As detailed in the review paper by Brodland [5], the lattice models, cell-centric models and the vertex dynamics modelshave successfully modeled a wide range of cell-cell interactions and cell aggregates phenomena. However, these modelshave no representation or a very simplified representation of cell geometries, cell-cell and cell-substrate interactions,cytoskeletal remodeling, cytoplasmic viscosity, etc. This greatly limits the ability to model realistic, smooth andanisotropic cell shape evolution, the mechanics of cell surface evolution due to cell-cell contact, the ability to controlcell volumes and to ensure proper dissipative dynamics. Modeling and understanding these processes is central tocharacterizing the process of growth and evolution of multi-cell aggregates that we refer to as a problem of soft packing.In this manuscript, we present a finite element method based phase field representation of cells and the resulting softpacking dynamics of cell aggregates. The phase field method is a popular numerical technique for simulating diffuseinterface kinetics at the mesoscale and has been widely used to model evolving interface problems such as crystalgrowth, solidification and phase transformations in alloys. In this method, the evolution of a species concentrationand/or phase is modelled using a set of conserved or non-conserved order parameters. The evolution of the orderparameters and the corresponding interface kinetics are governed by a system of parabolic partial differential equations,which are referred to as the Cahn-Hilliard formulation (for conserved order parameters) [6] and Allen-Cahn formulation(for non-conserved order parameters)[1]. The phase field representation of the cell geometries and the soft packing ofmulti-cell aggregates presented here allows for an improved representation of realistic cell shapes and the modelingof cell growth, division, and mechanical compaction processes intrinsic to the formation of multi-cell aggregates.An earlier attempt at modeling multi-cell aggregates using a phase field representation was outlined by Nonomura[26] using an Allen-Cahn representation of the cells; i.e. a non-conserved order parameter treatment. In contrast, theformulation presented in this paper treats cell mass as a conserved quantity and models the evolving cell clusters usinga Cahn-Hilliard representation. Furthermore, our treatment considers the mechanics of soft packing and allows for thenecessary anisotropic shape evolution of cells.In Section 2, we present the phase field formulation, its numerical implementation and simulations of cell growth anddivision. This is followed by the modeling of mechanics of soft packing and simulations of soft packing in Section 3.Closing remarks appear in Section 4.
Our formulation of the problem rests on a phase field representation with as many scalar fields as cells. The treatment iscentered on the definition of a free energy density, as a function of the scalar fields. In the non-mechanical version of theproblem, terms are constructed to model cell membranes by phase segregation of cell interiors from the extra-cellularmatrix, and contact inhibition, or intercellular repulsion, by penalizing overlapping scalar fields. We present thevariational treatment, the mechanisms that model cell division, numerical aspects, and an illustrative numerical example.
Let Ω ∈ R with a smooth boundary ∂ Ω . Scalar fields c k , k = 1 , . . . , N with c k ∈ [0 , serve to delineate the interiorand exterior of the cell numbered k . Here, the interior of cell k is ω k ⊂ Ω , where ω k = { X ∈ Ω | c k ( X ) = 1 } . Theexterior is Ω \ ω k . The free energy density function is built up beginning with the following form: ψ ( c k ) = αc k ( c k − + κ | ∇ c k | (1)where the double-well term, f ( c k ) = αc k ( c k − , enforces segregation into ω k and Ω \ ω k . In Equation (1), the secondterm enforces a diffuse cell-matrix interface (the cell membrane) of finite thickness, where κ controls the interfacethickness, and thereby the interfacial energy. For N cells in Ω , the above free energy density needs to be extended tomodel contact by adding a cell-cell repulsion term. The total free energy of the multi-cell aggregate is a functional Π[ c ] ,2 REPRINT - J
ANUARY
9, 2019defined as Π[ c ] := (cid:90) Ω ψ ( c , ∇ c ) d V = (cid:90) Ω N (cid:88) k =1 f ( c k ) + N (cid:88) k =1 κ | ∇ c k | + (cid:88) l (cid:54) = k N (cid:88) k =1 λc k c l d V. (2)Here, c = { c , . . . , c N } , and λ is a penalty coefficient that enforces repulsion between any two cells k, l thus modellingcell contact.Taking the variational derivative with respect to c k in Equation (2) yields δ Π k [ c ; w ] = dd (cid:15) (cid:90) Ω N (cid:88) k =1 f ( c k + (cid:15)w ) + κ | ∇ ( c k + (cid:15)w ) | + (cid:88) l (cid:54) = k λ ( c k + (cid:15)w ) c l d V (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = (cid:90) Ω w f (cid:48) ( c k ) − κ ∆ c k + (cid:88) l (cid:54) = k λc k c l dV + (cid:90) ∂ Ω wκ ∇ c k · n dS (3)where n is the unit outward normal vector to ∂ Ω . The chemical potential of the k th cell is identified as, µ k = f (cid:48) ( c k ) − κ ∆ c k + (cid:88) l (cid:54) = k λc k c l (4)At equilibrium, δ k Π[ c ; w ] = 0 for the k th cell, yielding µ k = 0 in Ω , and κ ∇ c k · n = 0 on ∂ Ω . The following parabolic partial differential equation, popularly known as the Cahn-Hilliard equation [6], imposes theconserved dynamics that governs the delineation and growth of the N − cell agglomerate, and of repulsion between cellpairs: ∂c k ∂t = − ∇ · ( − M ∇ µ k ) + s k (5)where the source term s k has been introduced, and M is the mobility, assumed to be constant. The dynamics of themulti-cell soft packing problem is governed by Equation (5) with the thermodynamics prescribed by Equation (4) andboundary condition κ ∇ c k · n = 0 on ∂ Ω for k = 1 , . . . N . Time discretization is carried out by the implicit, backward Euler method. Time instants are indexed by superscriptsin the following development, and the time step is denoted by ∆ t = t n +1 − t n . Starting with the initial conditions { c k , µ k } , and given the solution { c nk , µ nk } at time t n , the time-discrete versions of Equations (5) and (4) are, c n +1 k = c nk + ∆ t ( M ∇ · ( ∇ µ n +1 k ) + s k ) where µ n +1 k = f (cid:48) n +1 ( c k ) − κ ∆ c n +1 k + (cid:88) l (cid:54) = k λc n +1 k c n +1 l (6)The weak form, posed as a classical two-field mixed finite element formulation [4], is stated as follows: Find c n +1 k ∈ S = { c | c ∈ H (Ω) , ∇ c · n = 0 on ∂ Ω } and µ n +1 k ∈ S = { µ | µ ∈ H (Ω) } such that for all variations w ∈ V = { w | w ∈ H (Ω) } on c k and v ∈ V = { v | v ∈ H (Ω) } on µ k , respectively, the following residual equations In some simulations, a buffer zone, Ω (cid:48) , is needed around the simulation domain to inhibit unrealistic cell shapes resulting fromthe enforcement of κ ∇ c k · n = 0 on ∂ Ω . In addition, this buffer zone also acts as a membrane around the cell cluster. In the bufferzone, an additional term of the form (cid:80) Nk =1 λc k is added to the free energy density to penalize the movement of any cells from theactive simulation domain ( Ω ) to the buffer zone ( Ω (cid:48) ). In this initial boundary value problem, ∂ Ω is composed of only a Neumann (flux) boundary ( ∂ Ω h = ∂ Ω ) and hence the Dirichletboundary is a null set ( ∂ Ω g = ∅ ). REPRINT - J
ANUARY
9, 2019are satisfied: (cid:90) Ω wc n +1 k d V = (cid:90) Ω ( wc nk − ∇ w · ∆ tM ∇ µ n +1 k + w ∆ ts k ) d V (cid:90) Ω vµ n +1 k d V = (cid:90) Ω ( vf (cid:48) n +1 ( c k ) + ∇ v · κ ∇ c n +1 k ) dV + (cid:90) Ω v (cid:88) l (cid:54) = k λc n +1 k c n +1 l d V (7)Spatial discretization is implemented in a standard finite element framework and uses bilinear quadrilateral elementsleading to standard matrix-vector forms of the equations in (7). Two criteria of cell division are considered: Age-based and mass-based. In the age-based criterion, the cells divideafter reaching a specified physical age irrespective of their current mass. In the mass-based criterion, the cells divideafter reaching twice their mass at birth irrespective of the time it takes to reach that mass. If the cell division processis spatially symmetric and leads to bisection of the parent cell, and if the daughter cells have the same growth rate,then both these criteria can lead to a sequence of divisions at the same time instants. However, asymmetrical celldivision is common due to small perturbations, which, in our model, are introduced by small numerical differences.As a result, different sequences of evolution of the cell clusters can result with the two criteria. Furthermore, thedistinction introduced by absence/presence of the source term, s k , with the time-/mass-based criteria, respectively, inour implementation (see Section 2.7), contributes more significant differences. The implementation of these criteria isoutlined in Algorithm 1, and the process of cell division proceeds as follows: For an elliptical cell dividing along its minorprincipal axis, ω n +1 N +1 is defined such that ω nk is bisected into ω n +1 k and ω n +1 N +1 , ensuring meas ( ω n +1 k ) = 0 . meas ( ω nk ) and meas ( ω n +1 N +1 ) = 0 . meas ( ω nk ) , where meas ( ω k ) is the volume (area in two dimensions) of ω k . The regions ω n +1 k and ω n +1 N +1 are thus determined by the divisions of elliptical cells along their minor principal axes. This strategy ofbisecting an elliptical cell along its minor axes is motivated by observations of cell division in biological cells. However,recognizing that the shape of ω nk will, in general deviate from an ellipse, we define the division axis to lie along theaxis of the major principal moment of inertia through the center of mass. Note that, for an ellipse, the minor principalaxis, and the major principal moment of inertia axis through the center of mass coincide. At time t n +1 , a new interfaceforms between ω n +1 k and ω n +1 N +1 following a division of the k th cell at time t n that incremented the number of cells N (cid:55)→ N + 1 . This new interface is along the major principal axis of the moment of inertia tensor through the center ofmass of ω nk ; that is, of the k th cell’s interior at time t n .Numerical implementation of the division process involves introducing two new fields, c N +1 and µ N +1 , over theproblem domain to represent the new cell ω n +1 N +1 . This is followed by initializing c N +1 and resetting c k as follows: c n +1 N +1 ( X ) = c nk ( X ) ∀ X ∈ ω n +1 N +1 c n +1 N +1 ( X ) = 0 . ∀ X ∈ ω n +1 k c n +1 k ( X ) = 0 . ∀ X ∈ ω n +1 N +1 So computationally, during every cell division, two new fields are introduced, and the total number of degrees offreedom in the problem is increased by N nodes , where N nodes is the number of nodes/mesh-points in the discretizationof the problem domain. This would cause the number of degrees of freedom to linearly increase with the number ofcells, becoming computationally expensive for simulations with large clusters of cells. However, this limitation can beeliminated by active parameter tracking methods developed in the phase field community [31]. These methods leveragegraph theory to introduce “coloring schemes” that can represent arbitrary numbers of cells with a finite number of fields(typically less than 10). This is accomplished by restricting the support of each field to a neighbourhood that is onlyslightly larger than the corresponding cell, instead of the entire problem domain. A uniform time step, ∆ t = ∆ t , is chosen such that it ensures the stability and convergence of the Cahn-Hilliarddynamics. However, when cell divisions occur, adaptive time step control is necessary to equilibrate the large penaltyforces arising from sharply varying composition fields over the new daughter cells’ common boundary. To address theassociated, transient numerical stiffness, ∆ t is decreased by a factor of . × − m for a few time steps ( n div ) in orderto ensure convergence of the nonlinear system of equations. The time step is then sharply increased back to ∆ t = ∆ t until the next cell division process. Algorithm 2 details the implementation of this adaptive time step control in ourcode. Typical values used for the simulations presented in this paper are n div < and m < .4 REPRINT - J
ANUARY
9, 2019
Algorithm 1:
Cell division mechanisms, given < ε < ε (cid:28) Age based division criterion: T age k is the current age of the cell ω k , and T division k is the age at division. if ε < T division k /T age k − < ε then N (cid:55)→ N + 1 ;meas ( ω n +1 k ) = 0 . meas ( ω nk ) ;meas ( ω n +1 N +1 ) = 0 . meas ( ω nk ) ; end Mass based division criterion: m nk is the current mass of the cell ω k , given by: m nk = (cid:82) Ω c nk dV ; if ε < m nk / m k − < ε then N (cid:55)→ N + 1 ;meas ( ω n +1 k ) = 0 . meas ( ω nk ) ;meas ( ω n +1 N +1 ) = 0 . meas ( ω nk ) ; endAlgorithm 2: Adaptive time step control ∆ t = ∆ t ; for t ← to T ; t n +1 = t n + ∆ t do /* T is the total time of computation; ∆ t is the time step */ if cell is about to divide then ∆ t = ∆ t × . × − m ; a = 0 ; if ∆ t < ∆ t then a = a + 1 ; if a > n div then ∆ t = ∆ t × . × − m × ( a − n div ) m ; endendendend2.5 Adaptive mass source The growth of cells is controlled by the source term, s k , which can be determined from experimentally observed celldoubling time estimates. In this work, s k is a function of the mass ratio ν nk = m nk /m k : The default value of s k = s ,the average growth rate. Growth continues at this rate until the cell has doubled in mass and then the division processseparates the cell into two daughter cells. The two daughter cells then continue to grow at the rate s . However, the totalnumber of cells in the limit of optimal soft packing is given by N max = meas (Ω) /m k . If either the number of cellsapproaches N max or if no new scalar fields are available to initialize new daughter cells leading to ν nk > . , we turn offthe source, s n +1 k = 0 , so that cell k no longer grows. This adaptive mass source control appears in Algorithm 3. The two-dimensional, cell growth and soft packing formulation presented here has been implemented in the
C++ based deal.II open source parallel finite element library [3]. We use the
SuperLU direct solver [20] to solve the system oflinear equations obtained from the linearization of Equation (7). The linearization itself is obtained using the
Sacado algorithmic differentiation library of the open source
Trilinos project [14].5
REPRINT - J
ANUARY
9, 2019
Algorithm 3:
Adaptive mass source if X ∈ ω nk ; in cell k then s nk = s ; ν nk = m nk /m k ; /* ratio of current cell mass to its initial mass */ if N ≥ N max ; new phase field is not available thenif ν nk > . then s nk = 0 . ; endendend2.7 Illustrative examples for the progression of cell growth, division and compaction. We consider two cases for demonstrating the evolution of cell growth, division and contact leading to compaction: • Cell division at a constant volume , s k = 0 : Early stages of embryonic cell division in many organisms (e.g.Caenorhabditis elegans [8]) is known to occur at a constant embryonic volume. During these early stages,there is an increase in cell numbers, but the overall embryo volume stays fixed. This process is modeled inFigure 1, where a single circular cell divides into N max = 16 cells while maintaining the total cell volumefixed. Here we use the age based criterion for cell division as the mass of each daughter cell does not changeover time. • Cell division with a source , s k (cid:54) = 0 : In species where the early embryo grows in size, the total cytoplasmicvolume of all the cells increases with time. This process of growth of the embryonic cells can be modeledvia a positive source term in Equation (5). The corresponding simulation for N max = 12 is shown in Figure2. Here we use the mass based criterion for cell division as the positive source term leads to growth of thedaughter cells.The parameters used in these numerical examples appear in Table 1. Our computations begin with a single cell, ω ,of circular shape at the center of an elliptical “embryo” Ω . Each division axis is determined by the major principaldirection of the moment of inertia tensor through the center of mass, which distributes mass evenly in each daughter cell.This can be observed from the progression in Figures 1b–1f for division at constant total cell volume, and in Figures2a–2d for division and growth. Attention is also drawn to the delineation of daughter cells following each division, andof non-sibling cells from each other due to the repulsion terms in Equation (2).Figure 3 tracks the total free energy of the system for the case of cell division and growth with a source shown in Figure2. Note that the total free energy increases with time due to the mass supply, and that transient fluctuations occur at eachcell division due to the formation of a sharp boundary between the daughter cells and the transiently stronger repulsionbetween them. When the source term vanishes (after the last division and corresponding energy spike), the free energymonotonically decreases, as dictated by the dissipative nature of phase field dynamics. The initial conditions for thephase field for the computations in Figures 1 and 2 include a small random perturbation from a mean value to help drivethe Cahn-Hilliard dynamics. This perturbation results in a tilt in the major principal direction of the moment of inertiatensor, and asymmetric division of the cells ( ∼ . radians off the axes of symmetry).The computations in Figures 1 and 2 were performed on an elliptical domain (mesh shown in Figure 1a for the constantvolume case) with 5120 elements and 5185 nodes. Since each cell is represented by two scalar fields, c and µ , thedegrees of freedom associated with each cell are 10370. For the simulation shown in Figure 1 involving sixteen cells,this resulted in 165,920 total degrees of freedom.Table 1: Numerical values of parametersParameters α κ λ M s ∆ t Constant volume . − . − With source . − . . − REPRINT - J
ANUARY
9, 2019 (a) Spatial discretization. (b) Initial circular cell. (c) Division to two cells.(d) Division to four cells. (e) Division to eight cells. (f) Division to sixteen cells and com-paction.
Figure 1: A demonstration of the progression of cell division from one cell into sixteen cells at a constant total cellvolume. The underlying spatial mesh/discretization over the simulation domain is shown in (a). Cell interiors are shownin red and the cell membrane in cyan-yellow. (a) Initial circular cell. (b) Progression to four cells.(c) Progression to eight cells. (d) Progression to twelve cells.
Figure 2: A demonstration of the progression of cell division and growth from one cell into twelve cells driven by asource term. Cell interiors are shown in red and the cell membrane in cyan-yellow.
We now present an extension of the formulation to an elementary mechanical model that associates energy to globalshape changes of a cell. The model works by penalizing departures of the principal moments of inertia from theirinitial values, which are established when a specific cell is created. A more complete treatment of mechanics in the cellgrowth, division and soft packing problem would consider nonlinearly elastic interactions between cells. However,as we demonstrate, the simple shape change-based model does capture the essential mechanics of soft packing. Weinclude a study of the reduced moduli that govern shape change via a parametric study, before studying a sequence ofcell growth, divisions and soft packing with the added mechanical contributions.7
REPRINT - J
ANUARY
9, 2019Figure 3: Evolution of total free energy with time (normalized). Each spike in the energy curve corresponds to transientrepulsion between newly formed daughter cells following a cell division. The corresponding cell division events areshown in the inset sub-figures.
Consider the addition of the following term to the free energy density (2) in order to account for the principal momentsof inertia of a single cell: Π MI = dim (cid:88) i =1 δ i ( I i ref − I i ) (8)Here, δ i is a mechanical modulus penalizing variations in the i th principal moment of inertia, I i , from its reference state, I i ref , established at the birth of this cell. With this addition that accounts for changes in shape, the total free energyfunction takes on the form: Π[ c , ∇ c ] := (cid:90) Ω N (cid:88) k =1 f ( c k ) + N (cid:88) k =1 κ | ∇ c k | + (cid:88) l (cid:54) = k N (cid:88) k =1 λc k c l d V + N (cid:88) k =1 dim (cid:88) i =1 δ ik ( I ik ref − I ik ) (9)The moment of inertia tensor through the center of mass, I , is expressed in terms of its Cartesian components I ij , withthe phase field playing the role of the mass density in the traditional definition of this quantity: I [ c ] = (cid:20) I [ c ] I [ c ] I [ c ] I [ c ] (cid:21) = (cid:20) (cid:82) c ¯ X d V − (cid:82) c ¯ X ¯ X d V − (cid:82) c ¯ X ¯ X d V (cid:82) c ¯ X d V (cid:21) (10)where, with the center of mass X ci = (cid:82) Ω cX i d V (cid:82) Ω c d V , (11)the Cartesian coordinates relative to the center of mass are ¯ X i = X i − X ci . Then, the principal moments of inertiathrough the center of mass, I i , can be determined from the moment of inertia tensor in Equation (10) by an eigendecomposition governed by the Cayley-Hamilton Theorem: I k − I k tr I k + det I k = 0 (12)where I k denotes a principal value of the moment of inertia tensor, relative to the center of mass, of the k th cell, I k .8 REPRINT - J
ANUARY
9, 2019The variational machinery applied to computing the chemical potential now must account for the added, shape-dependent,moment of inertia terms: δ Π k [ c ; w ] = dd (cid:15) (cid:90) Ω N (cid:88) k =1 f ( c k + (cid:15)w ) + κ | ∇ ( c k + (cid:15)w ) | + (cid:88) l (cid:54) = k λ ( c k + (cid:15)w ) c l d V + dd (cid:15) N (cid:88) k =1 dim (cid:88) i =1 δ ik (cid:16) I ik ref − I ik [ c k + (cid:15)w ] (cid:17) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = N (cid:88) k =1 (cid:90) Ω w f (cid:48) ( c k ) − κ ∆ c k + (cid:88) l (cid:54) = k λc k c l d V − N (cid:88) k =1 dim (cid:88) i =1 δ ik (cid:16) I ik ref − I ik [ c k ] (cid:17) ˜ I ik [ c k ]+ (cid:90) ∂ Ω wκ ∇ c k · n dS (13)where, the term ˜ I ik arises from the variation of I ik , and is obtained from the Cayley-Hamilton Theorem to be: ˜ I ik = (cid:82) Ω w (cid:18) I ik tr ¯ I − ¯ X (cid:82) Ω c ¯ X d V − ¯ X (cid:82) Ω c ¯ X d V + 2 ¯ X ¯ X (cid:82) Ω c ¯ X ¯ X d V (cid:19) d V I ik − tr I , (14)with the tensor ¯ I = (cid:20) ¯ X − ¯ X ¯ X − ¯ X ¯ X ¯ X (cid:21) . (15)Extending Equation (4), the chemical potential is now defined as µ k = f (cid:48) ( c k ) − κ ∆ c k + (cid:88) l (cid:54) = k λc k c l − dim (cid:88) i =1 δ ik (cid:16) I ik ref − I ik (cid:17) ˆ I ik , (16)where ˆ I ik is ˆ I ik = I ik tr ¯ I − ¯ X (cid:82) Ω c ¯ X d V − ¯ X (cid:82) Ω c ¯ X d V + 2 ¯ X ¯ X (cid:82) Ω c ¯ X ¯ X d V I ik − tr I , (17)thus introducing a simple mechanical model that penalizes shape changes of the cell. When combined with thegoverning parabolic partial differential equation in conservation form (5), and the boundary condition κ ∇ c k · n = 0 on ∂ Ω for k = 1 , . . . N , we have a description for multi-cell growth, division and soft packing with the penalization ofshape change. The time discretized dynamics are now written as an explicit-implicit scheme. Given the initial conditions { c k , µ k } andthe solution { c nk , µ nk } , the time-discrete versions of Equations (5) and (16) are, c n +1 k = c nk + ∆ t ( M ∇ · ( ∇ µ n +1 k ) + s k ) µ n +1 k = f ,c n +1 k − κ ∆ c n +1 k + (cid:88) l (cid:54) = k λc n +1 k c n +1 l − dim (cid:88) i =1 δ ik (cid:16) I ik ref − I ik (cid:17) ˆ I ik . (18)Note that the mechanical terms exerting control over the shape do not carry time instant superscripts in Equation (18)for the sake of brevity. We have experimented with evaluating these terms at t n +1 in a fully implicit method, as well asat t n , in an explicit-implicit implementation. The complicated functional evaluations in Equations (14)-(17) make thefully implicit method notably more expensive than the explicit-implicit method, which has been used in the numericalexamples of Section 3.4. 9 REPRINT - J
ANUARY
9, 2019The weak form is stated as follows: Find c n +1 k ∈ S = { c | c ∈ H (Ω) , ∇ c · n = 0 on ∂ Ω } and µ n +1 k ∈ S = { µ | µ ∈ H (Ω) } such that for all variations w ∈ V = { w | w ∈ H (Ω) } on c k and v ∈ V = { v | v ∈ H (Ω) } on µ k ,respectively, the following residual equations are satisfied: (cid:90) Ω wc n +1 k dV = (cid:90) Ω ( wc nk − ∇ w · ∆ tM ∇ µ n +1 k + ws k ) dV (cid:90) Ω vµ n +1 k dV = (cid:90) Ω ( vf n,c k + ∇ v · κ ∇ c nk ) dV + (cid:90) Ω v (cid:88) k (cid:54) = l λc k c l dV − dim (cid:88) i =1 δ ik (cid:16) I ik ref − I ik (cid:17) (cid:90) Ω v ˆ I ik dV (19) In our numerical experiments with the above model for mechanical control of shape change in packing, we havefound that the effective values of the mechanical moduli δ ik that impose control on cell shape fall within the range δ ik ∈ [1 × , × ] . Thus, with δ k → and δ k in the above range, the model produces elliptical single cells ofincreasing aspect ratio with the major axis corresponding to the smaller principal moment of inertia I k < I k . Thisparametric study appears in Figure 4. Much larger values δ ik > × make the problem stiff, and convergence of thenonlinear solver becomes difficult. Other aspects of the numerical implentation remain the same as in Sections 2.2-2.5. In order to illustrate the control of cell shape by the mechanical moduli δ ik , we set δ k → and δ k ∈ [1 × , × ] .The results appear in Figure 4 for a single cell. The modulus δ takes on values × , × , and × in Figures4a, 4b and 4c, respectively. It can be seen from these three figures that the larger mechanical modulus, δ constrainsgrowth along the minor principal axis of the elliptical cell that forms. As δ becomes larger, it tends to produce anoblate shape, as the final equilibrium state in Figure 4c. In comparison, at the lower value of δ = 5 × , this modulusdoes not much affect the ellipticity of cell shape, as seen in Figure 4a.The progression of cell division from a single mother cell into twelve daughter cells is demonstrated in Figure 5. Notethe tight packing and shape changes attained by the cells. We also draw attention to the differences in cell shapesbetween Figures 2 and 5. This difference is especially notable at the 12-cell stage, and is due to the added penalizationof cell shape change already demonstrated in Figure 4. We expect, also, that the shapes in Figure 5 are more physicallyaccurate because they account for shape change, albeit by a simple mechanical model.The accompanying evolution of the total free energy is shown in Figure 6. As in Figure 3, transient fluctuations occurwith each cell division due to the formation of a sharp boundary between the daughter cells, and the transiently strongerrepulsion between the daughter cells. However, in this case, the transient fluctuations are more prominent and spreadout due to the rapid changes in daughter cell shapes that themselves result from the elastic repulsion following division.Compared to Figure 3, the free energy values in Figure 6 are much higher. This is due to the additional penalization inthe form of the mechanics term whose relative magnitude has been made higher than the regular Cahn-Hilliard andoverlap penalization terms. As observed before, the height of a spike and its width on the time axis increase with thenumber of cell divisions occurring in that time interval and the mass source causes the gradual increase in the freeenergy over time. In this paper we have presented a phase field-based diffuse interface framework for modeling the growth, divisionand packing of multi-cell aggregates. The model allows for high fidelity representation of smooth, anisotropic andunsymmetric cell geometries, cell-cell contact and the resulting mechanical compaction processes intrinsic to softpacking in multi-cell aggregates. The cells are delineated by conserved scalar phase fields driven by the growth,division and compaction processes. The zero level set of the phase field representing a given cell also identifies thecell membrane. The driving force is the minimization of a free energy functional, and the governing equations arevariationally derived to yield a system of parabolic partial differential equations. The salient features of the formulationare: 10
REPRINT - J
ANUARY
9, 2019 (a) δ = 50000 (b) δ = 100000 (c) δ = 200000 Figure 4: Fully developed single cell shape change due to anisotropic mechanical moduli δ → and δ ∈ [1 × , × ] . (a) Mother cell (b) Four daughter cells(c) Eight daughter cells (d) Twelve daughter cells Figure 5: Cell division with penalization of cell shape change. The shapes of the 12-cell cluster differ from the 12-cellcluster in Figure 2 due to the additional effect of mechanics. Cell interiors are shown in red and the cell membrane incyan-yellow. • Being a field formulation, no discreet interface evolving mechanisms are needed, such as those employed inlattice, cell-centric and vertex dynamics models. • The dynamics occurs at time scales controlled by physically meaningful mechanisms: the growth rate and thedoubling time of cells. Absent are severely reduced (by orders of magnitude) time steps needed to equilibratethe very high frequency response of vertices in, e.g., cellular automata models. We note also that such, veryhigh frequency response, is not physically realistic in the highly dissipative setting of cell biophysics. On theother hand, the first-order dynamics of phase field models does impose dissipation that can be tuned to therelevant timescale. • The cell boundary is represented by a continuously differentiable zero level set of the phase field and can thusrepresent general, smooth cell shapes without being limited to polygons or a jagged representation of the cellboundary. 11
REPRINT - J
ANUARY
9, 2019Figure 6: Evolution of the total free energy with time (normalized). Each spike in the energy curve corresponds totransient repulsion between newly formed daughter cells following a cell division, and the corresponding cell divisionevents are shown in the inset sub-figures. • The mechanics of soft packing is modelled by penalizing departures of the principal moments of inertia of eachcell from their initial values. This is a very simplified representation of the anisotropic mechanical response ofthe underlying cytoskeleton and cell membrane.These advantages significantly differentiate our model from cellular automata based cell-centric and vertexdynamics models used in the popular frameworks for modelling cell populations such as Cell-based Chaste [22]and CompuCell3D [30]. However, the increased accuracy of cell shapes and mechanics does come at a highercomputational cost needed to solve coupled partial differential equations over the problem domain. The targets ofthis numerical framework are applications requiring an accurate representation of cell shapes and modeling of longrange mechanical interactions. In biological processes involving packing of cells in a confined space, interestingspatial differentiation/patterning of cell clusters are observed (e.g. early stages of Zygote-Morula-Balstocyst-Gastrulaevolution in embryonic development [18], stem cell heterogeneity in tumors [27] and tumor shape evolution [21, 28]).Intrinsic to the process of packing is the effect of mechanical interactions of the cells among themselves and with theextra cellular environment. These phenomena operate on length scales comparable to the dimensions of the cells, andalso introduce the effect of cell shape in the packing dynamics. A case in point is the early stage of embryogenesis(Zygote to the Morula and the Blastocyst) that involves the growth of a cell cluster beginning with a single cell thatdivides and grows to a few hundred. An attempt at modelling some of these processes (cleavage and compaction) in theearly stage of embryogenesis using this framework is shown in Figure 7. The cell packing morphology predicted by ourcomputational framework is compared to experimentally observed embryo morphology patterns in nematodes. As canbe seen, the close correlation of cell volumes, shapes and their spatial distribution is encouraging.Currently, the framework does not include a model for cell migration and only has a very simplified repre-sentation of cell-cell mechanical interactions through the penalization of cell shape. However, migration can bemodelled by incorporating an advection term. The authors have demonstrated such models previously for stress-drivenmovement of tumor stem cells [28], and recently for neuronal migration in the developing brain [32]. We note, also,that the simplified mechanics model, when applied to the shape changes of a sphere into an arbitrary ellipsoid, hasreproduced the essential features of the neo-Hookean strain energy function in our studies (to be reported in a futurecommunication). Incorporation of pointwise constitutive models (e.g. general hyperelastic models, viscoelasticity andporoelasticity), required for a physically complete treatment of mechanics, also needs the formulation of the governingequations of mechanics in a fully Eulerian setting [19]. Such extensions will be incorporated into future developmentsof this framework. Lastly, the computational cost of the current model can be significantly reduced by removing thelinear dependence of the number of degrees of freedom on the number of cells using the active parameter trackingmethod, as described in Section 2.3. This, also, will be presented in a future communication.12
REPRINT - J
ANUARY
9, 2019 (a) Early embryogenesis of nematodes up to the eight-cell stage observed in
C. elegans , (a-c), and
Prionchulus sp. , (g-i). Some of the cells are either partially or completely hidden from view in (c) and(i). Figure reproduced from Schierenberg [29] (Original figure distributed under Creative CommonsAttribution License).(b) Computations showing up to the third generation of divisions resulting in eight cells at a constant totalcell volume.
Figure 7: Comparison of cell morphologies during early embryogenesis observed in (a) nematodes (
C. elegans and
Prionchulus sp. ), and (b) the computational model.
References [1] Samuel M. Allen and John W. Cahn. A microscopic theory for antiphase boundary motion and its application toantiphase domain coarsening.
Acta Metallurgica , 27(6):1085 – 1095, 1979.[2] Silvanus Alt, Poulami Ganguly, and Guillaume Salbreux. Vertex models: from cell mechanics to tissue mor-phogenesis.
Philosophical Transactions of the Royal Society of London B: Biological Sciences , 372(1720),2017.[3] W. Bangerth, R. Hartmann, and G. Kanschat. deal.II – a general purpose object oriented finite element library.
ACM Trans. Math. Softw. , 33(4):24/1–24/27, 2007.[4] F. Brezzi and M. Fortin.
Mixed and Hybrid Finite Element Methods . Springer-Verlag, 1991.[5] G. W. Brodland. Computational modeling of cell sorting, tissue engulfment, and related phenomena: A review.
Applied Mechanics Reviews , 57:47–76, 2004.[6] John W. Cahn and John E. Hilliard. Free energy of a nonuniform system. i. interfacial free energy.
The Journal ofChemical Physics , 28(2):258–267, 1958.[7] Alexander Fletcher, Miriam Osterfield, Ruth E. Baker, and Stanislav Y. Shvartsman. Vertex models of epithelialmorphogenesis.
Biophysical Journal , 106(11):2291 – 2304, 2014.[8] S.F. Gilbert.
Developmental Biology. th edition. Sinauer Associates, 2000.[9] James A. Glazier and Francois Graner. Simulation of the differential adhesion driven rearrangement of biologicalcells.
Phys. Rev. E , 47:2128–2154, Mar 1993.[10] Narendra Goel, Richard D. Campbell, Richard Gordon, Robert Rosen, Hugo Martinez, and Martynas Yaas.Self-sorting of isotropic cells.
Journal of Theoretical Biology , 28(3):423 – 468, 1970.13
REPRINT - J
ANUARY
9, 2019[11] Narendra S. Goel and Gary Rogers. Computer simulation of engulfment and other movements of embryonictissues.
Journal of Theoretical Biology , 71(1):103 – 140, 1978.[12] Francois Graner. Can surface adhesion drive cell-rearrangement? part i: Biological cell-sorting.
Journal ofTheoretical Biology , 164(4):455 – 476, 1993.[13] Francois Graner and James A. Glazier. Simulation of biological cell sorting using a two-dimensional extendedpotts model.
Phys. Rev. Lett. , 69:2013–2016, Sep 1992.[14] Michael Heroux, Roscoe Bartlett, Vicki Howle Robert Hoekstra, Jonathan Hu, Tamara Kolda, Richard Lehoucq,Kevin Long, Roger Pawlowski, Eric Phipps, Andrew Salinger, Heidi Thornquist, Ray Tuminaro, James Wil-lenbring, and Alan Williams. An Overview of Trilinos. Technical Report SAND2003-2927, Sandia NationalLaboratories, 2003.[15] Hisao Honda. Description of cellular patterns by dirichlet domains: The two-dimensional case.
Journal ofTheoretical Biology , 72(3):523 – 543, 1978.[16] Hisao Honda. Geometrical models for cells in tissues.
International Review of Cytology , 81:191 – 248, 1983.[17] Hisao Honda, Hachiro Yamanaka, and Goro Eguchi. Transformation of a polygonal cellular pattern during sexualmaturation of the avian oviduct epithelium: computer simulation.
Development , 98(1):1–19, 1986.[18] Joseph Itskovitz-Eldor, Maya Schuldiner, Dorit Karsenti, Amir Eden, Ofra Yanuka, Michal Amit, Hermona Soreq,and Nissim Benvenisty. Differentiation of human embryonic stem cells into embryoid bodies compromising thethree embryonic germ layers.
Molecular medicine , 6(2):88, 2000.[19] K. Kamrin, C. H. Rycroft, and J. C. Nave. Reference map technique for finite-strain elasticity and fluid–solidinteraction.
Journal of the Mechanics and Physics of Solids , 60(11):1952–1969, 2012.[20] Xiaoye S. Li. An overview of SuperLU: Algorithms, implementation, and user interface.
ACM Transactions onMathematical Software , 31(3):302–325, September 2005.[21] K. L. Mills, R. Kemkemer, S. Rudraraju, and K. Garikipati. Elastic free energy drives the shape of prevascularsolid tumors.
PloS one , 9(7):e103245, 2014.[22] Gary R Mirams, Christopher J Arthurs, Miguel O Bernabeu, Rafel Bordas, Jonathan Cooper, Alberto Corrias,Yohan Davit, Sara-Jane Dunn, Alexander G Fletcher, Daniel G Harvey, et al. Chaste: an open source c++ libraryfor computational physiology and biology.
PLoS computational biology , 9(3):e1002970, 2013.[23] Atsushi Mochizuki, Naoyuki Wada, Hiroyuki Ide, and Yoh Iwasa. Cell-cell adhesion in limb formation, estimatedfrom photographs of cell sorting experiments based on a spatial stochastic model.
Developmental Dynamics ,211(3):204–214, 1998.[24] Payman Mosaffa, Nina Asadipour, Daniel Millán, Antonio Rodríguez-Ferran, and Jose J Muñoz. Cell-centredmodel for the simulation of curved cellular monolayers.
Computational Particle Mechanics , 2(4):359–370, Dec2015.[25] Payman Mosaffa, Antonio Rodríguez-Ferran, and José J. Muñoz. Hybrid cell-centred/vertex model for multicellu-lar systems with equilibrium-preserving remodelling.
International Journal for Numerical Methods in BiomedicalEngineering , 34(3):e2928, 2017.[26] Makiko Nonomura. Study on multicellular systems using a phase field model.
PLOS ONE , 7(4):1–9, 04 2012.[27] T. Reya, S. J. Morrison, M. F. Clarke, and I. L. Weissman. Stem cells, cancer, and cancer stem cells. nature ,414(6859):105, 2001.[28] S. Rudraraju, K. L. Mills, R. Kemkemer, and K. Garikipati. Multiphysics modeling of reactions, mass transportand mechanics of tumor growth. In
Computer Models in Biomechanics , pages 293–303. Springer, 2013.[29] E. Schierenberg.
Embryological variation during nematode development (January 02, 2006), WormBook, ed. TheC. elegans Research Community.
WormBook, 2006.[30] Maciej H Swat, Gilberto L Thomas, Julio M Belmonte, Abbas Shirinifard, Dimitrij Hmeljak, and James A Glazier.Multi-scale modeling of tissues using compucell3d. In
Methods in cell biology , volume 110, pages 325–366.Elsevier, 2012.[31] Srikanth Vedantam and BSV Patnaik. Efficient numerical algorithm for multiphase field simulations.
PhysicalReview E , 73(1):016703, 2006.[32] S.N. Verner and K. Garikipati. A computational study of the mechanisms growth-driven folding patterns on shells,with application to the developing brain.