A microtubule-based minimal model for spontaneous and persistent spherical cell polarity

We propose a minimal model for the spontaneous and persistent generation of polarity in a spherical cell based on dynamic microtubules and a single mobile molecular component. This component, dubbed the polarity factor, binds to microtubules nucleated from a centrosome located in the center of the cell, is subsequently delivered to the cell membrane, where it diffuses until it unbinds. The only feedback mechanism we impose is that the residence time of the microtubules at the membrane increases with the local density of the polarity factor. We show analytically that this system supports a stable unipolar symmetry-broken state for a wide range of parameters. We validate the predictions of the model by 2D particle-based simulations. Our model provides a route towards the creation of polarity in a minimal cell-like environment using a biochemical reconstitution approach.

and the MT conservation law The PF can either be adsorbed to the membrane, on which it diffuses with (angular) diffusion constant D = D b /R 2 and from which it unbinds with rate k u , freely diffuse within the interior of the cell, where we assume it diffuses with effective diffusion constant D i = ∞ or bind to MTs. Its conservation law reads where C f is the number of free PFs in the interior, C m the number bound to MTs and c b the density on the membrane. We assume that the binding kinetics of the PF to the MTs is fast with respect to the other processes, so it can be described as an equilibrium, with the density (per unit length) of PFs on MTs given by where l tot (t) = dω is the total length of all MTs in the system, and l1 2 is a parameter that sets the affinity of the PFs for binding to the MTs. The PFs bound to the MTs are transported towards the MT plus end with speed v m , where they are delivered to the membrane, if the MT they are bound to is attached to the membrane, or they simply fall off. The time evolution of the PF density on the membrane is thus described by where ∆ω is the Laplacian on the unit sphere. The final ingredient of the model we need to specify is the dependency of the residence time of the MTs at the boundary on the local density of PFs, which we parametrize as Here σ is a dose-response function, which depends on the reduced density γ ≡ c b /c * , where c * is defines the relevant density scale parameter. In principle the only properties we require of σ (γ) is to monotonically decrease from its maximal value σ (0) = 1 and to vanish sufficiently fast for large arguments, i.e. lim γ→∞ γσ (γ) = 0 (see Section V). For concreteness sake, however, we will in the following adopt a standard sigmoidal Hill type function σ (γ) = 1 1 + γ p , p > 1. (13)

II. STEADY STATE
In this section we derive the steady state solutions to the MT dynamical equations. The steady state equations for the MTs are given by with boundary conditions Combining Eqs. (14), (15), (17) and (19) for all l ∈ [0, R]. This allows the elimination of m − i (l,ω) in Eq. (14) and its solution with the aid of Eq. (16): where is the mean MT length in the absence of boundaries. We can now use the MT conservation law and the boundary condition Eq. (19) which here reads to determine This shows that in steady state all MT densities are fully determined by the density of PFs at the membrane. For convenience, we now introduce the following dimensionless quantities which respectively are proportional to the number of MTs of length < R, corresponding to the cell interior, and of length R, thus reaching the boundary, in the absence of boundary effects, i.e. dependent only on the MT dynamical parameters. We now turn to the explicit calculation of the total length of the MTs (Eq. (10)), which works out as We introduce the dimensionless quantities proportional to the mean length stored in MTs of length < R and length equal to R respectively, again in the absence of boundary effects. These definitions allow us to concisely and So, in steady state both the density c m of PFs on the MTs and the density m b (ω) of MTs at the boundary functionally depend on the density c b (ω) of PFs in the membrane. This allows us to formulate a single autonomous equation for the latter density where, combining Eq. (9) and the steady state solutions, the effective binding rate K b of PFs to membrane is given by We now note that since l t ot is proportional to the overall MT density (Eq. (31)), K b is invariant under the scaling m → αm and l1 2 → αl 1 2 . This implies that any change in overall MT density can be exactly compensated by an correspondingly increased affinity of PFs to bind to MTs, which leaves the density of PFs in the membrane, and hence the behaviour of the system as a whole, unchanged.
Clearly, Eq. (33) admits an isotropic solution c b (ω) =c b (C) for any value of the number of PFs. This solution obeys While not soluble on closed form, we can prove the property thatc b (C) is monotonically increasing with C and is asymptotically linear,c b (C) ∝ C, C → ∞ (neglecting saturation of binding), so thatc b (C) can take on any positive value (see SI Section IV), a property which we will need below.

III. BIFURCATION ANALYSIS
We now ask whether Eq. (33) can also support an anisotropic solution. To that end we adopt a standard bifurcation approach and consider the first order response of the isotropic solution to a perturbation of the type where the square of the dimensionless wavenumber is given by and As the spherical harmonics Y m n (ω) are eigenfunctions of the spherical Laplacian, they are also solutions to the bifurcation equation (36) provided Ω 2 (γ) = n (n + 1) δ. Provided there is a value of the relative membrane PF density for which Ω 2 (γ) > 0, a suitably small value of the angular displacement parameter δ can be found, such that the unipolar (n = 1) solution Y 0 1 (ω) ∝ cos θ is accessible. Given that Ω 2 (0) = −1 by definition, we thus need to ensure that Ω 2 (γ) has zeros. Explicitly so that Ω 2 (γ) = 0 leads to quadratic in γ p which has two real solutions provided Importantly, this shows that there is an upper bound condition on the parameter η. Note that η can be made arbitrarily small by both increasing the residence time of the MTs at the boundary, i.e. reducing r u (∞), while simultaneously decreasing the probability of MTs reaching the boundary, either through reducing the nucleation rate r n or decreasing the mean MT length l, which exponentially suppresses the factor µ b (cf. Eq. (27)). This analysis suggests that the smaller the parameter η, the more competitive advantage a MT stabilized at the boundary has; either because it stays much longer at the membrane allowing it to deposit yet more PFs, and/or because it has fewer MTs in other orientations competing for stabilization. Using Eq. (39) we can also determine the maximum value of Ω 2 (γ) , which is given by Unipolar solutions to the bifurcation equation can therefore only exist provided (recall that This requirement can readily be interpreted: if a PF released by a stabilized MT travels too far along the membrane before unbinding, it may contribute to the stabilization of competing MTs in other directions, thus decreasing the propensity for polarization. Finally, explicitly solving Ω 2 (γ) = 0 in the limit η ↓ 0, yields an absolute lower bound on the number of PFs necessary to achieve polarization, which is given by γ min = (p − 1) −1/p . Such a bound is expected, since a minimal density of PFs is required to achieve stabilization of MTs in the first place. Note that the bounds on all these three parameters provide an a posteriori explanation of the requirement on the Hill coefficient p > p min = 1.

IV. PROPERTIES OF THE ISOTROPIC SOLUTION
Here we derive the fact that the density of PFs on the membrane in the isotropic state increases monotonically with the total number of PFs C, which is a minimal requirement for the bifurcation to the polarized state to occur.
The explicit form for the effective binding rate in the isotropic state that appears in Eq.
(35) is with Trivially,c b (C = 0) = 0, and an explicit calculation shows that Taking the derivative of Eq. (35) with respect to the total number of PFs C and rearranging, .
We now note that, from Eq. (43) Thus d dCc b (C) can only change sign through a pole (ruled out as K b (C,c b ) is bounded) or a cusp (ruled out by smoothness of K b (C,c b )), showing thatc b (C) is monotonically increasing.
Finally, the form of Eq. (43) and the asymptotic behaviour of r u (c b ) = r u (∞) + O (c b ) −1−ε asc b → ∞, which follows from the Ansatz lim x→∞ xσ (x) = 0 on the doseresponse function, dictate thatc b (C) is asymptotically linear. An explicit calculation gives Together these results imply thatc b (C) can take on any positive value by suitably choosing C.

V. REQUIREMENTS ON THE DOSE-RESPONSE FUNCTION σ (γ)
Here we derive the minimal requirements the dose-response function σ (γ =c b /c * ) must meet in order for polarization to be possible.
A necessary requirement for the bifurcation to the polar state to occur is that Ω 2 (γ) > 0.
This implies the requirement This criterion is met (for suitably small η) if S (γ) has at least one zero for a finite value of γ. We first note that by definition S (0) = −1. Next consider By imposing the requirement I ∞ = 0, i.e. lim γ→∞ γσ (γ) = 0, we guarantee that S (γ) > 0 for some γ ∈ (0, ∞) , so that it has at least one zero. If we furthermore assume that σ (γ) has no more than a single inflection point on (0, ∞), the zero of S (γ) is moreover unique.

VI. EFFECT OF A FINITE TUBULIN POOL
In our model we have so far assumed that the amount of tubulin is not a limiting factor.
Since, however, we are working in a finite volume it is a reasonable question to ask to what extent our results are robust against possible finite tubulin pool size limitations. Here we address this question, by explicitly modelling the effect a finite pool has on the MT dynamics, specifically the growth speed and the nucleation rate. For simplicity, we disregard the effects of the capping of lengths due to cell boundary, focussing on the first order effects.

A. Dynamical equations
We assume the cell contains a finite amount of tubulin, which is represented in terms of a total microtubule length L tot . There are M nucleation sites from which MTs can be nucleated with a rate given by where L f ree is the available free (non-polymerized) tubulin length, L1 2 a cross-over parameter, which distinguishes the regime of low availability, in which the rate is strongly limited by the available free pool size, and the regime of high availability, in which the rate ultimately becomes independent of the free pool size. The Hill coefficient q > 1, describes the potential cooperativity necessary in the nucleation a new microtubule. Once nucleated, microtubules which at low availability captures the linear dependence of growth on the free tubulin density.
The other dynamical parameters v − , r + and r − are considered to be independent of the available amount of free tubulin. Denoting by M 0 (t) the total number of free nucleation sites (orientation is unimportant in this context), the dynamics of the system is described by the equations with boundary condition As a check on these equations we note that they should obey the conservation laws Now note that where we have assumed lim l→∞ m + (l, t) = lim l→∞ m − (l, t) = 0.

B. Steady state solution
We now consider the steady state, for which with boundary condition We note that Eqs. (65) and (66) lead to wherel We can determine M 0 , m 0 and L f ree through the identities so that We now take as our length scale l 0 = L1 2 /M and introduce λ = L/L1 2 , λ f ree = L f ree /L1 2 and λ (λ f ree ) =l (L f ree ) /l 0 so that and note that Below, we now prove the following two inequalities: The first inequality implies that, due to the finite tubulin pool, the MTs are on average shorter than in the saturated case. This decreases the fraction of MTs reaching the boundary, and hence decreases the parameter η (Main text, Eq. (40)), which in turn enhances the propensity to polarize. The second inequality implies that the average number of active MTs decreases due to the finite tubulin pool. However, as we show in the Main text (Section Simulations) the model can in fact be made robust against this decrease.

VII. SIMULATION STATISTICS
To calculate the order parameter S 1 for each value of the control parameter C, we performed several runs and then we averaged over them. The number of runs that was executed for each C is included in Tables I through VI.