Supporting Text S1: Inﬂuence of cell shape, inhomogeneities and diﬀusion barriers in cell polarization models

Supporting Text S1: Influence of cell shape, inhomogeneities and diffusion barriers in cell polarization models Wolfgang Giese, Martin Eigel, Sebastian Westerheide, Christian Engwer, Edda Klipp1,∗ 1 Theoretical Biophysics, Humboldt-Universität zu Berlin, Berlin, Germany 2 Weierstrass Institute, Berlin, Germany 3 Institute for Computational and Applied Mathematics, University of Münster, Münster, Germany 4 Cells in Motion (CiM) Cluster of Excellence EXC 1003, University of Münster, Münster, Germany ∗ Corresponding author: Edda Klipp, e-mail: edda.klipp@rz.hu-berlin.de


Modeling framework
In our models, the computational domain comprises the cytosolic volume V cyt ⊂ R 2 , the outer cell membrane M cell ⊂ ∂V cyt and membranes M org = ∂V cyt \ M cell that separate the organelles from the cytosol. A sketch of the computational domain can be seen in Figure S1. Furthermore, we equip our models with diffusion rates which are not necessarily constant along the cell compartments M cell and V cyt . Given some initial values u(·, 0) on M cell and v(·, 0) in V cyt , the model equations incorporating the described cell geometry and inhomogeneous diffusion read with boundary conditions −D c ∇v · n = 0 on M org × [0, T ] for the cytosolic equation (2). Here, n denotes the vector field of outer unit normals on M cell and M org , respectively. For this system of equations, we have mass conservation where K ∈ R. This conservation property directly follows from the model equations' weak formulation, see Section 2.1.

Numerical methods used for simulations
In the following, the numerical approach is described that is used for simulating the class of mathematical models from Section 1. As for all numerical approaches, the fundamental idea is to discretize the mathematical model. This process yields a system of algebraic equations which can be solved using a computer.
To separate discretization in space and time, we use the well-known method of lines, see e.g. [1], and start with a semidiscretization in space using conforming first-order finite element methods based on triangular meshes. Given a triangular mesh which describes the geometrical setup introduced in Section 1, the equation for the volume species is treated by the standard conforming finite element approach, employing Lagrange basis functions of polynomial degree one. To treat the equation describing the membrane-bound species, we apply a surface finite element method on the boundary of the mesh using a restriction of the same volumetric basis functions. This enables implementing the numerical approach with tools provided by standard software frameworks for scientific computing. The idea of performing spatial discretization by combining the conforming finite element method and surface finite elements on the same mesh is related to the procedure presented in [2], where a similar approach is used for the discretization of a coupled elliptic model problem.
With the method of lines, different schemes can be employed for the discretization in time. In accordance with our spatial discretization, we restrict ourselves to first-order schemes. We use both a fully implicit scheme and a semi-implicit scheme inspired by the implicit-explicit (IMEX) Euler method presented e.g. in [3]. The idea of the IMEX method is to treat the spatially discretized reaction part of the system explicitly. Therefore, on the one hand the membrane equation and the cytosolic equation are decoupled and thus can be treated separately, and on the other hand the non-linearities are treated explicitly which enables using a linear solver. On the downside, an explicit treatment of the reaction part can affect the stability of the scheme in the reaction-dominated case and for a stiff reaction part in gereral [3]. Our semi-implicit scheme decouples membrane and cytosolic equations while still treating the reaction part of each separate equation and its non-linearities implicitly. For each equation, this is done by treating only the unknowns of the other equation explicitly.

Weak formulation
First, we derive a weak formulation of model equations (1) -(4) in order to apply finite element methods for semidiscretization in space. Let V bulk := H 1 (V cyt ) denote the usual Sobolev space containing weak solutions of elliptic equations in the bulk domain V cyt . A natural counterpart containing weak solutions of elliptic equations on hypersurfaces are surface Sobolev spaces [4][5][6]. To treat equation (1) on the closed hypersurface M cell , we therefore define the surface Sobolev space V sur := H 1 (M cell ).
Multiplication of model equations (1) and (2) with some test functions ϕ m ∈ H 1 (M cell ) respectively The weak formulation of model equations (1) -(4) now is to look for a solution (u, v) Note that also the constant test functions ϕ m ≡ 1 respectively ϕ c ≡ 1 are permitted which yields mass conservation (5).

Semidiscretization in space
To obtain a semidiscretized system, we combine the conforming finite element approach (FEM) and a surface finite element method (SFEM). The FEM is a standard approach which is well-known to literature. See e.g. [7,8] to gain an insight into the methodology. An SFEM developed in [5] can be seen as a natural generalization, as the idea of FEM is transfered to elliptic equations on hypersurfaces. Its extension [6] to treating parabolic equations, like membrane equation (1), provides the basis for the SFEM that we use.
Both approaches are based on an approximation of the bulk domain V cyt and the hypersurface M cell , each by a triangulable geometrical object, and corresponding meshes. For simplicity, we assume V cyt to be a polyhedral domain that can be exactly represented by a triangular mesh T bulk h . With M cell being part of the boundary of V cyt , it corresponds to a set of boundary entities of T bulk h which make up a surface mesh T sur h ⊂ T bulk h . Each method uses its corresponding mesh to set up a finite-dimensional function space usable for spatial discretization of the model equations. In particular, we replace the function spaces V sur and V bulk by finite-dimensional conforming function spaces V sur As discrete function spaces, we employ the node-based Lagrange spaces of polynomial degree one on The semidiscrete solution (u h , v h ) can be represented as . Therefore, we get a system of ordinary differential equations with mass and stiffness matrices and vector valued right hand side functions f m , f c ,

SFEM matrix/vector assembly via volumetric FEM basis functions
As we use node-based Lagrange spaces of the same polynomial degree based on the same triangular mesh, the set of Lagrange nodes of V sur h can be described as X sur

Moreover, the employed Lagrange basis functions have the property, that
Thus, instead of directly implementing the SFEM ansatz space V sur h , the constrained function space can be utilized to perform the assembly of M m , S m , f m and f c . In particular, the restriction of its basis can be calculated as

Fully discretized systems
System (19) can be equivalently written as using the notation Employing the backward Euler method in (29), we end up with a system of nonlinear algebraic equations which can be solved e.g. using a multidimensional Newton's method. However, solving this fully discretized system can get very time-consuming, especially when a fine finite element mesh is used, as this results in high-dimensional coefficient vectors. For this reason, we employ a semi-implicit scheme for temporal discretization which decouples equations (27) and (28). In particular, we use the backward Euler method separately for each equation, while treating the unknowns of the other equation explicitly.
In contrast to the IMEX Euler method (see e.g. [3]), the unknowns of each separate equation are thus treated fully implicit. The fully discretized system then reads and can be solved e.g. using a multidimensional Newton's method for both equations in parallel.
Our experiments have shown no obvious qualitative differences between the solutions computed with both schemes (31) and (32).

Simulation framework
The presented numerical approach can be implemented with tools provided by standard PDE software frameworks. For the assembly of the required matrices, the framework has to provide the space of simplicial Lagrange finite elements of order one on a triangular mesh. Furthermore, it either has to provide the space of simplicial Lagrange finite elements of order one on the boundary of the same mesh or a mechanism for constraining the degrees of freedom of the volumetric Lagrange space which do not lie on the boundary of the mesh. The latter mechanism usually is available in those frameworks, since constrained degrees of freedom are frequently used to implement Dirichlet boundary conditions. All of our simulations were performed using the Distributed and Unified Numerics Environment (DUNE) [9,10]. The numerical discretization schemes were implemented using the discretization module DUNE-PDELab which is based on DUNE. It provides the volumetric finite element space V bulk h which in addition can be constrained for SFEM matrix/vector assembly using the space ∂V bulk h . Furthermore, it features an easy to use assembly infrastructure, as well as the linear solver, nonlinear solver and time-stepping schemes that were used.
The finite element meshes were generated with Gmsh [11], examples of which are depicted in Figure  S4. All of these meshes comprise several thousand elements to guarantee a high precision of the method. Additionally, meshes are refined near the outer cell membrane M cell in our simulations to enable an accurate computation of the coupling process on the surface.
All code and meshes used in this work are freely available from the authors on request.

Alternative approaches
An alternative approach usable for simulating the mathematical models is presented in [12]. It uses an implicit description of the cellular geometry via level set functions and seems to be promising especially in the context of geometrical setups that are more complex than those investigated in this work, like those arising from real microscopy data. Particularly, it might be useful for coping with a cellular geometry that even evolves in time. Furthermore, already for simple geometries it can be advantageous as there is no need for a customized mesh generation. This can ease the simulation workflow. A similar approach can be found in [13], which also has the advantage of using an implicit geometry description, but uses a diffuse-interface representation via a phase-field function.
Another alternative approach is presented in [14] and implemented in the cell modeling and simulation software Virtual Cell [15]. Though employing a cell-centered finite volume scheme on a regular mesh and thus approximating the membrane M cell in a staircase manner, convergence results for static spherical geometries are presented in [14] that render the method well-suited in the considered case.

Linear stability analysis
We used linear stability analysis as in [16] to capture the essential behavior and parameter dependencies of the GOR and WP model. For the analysis we consider the following simplifications. First, we assume the computational domain to have a circular shape with radius R ∈ R, centered at the origin of ordinates. In polar coordinates (r, φ), the Laplace-Beltrami operator on the boundary can then be expressed as Note that ∆ = ∆ r + ∆ φ in polar coordinates, with ∆ r := ∂ 2 ∂r 2 + 1 r ∂ ∂r . Second, we assume that there are no inhomogeneities on the membrane as well as in the cytosol. Expressed in polar coordinates, the equations from Section 1 then read with only one boundary condition for the cytosolic equation (34). Let (ū 0 ,v 0 ) be a steady state of our system such that f (ū 0 ,v 0 ) = 0. We consider small perturbations from this steady state which are described by functions ξ m = u −ū 0 and ξ c = v −v 0 . For these perturbations, we have in the linear regime, and boundary condition

Power series expansion using two-dimensional polar harmonics as ansatz
Since R 2 ∆ Γ h = ∂ 2 h ∂φ 2 = −k 2 h for h(φ) = cos(kφ) and h(φ) = sin(kφ), the eigenfunctions of ∆ φ | r=R are the sine and cosine functions. Therefore, we expand the functions ξ m and ξ c as follows: Due to linearity of the linearized system (36) -(38) with respect to solutions (ξ m , ξ c ), it is sufficient to identify solutions (ξ k m , ξ k c ), as defined in equations (39) -(40), via their corresponding parameters first. We thus determine the parameters a k m , b k m , a k c , b k c and λ k separately for each k = 0, 1, 2, . . . . Inserting the ansatz for ξ k c into cytosolic equation (37), we get After canceling, this takes the form which leads to Mathematical solutions of this equation are modified Bessel functions of the first kind. We thus obtain A k (r) = I k (r λ k /D c ). Inserting ξ k m and ξ k c into cytosolic boundary condition (38), we get Since the choice of φ is arbitrary, this leads to Therefore, a k c as well as b k c can be expressed in terms of the coefficients a k m and b k m , respectively. This yields Putting ξ k m and ξ k c into equation (36) for processes on the membrane, we can deduce .
After canceling the nonzero terms, we finally get . (43) Note that initial conditions for system (36) -(38) would be required in order to determine the free parameters a k m and b k m . Having an equation for the parameters λ k , we now want to identify conditions which lead to growing eigenmodes. Following [17], the derivative of I k can be expressed as Since I k (x) > 0 for x > 0 and for all k = 0, 1, 2, . . . , we have the inequality We use this inequality for the derivative ∂ ∂r I k and furthermore assume f u (ū 0 ,v 0 | r=R ), f v (ū 0 ,v 0 | r=R ) ≥ 0. This holds true for the GOR model in the stable steady state and for their WP model in the transient steady state (v 0 ,ū T 0 ), where phase seperation occurs. Therfore, we obtain: . This finally yields the following inequality for λ k : In a physiological range, the term fu(ū0,v0| r=R )fv(ū0,v0| r=R ) is very small and we get For k = 0, the fundamental solutions ξ 0 m (φ, t) = a 0 m exp(λ 0 t) and ξ 0 c (r, φ, t) = A 0 (r)a 0 c exp(λ 0 t) do not depend on the polar angle φ and, therefore, do not contribute to cellular asymmetry. Hence, we only consider wave numbers k = 1, 2, 3, . . . in the following. For cell sizes R small with we have no growing eigenmodes since λ 1 λ 2 . . . . For cell sizes R large with we have more than one growing eigenmode and therefore several clusters.

Cluster width approximation for the WP model
In the following, we want to examine the size of the cluster for the WP model analytically. Stimulation of the system, as performed in the simulations in the main text, causes a traveling wave on the cell surface that is pinned due to mass conservation. The cluster grows mainly in width but only slightly in height.
In addition to the POL measure from the main text which describes the relative height of the cluster, we therefore introduce a new measure for the relative cluster width as follows: where M + is the excited region of the cell surface. The excited region of the cell surface is defined as that part of the surface, where the surface concentration u is higher than the average surface concentrationū: Therefore, the POL W measure compares the size of the excited part and the total surface area of the cell.
Due to the structure of the reaction kinetics of the WP model, the concentration level of active molecules u on the cell surface separates the cell surface into two parts. The excited part M + ⊂ M cell of the cell surface is at concentration level u + and the part M − = M cell \ M + which is not excited is at concentration level u − . At steady state, the concentration in the volume is almost equally distributed and takes the value v c . The cytosolic concentration v c is in the range of parameters where f (·, v c ) has three roots u − , u T and u + , with u − < u T < u + . The critical steady state concentration v c can be calculated solving the Maxwell condition [18]: Note that the value for v c is independent of the cell size. We obtained the values v c ≈ 0.179 µM , u − ≈ 0.017 µm µM and u + ≈ 0.11 µm µM .
Due to mass conservation, in steady state the approximation holds. Solving this equation for the excited area of the cell surface |M + |, we get Using this relation, we can easily approximate the POL W measure: In the case of a circular cell, this becomes where d = 2 is the dimension of the computational representation of the cell. Therefore, the relative cluster width grows linearly with the cell size. Note that for the corresponding 1D model, the ratio |V cell | |M cell | is constant and the approximated relative cluster width in (46) becomes independent of the cell size. A comparison of the simulated and the approximated relative cluster width is presented in Figure  S6.

Figure S2
The right hand side of both models, plotted for a fixed volume concentration v. In case of the GOR model, we obtain two homogeneous steady states. The first steady state at u 1 = 0 is stable while the second at u 2 > 0 is transient; e.g. we have a membrane/cytosolic flux on the membrane for u > u 2 (region II in the plots) and a dissociation to the cytosol for u < u 2 (region I in the plots). In case of the WP model, we observe a bimodal behavior. For the shown range of values for v, we can either have one, two or three steady states. If there are three steady states, say u 1 , u 2 and u 3 with u 1 < u 2 < u 3 , the steady states u 1 and u 3 are stable while u 2 is transient. If the membrane concentration u is in the regions I or III of the plots, there is a flux from the cytosol to the membrane. If the membrane concentration u is in the regions II or IV, there is a dissociation from the membrane to the cytosol.  t=15s t=120s S1 S2 S1 S2 S1 S2 S1 S2 S1 S2 S1 S2 S1 S2 S1 S2  Figure S4 The influence of cytosolic diffusion barriers on the GOR and the WP model is investigated, see results section (D). Organelles (represented by circles or ellipses) which serve as diffusion barriers are placed at different positions. Two stimuli S1 and S2 of different strengths are applied. Stimulus S1 at the top prevails only if the organelle is sufficiently far away, i.e., if it has a bottom position. With the organelles placed close to the stronger stimulus S1, S1 is supressed and S2 dominates.   Figure S7 Excitation of the system with one weakly graded and one strongly graded signal. Both signals are persistent and defined by k S (x, t) = c · x 2 − (x min ) 2 , where x min is the lowest point along the x 2 -direction of the two dimensional cell (opposite side of the protrusion). For the weak gradient, we choose c = 1.0 · 10 −4 and for the strong gradient, we choose the ten-fold higher value c = 1.0 · 10 −3 . The signals are plotted along the surface of each cell shape (see top figures of A and B). In case of the WP model, the simulations show that the system is strongly shape-dependend. The circular cell polarizes toward the gradient as expected. On the contrary, for the cell with a protrusion, polarization is either strongly delayed (strong gradient) or occurs on the opposite side of the protrusion (weak gradient).