Pom1 gradient buffering through intermolecular auto-phosphorylation

Concentration gradients provide spatial information for tissue patterning and cell organization, and their robustness under natural fluctuations is an evolutionary advantage. In rod-shaped Schizosaccharomyces pombe cells, the DYRK-family kinase Pom1 gradients control cell division timing and placement. Upon dephosphorylation by a Tea4-phosphatase complex, Pom1 associates with the plasma membrane at cell poles, where it diffuses and detaches upon auto-phosphorylation. Here, we demonstrate that Pom1 auto-phosphorylates intermolecularly, both in vitro and in vivo, which confers robustness to the gradient. Quantitative imaging reveals this robustness through two system’s properties: The Pom1 gradient amplitude is inversely correlated with its decay length and is buffered against fluctuations in Tea4 levels. A theoretical model of Pom1 gradient formation through intermolecular auto-phosphorylation predicts both properties qualitatively and quantitatively. This provides a telling example where gradient robustness through super-linear decay, a principle hypothesized a decade ago, is achieved through autocatalysis. Concentration-dependent autocatalysis may be a widely used simple feedback to buffer biological activities.

Our model of Pom1 gradient formation assumes intermolecular phosphorylation of Pom1. We also assume a large number n of phosphorylation sites (at least 6 have been documented to be relevant for gradient formation [1]). We simplify the geometrical aspect of the problem by assuming a rotational symmetry of the cell. The cell shape is well described by a capsule, i.e., a cylinder with half spheres at the poles. However, we start by approximating the cell as a cylinder, allowing us to assume that the gradient takes place on a straight line. The impact of this assumption on the profile is negligible, as documented in Section 7.
Let P (x, t) be the total amount of Pom1 at position x and time t, and P i (x, t) the amount of Pom1 phosphorylated i times. We thus have n 0 P i = P . We further assume that all P i have the same diffusion and phosphorylation constants (respectively D and β), but different detachment rates κ i that increase with the phosphorylation state i, with κ 0 = 0 [1]. We also assume, consistent with the data presented in [1], that the amount of phosphorylated Pom1 that reattaches * The author order is not the same as in the main text to highlight the contributions specific to the theoretical analyses  to the membrane is negligible. We therefore assume that Pom1 molecules are brought in non-phosphorylated state to the cortex at a rate S(x, t) (which may depend on other factors such as the Pom1 concentration in the cytosol or the availabily of the Tea4-Dis2 complex). Pom1 dynamics can then be described as follows Summing those equations yields In order to estimate the P i , we assume for a moment that diffusion is small compared to detachment and phosphorylation rates, such that we can make the following adiabatic approximation: the various phosphorylation states of Pom1 at a given position are always in quasi steady-state neglecting the effect of diffusion. We can thus represent the phosphorylation and detachment dynamics with the state transitions depicted in Fig. 1. At steady state, the flows of molecules in and out of each state are equal such that we have: It follows that and we assume that the amount of Pom1 in the last state P n is negligible as the κ i increase with i. The overall detachment rate given by the sum term in (1) is then given by The last equation can be deduced by considering that at steady state, the rate of molecules entering and exiting the meta-state defined by the union of all phosphorylated states (see rectangle in Fig. 1) are equal. 1 Now, using (3), we have Given that κ i increases with i and assuming κ 0 = 0, we can model the following relationship between κ i and i where c is a parameter quantifying how fast κ i increases with i. It can be shown (see Appendix A), that if βP κ is large enough, we have where Γ is the standard Gamma function. Inserting this into (5) leads to indicating that P 0 ∝ P η , with η ≤ 1. Due the mixing caused by the actual diffusion, the relative amounts between P 0 and P tend to homogenize across the profile. A perfect mixing would indeed imply that P 0 ∝ P , that is η = 1. Hence, c/(c + 1) must be considered a lower bound on the effective value of η. We can thus rewrite (1) as where 1 < γ = 1 + η ≤ 2 and α = βP 0 /P η can be seen as the effective detachment rate. 1 It can also be shown more formally by the following formula n i=0 a i i j=1 , which can be easily proven by induction on n and considering the product on the RHS being negligible for large enough n. We thank anonymous contributors to the math.stackexchange forum for pointing this out.  Figure 2: Convergence of P 0 to a power of P in (8) as the ratio βP/κ increases for various levels of cooperativity indicated by c in (6). If βP is at least somewhat bigger than κ the approximation P 0 ≈ P c c+1 is accurate.
A simple model where Pom1 detachment is caused by the repulsive Colomb forces between phosphate groups and the cell membrane (both negatively charged), would lead to c ≥ 1 or η ≥ 1/2 as the force increase linearly with the number of phosphorylated sites. Any kind of cooperativity between the phosphorylated sites for the membrane detachment would increase c and bring η closer to 1.

Gradient shape
In this section we describe an analytical solution to the simplified model given by (9), assuming that the source is localized at the cell pole. Discarding the domain of Pom1 attachment around the cell pole, the gradient is given by the solution of the following equation at steady-state Making the Ansatz that P (x) = (ax + b) c , defining A as the amplitude of the gradient at the pole A = P (0) and imposing lim x→∞ P (x) = 0 because the cell length is large compared to the decay length, we get to the solution (see also [2]) with If we assume γ = 2, this becomes where x 0 can be thought of the natural length scale of the gradient (1.97 microns on average in our data).

Source quantification
We now turn to the source term S that was left unspecified in the above sections. It describes the attachment of non-phosphorylated Pom1 to the membrane at the cell poles. It is believed that the Tea4-Dis2 phosphatase complex is actively brought to the cell poles by the microtubules. It binds and de-phosphorylates cytoplasmic Pom1, which then binds to the membrane at the cell pole. As a first approximation, S can thus be assumed to be proportional to both the Tea4 concentration at the pole and the cytoplasmic Pom1 concentration.
Since the cell is a closed system and we assume a steady state gradient, the amount of Pom1 that attaches to the membrane from the cytoplasm is equal to the overall amount of Pom1 that detaches from the membrane into the cytoplasm. In other words, the rates of total attachment and detachment over the whole profile are equal. The total attachment S is thus given by If the number of Pom1 molecules attaching to the membrane is proportional to the Tea4 concentration, this relationship predicts the following power law between Tea4 and Pom1 at the pole: For γ = 2, this implies a 2/3 power law. Precisely this has been observed experimentally (see Figure 1 of the main text), confirming that γ is equal or close to two.

Buffering
Knowing the shape of the gradient profile, we can now quantify the buffering such a mechanism provides. Defining the decay length λ as the position at which the profile divides its amplitude by e = 2.7 . . . , we have: leading to The model thus predicts that the decay length decreases with the inverse square root of the amplitude at the pole. This has also been verified experimentally (see Figure 1 of the main text).
Following [3], we can compute the effect dx of a small variation dS in the source term S on the positional information Inserting (18) into (9) and defining a = 1 Local sensitivity to S is thus given by This means that the absolute imprecision on the position x is independent of x (as can be expected from a diffusive gradient) and increases slower than the relative imprecision on the Pom1 attachment rate ∆S S . The global sensitivity can then be computed by considering the difference in positional information ∆x induced by multiplying the source term S 0 by a factor k > 1.
This shows that x 0 constitutes an upper bound on ∆x. Moreover, it implies that an overall higher attachment rate will translate into a more precise positional information, as x 0 can be reduced by bringing more Pom1 to the pole. As expected and illustrated in Fig. 3, ∆x is smaller than for a non-buffered exponential profile, given by (1 − √ e)log k [4]. A realistic upper value of k = 8 corresponds to a positional shift of about half of the length scale, which fits well with the observed variability of the Cdr2 domain boundary which has a standard deviation of about 30% of the largest length scale [5].
We can also consider two profiles of different amplitudes at the pole A 1 and A 2 . Assuming γ = 2, at a distance x L away from the pole the ratio between the concentrations of the two profiles P 1 and P 2 is given by This indicates that high amplitudes gradients are buffered at a shorter range than low amplitude gradients and that reducing the diffusion constant (for example with a cluster formation mechanism) increases the phosphorylation-driven buffering.

Total Pom1 in gradient
The total amount of Pom1 in a gradient can be computed as the area under the gradient curve. This is given by If the gradient extends to a large enough range such that it can be approximated by infinity, a square power law between the total amount of Pom1 and the amplitude at the pole is expected and consistent with our experimental data (see Supplemental Fig. S1F). Simulations indicate that such an approximation can be made only for ranges bigger than λ (or x 0 ) by at least an order of magnitude. For smaller ranges the power law is between 1 and 2.

Numerical validation
In addition to experimental evidence, numerical simulations were performed to verify that violations of the assumptions described above do not hamper the validity of the model. Using the standard Matlab numerical solver and assuming six phosphorylation sites, the gradients were generated using (1), for various values of α/κ. The diffusion constant D is not a free parameter of the model, and was set to one. Using a Gaussian Adaptation optimizer [6], the α parameter was then optimized to match the simple trans-phosphorylation model, i.e., (9) with γ = 2. The results, shown in Fig 4 indicate that the simple model provides a good approximation of the detailed model, for the tested gradient shapes.

Cell geometry
In the above model, the cell was assumed to be a cylinder to allow an analytical solutions to the model. We show here that this assumption has a minimal impact on the results. In (9), the only term that is affected by the cell geometry is the diffusion term D ∂ 2 P ∂x 2 . We thus consider diffusion at the cell tip as diffusion on a sphere. The diffusion term is given by the divergence of the flux J of P . Considering the standard spherical coordinates defined by (r, θ, φ), we have where J θ = D∂P r∂θ , and J r = J φ = 0 because the diffusion is limited to the cortex and because of rotational symmetry. Hence  Table 1: Various hypotheses for a negative correlation between Pom1 gradient amplitude A and decay length λ. The top three rows correspond to cases of linear detachment and the last row corresponds to our supra-linear detachment hypothesis. Unlike the total cortical Pom1, which varies across cells, the diffusivity and affinity of a Pom1 molecule to the membraine are very likely to be constant across cells. The top three hypotheses can thus be safely rejected. Furthermore, as shown in line 3, the coordination between the varying parameters required to account for a -1/2 power law seems not realistic. because x = rθ. As can shown in Fig. 5, using this last expression for the diffusion in the spherical part of the cortex produces very similar profiles than the ones obtained with the cylindrical approximation. In those simulations and consistently with the experimental observations, a Gaussian source was considered with a standard deviation set equal to the radius of the sphere. This confirms that our cylindrical approximation, which is necessary for the derivation of the analytical solutions, is justified. As for the point source approximation, it has been previously shown to only impact the solution close to the source [7].

Aternative hypothesis
It could be argued that the decreasing relationship between the Pom1 gradient amplitude and its decay length is not due to a buffering mechanism but rather to a very tight control of cortical Pom1 that ensures it remains constant across cells.
Assuming a linear decay, we have: where S is likely to depend on the cytoplasmic Pom1 concentration P c . This assumes that the production and degradation of Pom1 are negligible and Pom1 only cycles between the cortex and the cytoplasm. Assuming a point source S(x, t) = Sδ(x) and solving (33) for steady state results in the classical exponential gradient The relationship between A and S can be found either by considering (34) at steady state or equivalently by integrating (33) over x at steady state. The result for a point source is: where P tot is the total cortical Pom1 over the profile. A reasonable hypothesis for S would assume that it is proportional to both the Tea4 concentration T at the pole and the cytoplasmic Pom1 S = kP c T . However, our data (see Supplementary Figure S1) indicate that A is highly correlated with T but not with P c , maybe because variations in T are much larger than variations in P c or the reaction takes place at P c saturation level. So we may remain general and write that S = k(P c )T and restrain from any assumption on the potential role of P c in the Pom1 attachment rate. Table 1 first lists three hypotheses explaining the observed buffering within the framework of a linear detachment and homogeneous diffusion. Those hypotheses are very unlikely as they assume that the variability in gradient shapes originate in the variability in Pom1 diffusion or detachment rate while maintaining a tight control of the total cortical Pom1 across cells. This is not consistent with our data, which display a large variability in total cortical Pom1 across cells. On the contrary our hypothesis assumes constant Pom1 properties across cells and sees the Pom1 attachment rate to the membrane as the main source of variability across cell (hypothesis 4 in Table 1). This is very likely to be the case given the known stochasticity of microtubule dynamics even within the same cell.
where Γ is the standard Gamma function.
Proof Let us define a = u c+1 and call the above sum s. Then we get The main idea of the proof is the observation that the last infinite series above behaves like a Riemann sum as u → ∞. To this end, we fix M > 0 and divide the sum into two parts s 1 and s 2 as defined below:  Taking M → ∞ we finally obtain that f (c) = lim s exists and is equal to ).