## Abstract

The polarization of cells and tissues is fundamental for tissue morphogenesis during biological development and regeneration. A deeper understanding of biological polarity pattern formation can be gained from the consideration of pattern reorganization in response to an opposing instructive cue, which we here consider using the example of experimentally inducible body axis inversions in planarian flatworms. We define a dynamically diluted alignment model linking three processes: entrainment of cell polarity by a global signal, local cell–cell coupling aligning polarity among neighbours, and cell turnover replacing polarized cells by initially unpolarized cells. We show that a persistent global orienting signal determines the final mean polarity orientation in this stochastic model. Combining numerical and analytical approaches, we find that neighbour coupling retards polarity pattern reorganization, whereas cell turnover accelerates it. We derive a formula for an effective neighbour coupling strength integrating both effects and find that the time of polarity reorganization depends linearly on this effective parameter and no abrupt transitions are observed. This allows us to determine neighbour coupling strengths from experimental observations. Our model is related to a dynamic 8-Potts model with annealed site-dilution and makes testable predictions regarding the polarization of dynamic systems, such as the planarian epithelium.

## 1. Introduction

Epithelial tissues can be considered as two-dimensional sheets of densely packed cells. The properties of epithelia are highly regulated and instrumental for morphogenesis during biological development and regeneration [1]. A key property of epithelia is the establishment of different membrane domains on either side of the plane, termed apical and basal. This process polarizes epithelial cells perpendicular to the plane [2]. Molecules responsible for the establishment of apico-basal polarity include phosphoinositides, various GTPases, and the Crumbs and PAR complexes [3].

By asymmetrically localizing an independent set of molecules including Frizzled/Flamingo and Fat/Dachsous along an axis perpendicular to the apical–basal axis, cells of many epithelia superimpose a second polarity pattern *within* the plane, termed planar cell polarity (PCP) [4]. PCP controls fundamental processes during embryonic development and tissue regeneration in many species including actin filament orientation, convergence–extension, tissue reshaping, sensory organ formation, wing hair orientation, directional tissue growth and animal locomotion [4–6].

Mechanistically, PCP and the resulting planar tissue polarity integrate two general classes of inputs. (i) Global cues provided by the slope of tissue-scale gradients. These can consist of ligand concentration profiles [7,8], gene expression gradients [9] or mechanical shear stress [10]. (ii) Local cues are provided by cell–cell coupling. The alignment of cell polarization vectors among neighbouring cells propagates anisotropies from tissue boundaries or mutant clones and is mediated by the differential distribution of PCP and/or Fat/Dachsous components across cell/cell interfaces [4,11–13]. These mechanisms are universally found across many species and tissues. In most contexts, both inputs act synergistically to establish and maintain planar tissue polarity [14].

Theoretical studies of the collective phenomena of PCP confirmed that cell–cell neighbour coupling fosters a uniform polarity response of all cells to noisy and non-monotonous tissue-scale signals [6,15–23]. In particular, weak and even transient biases stemming from a polarized boundary or graded signal suffice to orient an entire epithelium when present from the *onset* of PCP dynamics in initially unpolarized cells [20,21]. Understanding of the underlying principle can be gained from statistical physics: the *q*-Potts model studies two-dimensional lattices that allow discrete polarization vectors in [24].

Indeed, the emergence of long-range order in PCP bears analogy to ferromagnetism. In the above models, each lattice node (or cell) carries a vectorial magnetic moment, analogous to a cell's PCP vector. That system's energy decreases by favouring configurations where individual magnetic moments align among neighbours and with the vector of an external magnetic field [25]. When fluctuations that tend to randomize individual magnetic moments are below a critical value, then long-range order and a system-wide net magnetization emerge spontaneously also in the absence of an external bias [26]. Analogously, PCP patterns in mutant tissue of fly wings and in model simulations, which abolished or decoupled the external bias, show spontaneously emerging order [20].

Contrary to the fixed arrangement of spins in ferromagnetic matter, however, biological tissues are composed of living cells that are born, age and become eliminated from the tissue. Tissues often exist much longer than their constituting individual cells and many tissues maintain their polarized state despite continuous cell turnover. In general, there are two scenarios how new cells can establish their PCP: they either inherit PCP from their polarized mother cells or polarize de novo. Since PCP signalling depends on the state of neighbouring cells, such cell turnover not only modulates the PCP state locally but also constitutes a topological perturbation of the cell arrangement, modulating the number of signalling neighbours. This may fundamentally alter the system dynamics beyond that of the classical Potts model.

In statistical physics, variants of the Potts model with zero magnetic moments for a subset of nodes have been studied, called diluted Potts models. The positions of the zero states are either fixed (quenched site-dilution) or in thermodynamic equilibrium with the other states (annealed site-dilution). The main interest there is in the equilibrium distribution and properties of the asymptotic state [24]. However, less is known about the duration and trajectories of transient dynamics approaching the asymptotic state, which are more relevant to biology, and on the impact of site-dilution on them. Since transient dynamics depend on the specifics of site-dilution, or cell turnover in biological terms, it is important to model the process of cell turnover directly according to a specific experimental system if one wants to exploit the physical models in a biological context.

Our work has two objectives. First, it bridges the gap between existing models with/without static site-dilution and the dynamics of polarity in tissues with cell turnover. We propose a dynamic model, similar to an 8-Potts model with annealed site dilution, termed *dynamically diluted alignment model* in the following, for the study of planar polarity formation and maintenance in biological tissues. Second, it elucidates the transient dynamics approaching the asymptotic state. We propose that new insight into polarity pattern formation can be gained from analysing the particular transient dynamics of polarity reorganization when an initially coherent polarity pattern is confronted with an opposing instructive signal. We therefore ask, how the contradiction between inputs is resolved and how the time requirement for conflict resolution depends on parameters, especially the cell birth and death rates.

The biological inspiration for our approach is the experimentally inducible inversion of global body plan polarity in the planarian *Schmidtea mediterranea* [27,28]. The regeneration of a second head instead of a tail (figure 1*a*–*c*) can be assumed to constitute a conflicting cue for the polarization pattern in pre-existing tissues. The multi-ciliated ventral epithelium is likely to be one such planarly polarized tissue [29,30]. Its cilia drive the gliding locomotion of planarians, implying consistent polarization of individual cilia and thus of the ventral epithelium as a whole [31,32]. Consequently, the movement of the animals may also inform us about polarization phenomena within the epithelium. In experimentally generated double-headed animals, each head moves in the opposite direction, thus giving rise to a continuous tug-of-war between the two heads with little net movement but stretching and thereby thinning the bulk tissue [27]. We interpret the balanced bi-polar movement as evidence for a re-polarization of the pre-existing epithelium (grey area in figure 1*e*,*f*) in response to an instructive cue provided by the new head. Moreover, we hypothesize that the cue constitutes a gradient of a signalling molecule, analogous to the Wnt gradient that patterns the planarian tail [33–35]. This interpretation is further supported by the observed symmetric inward motion of both body halves in double-tailed planaria (see electronic supplementary material, movie S3, of reference [27]; and electronic supplementary material, movie S2, of reference [31]).

Our work explores the question of how the polarization of a cell field responds to the superposition of a conflicting long-range signal. We assume that constituent cells undergo continuous turnover via the integration of new, initially unpolarized cells born outside the tissue (the equivalent of planarian neoblast progeny [36,37]) and the balanced extrusion of old, polarized cells in a dynamic steady state. A transiently naive cell presents no polarity information to its neighbours, and with a certain rate turns into a polarized state itself [38]. Such dynamic loss and re-establishment of polarization is here modelled as dynamic (annealed) site dilution of an 8-Potts model. Each cell in our model occupies exactly one node of a lattice. The introduced dynamics extends the classical Potts model in several ways described below. Note that other extensions of the Potts formalism exist that must not be confused, in particular, the ‘Cellular Potts Model’, in which an individual cell occupies multiple lattice nodes [39,40].

This article has the following structure. We first develop a dynamically diluted alignment model in the framework of interacting particle systems (IPS) which accounts for a global orienting signal, local coupling and cell turnover. This dynamically diluted alignment model allows us to study the effects of cell turnover on polarity patterns. Specifically, we ask whether and how polarity patterns with coherent initial polarization counter-directional to the global signal reorganize, to resolve the conflict between local and global directional cues, and what the time requirement is if they do so. We consider a polar alignment order parameter and identify the corresponding time of minimal order as the key characteristic of transient dynamics. Simulating the full model and by theoretical as well as numerical analysis of a mean-field approximation, we then show that cell–cell neighbour coupling in addition to its synergistic and noise-filtering role mentioned above *retards* the response of planar tissue polarity to dynamically changing global inputs, whereas cell turnover *accelerates* it. Finally, we establish a relation of the system parameters that determines the time requirement for polarity reorganization. We close with a discussion of these results.

## 2. Mathematical model of cell polarity and turnover

### 2.1. Model definition

We define an IPS [41–43] model for tissue polarity dynamics at the cellular level, which incorporates polarity alignment with respect to a global signal and to neighbours' polarity vectors as well as cell turnover. The model cells occupy the nodes of a finite two-dimensional lattice *S* that represents the epithelial tissue subjected to initially conflicting signals, as for instance the grey-shaded area in figure 1*e*,*f*. The cellular scale of granularity allows us to describe the essential interactions yet keeps the model analytically tractable, in analogy to the variants of the Potts model studying ferromagnetism.

Each cell is equipped with one of nine polarization states, see figure 2*a*, as follows. Thereby, the highly asymmetric concentration profile of PCP complexes along the cell membrane of a *polarized* cell, which determines the cell's polarity orientation, is abstracted as one unit vector per cell pointing towards the highest membrane accumulation of a selected PCP component. The directions of the unit vectors are discretized yielding the eight states
2.1Naive cells, in planaria resulting from progeny of stem cells immigrating into the epithelium and presenting no polarity information for an initial period of time, are represented by a ninth, *unpolarized*, state **e**_{0}≔(0, 0). Thus, a cell at node *z* ∈ *S* has a polarization state *η*_{z} which is an element of *W*≔{**e**_{i}, *i* = 0, …, 8}. The state space of the whole system is *W*^{S}. An element ** η** = (

*η*_{z})

_{z ∈S}∈

*W*

^{S}of the state space is called a configuration and describes the global state of the system. The model dynamics comprise two processes acting on individual cells: polarity alignment and cell turnover. Polarity alignment, in turn, is directed by two signals: local neighbours' polarities and a global orienting signal. As asymmetric protein complexes bridge adjacent cell membranes, polarization of each cell tends to align with neighbours' polarization vectors (see Introduction). In the model, a cell's neighbourhood is defined as those cells sharing a cell–cell interface with that cell. This is implemented by considering von Neumann neighbourhood in a square lattice which we complete with periodic boundaries. Assuming approximately equal lengths of cell–cell interfaces, we use the equally weighted average polarization vector 2.2of neighbours

*x*to node

*z*as the local director of polarity alignment. Here, #

*A*denotes the number of elements of any set

*A*. Our analysis below indicates that the specific choice of the lattice geometry, the neighbourhood template and spatial asymmetries with respect to the neighbourhood information do not affect our main results.

Additionally, a global vector of polarity alignment is considered, representing the slope orientation ** s** = (

*s*

_{x},

*s*

_{y}) of a tissue-scale gradient. The local director

*ν*_{z}(

**) and the global director**

*η***can be differently weighted by a neighbour coupling strength**

*s**ε*

_{n}≥ 0 and a coupling strength to the global signal

*ε*

_{s}≥ 0, respectively. Both weighted vectors are then summed vectorially to yield the reference orientation 2.3for node

*z*(figure 2

*b*). Considering a

*polarized*cell over time, a change of polarity to any new direction is modelled as more probable the more the new direction is aligned with the reference orientation

**w**, but it is assumed to be

*independent*of the current polarization direction of the considered cell. We deliberately consider abrupt changes in polarization direction, because protein complexes bridging pairs of membranes from neighbouring cells cannot shift across cell vertices, but disassemble at a given cell interface and assemble anew at another interface. The degree of alignment is measured by the standard scalar product 〈 · , · 〉 in . The rate

*c*

_{z}for changing polarization direction in node

*z*∈

*S*from polarized state

*η*_{z}∈

*W*\{

**e**

_{0}} to another polarized state

**e**

_{i},

*i*= 1, …, 8, while keeping all other nodes unchanged is then defined as 2.4where parameter

*γ*gives the overall pace of polarity reorientation, in analogy to previous alignment models [44–46]. Note that these rates minimize the energy functional 2.5where the first term measures alignment between neighbouring cells and the second term alignment with the global signal.

Additionally, cell turnover (ageing and replacement by naive cells) occurs independently of polarization direction. In the model, we let a polarized cell *z* ∈ *S* with *η*_{z} ∈ *W*\{**e**_{0}} change into the unpolarized state **e**_{0} with death rate *δ* ≥ 0,
2.6The establishment of any polarization direction from scratch in an unpolarized cell is modelled with de novo polarization rate *β* ≥ 0. By setting
2.7the polarization directions after de novo polarization are distributed as those in a re-alignment step, cf. equation (2.4) and note that *γ* cancels. Different from the diluted models in statistical physics, this type of cell turnover with a constant, state-independent death rate cannot be described in terms of an energy functional.

Taken together, the model equations (2.1)–(2.4), (2.6) and (2.7) define the transition rates of a continuous time Markov chain (** η**(

*t*))

_{t≥0}or more specifically an IPS [41–43], which we call the

*dynamically diluted alignment model*. See figure 2 for an illustration of the model dynamics. This dynamically diluted alignment model allows us to study the effects of cell turnover on polarity patterns. Specifically, we ask whether and how polarity patterns with coherent initial polarization counter-directional to the global signal

**reorganize, to resolve the conflict between local and global directional cues, and what the time requirement is if they do so.**

*s*### 2.2. Parametrization

The model parameters are *δ*, *β*, *γ*, *ε*_{n}, *ε*_{s}, ** s**, see §2.1 and figure 2, and the end time

*t*

_{max}of a simulation. The complete list of symbols is given in electronic supplementary material, table S1. We non-dimensionalize by choosing 1/

*γ*as the time unit of the model. Hence, the rates

*β*and

*δ*become dimensionless parameters, and we omit

*γ*whenever possible. The cell death rate

*δ*and de novo polarization rate

*β*may not be directly experimentally accessible. However, there are two quantities that are measurable, the average duration

*τ*of a cell's complete cycle through unpolarized and polarized states and the fraction of polarized cells, which can together be used to determine

*δ*and

*β*as follows. First, since a single polarized cell loses polarization at rate

*δ*, it remains polarized on average for 1/

*δ*time units. Analogously, unpolarized cells remain so for 1/

*β*time units on average. Hence, the average duration

*τ*of a complete cycle through unpolarized and polarized states per cell is related to

*β*and

*δ*by 2.8where the denominator

*γ*accounts for the unit model time.

Second, the fraction of polarized cells *p*_{eq} among all cells equals, after any initial transients have decayed, the ratio of the expected duration of the polarized state over the expected duration of the whole cycle of a cell,
2.9
Solving equations (2.8) and (2.9) for (*β*, *δ*) yields
2.10

The timescale *δ*^{−1} for cell death is typically of the order of days, while the timescale *β*^{−1} for de novo polarization is of the order of hours. The timescale of PCP alignment, given by the unit time *γ*^{−1} = 1, is also of the order of hours, since reorganization of transmembrane proteins from one to another side of the cell membrane requires endocytosis, vesicular transport and exocytosis. These processes depend on the polarized organization of the cytoskeleton, which has to be remodelled simultaneously. Hence, we explore the model behaviour for *β* ∈ {0.1, 1, 10} and vary *δ* = 0, 0.2*β*, …, 10*β* for each value of *β*. By equation (2.9), this corresponds to *p*_{eq} ranging from 1, where all cells are polarized (the case of dynamic Potts model without dilution), to a fraction of polarized cells approximately 0.099, effectively covering the fraction of polarized cells in all known polarized epithelia. The resulting values of *τ* range in [0.11, ∞]. For planarians, an experimental observation time *T* = 14 d [37] then corresponds to a cycle length from immigration into the tissue via polarization until the cell's death from 1.5 d upwards, which is realistic. Below, the impact of the dimensionless weights *ε*_{n}, *ε*_{s} that represent cellular sensitivities to neighbours' polarity and global signal, respectively, will be studied in detail. Since the reference orientation **w** ≔ *ε*_{n}** ν**(

**) +**

*η**ε*

_{s}

**is a weighted sum of the two directional cues, it suffices to keep**

*s**ε*

_{s}= 1 fixed and vary

*ε*

_{n}∈ [0, 5]. We also checked that no abrupt changes of model behaviour occur for

*ε*

_{n}= 1 and variation of

*ε*

_{s}∈ [0, 5] (electronic supplementary material, figure S7). The model is symmetric w.r.t. the discretized directions

**e**

_{1}, …,

**e**

_{8}, so the

**e**

_{8}direction can be chosen to coincide with the direction of the global signal, i.e.

**s**= (1, 0) =

**e**

_{8}.

The initial configuration represents a homeostatic tissue with fraction *p*_{eq} of polarized cells, where the dominant polarization direction is opposite to the global signal ** s**. Therefore, we assign the initial state

**e**

_{0}with probability 1 −

*p*

_{eq}=

*δ*/(

*β*+

*δ*) to each cell independently and set the initial states of polarized cells as if each cell had experienced prior to simulation start coherent global and local directors

**s**

_{init}=

*ν*_{init}= (−1, 0) =

**e**

_{4}= −

**s**, opposite to

**. The latter means that each polarized cell is assigned**

*s***e**

_{k}with probability 2.11independently, where is a normalization constant. This way, cells are initially coherently polarized with main polarization direction

**e**

_{4}= (−1, 0), which indeed conflicts with

**s**= (1, 0) (figure 3

*A*).

### 2.3. Observables

To decide whether the polarity pattern adapts to the global signal, and if so to quantify the time requirement for reorientation, we track the mean polarization
2.12of an evolving configuration ** η**(

*t*) during the simulation. The behaviour of this polar alignment order parameter is best described in terms of modulus ∥

**p**∥ and angle to the positive

*x*-axis ang(

**p**). Owing to the stochasticity of the model, both modulus and angle fluctuate, but these variations decrease with increasing lattice size (electronic supplementary material, figure S1). The modulus ∥

**p**∥ ∈ [0, 1] characterizes the degree of alignment among all cells' polarization directions. If ∥

**p**∥ = 0 then no globally dominant polarization direction exists, whereas in case of ∥

**p**∥ = 1 all cells are polarized into the same direction. However, ∥

**p**∥ = 1 requires that all cells are polarized. Owing to cell turnover, the actual fraction of polarized cells 2.13fluctuates around

*p*

_{eq}and obeys ∥

**p**∥ ≤

*p*

_{p}≤ 1. Hence, ∥

**p**∥/

*p*

_{p}≤ 1 characterizes the degree of alignment among the

*polarized*cells where equality holds if and only if all polarized cells share one direction. The angle ang(

**p**) indicates the predominant polarization direction, and exhibits switching behaviour if the polarity pattern reorients. Hence, a high value of ∥

**p**∥ together with an angle ang(

**p**) approximately oriented parallel to the global vectorial signal

**will inform us that polarity reorientation has occurred. Then local and global signals are coherent and the configurations are in a stochastic dynamic equilibrium.**

*s*To quantify the time requirement of reorientation, we monitor the time needed to reach the dynamic equilibrium. We observe that, as a prerequisite for polarity pattern reorientation, the degree of alignment measured by ∥**p**∥ diminishes until ∥**p**∥ attains a distinct minimum, and that ang(**p**) undergoes the fastest change around that time, cf. figure 3*B*. In addition, the first phase in which the initially coherent polarity pattern resolves into a minimally ordered transient state, is of particular interest to study the influence of conflicting signals, while the subsequent evolution of a disordered system towards an ordered state due to the global signal has been studied before [20]. Therefore, we use the time of minimal order
2.14as the characteristic, statistically robust time for conflict resolution in tissue polarity reorganization, cf. figure 3*B*,*C*.

### 2.4. Model analysis

#### 2.4.1. Simulation

To sample from the trajectories of our stochastic model, we employ the exact stochastic simulation algorithm by Gillespie [47] in an efficient implementation for the IPS [43]. Three of six parameters are held fixed as described in §2.2: *γ* = 1 by non-dimensionalization, **s** = (1, 0), *ε*_{s} = 1. The initial configuration at *t* = 0 is specified by equation (2.11). The other three parameters are varied as *ε*_{n} = 0, 0.5, 1.0, …, 5.0, *β* = 0.1, 1, 10 and *δ* = 0, 0.2*β*, 0.4*β*, …, 10.0*β*. Simulations are carried out on a 100 × 100 lattice with periodic boundary conditions until the simulated time exceeds . See figure 3*A*; electronic supplementary material, movie for an example simulation, and electronic supplementary material, figure S1 for a justification that the lattice size is sufficient.

#### 2.4.2. Mean-field analysis

Using mean-field approximation, we derive an ODE which approximates the temporal evolution of **p**(*t*) = **p**(** η**(

*t*)). See electronic supplementary material, S2 for more details of the following derivations. Denote the fractions

**(**

*a**t*) = (

*a*

_{0}(

*t*), …,

*a*

_{8}(

*t*)) of nodes in states

**e**

_{0}, …,

**e**

_{8}at time

*t*in the IPS model by 2.15Then the mean polarization vector can be expressed as 2.16and the fraction of polarized cells is .

The mean-field approximation (MFA) simplifies the rates c_{z}(** η**,

**e**

_{i}) defined in equations (2.4) and (2.7) to (cf. electronic supplementary material, S2) 2.17

In the limit for increasing lattice size, the *a*_{i}'s become continuous quantities and their dynamic behaviour can be described by an ODE system [48–50]
2.18Here and denote the counterparts of ** a** and

*a*

_{i}under MFA. We call equation (2.18) the

*mean-field model*and note that the overall error of approximation introduced in equation (2.17) increases with

*ε*

_{n},

*ε*

_{s}and the fraction of polarized cells

*p*

_{p}. Note that other, even irregular, lattice geometries, different neighbourhood templates and spatially asymmetrically weighted neighbour polarity information yield the same MFA equation (2.18) as long as neighbour polarity information is weighted independently of the considered cell's polarity state.

The fraction of unpolarized cells tends to the unique, globally attracting equilibrium , which is in perfect agreement with the dynamic equilibrium of death and de novo polarization in the original IPS (equation (2.9)). This equilibrium simplifies the ODE system (2.18) to
2.19which can be summed to an ODE for , see electronic supplementary material, eq. (S13). Solutions of (2.19) preserve the symmetry of the initial condition with respect to the *x*-axis, which was imposed by setting **s** = (1, 0), *a*_{1}(0) = *a*_{7}(0), *a*_{2}(0) = *a*_{6}(0) and *a*_{3}(0) = *a*_{5}(0) (cf. equation (2.11)). For such a symmetric initial condition, it holds for all times *t* > 0. Numerical solutions of equation (2.19) are shown in §3.2.

#### 2.4.3. Linearization of mean-field model

One can approximate the nonlinear mean-field model (2.19) further, see electronic supplementary material, S3, to obtain an analytically tractable ODE
2.20where an overline over symbols indicates this second approximation. Again, an ODE for (electronic supplementary material, eq. (S15)) follows. The analytical solution (electronic supplementary material, eq. (S16)) implies that polarity reverses in the linearized mean-field model only for *ε*_{n}(*β*/(*β* + *δ*)) < 2. In this case, the time of minimal order is uniquely determined as (see electronic supplementary material, eq. (S17))
2.21

## 3. Results

### 3.1. Simulation of the dynamically diluted alignment model

#### 3.1.1. Time course in simulations

For all parameter combinations studied, the polarity pattern resolved the initial conflict within the simulated time window by reversing the main polarization direction and adapting to the global signal. Even for the largest chosen neighbour coupling strengths *ε*_{n} = 5.0, we never observed the frustrated initial condition to persist. The time courses of ∥**p**∥ and ang(**p**) exhibit several common characteristics independent of the specific parameter sets, described as follows.

Figure 3*A* and electronic supplementary material, movie show an example trajectory of the dynamically diluted alignment model together with the time course of the observable **p** derived from equation (2.12) in figure 3*B*. According to the initialization (cf. equation (2.11)), the system starts globally ordered where the dominant orientation **e**_{4} = (−1, 0) is opposed to the global signal **s** = (1, 0). Single cells or small cell patches deviate from their initial polarization direction (figure 3*A**a*) and others follow until there is hardly any predominant polarization direction around the time of minimal order (figure 3*Ab*). This happens independent of lattice size (see electronic supplementary material, figure S1). In succession, the polarity pattern approaches a state of well-aligned polarization directions, dominated by **e**_{8} in coherence with the global signal (figure 3*A**c*). These configurations form a dynamic equilibrium that persists (cf. figure 3*Ac*, *Ad*). The transition from alignment among cells conflicting with the global signal to alignment with the global signal happens via a disordered state when each polarization direction is approximately equally abundant (see electronic supplementary material, figure S3A).

These dynamics are well recapitulated in the time course of the polarity alignment order parameter **p**, see figure 3*B* and §2.3. Its modulus ∥**p**∥ starts at a high value not above *β*/(*β* + *δ*) and decreases to a distinct minimum at the time of minimal order *T*_{mo} as more and more cells leave their initial directions. The modulus ∥**p**∥ increases again immediately after attaining the minimum as a growing majority of cells adopts states **e**_{1}, **e**_{7} and finally **e**_{8} when ∥**p**∥ reaches a plateau. The angle ang(**p**) remains almost constant ≈−*π* up to shortly before *T*_{mo}, increases steeply around *T*_{mo} to a level of ≈0 around which it then fluctuates. Note that the switching in ang(**p**) is also possible as a decrease from ≈*π* to ≈0 for symmetry reasons, cf. electronic supplementary material, figure S3B. Hence, the time of minimal order *T*_{mo} is distinguished not only by the least degree of alignment ∥**p**∥ in the polarity pattern, but also by the fast, switch-like change of the angle ang(**p**). This indicates that pattern reorientation takes place during the short phase of low polarity alignment around *T*_{mo}, and that the previous decrease in order is a precondition. Therefore, *T*_{mo} is appropriate to study the time requirement of polarity pattern reorientation. Moreover, the time of minimal order *T*_{mo} is statistically robust. It differs only slightly between sampled trajectories of **p**(*t*) from a fixed parameter set, see figure 3*B*, lower panel, and electronic supplementary material, figure S3B.

#### 3.1.2. Parameter dependence of time of minimal order *T*_{mo}

As described in §2.2, it suffices to vary replacement rate *δ* and neighbour coupling strength *ε*_{n} to explore the model behaviour while keeping de novo polarization rate *β* and coupling strength to global signal *ε*_{s} fixed. Figure 3*C*,*D* shows the parameter dependence of the time of minimal order for empirical standard deviations see electronic supplementary material, figure S2B. First, we find that the time of minimal order grows with *ε*_{n}. This is plausible as neighbour coupling hinders cells from breaking free from their initially frustrated polarization direction. Second, we observe a decline in the time of minimal order with growing replacement rate *δ*. This is plausible as well, since higher values of *δ* for fixed *β* imply an increased fraction of unpolarized cells, cf. equation (2.9), and the average neighbour polarization shrinks in absolute value. Therefore, the reference orientation ** w**, which is the vectorial sum of global and local directors (see equation (2.3) and figure 2

*b*) shifts towards the global signal, which reduces the decelerating effect of neighbour coupling and accelerates reorientation. For

*δ*≫

*β*,

*T*

_{mo}depends less strongly on

*ε*

_{n}as evident from the smaller slope of the corresponding red data points compared to all others in figure 3

*D*. Actually,

*T*

_{mo}might approach some plateau for

*δ*→ ∞, while

*β*= const since polarized cells become virtually isolated, cf. equation (2.9). The isotemporales in figure 3

*D*connecting parameter sets with equal

*T*

_{mo}highlight the decelerating effect of increased neighbour coupling and the acceleration by cell turnover. The results for

*β*= 0.1 and

*β*= 10 agree qualitatively with the case

*β*= 1, for details see electronic supplementary material, figure S2A–C. We checked that the model behaviour is qualitatively unchanged for a different lattice size, see electronic supplementary material, figure S1 for results on a 20 × 20 lattice.

### 3.2. Numerical solution of the mean-field model

To unravel the quantitative dependence between the accelerating impact of *δ*, the decelerating effect of *ε*_{n} and the role of the de novo polarization rate *β* for the time of minimal order, we analyse the mean-field approximation. We numerically solve the mean-field model (2.19) in Morpheus [51,52] using the Runge–Kutta discretization scheme with time step 10^{−6} and the initial configuration specified by equation (2.11). The temporal evolution of the mean-field fractions of cells in state **e**_{i} is sketched exemplarily in figure 4*A*, and shown together with the observable derived from electronic supplementary material, eq. (S12) in panels B and C, respectively. For all parameter combinations studied in the IPS model, the polarity patterns described by the mean-field model, equation (2.19), reorganize within the simulated time window as well, reversing the main polarization direction and adapting to the global signal, except for the border cases of *δ* = 0, *ε*_{n} ∈ {4.5, 5}. The time courses of and for the turning cases exhibit several common characteristics independent of the specific parameter sets and with those for the original model, described in electronic supplementary material, S6 in detail. In particular, there is a distinct time of minimal order which characterizes the time requirement of polarity pattern reorientation in the mean-field model. For the non-turning cases, the neighbour coupling strength *ε*_{n} is too high such that the mean-field approximation is no longer usable as discussed in §2.4.2 and exemplarily shown in electronic supplementary material, figure S4.

#### 3.2.1. Parameter dependence of time of minimal order

The measured time of minimal order from our numerical simulations of the nonlinear mean-field model (2.19) is reported in figure 4*D* as isotemporales alongside data for *T*_{mo} from the IPS simulations. The two datasets coincide very well qualitatively and quantitatively. We conclude that the mean-field approximation is a suitable tool to study the polarity reorientation in the dynamically diluted alignment model throughout a wide parameter range, and that the qualitative interpretation of parameter dependencies extends from the IPS (cf. §3.1) to the mean-field approximation. For border cases *ε*_{n} > 4.5 and *δ*↘*δ*_{crit} for critical value *δ*_{crit} = *δ*_{crit}(*ε*_{n}) > 0, the time of minimal order grows faster than exponentially as evidence of the transition to non-turning polarity patterns (figure 4*E*).

For *ε*_{n} = 0, cells evolve independently and the time of minimal order can be determined analytically from an ODE as (electronic supplementary material, S5)
3.1Moreover, coincides with time of minimal order for the averaged order parameter 〈**p**〉, see electronic supplementary material, S5 and figure S5.

### 3.3. Linearized mean-field ODE

Because of the good quantitative agreement in the time of minimal order between the dynamically diluted alignment model (*T*_{mo}) and its mean-field approximation (), we investigate the latter model further in the linearized mean-field form of equation (2.20), an ODE system that allows analytical treatment. The stability behaviour of the solution for the order parameter ODE depends on *A* ≔ (*δ* + 8)(−1 + (*ε*_{n}/2)(*β*/(*β* + *δ*)) (see electronic supplementary material, eqs. (S15) and (S16)). For *A* < 0 or equivalently *ε*_{n}(*β*/(*β* + *δ*)) < 2, the solution converges to a unique stable equilibrium. For *A* > 0 or equivalently *ε*_{n}(*β*/(*β* + *δ*)) > 2, there is a unique but unstable equilibrium and the solution diverges. The order parameter in steady state, whether stable or not,
3.2is aligned with the global signal ** s** in the stable, therefore convergent, case, but counter-directional in the unstable, therefore divergent, case. This together with the fact that for 2 −

*ε*

_{s}∥

**s**∥

*β*/(

*β*+

*δ*) <

*ε*

_{n}(

*β*/(

*β*+

*δ*)) < 2 +

*ε*

_{s}∥

**s**∥

*β*/(

*β*+

*δ*) in contrast with ∥

**p**(

*t*)∥ ≤ 1 in the IPS indicates that the linearized ODE system (2.20) and derived equations ((2.21), electronic supplementary material, eqs. (S15) and (S16)) are an appropriate approximation of the dynamically diluted alignment model only in the convergent case

*ε*

_{n}(

*β*/(

*β*+

*δ*)) < 2. The divergent case is a spurious solution introduced by the errors of MFA and linearization that both grow with

*ε*

_{n}(

*β*/(

*β*+

*δ*)). Figure 5

*b*visualizes the parameter dependence of for

*ε*

_{s}= 1,

*s*

_{x}= 1 (cf. equation (2.21)) and the initial condition given by equation (2.11). The results coincide with

*T*

_{mo}in the IPS description qualitatively and even quantitatively, except for divergence of for

*ε*

_{n}(

*β*/(

*β*+

*δ*))↗2 as evidence for the transition to the divergent case (compare figure 3

*D*to figure 5

*b*). Note that in equations (2.21), (3.2) and electronic supplementary material, eq. (S16), the neighbour coupling strength

*ε*

_{n}is rescaled with a factor

*β*/(

*β*+

*δ*) that equals the fraction of polarized cells

*p*

_{eq}, cf. equation (2.9), thereby modulating neighbour influences with cell turnover. This suggests revisiting the empirical results from the IPS simulations and nonlinear MFA as a function of this parameter combination, which we call the

*effective*neighbour coupling strength, 3.3

### 3.4. Analysis of simulation results through parameter rescaling

The identification of the effective neighbour coupling strength *ε*^{eff}_{n} in the analysis of the linearized mean-field model above suggests rescaling the simulation results according to this parameter combination. Figure 6 shows the data from our stochastic IPS model (from figure 3*C*,*D*) and from numerical solutions of the nonlinear mean-field model (from figure 4*D*) versus the rescaled parameter *ε*^{eff}_{n} = (*β*/(*β* + *δ*))*ε*_{n}. We find that all data of figures 3*C*,*D* and 4*D* collapse onto approximately linear relationships. For the IPS model with *ε*_{s} = 1,
3.4where values in brackets denote (mean ± s.d.) an orthogonal distance regression. The inverse relationship reads
3.5with experimentally accessible variables *T*_{mo} and *β*/(*β* + *δ*) = *p*_{eq} on the right-hand side. The small remaining scatter of the data points around the black line in figure 6*b* is largely given by a *δ*-dependent offset as the ordered colours of the data points indicate. This *δ*-dependent offset is analytically known for *ε*_{n} = 0 by equation (3.1). Subtracting this offset reduces the scatter further, as shown in electronic supplementary material, figure S6.

In the mean-field model, the critical value *δ*_{crit}(*ε*_{n}) for the transition to non-turning dynamics is transformed into a critical effective neighbour coupling strength *ε*^{eff}_{n,crit} ≈ 4.60 valid throughout the considered parameter space. According to our simulation data and model analysis, the impact of cell replacement on tissue polarity reorganization can be well described as weakening of the neighbours' influence on each single cell's polarization by the fraction of unpolarized cells. Altogether, neighbour coupling retards polarity pattern reorganization, whereas cell turnover accelerates it. The time of minimal order and effective neighbour coupling strength are related by an approximately linear function.

### 3.5. Determination of parameter values

Our model analysis allows estimation of the actual parameter value of *ε*_{n} from experimental observations. Equation (3.5) only requires *p*_{eq} and *T*_{mo}, that can be measured, and indirectly the time unit 1/*γ* which relates to the frequency of polarity changes, see equation (2.4). To determine the latter, we see different options. First, it might be obtained using time-lapse imaging of sub-cellular polarity markers. Second, the analysis of calibrated models of the molecular processes underlying tissue polarity may provide estimates for *γ*. Third and most realistically, the dimensionless equation (3.4) can be written as
3.6where denotes the value *with* physical time unit. Measuring in experiments for constant neighbour coupling strength *ε*_{n} and at least two different values of *p*_{eq} = *β*/(*β* + *δ*), the instances of equation (3.6) form a linear equation system that can be solved for (*γ*, *ε*_{n}). Variations of cell death rate *δ* and/or de novo polarization rate *β* are feasible using drugs, RNAi and other techniques that interfere with pathways of cell proliferation and apoptosis. The absolute value of the neighbour coupling strength *ε*_{n}, obtained through any of the above ways, can then be interpreted in relation to the coupling strength to the global cue, here chosen as *ε*_{s} = 1.

Beyond such interpretation of absolute values, the comparison of inferred values for *ε*_{n} among multiple perturbed conditions within a screening approach allows us to disentangle the mechanistic effects of perturbations on neighbour coupling versus cell turnover versus intracellular polarity dynamics. Our theoretical results have resolved how all three contributions jointly determine the timescale of polarization reorientation. As our proposed protocol only requires measurement of quantities that are accessible from still images, this theory alleviates the need for life imaging of the same specimen.

## 4. Discussion

We have here posed the question how local and global instructing signals are integrated in a PCP system. In particular, conflict resolution, as observed in double-headed planaria, is proposed as a valuable source of information about which signals dominate the temporal evolution which are not accessible from studying random initial conditions. To study this problem theoretically, we propose a cell-based IPS model accounting for local cell–cell coupling of polarity, for sensitivity to global ligand gradients and for the impact of cell turnover as present in planaria, termed dynamically diluted alignment model.

Analysing the model numerically and analytically, we find that the global signal dictates the final tissue polarization orientation independent of the strength of neighbour coupling between cells and the amount of cell turnover. The temporal evolution from an initially polarized tissue conflicting with the global signal to the final polarized state in accordance with the global signal occurs via a disordered state in which each polarization direction is approximately equally abundant. We introduce this time of minimal order *T*_{mo} as an observable to measure the timescale of conflict resolution and study its dependency on cell turnover rate and neighbour coupling strength. It turns out that neighbour coupling retards polarity pattern reorganization, whereas cell turnover accelerates it, and that dependencies are gradual without abrupt transitions. We employ mean-field analysis of the IPS model to derive an ODE for the temporal evolution of the average polarization and an equation for the time of minimal order which well approximate the IPS data. As a result, we identify an effective neighbour coupling strength which integrates the parameters of cell turnover and neighbour coupling, and demonstrate that the time of minimal order in the dynamically diluted alignment model depends linearly on the effective neighbour coupling strength.

The dynamically diluted alignment model developed here extends former approaches where the effects of neighbour coupling on PCP polarity establishment were studied by means of a cell-grained tissue polarity model [20]. Our model accounts for cell turnover which is present in most biological tissues. From a theoretical point of view, dynamical site dilution is a generalization of static site dilution as studied in statistical physics in the context of ferromagnetism. Annealed site dilution of Potts models has been formulated previously, but studies so far have focussed on steady-state properties rather than temporal dynamics. Moreover, PCP is a widespread phenomenon in live matter and provides an experimental realization of an annealed site-dilution Potts model.

In the model, we neglect the molecular details of specific pathways underlying cell polarity reorganization, because we are interested in polarity conflict resolution between cellular and tissue scales and aim for a model that is still analytically tractable. In planaria, the molecular details of the PCP pathway and its upstream global signals are slowly unravelling, but their role for planar tissue polarity needs to be studied further [30,33–35]. For our model, all contributing tissue scale signals were subsumed into the abstract global director ** s**. In the same spirit, the orienting cues of the local cell–cell coupling are summarized as an average neighbour polarization vector. Independent of the specific molecular pathways, cell–cell bridging complexes do not move across three-cell junctions, but are degraded at one-cell interface and assembled anew at another interface. Therefore, we deliberately do not consider gradual changes in polarization direction, but model the change of cell polarity direction as independent of the current polarization direction of the considered cell. The resulting model is equally well applicable to study Frizzled/Flamingo- or Fat/Dachsous-based patterning or other mechanisms of tissue polarity at cellular resolution in biological tissues with and without cell turnover.

We have focussed on planaria as a biological reference system. Another classical system to study tissue polarity is *Drosophila* wing development. While the polarity alignment mechanisms of *Drosophila* and planaria can be considered similar in our model abstraction, the process of cell turnover in *Drosophila* differs as polarized epithelial cells divide, polarity can be inherited and naive unpolarized cells do not exist. Also, in *Drosophila* wing development, mechanical processes including hinge contraction, cell elongation and cell rearrangements cause local tissue shear that contributes to polarity reorientation. Interestingly, this remodelling of the *Drosophila* wing is also characterized by a transient phase of reduced global polarity order [10] as observed in our model for planaria.

Our model can be easily extended to include other variants of cell turnover, cell rearrangements or oriented cell division as long as the total number of cells is maintained. We speculate that the inheritance of polarity cancels the diluting effect of cell death and therefore decelerates tissue polarity conflict resolution in *Drosophila* wing development. On the other hand, to explicitly incorporate cell elongation, cell packing changes or tissue growth, other modelling frameworks appear more suitable, for example, the Cellular Potts Model, although the refined spatial resolution hinders analytical tractability [39,40].

The IPS model follows the inherent discretization of tissue into cells and describes the dynamics in continuous time as a stochastic process, which reflects noise and randomness in molecular interactions underlying the polarity patterns. It was described here for a square lattice with von Neumann neighbourhood but the model definition is valid as well for other lattices, like hexagonal or even image-derived lattices, and for more general neighbourhood templates. We expect similar results for other lattice geometries or inhomogeneously weighted neighbour polarity information since the mean-field approximation, which is independent of the actual spatial arrangement, closely agrees with the original model. A model extension that includes cell migration and cell division is straightforward (e.g. [53–55]), but their effects on tissue polarity patterns have not been in the focus of our study. The number of polarized states is an implicit model parameter.

Our choice of eight polarization vectors is the smallest number that conforms with the four-fold rotational symmetry of the lattice and allows, besides perpendicular and opposing polarity directions, also partial alignment with a directing polarity signal. Also, it is known that the planar *q*-Potts model without external field exhibits a first-order phase transition only for *q* > 4, whereas the case *q* = 4 can be reduced to *q* = 2, i.e. the Ising model [24]. Note that by considering an ordered initial condition and a conflicting external field in this work, the emergence of spontaneous order and corresponding phase transitions cannot be studied. From the analytical point of view, the choice of more than four-cell polarization states also contributes to the good agreement between the IPS data and mean-field approximation. In the limit *q* → ∞, the *q*-Potts model yields the XY-model with continuous angle space and the Beresinkii–Kosterlitz–Thouless transition [56]. However, we do not consider this limit an appropriate model for polarity patterns in epithelia since there are a number of three-cell junctions around each cell where transmembrane protein complexes cannot form. Hence, a continuous polarity angle is only possible within finite angle ranges that are separated by narrow excluded angle ranges. The transient dynamics of such a piecewise continuous model may, on short timescales, resemble that of the XY-model but on long timescales that of the Potts model. Concerning the results of our analysis on the asymptotic steady state, we remark that the observed universal dominance of the global signal on the long-term tissue polarization direction is in agreement with the known behaviour of Potts models [24]. However, the detailed temporal dynamics and the question of which signal dominates in a conflicting situation in the context of planar tissue polarity have not been studied before. This indicates planar tissue polarity as a versatile experimental framework that represents Potts models with dynamic site-dilution.

Given the knowledge that the global signal determines the long-term behaviour of planar tissue polarity, it is plausible that neighbour coupling retards polarity pattern reorganization since it enforces the maintenance of the initial, locally coherent but globally conflicting direction. For the same reason, cell replacement has an accelerating effect since it facilitates resolution of contradictory signals by reducing the number of polarized neighbours. By quantifying that *T*_{mo} depends linearly on the effective neighbour coupling strength *ε*^{eff}_{n} = (*β*/(*β* + *δ*))*ε*_{n} (figure 6), our theory enables the estimation of the neighbour coupling strength from equation (3.5) or (3.6) by measuring *T*_{mo} and the fraction of polarized cells *p*_{eq} = *β*/(*β* + *δ*). This allows us to infer the relative importance of global versus local directing signals, and to predict the effects of altering cell turnover on the timescale of tissue polarity reorganization. In particular, in the finite time frame of an experiment, it is possible that modified cell turnover prolongs the time required for reorganization beyond the observation window.

Mean-field approximation of each cell's local director field by the average of all fields allows further analytical results, notably the derivation of *ε*^{eff}_{n}, but neglects spatial correlations that become important for high neighbour coupling strength *ε*_{n}. This discrepancy induces a phase transition in the mean-field results that is not observed in our dynamically diluted alignment model itself. However, the artificially introduced phase transition occurs not until *ε*^{eff}_{n} ≈ 4.6, making mean-field approximation a valuable tool that provides qualitative and quantitative match for a wide range of parameters. The linearized mean-field model as a further simplification (equation (2.20)) provides an analytical expression for the polarity state and the time of minimal order (equation (2.21)), but shifts the spurious phase transition down to *ε*^{eff}_{n} = 2. The mean-field analysis can be improved when the independence assumption is replaced by a kind of pair approximation [57]. This can reduce approximation errors and extend the range of applicability towards higher neighbour coupling strength *ε*_{n}, for which the mean-field results under the independence assumption so far deviate from the IPS simulation results.

The proposed model and its analysis performed here are ready to be applied to quantitative data. Direct or indirect measurements of tissue polarity with cellular resolution in *S. mediterranea* would allow us to determine the model parameters, since quantified local alignment of polarity directions as a time- and space-dependent order parameter or decay of correlations directly link experimental data to observables of the model. The determined parameter set then implies further model predictions that could be tested experimentally. The model is also applicable to discriminate between several hypothetical ligands that might provide the global signal. If they have different spatio-temporal concentration profiles, these can be tested *in silico* to reproduce the observed tissue polarity reorientation patterns.

## Data accessibility

Electronic supplementary material accompanies this study and provides electronic supplementary material text §§1–6, a table summarising all variable names and notations, electronic supplementary material, figures S1–S7, a movie of the typical dynamics of the IPS model, and data files with all measured times of minimal order for the IPS and MF models.

## Authors' contributions

K.B.H., A.V.-B. and L.B. designed the study and developed the model. J.C.R. designed the experiments. K.B.H. performed model simulations and analysis and generated figures. K.B.H., A.V.-B., J.C.R. and L.B. interpreted the data, wrote and revised the manuscript. All authors gave their final approval for publication.

## Competing interests

We have no competing interests.

## Funding

This work was supported by the German Federal Ministry of Education and Research (BMBF) under funding codes 0316169 and 031L0033. A.V.-B. acknowledges the support by Sächsisches Staatsministerium für Wissenschaft und Kunst (SMWK) in the framework of INTERDIS-2. The simulations were performed on HPC resources granted by the ZIH at TU Dresden.

## Acknowledgements

The authors are grateful to Hanh Thi-Kim Vu for the help with imaging planaria and to Walter de Back, Andreas Deutsch, Benjamin Friedrich, Michael Kücken, Fabian Rost, Jörn Starruß and Carsten Timm for fruitful discussions. We thank anonymous referees for their constructive reviews.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3887818.

- Received June 23, 2017.
- Accepted September 11, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.