Turing pattern formation on the sphere is robust to the removal of a hole

The formation of buds on the cell membrane of budding yeast cells is thought to be driven by reactions and diffusion involving the protein Cdc42. These processes can be described by a coupled system of partial differential equations known as the Schnakenberg system. The Schnakenberg system is known to exhibit diffusion-driven pattern formation, thus providing a mechanism for bud formation. However, it is not known how the accumulation of bud scars on the cell membrane affect the ability of the Schnakenberg system to form patterns. We have approached this problem by modelling a bud scar on the cell membrane with a hole on the sphere. We have studied how the spectrum of the Laplace–Beltrami operator, which determines the resulting pattern, is affected by the size of the hole, and by numerically solving the Schnakenberg system on a sphere with a hole using the finite element method. Both theoretical predictions and numerical solutions show that pattern formation is robust to the introduction of a bud scar of considerable size, which lends credence to the hypothesis that bud formation is driven by diffusion-driven instability. Supplementary Information The online version contains supplementary material available at 10.1007/s00285-023-02034-z.

S1 FEM-meshes of the unit sphere with a single hole of growing radius ε = 0.3 ε = 0.4 ε = 0.5 ε = 0.6 Increasing geodesic hole radius, ε Fig. S1 FEM-meshes with an increasing geodesic hole radius ε.The meshes are presented with a growing area of the holes from left to right where the geodesic radii of the holes are ε = 0.3, ε = 0.4, ε = 0.5 and ε = 0.6, respectively.In our experimental design, we used 15 meshes with a single hole in them corresponding to the geodesic radii ε = 0, 0.05, 0.10, . . ., 0.70 where ε = 0 corresponds to the mesh without any hole, i.e., the standard unit sphere S 2 .On all of these meshes, we simulated the Schnakenberg model with homogeneous Neumann boundary conditions and stochastic initial conditions, and specifically these initial conditions corresponded to small perturbations around the steady states of the Schnakenberg model.All meshes where generated using Gmsh [1] and the FEM simulations were run in FEniCS [2,3].

S2 The eigenfunctions of the Laplace-Beltrami operator
The eigenfunctions of the Laplace-Beltrami operator on the sphere are the spherical harmonics.They come in various forms, e.g., complex and real, with and without the Condon-Shortley phase factor.Since we set out to replicate the results in [4] to validate our numerical implementation, we chose to use spherical harmonics in accordance to what was used in [4].This means that we have used the real part of the complex spherical harmonics and taken The utilized eigenfunctions of the Laplace-Beltrami operator are thus where Y n,m are the real spherical harmonics.We list some of the real spherical harmonics in Table S1 and Table S2.Also, we present a visualisation of the eigenfunctions Y m n corresponding to the eigenmodes n = 0, 1, 2, 3, 4 when they are projected onto our mesh of the unit sphere (Fig. S2).
Eigenfunctions Y m n of the Laplace-Beltrami operator Table S1 The real spherical harmonics corresponding to the eigenmodes n = 0, 1, 2, 3.

Eigenmode, n
Real sph.harm., Yn,m(x, y, z), (x, y, z) ∈ S 2 , r = S3 The perturbed eigenvalues in the cases n = 3, 4 In addition to the results presented in the article, we visualised the perturbed eigenvalues as functions of the hole radius in the cases when n = 3 and n = 4 (Fig. S3).As can be seen, the parameter choices d = 18 and γ = γ c (n), where γ c (n) is the critical reaction strength parameter for the eigenmode n, trap the perturbed eigenvalues corresponding to the particular value of n, when (a, b) = (0.20, 1.00).

S4 More results from investigating the effect of holes on pattern formation on the sphere
Here, we present more results from the simulations that were run.We used the following parameters of the Schnakenberg model: (a, b) = (0.20, 1.00), (S23) and γ was set to its critical value, i.e., γ = γ c (n) where n = 1, 2, 3, 4 is the particular eigenmode of interest.In this section, we present three types of plots for each parameter set: 1.The spectral coefficients of a single numerical solution u at time t = 50 computed on the mesh without a hole (for n = 1, 2, 3, 4), 2. The spectral coefficients for all numerical solutions u at time t = 50 computed on all meshes with a single hole of increasing radius (for n = 3, 4), 3. The quantitative properties for all numerical solutions u at time t = 50 computed on all meshes with a single hole of increasing radius (for n = 3, 4).
For each parameter set, we run the simulations on 15 meshes, and on each mesh we run 20 repetitions to account for the stochasticity in the initial conditions.In total, this means that we run: When we run the simulations on the meshes with holes with multiple repetitions, we present the 95%, 50% and 5% quantiles of the quantitative measures.
S4.1 (n, d) = (1, 20.00) Spectral coefficients, U m n Geodesic radius of the hole, ε To make sure that the end time t = 50 that has been used throughout our simulations actually resulted in a stable pattern, we run a ten times longer simulation until t = 500 (Fig. S12).As can be seen, the time scale t = 50 is indeed long enough to capture the long term behaviour of the system..00,20.00, 6.87), we run a simulations until t = 500 which is 10 times longer than all other simulations presented in this work.As can be seen, the concentration profile for these parameters is identical to the concentration profile obtained at t = 50 using the same exact parameters S4B.

Fig. S2
Fig. S2 The eigenfunctions Y m n of the Laplace-Beltrami operator.The eigenfunctions Y m n corresponding to the eigenmodes n = 0, 1, 2, 3, 4 are projected onto the mesh approximating the unit sphere.

Fig. S3
Fig. S3 Perturbed eigenvalues as a function of the geodesic radius of the hole on the sphere.The perturbed eigenvalues λ m n (ε) are plotted against the hole radius ε when n = 1, 2, 3, 4 and m = 0, 1, . . ., n.Also, the upper bound γM (a, b, d) and the lower bound γL(a, b, d) in the Turing condition involving the eigenvalues giving rise to patterns are illustrated in the dashed lines.The parameters defining these boundaries are chosen to (a, b) = (0.20, 1.00) and the value of γ is set to the critical value, i.e., γ = γc(n) for a particular eigenmode n.The upper bound γM (a, b, d) and the lower bound γL(a, b, d) are illustrated in two cases: (A) (n, γ, d) = (3, 41.24, 18.00) and (B) (n, γ, d) = (4, 68.73, 18.00).
Fig. S4 Spectral coefficients of the concentration profile of the active component at time t = 50 when (n, d) = (1, 20).The concentration profile u(x, t = 50), x ∈ S 2 resulting from the rate parameters (a, b, d, γ) = (0.20, 1.00, 20.00, 6.87), where γ = γc(n = 1), is illustrated in two ways.(A) The spectral coefficients show that these rate parameters isolate the excited eigenmodes n = 0 and n = 1.(B) The concentration profile of the active component at time t = 50 consists of one pole, i.e., high concentration region, and it is bounded between the values u min = 0.69 and umax = 1.85.

S4. 2
Fig. S5 Spectral coefficients of the concentration profile of the active component at time t = 50 when (n, d) = (2, 18).The concentration profile u(x, t = 50), x ∈ S 2 resulting from the rate parameters (a, b, d, γ) = (0.20, 1.00, 18.00, 20.62), where γ = γc(n = 2), is illustrated in two ways.(A) The spectral coefficients show that these rate parameters isolate the excited eigenmodes n = 0 and n = 2. (B) The concentration profile of the active component at time t = 50 consists of two poles, i.e. high concentration regions, and it is bounded between the values u min = 0.93 and umax = 1.70.

S4. 3
Fig. S6 Spectral coefficients of the concentration profile of the active component at time t = 50 when (n, d) = (3, 18).The concentration profile u(x, t = 50), x ∈ S 2 resulting from the rate parameters (a, b, d, γ) = (0.20, 1.00, 18.00, 41.24), where γ = γc(n = 3), is illustrated in two ways.(A) The spectral coefficients show that these rate parameters isolate the excited eigenmodes n = 0 and n = 3. (B) The concentration profile of the active component at time t = 50 consists of several poles, i.e. high concentration regions, and it is bounded between the values u min = 0.83 and umax = 1.70.

u
Fig. S12 A single simulation with (n, d) = (1, 20) until t = 500.Using the parameters (a, b, d, γ) = (0.20, 1.00, 20.00, 6.87), we run a simulations until t = 500 which is 10 times longer than all other simulations presented in this work.As can be seen, the concentration profile for these parameters is identical to the concentration profile obtained at t = 50 using the same exact parameters S4B.