{'fraction_non_alphanumeric': 0.04802905959908516, 'fraction_numerical': 0.01690658773936051, 'mean_word_length': 4.2040840140023334, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Generation and motion of interfaces in a mass-conserving reaction-diffusion system

Reaction-diffusion models with nonlocal constraints naturally arise as limiting cases of coupled bulk-surface models of intracellular signalling. In this paper, a minimal, mass-conserving model of cell-polarization on a curved membrane is analyzed in the limit of slow surface diffusion. Using the tools of formal asymptotics and calculus of variations, we study the characteristic wave-pinning behavior of this system on three dynamical timescales. On the short timescale, generation of an interface separating high- and low-concentration domains is established under suitable conditions. Intermediate timescale dynamics is shown to lead to a uniform growth or shrinking of these domains to sizes which are fixed by global parameters. Finally, the long time dynamics reduces to area-preserving geodesic curvature flow that may lead to multi-interface steady state solutions. These results provide a foundation for studying cell polarization and related phenomena in biologically relevant geometries.

Reaction-diffusion models with nonlocal constraints naturally arise as limiting cases of coupled bulk-surface models of intracellular signalling. In this paper, a minimal, mass-conserving model of cell-polarization on a curved membrane is analyzed in the limit of slow surface diffusion. Using the tools of formal asymptotics and calculus of variations, we study the characteristic wave-pinning behavior of this system on three dynamical timescales. On the short timescale, generation of an interface separating high-and low-concentration domains is established under suitable conditions. Intermediate timescale dynamics is shown to lead to a uniform growth or shrinking of these domains to sizes which are fixed by global parameters. Finally, the long time dynamics reduces to area-preserving geodesic curvature flow that may lead to multi-interface steady state solutions. These results provide a foundation for studying cell polarization and related phenomena in biologically relevant geometries.

Introduction.
Reaction-diffusion models are essential tools for understanding the spatial selforganization of chemical patterns inside the cell. Over the past decade, research interest has coalesced around two key features which set intracellular dynamics apart from other classes of models. First, in contrast to traditional models which exist on a single domain, these models are often bulk-surface models, in that they feature distinct diffusion processes within the 3D cytosolic volume of the cell and on the 2D cell membrane coupled by a nonlinear boundary condition [46]. Second, many recent studies have emphasized systems subject to mass-conservation [24,42]. These properties are particularly prevalent in models of cell polarization, a crucial process by which the spatial distribution of proteins within a cell becomes highly localized to a region of the membrane as a result of spontaneous symmetry breaking or an external guiding cue [16]. Polarization is essential for a great variety of biological phenomena, including guiding developmental outcomes, establishing axes for cell division and guiding locomotion in motile cells, and so has motivated significant interest from the mathematical biology community [16].
The attention of applied mathematicians is increasingly turning to understanding how various forms of spatial heterogeneity influence dynamics and steady state behavior. Research that has emphasized the effect of spatially varying kinetic parameters are the most obvious examples of this trend [21,41]. Over the last few years, the question of what role cell geometry plays in guiding polarization has drawn considerable interest. While initial polarization models were reduced to 1D systems, computational advances have enabled numerical studies on 2D and fully 3D domains [6,15]. Early steps in addressing this problem have produced simulations that are highly suggestive that localization is closely tied to surface curvature via minimization of interfacial length, but this principle has yet to be demonstrated by formal analysis [14,20]. Very recently, a novel numerical framework was introduced allowing for efficient simulation of cell polarization models on surfaces of revolution, which in particular motivates this particular study [37].
In this paper, we perform a formal asymptotic analysis of the surface-bound version of the wavepinning model first proposed in Ref. [39], which features both mass-conservation and bulk-surface coupling. This model has attracted considerable interest as a minimal theory for polarization, and has been successfully applied to model systems such as the Rho-GTP pathway [56] or Ezrin polarization in embryonic mouse cells [58]. The emergence of spontaneous symmetry breaking in this system has been intensely studied, and the dependence of the initial instability of the spatially uniform solution on various parameters is well-characterized [33,43,48,56]. While previous studies considered the asymptotics of wave-pinning as surface diffusion becomes small, they have generally been concerned with 1D domains or lower orders of perturbation than we discover are needed to capture the full effect of 3D-embedded domain geometry [14,40]. Other works which have considered fully 3D domains have limited analysis to linear stability studies [20,47]. Our research here uses the timescale separation techniques to probe the generation and propagation of interfaces in two-species reaction-diffusion systems to the specific case of the wave-pinning model (for related rigorous studies, see [12,27,51]). We note that for single bistable reaction-diffusion equations in the Euclidean setting the generation and propagation of interfaces is by now well understood mathematically [1,5,11]. In the context closely related to our problem, results on propagation of interfaces for bistable reaction-diffusion equations on Riemannian manifolds were recently obtained in Refs. [44,45].
In this work, we demonstrate that understanding the patterning outcome of polar domains on a surface of arbitrary shape requires a thorough analysis of the long-timescale behavior of the associated mass-conserving reaction-diffusion equation. Working in the limit of slow surface diffusion, we establish a separation of the dynamics into three distinct time scales: first, the initial generation of the interface, then, uniform growth or shrinking of the domains until their areas converge to the steady state values dependent only on the global parameters, and finally, evolution of the interface by area-preserving geodesic curvature flow that results in a finite union of geodesic disks as time goes to infinity. In doing so, we demonstrate that stable steady states with multiple disjoint interfaces exist on biologically plausible domain shapes, something impossible without geometric effects. This paper is organized as follows. In Sec. 2, we formulate our model of polarization on a closed surface and define three sub-problems whose limiting behavior should describe the effective dynamics on different asymptotic timescales. In Sec. 3, we carry out a preliminary analysis of existence and stability of the uniform states and, in particular, identify the parameter regimes in which only nonuniform steady states can be stable. In Secs. 4 through 6, we asymptotically derive the limiting sub-problems as the surface diffusion coefficient tends to zero and examine their dynamical behavior. Finally, in Sec. 7 we perform a number of numerical tests to corroborate the predictions of the asymptotic theory, and in Sec. 8 we make our concluding remarks.

Model Summary.
Let Ω ⊂ R 3 be an open, bounded, connected set with a sufficiently regular boundary ∂ Ω. A basic model of cell polarization can be written as (see [15] for more details): Here, (2.1) describes a reaction-diffusion process on the surface ∂ Ω of the surface-bound protein concentration B, while (2.2) gives bulk diffusion of the concentration C of the same protein in the bounded volume Ω, and (2.3) is the boundary condition coupling the two. The operator ∇ 2 ∂ Ω is the Laplace-Beltrami operator on the surface, while ∇ 2 Ω refers to the standard 3D Laplacian. While we start with a model where membrane-bound dynamics is already purely two-dimensional, the validity of this class of model as a limiting case of a membrane of finite thickness was demonstrated in Ref. [31].
We are principally interested in the regime where D C D B , and in this regime it is reasonable to treat the bulk concentration as spatially uniform [15]. As mass is conserved globally, we define the quantity C tot as the total number of protein molecules divided by the bulk volume V of Ω, so that the bulk concentration can be expressed as Substituting this back into (2.1) yields a non-local equation for surface-bound species A is a rescaling of Ω that ensures that the rescaled domainΩ has boundary of unit area. Dropping the tildes and the subscript ∂ Ω for simplicity of notation from now on, we arrive at the dimensionless form for x ∈ ∂ Ω, (2.10) with the dimensionless parameters 1 and where we defined f (u) = β + u ν γ ν +u ν for further notational convenience. Together with the already dimensionless parameters β and ν, the parameters α, γ and δ define the parameter space for our problem. As was already noted, we have |∂ Ω| = 1 now.
The long timescale dynamics of this problem for a purely 1D model was previously studied via asymptotic expansions in [40], and later on a disc in [14]. In Ref. [15], exact solutions were constructed for the special case of a spherical domain and ν = ∞. Here we confine ourselves to the specific case of ν = 2 for the sake of analytical tractability, but treat general spatial domains. A few useful remarks on this system can be made based upon existing results. First, (2.1)-(2.5) represents a special case of a more general class described in Ref. [52], and per the results therein, there exists a unique classical solution (B,C) such that the functions B and C are smooth and uniformly bounded. Further, Theorem 2.6 of Ref. [26] demonstrates that our shadow system (2.8)-(2.10) likewise has a unique weak solution, and that solution is the limit of the solutions to (2.1)-(2.5) as D B → ∞, justifying the use of this reduced approach.
In the asymptotic analysis described in the following sections, we find it convenient to frame our problem as three equivalent problems parametrized by δ 1: Each problem describes the dynamics at the timescale of a different order in δ : (P δ 0 ) is the problem on the original O(1) timescale, (P δ 1 ) is the problem on the O(δ −1 ) long timescale and (P δ 2 ) is the problem on the O(δ −2 ) longer timescale. Crucially, we will show that as δ → 0 the initial condition of the second and third problem may be taken to be the infinite time limit of the solution at the previous stage, thus allowing to connect the solutions at different timescales. In our analysis below, we will formally obtain the limit behavior of each of these problems when δ → 0.

Preliminaries.
In this section we introduce some basic facts about our system, and establish some notation. We begin by examining the behavior of the spatially uniform steady states of (2.8). These states satisfy and with our choice of the nonlinearity (recall that ν = 2) this reduces to a cubic equation For α, β , γ > 0, (3.2) will have at least one real positive root u 0 , with two additional positive roots u ± emerging when its discriminant We adopt the convention u − < u 0 < u + , with u ±,0 denoting any one of these roots for shorthand. Setting ∆ = 0 allows one to derive the conditions for the existence of three roots. We find the requirements for three real roots to be Slices of the boundary between parameter space regions are plotted in the γ − α plane in Fig. 3.1 for various values of β .
We next examine the linear stability of the uniformly stable states in the usual manner, examining the behavior of perturbations of the form u ±,0 (x,t) = u ±,0 + εη ±,0 (x,t) for ε 1. Denoting the orthonormal eigenfunctions of the Laplace-Beltrami operator as where λ k are the corresponding eigenvalues, we write η ±,0 (x,t) = ∑ ∞ k=0 a ±,0 k e σ ±,0 k t v k (x). Note that for a compact surface, the eigenvalues λ k are real with 0 The linearized equation about the root u 0,± takes the form where δ k,0 is the Kronecker delta symbol, and we took into account that v 0 = 1 and ∂ Ω v k dS = 0 for all k > 0. Then the following dispersion relation holds: In the three-root parameter region, σ ± 0 < 0 and σ 0 0 > 0, while outside this region, σ 0 0 < 0, as can be verified numerically. Furthermore, one can see that if δ is sufficiently large then σ ±,0 k < 0 for all k > 0. So, as expected, for sufficiently large diffusion the stability of each fixed point can be determined entirely by its stability with respect to uniform perturbation corresponding to k = 0.
Conversely, for sufficiently small values of δ the steady states u ±,0 may lose stability with respect to non-uniform perturbations corresponding to k > 0. In this case the fastest growing mode corresponds to Then for δ < δ 0,± it follows that σ 0,± 1 > 0. As our subsequent analysis emphasizes the limit where δ → 0, it is useful to examine the boundaries of the region where f (u ±,0 )(1 − αu ±,0 ) = 1 for each fixed point. Let b 1 (α, γ) < b 2 (α, γ) be the two real roots of the depressed quartic We plot a slice of the parameter space at β = 0.025 in Fig. 3.2 in order to illustrate how the parameter space is partitioned. One can see that σ 0 1 > 0 for all points within the three root region, whereas σ ± 1 > 0 in distinct but overlapping regions. Importantly, there is a large parameter region in which no uniform state is linearly stable and, therefore, pattern formation must occur from a generic initial condition. Note that for finite δ the diffusion exerts a stabilizing effect and shrinks the regions with positive growth rates, but even for a relatively large value of δ (e.g. δ ≈ 0.1) the stability regions look nearly indistinguishable from those in Fig. 3.2. As a final remark, previous examinations of wave-pinning models have observed that even in the absence of linear instability, such systems tend to be easily destabilized by small finite perturbations -it is hypothesized that many real systems exist entirely within the stable regime but nonetheless polarize due to large scale fluctuations [40]. Regions of the parameter space at β = 0.025 colored according to the stability of the uniform steady state in the limit δ → 0. Regions i and ii denote where there is a unique uniform steady state u 0 . In the former, σ 0 k < 0 for all k, while in the latter σ 0 k > 0 for some k > 0. In iii, iv and v, three uniform steady states A final region exists on the boundary of v and i where u − becomes destabilized but u + remains stable, but this region is too narrow to reasonably visualize on this plot.

Generation of Interface.
At the zeroth order in δ , the diffusion term in (P δ 0 ) may be dropped to the leading order to obtain a closed system of equations for u 0 = lim δ →0 u δ 0 : provided g = lim δ →0 g δ 0 is smooth and the limit is uniform. All spatial coupling now occurs via the parameter U 0 . In essence, this problem reduces to an infinite system of ODEs connected via a mean-field term. Our goal in this section is to demonstrate that at almost every point x ∈ ∂ Ω, the value of u(x,t) is expected to converge to one of two specific steady state values as t → ∞ for generic choices of g. This problem has been considered before in previous studies of wave-pinning, and conventionally the local bistability of (4.1) for fixed U 0 was used to justify the assumption that u 0 indeed converges [14,40]. Nevertheless, for variable U 0 , an additional non-trivial argument is necessary to make such a conclusion. Here, we demonstrate that problem (P 0 0 ) is in fact a gradient flow generated by the following energy functional: This function is bounded from below for all β > 0, and the second term in (4.2) is nonnegative, so we have that H[u] is bounded from below. Taking the variational derivative of H yields which is a gradient flow generated by the energy functional H[u] with the mobility f (u). As f (u 0 ) > 0 uniformly, we have Here equality only occurs when u 0 is one of the steady state solutions to (4.1), which due to the dissipative nature of the dynamics are the only elements of the attractor of the long-time dynamics of the solutions of (P 0 0 ) [25]. Away from the steady states, H[u 0 (·,t)] is both decreasing in time and bounded below, so as t → ∞, H(u 0 (·,t)) → H ∞ . We note that existence of solutions for problem (P 0 0 ) with positive bounded initial data follows from the standard theory of ODEs in Banach spaces, yielding a unique solution [10]. In particular, this gives U 0 ∈ C ∞ ([0, ∞)). Therefore, by the standard existence theory for ODEs we have that u(x, ·) ∈ C ∞ ([0, ∞)) for each x ∈ ∂ Ω. In fact, it is not difficult to see that the solutions of (P 0 0 ) remain positive and bounded independently of t, as they should. Thus, the arguments leading to (4.5) are justified. Furthermore, the boundedness of u 0 yields boundedness of ∂ 2 t u 0 as well, which together with (4.5) implies that ∂ t u 0 (·,t) → 0 in L 2 (∂ Ω) as t → ∞.
Consider a sequence of t n > 0 such that t n → ∞ as n → ∞. Up to an extraction of a subsequence (not relabeled), we then have ∂ t u 0 (·,t n ) → 0 a.e. in ∂ Ω. Since U 0 (t n ) is bounded, upon a further extraction of a subsequence we have U 0 (t n ) → U ∞ 0 as n → ∞. Therefore, for a.e. x ∈ ∂ Ω we have that u(x,t n ) converges to one of the roots of the equation This equation has at most three real roots, which we label h − (U) < h 0 (U) < h + (U), and plot their dependence on U in Fig. 4.1. The conditions for all three roots to be real are 0 < β < 1/8 and correspond to saddle node bifurcations. We calculate the explicit dependence of these on the kinetic parameters as Notably, the form of Φ (h) implies that dh ± (U) dU < 0. Therefore, the roots u 0 = h ± (U ∞ 0 ) represent the two linearly stable equilibria of (4.1), while u 0 = h 0 (U ∞ 0 ) represents an unstable equilibrium for U 0 = U ∞ fixed. Thus, generically we would expect that for almost every x ∈ ∂ Ω we have , the following system of equations arises in the limit as n → ∞: We now consider the above equations as an algebraic system with U ∞ 0 fixed and |∂ Ω ± 0 | and H ∞ as variables. The conservation of area (4.9) ensures that d |∂ Ω − 0 |, and differentiating (4.10) yields (4.12) dh Rearranging the terms then establishes that |∂ Ω + 0 | is an increasing function of U ∞ 0 , since As the final step, taking the derivative of (4.11) and using the identity Φ (h ± (U ∞ 0 )) = −αU ∞ 0 yields the relationship The term inside the brackets only equals zero at a single value of U ∞ 0 , corresponding to the double tangent construction (a straight line that touches the graph of Φ(u) at two points). For all other values of U ∞ 0 , we thus have dH ∞ /dU ∞ 0 = 0, as shown in Fig. 4.2 for several choices of α. It follows that there are at most two isolated solutions U ∞ 0 to (4.9)-(4.11) for a given value of H ∞ , i.e., the set of possible values of U ∞ 0 for different choices of the sequence (t n ) is discrete. Thus, since the ω-limit set of u 0 is connected [25], convergence of H(u 0 (·,t n )) → H ∞ as n → ∞ fixes U ∞ 0 independently of the subsequence. So we have lim t→∞ U 0 (t) = U ∞ 0 , and by the ODE stability argument we then have lim t→∞ u(x,t) = h ± (U ∞ 0 ) as well for each x ∈ ∂ Ω ± 0 , respectively. This arguments also shows that up to sets of measure zero the sets ∂ Ω ± 0 are independent of the sequence, and so we finally obtain a full limit: for some ∂ Ω ± 0 ⊂ ∂ Ω depending only on the parameters and the initial condition and satisfying |∂ Ω − 0 | + |∂ Ω + 0 | = 1. In other words, the solution of problem (P 0 0 ) is expected to generically converge to a piecewiseconstant function corresponding to the two stable equilibrium branches of the bistable nonlinearity for a certain long time limit value of U 0 . This result can be interpreted as follows: for small but finite value of δ the solution of problem (P δ 0 ) will become nearly piecewise constant on the timescale 1 t δ −1 almost everywhere in ∂ Ω, except in a small transition region of width of order δ . Thus, an interface is expected to form on this timescale, which will serve as an approximate initial condition for problem (P δ 1 ) at the corresponding rescaled timescale δ t 1 for the latter. We note that while one should expect u = h ± (U ∞ 0 ) to be generically selected as the limiting values of the solution away from the interfaces, for arbitrary initial conditions the possibility of h 0 (U ∞ 0 ) on a set of positive measure may not be a priori excluded. This prevents us to make (4.15) into a rigorous conclusion about the long-time limit behavior for problem (P 0 0 ). This sort of issue is well known in the studies of gradient flows and its further treatment would require a more delicate analysis that goes beyond the scope of the present paper (for a treatment of a closely related problem, see [3]).
Lastly, we point out that the arguments above do not a priori exclude the possibility of |∂ Ω + 0 | = 0 or |∂ Ω − 0 | = 0, depending on the initial conditions. In these cases the system ends up in a uniform state and no further non-trivial dynamics governed by (2.8) -(2.10) is expected to occur. This outcome of the dynamics may be interpreted as a failure of cell polarization. However, such an outcome would be non-generic for the kinetic parameters corresponding to the regions ii and iii in Fig. 3.2, in which all uniform states are linearly unstable, consistently with the conclusion at the end of Sec. 3. In contrast, the case 0 < |∂ Ω + 0 | < 1 corresponds to a patterned state as the outcome of the dynamics of (P 0 0 ), in which the boundary between ∂ Ω + 0 and ∂ Ω − 0 represents an interface between the high-and low-concentration domains, respectively.

Convergence of surface concentration.
Having established the conditions for generation of the interface bounding uniform domains via non-local coupling, the next step is an examination of the evolution of these interfaces. At O(δ −1 ) timescale (in the original variables), the dynamics of our system is described by problem (P δ 1 ). We claim that in the limit δ → 0 the solution u δ 1 of this problem converges to that of (5.1) Here ∂ Ω ± 1 (t) are the time-dependent subsets of ∂ Ω corresponding to the quasi-steady states h ± (U 1 (t)) and Γ(t) ⊂ ∂ Ω is a smooth closed curve or a collection of curves representing the boundary between ∂ Ω + 1 (t) and ∂ Ω − 1 (t), i.e., Γ(t) is the interface. With some abuse of notation, we denote by ∂ Γ(x,t)/∂t the normal velocity of Γ(t) at each point x ∈ Γ(t) in the direction of ∂ Ω − 1 (t), with the function c(U 1 ) to be specified. As was already discussed at the end of Sec. 4, the sets ∂ Ω ± 1 (t) at t = 0 coincide with the sets ∂ Ω ± 0 obtained in the long time limit of the solution of problem (P 0 0 ). To aid in the derivation of (P 0 1 ), it is convenient to define Ψ = Ψ(x,t) which is the signed distance function from Γ(t), with Ψ > 0 corresponding to u 1 = h + (U 1 (t)). Defining the stretched distance z = Ψ(x,t)/δ to Γ(t), we seek the solution of (P δ 1 ) in the neighborhood of Γ(t) in the form u δ 1 (x,t) v(z,t), where v is an unknown function [18]. Then to the leading order in δ we obtain in some tubular neighborhood of Γ(t): where the function c(U) is obtained by solving the traveling wave equation The interface Γ(t) is then reconstructed from the condition Ψ(x,t) = 0 for x ∈ Γ(t). It is well known that for all values of U for which the nonlinearity in (5.3) is of bistable type, this equation has a unique solution (up to translations) for a unique value of c(U) [19]. In particular, there is a unique value U = U ∞ 1 such that c(U ∞ 1 ) = 0, which is obtained by solving the equation and the sign of c(U) coincides with the sign of the integral in (5.4) for U = U ∞ 1 . It is also possible to show that dc(U)/dU < 0. Note that (5.2) is equivalent to the last equation in (5.1), with the function c(U) given implicitly by the solution of (5.3). The rest of the equations in (5.1) are obtained by passing to the limit δ → 0 in the ansatz for u δ 1 . The dependence c(U) obtained from the numerical solution of (5.3) for a particular choice of the parameters is illustrated in Fig. 5

.1. This calculation entailed time integration of
for fixed values of U, the initial condition and the boundary conditions ∂ x u(±L, 0) = 0. Integration was performed using the Dedalus spectral package on a domain [−L, L] with L = 100, discretized into 1024 modes [8]. Each simulation was allowed to relax until t = 100 with dt = 0.1, at which point the front velocity was calculated as To study the long time behavior of the solutions of problem (P 0 1 ), we argue as in [27] and note that the growth rate of |∂ Ω ± 1 (t)| is related to c(U 1 ) as From this, we have .
As 0 < k(t) < ∞ uniformly so long as the interface persists (recall that dc(U)/dU < 0), as t → ∞ we find that U 1 (t) → U ∞ 1 exponentially fast and, therefore, Γ(t) → Γ ∞ 1 , with Γ ∞ 1 separating the two limiting domains ∂ Ω ± 1,∞ in ∂ Ω such that lim t→∞ u 1 (x,t) = h ± (U ∞ 1 ) for all x ∈ ∂ Ω ± 1,∞ . Again, this behavior can be interpreted as follows: for times 1 t δ −1 the solution of problem (P δ 1 ) will approach the long time limit of the solution of problem (P 0 1 ), except in a transition layer of width δ . This solution corresponds to the times δ t 1 for problem (P δ 2 ) and will, therefore, serve as the approximate initial condition for the latter.
An interesting feature of the long time behavior of problem (P 0 1 ) is that as t → ∞, the total amount U of material on the surface converges to a value that depends only on α, β , and γ, and is independent of the initial condition. Furthermore, the long time value of (5.11) |∂ carries no explicit dependence of diffusion or domain geometry, other than through the dependence of the dimensionless parameter α on the surface to volume ratio A/V . The dependence of U ∞ 1 on α and γ can be expressed as U ∞ 1 = 1−b(β )γ α for some function b(β ). We illustrate this fact in Fig. 7.2, which demonstrates the convergence of U(t) → U ∞ from three distinct initial conditions. Some consideration is merited on the possible elimination of fronts: arbitrary initial conditions imply an arbitrary number of closed interfaces might be created during the initial dynamics, and some of these interfaces may vanish under the dynamics described above. However, provided we restrict ourselves to a parameter choice where β , γ)), as well as U − < U ∞ 1 < U + as per the previous section, both h + and h − will exist as real roots and conservation of mass will require |∂ Ω + 1 (t)|, |∂ Ω − 1 (t)| > 0 for all t. Thus, in this regime, at least one interface exists as t → ∞. 6. Area-preserving geodesic curvature flow. In the previous section, we demonstrated that on the O(δ −1 ) timescale (in the original variables) the dynamics results in a stationary interface with a prescribed quantity of matter on the surface and a fixed area of the domains bounded by the interface. We now show that on longer timescales there is no longer a net transport of mass between the boundary ∂ Ω and the enclosed volume Ω, and the dynamics may be reduced to a gradient flow driven by an interface length subject to domain area conservation [34].
The dynamics on the O(δ −2 ) timescale (again, in the original variables) is described by problem (P δ 2 ). We now claim that in the limit δ → 0 the solution u δ 2 of this problem converges to that of (6.1) Here, as before, ∂ Ω ± 2 (t) are the time-dependent subsets of ∂ Ω corresponding to the quasi-steady states h ± (U ∞ 1 ) and Γ(t) ⊂ ∂ Ω is a smooth closed curve or a collection of curves representing the boundary between ∂ Ω + 2 (t) and ∂ Ω − 2 (t). The function k g (x,t) refers to the geodesic curvature of Γ(t) at x ∈ Γ(t), with the sign convention that k g < 0 if ∂ Ω + 2 is a small geodesic disk, and ξ (t) is a Lagrange multiplier ensuring conservation of |∂ Ω ± 2 (t)|. The sets ∂ Ω ± 2 (t) at t = 0 now coincide with the sets ∂ Ω ± 1,∞ obtained in the long time limit of the solution of problem (P 0 1 ). To derive this problem, we again denote by z = Ψ(x,t)/δ the stretched distance to the interface and seek the solution in the form [18] where v, w and ζ are to be found. Substituting this ansatz into (P δ 2 ) yields (6.3) To eliminate the O(δ −1 ) term from the equation, we choose v to be a solution of (5.3) with c(U) = 0, which exists if and only if ∂ Ω vdS = U ∞ 1 . By continuity, the latter can always be fixed by a suitable choice of ζ (t) = O(1) for a given Γ(t) for which the condition Having now eliminated the O(δ −1 ) term from (6.3), we can write the solvability condition for w, which is obtained by multiplying the equation by ∂ z v and integrating over z [18]: We note that if y ∈ Γ(t) is the projection of x ∈ ∂ Ω from a small tubular neighborhood of Γ(t) onto Γ(t), then ∇ 2 Ψ(x,t) = ∇ 2 Ψ(y,t) + o(1) for |x − y| 1. Thus, after dropping the o(1) terms and rearranging the formula, we can rewrite (6.4) as The obtained equation (6.5) is nothing but the last equation in (6.1) written in terms of the distance function [2], recalling the definition (with our sign convention) of geodesic curvature k g (x,t) = ∇ · n(x,t), where n(x,t) = ∇Ψ(x,t) is the unit normal to the curve Γ(t) at x ∈ Γ(t) that lies in the tangent plane to ∂ Ω at x and points towards ∂ Ω + 2 (t). Finally, conservation of |∂ Ω ± 2 (t)| follows by passing to the limit δ → 0 in our ansatz, and the initial condition is given by the long time limit of the solution of problem (P 0 1 ), i.e., we have Γ(t) = Γ ∞ 1 for t = 0. Note that the value of |∂ Ω + 2 | is given explicitly by the right-hand side of (5.11).
Since the area of |∂ Ω ± 2 (t)| is constant for the solutions of problem (P 0 2 ), we have ∂ t |∂ Ω ± 2 (t)| = 0. Furthermore, since the growth rate of a region's area is equal to the integral of its normal velocity, we can eliminate ξ (t) to obtain a purely geometric expression for the normal velocity: where k g (·,t) Γ(t) = 1 |Γ(t)| Γ(t) k g (x,t)ds and |Γ(t)| is the length of Γ(t). This equation describes an areapreserving geodesic curvature flow on a surface. A geometric identity relates the evolution of interface length to curvature as d dt |Γ(t)| = − Γ(t) k g (x,t) ∂ Γ(x,t) ∂t ds [30]. From this, we obtain an inequality which demonstrates that the interface length |Γ(t)| is strictly decreasing except for when k g is constant along Γ(t). The interpretation of this result is that for small δ the long time limits of the solutions of problem (P δ 2 ) will be close to domains of prescribed area and locally minimal interfaces of constant curvature (geodesic disks).
The problem of identifying locally minimal interfaces on a given surface has been studied at length over the years, and while no comprehensive solution exists, there is a considerable literature of results which can inform our understanding [28,38,50]. In particular, it is known that for every 0 < a < 1 there exists a set ∂ Ω + ⊂ ∂ Ω whose boundary Γ globally minimizes the perimeter among all sets with |∂ Ω + | = a [50]. In this case Γ is a smooth curve of constant geodesic curvature. Framed in terms of our reaction-diffusion problem, this means that for appropriate choices of α, β , and γ, a non-uniform steady state for (P 0 2 ) should always exist regardless of the particular shape of Ω. Whether the solution of problem (P 0 2 ) converges as t → ∞ to a global minimizer, or, indeed, converges, however, is not known a priori. We note that due to the area conservation the interface Γ(t) associated with the solutions of problem (P 0 2 ) cannot vanish. Furthermore, due to the gradient flow nature of the dynamics the ω-limit set of these solutions consists only of steady states, and convergence to one of the steady states is guaranteed if the set of the steady states is discrete. In particular, this would be the case if the global minimizer of the perimeter were the unique constant curvature solution with the prescribed area. We note that the latter would represent the most robust and reproducible scenario for cell polarization.
We can also use existing variational results to construct surfaces which yield some surprising solutions. It has been widely suggested that strictly mass-conserving reaction-diffusion equations will exhibit uninterrupted coarsening and admit steady states with only a single domain, but this phenomenon has only been studied in the Euclidean setting [7,55]. Multi-domain polarization has in fact been observed in real cells, and failure to account for this is seen as a failure point of minimal models such as the one in this paper [13]. However, we find that when geometric effects are properly considered, multi-cap solutions are possible on some surfaces, which we now demonstrate by constructing examples. Per the theorems in Ref. [49], if ∂ Ω is a surface of revolution about the z-axis that is symmetric with respect to reflections around the xy-plane, and if ∂ Ω has Gaussian curvature K(z) that is strictly increasing in |z|, then the locally minimal interfaces bounding a fixed area can be fully classified [49], and are of the types depicted in Fig. 7.1. All stable interfaces lie in the planes z = const, which means k g can be written as a function of z, and it is convenient to define a function on ∂ Ω of the form L(z) = k 2 g (z) + K(z). In addition to single-interface solutions, there are two possible solutions with disjoint interfaces. Case 1 features two interfaces at ±z 0 for some z 0 > 0, respectively, which are stable provided they lie within the region where L(z 0 ) < 0. Case 2 has one interface at z 1 (s.t. L(z 1 ) < 0) and the other at some z 2 such that z 1 z 2 < 0 and L(z 2 ) > 0. As the positions of ∂ Ω + and ∂ Ω − are interchangeable under geodesic curvature flow, each case can either represent a two cap solution, or a single high-concentration band containing z = 0. In essence, these example solutions show the existence of a geometry-induced interruption to coarsening. As an aside, results such as those in Refs. [38,49,50] do rule out multi-cap solutions on symmetric convex surfaces such as spheroids.

Numerical Study.
We now verify numerically the existence of nontrivial steady state solutions predicted by the asymptotic analysis of (P δ 2 ). We simulated (2.8) -(2.10) on the surface shown in Fig. 7.1 for each of the possible configurations, using a recently developed numerical scheme [36]. (for a = 0.2728 and b = 0.3896) about the vertical axis. The simulations were performed on a spatial discretization of 128 × 128 modes and integrated with dt = 0.1 until t = 500, at which point a numerical steady state was achieved. In each case, the initial conditions were step functions approximating the desired steady state. The steady states were verified using a GMRES-based Newton-Krylov method [29]. This method further allowed us to compute approximate values of the dominant eigenvalues, and in each case these were negative, consistent with the proposed stability of these states. While the constraint on U ∞ for the single interface case is as described in the previous section, both case 1 and 2 of the double interface solutions have a narrower range of acceptable U ∞ , owing to the constraints on the locations of the interfaces. As such, we chose a different α for each simulation to ensure a viable solution.
To corroborate the multiple stages of the dynamics predicted by our asymptotic analysis, we performed two additional numerical studies of (2.8) -(2.10), again using the integration methods introduced in Ref. [36]. Recall that the results of Sec. 4 indicate convergence of U(t) to a constant value U 1 depending on the initial condition on a large O(1) time scale. Furthermore, by the results of Sec. 5 the function U(t) converges to the final value U ∞ 1 (α, β , γ) independent of initial conditions on a large O(δ −1 ) time scale. As the first test of this, we examine the convergence of U(t) in three simulations. These were initialized at u 0 + γσ (x)/2, where u 0 is as defined in Sec. 3, and σ (x) is a random function taking on values between −1 and 1, with the spatial correlation length 0.6. Further details on how this distribution was generated can be found in the supplementary materials of Ref. [36]. Simulations are performed on a sphere of unit area, discretized to 256 × 256 modes, and integrated with a timestep of dt = 0.1. The results of the simulations for δ = 10 −3 are plotted in Fig. 7.2 (see the figure caption for the values of other parameters). The trajectories of all three examples can be split into two types of behavior. Dynamics from t 0 = 0 to t 1 ≈ 10 1 exhibit a rapid change in U(t) towards an approximately constant value: this is the timescale of the interface generation. On the other hand, on the interval from t 1 ≈ 10 1 through t 2 ≈ 10 4 , the different instances of U(t) converge to the predicted steady state at U ∞ 1 from Sec. 5. A second numerical experiment is to study the behavior of the interface length over time. For δ 1, an effective proxy for the total interface length is the quantity ∂ Ω |∇u(·,t)|dS, as in the last stage of the evolution the profile u(x,t) will be close to two fixed values everywhere outside the thin transition layers approximating the interfaces. We initialize a simulation at a small perturbation about the steady state as u 0 + 0.01σ (x) and plot the evolution of the interface in Fig. 7.3. The simulation domain was a unit-area prolate spheroid with an aspect ratio of 2, discretized into 128 × 128 modes. As before, our timestep was dt = 0.1.
For our choice of δ = 10 −2 here, the distinct regimes corresponding to the limit behaviors of (P δ 0 ), (P δ 1 ) and (P δ 2 ) as δ → 0 are very clearly expressed. Much as in the previous example, the timescale of interface formation is on the order of t ≈ 10 1 : the order of magnitude increase in δ between the two figures has a minimal effect on the zeroth order behavior of the system. At the end of this interval, surface concentration is fully partitioned into approximately uniform domains ∂ Ω ± . Once the quantity |∇u| peaks, it begins a period of slower evolution consistent with the predicted front motion of (P 0 1 ), which then plateaus once U ∞ 1 is reached. A slower evolution of |∇u| , which is now monotonically decreasing, subsequently occurs under the approximate (P 0 2 ) dynamics. A brief note should be made of the step-like shape of the plot in this final interval. Contraction under (P 0 1 ) front motion reduces this simulation to two spots, and the rapid drop in interface length which occurs around t = 2000 represents the annihilation of one of these. After this, the remaining spot persists for a long period near midway between the two poles before finally migrating to one end of the spheroid. The final steady state is at the point of maximum Gaussian curvature, consistent with the geometric results noted in the previous section. Taken together, these numerical tests represent strong evidence that the formal asymptotics developed here accurately reproduces the behavior of (2.8) -(2.10) for small but finite choices of δ .
8. Discussion. In summary, we performed an extensive analysis of the asymptotic behavior of a minimal model of cell polarization on a 2D manifold. Our results establish a hierarchy of three timescales of the dynamics in the case of small surface diffusivity. On the short timescale, we derived an evolution equation in the form of a globally coupled system of ODEs. While our initial problem itself appears to lack a variational structure, the reduced dynamics is shown to be a gradient flow generated by a nonlocal energy functional. From this, we were able to demonstrate generation of interfaces between high-and low-concentration regions for generic initial data. Our argument relied on the assumption that the solution for almost all points of the surface approached one of the two stable branches of the bistable nonlinearity, which is expected to be generically true but cannot be established for all initial data to make the result fully rigorous. Let us mention that in Ref. [3], in a closely related nonlocal Allen-Cahn problem, the authors were able to avoid this kind of assumption by constructing a second Lyapunov functional (in fact, an infinite family of Lyapunov functionals) to prove full convergence in time at the expense of allowing unstable solutions in the limit as well. However, the direct coupling between the nonlocal and nonlinear terms in problem (P 0 0 ) greatly complicates this sort of argument, so it remains unclear whether this approach could be adapted to our problem.
On the intermediate timescale, we demonstrated in the sharp-interface limit that the interfacial velocity is uniform and depends only on the enclosed area, and that this enclosed area approaches a final value determined only by the kinetic rates. This area is determined strictly by the dimensionless parameters of our model, and thus is independent of the surface geometry. This fact raises a question as to how robust cell polarization can be to perturbations of the cell shape. This question could be particularly interesting in light of recent interest in coupled chemo-mechanical models of deforming cell membranes [17,22,35,37].
Finally, on the long timescale, we demonstrated that the dynamics reduce to area-preserving interface motion by geodesic curvature. This result formally establishes the intuitive notion that on curved surfaces wave-pinning phenomena is driven by a gradient flow generated by the interface length. The steady state of such a flow has a highly nontrivial dependence on surface shape -the ability of some surface geometries to support several distinct solutions suggests important implications for cellular decision-making and addresses a perceived failure of minimal mass-conserving reaction-diffusion models. One particularly interesting question for future work characterizing mechanisms to select between different steady states. Very recently theoretical work has shown that chemical gradients can enforce state transitions in simple geometries [9]. Examining similar effects on geometries such as those in Fig. 7.1 could provide useful insight into biological control mechanisms.
There is merit in some further comparison between our results and some deceptively similar results from flat surfaces. In the 1D version of this system, multi-domain solutions can also form, but these are strictly metastable and vanish on exponentially slow timescales [40,53]. By contrast, the multi-domain solutions discussed here can be truly stable. Previous works such as Ref. [57] have also showed that domain geometry can establish multi-domain solutions in wave-pinning models. One key difference is that these studies examined flat 2D domains bounded by no-flux boundaries of various shapes, rather than curved surfaces. Both versions can produce multi-domain solutions, but a key difference lies in the time needed to reach steady states. Specifically, a circular domain in the flat scenario will be metastable if it is sufficiently far from the boundary, and like in the 1D case discussed above it will only approach the steady state location exponentially slowly. This slowness stems from the approximate translational invariance of the interfaces away from the boundaries. Non-uniform local curvature breaks this invariance, making the convergence to steady state domain locations considerably faster, on the algebraic rather than exponentially long timescale in terms of the inverse diffusivity.
The wave-pinning model studied here represents one of the simplest possible characterizations of coupled bulk-surface dynamics in a single cell. A great number of more complex models have been proposed over the past decade, including ones featuring multiple chemical species [43] and advection [4,23]. Additionally, recent experiments have demonstrated that bulk-surface systems are capable of much more complex pattern formation than discussed here, such as the spiral wave dynamics observed in starfish embryo [32,54], establishing a long-term need for further analysis of this class of system.