Stability Analysis of a Bulk-Surface Reaction Model for Membrane-Protein Clustering

TL;DR

This paper investigates the stability and pattern formation of protein clusters on cell membranes using a reaction-diffusion model, revealing how ligand recruitment can lead to stable, heterogeneous protein distributions.

Contribution

It provides a theoretical analysis of membrane protein clustering through a reaction-diffusion model, including diffusion-driven instability estimates and a threshold for pattern formation.

Findings

Ligand recruitment promotes stable protein cluster formation.

A threshold exists for ligand recruitment to induce pattern formation.

Single-patch steady states emerge under certain conditions.

Abstract

Protein aggregation on the plasma membrane (PM) is of critical importance to many cellular processes such as cell adhesion, endocytosis, fibrillar conformation, and vesicle transport. Lateral diffusion of protein aggregates or clusters on the surface of the PM plays an important role in governing their heterogeneous surface distribution. However, the stability behavior of the surface distribution of protein aggregates remains poorly understood. Therefore, understanding the spatial patterns that can emerge on the PM solely through protein-protein interaction, lateral diffusion, and feedback is an important step towards a complete description of the mechanisms behind protein clustering on the cell surface. In this work, we investigate the pattern formation of a reaction-diffusion model that describes the dynamics of a system of ligand-receptor complexes. The purely diffusive ligand in the cytosol can bind receptors in the PM, and the resultant ligand-receptor complexes not only diffuse laterally but can also form clusters resulting in different oligomers. Finally, the largest oligomers recruit ligands from the cytosol in a positive feedback. From a methodological viewpoint, we provide theoretical estimates for diffusion-driven instabilities of the protein aggregates based on the Turing mechanism. Our main result is a threshold phenomenon, in which a sufficiently high recruitment of ligands promotes the input of new monomeric components and consequently drives the formation of a single-patch spatially heterogeneous steady-state.

Equations

$$\mathrm{ch}{\text{A}}_{j-1}+{\text{A}}_1<=>{\text{A}}_j\text{for}j=2,3,...,N.$$

$$f(u,a_1,a_N)=(k_0+k_b{a}_N)u-k_d{a}_1$$

$${\partial }_{t}u$$

$${\partial }_t{a}_1$$

$$-k_g{a}_1\left(\sum _{l=2}^{N-1}{a}_l\right)+\sum _{j=3}^N{k}_j{a}_j$$

$${\partial }_t{a}_2$$

$${\partial }_t{a}_j$$

$${\partial }_t{a}_N$$

$$-D_u\left(\mathbf{n}\cdot \nabla u\right)=(k_0+k_b{a}_N)u-k_d{a}_1.$$

$$n_U+n_{A_1}+2n_{A_2}+...+N{n}_{A_N}$$

$$M(t):={\int }_{\Omega }u(x,t)dx+\sum _{j=1}^N\left\{j\cdot {\int }_{\Gamma }{a}_j(x,t)ds\right\}$$

$$\frac{d}{dt}M$$

$${\int }_{\Omega }{\partial }_{t}udx$$

$$=D_u{\int }_{\Gamma }(\nabla u\cdot \text{n})ds$$

$$=-{\int }_{\Gamma }[(k_0+k_b{a}_N)u-k_d{a}_1]ds$$

$$\sum _{j=1}^N\left\{j\cdot {\int }_{\Gamma }{\partial }_t{a}_{j}ds\right\}$$

$$={\int }_{\Gamma }\left[D_1\Delta a_1+(k_0+k_b{a}_N)u-k_d{a}_1-2k_m{a}_1^2+2k_2{a}_2-k_g{a}_1\left(\sum _{l=2}^{N-1}{a}_l\right)+\sum _{j=3}^N{k}_j{a}_j\right]ds$$

$$+{\int }_{\Gamma }\left[2\cdot D_2\Delta a_2+2\cdot \left\{k_m{a}_1^2-k_g{a}_1{a}_2-k_2{a}_2+k_3{a}_3\right\}\right]ds$$

$$+\sum _{j=3}^{N-1}{\int }_{\Gamma }\left[j\cdot D_j\Delta a_j+j\cdot \left\{k_g{a}_1{a}_{j-1}-k_g{a}_1{a}_j-k_j{a}_j+k_{j+1}{a}_{j+1}\right\}\right]ds$$

$$+{\int }_{\Gamma }N\cdot \left[D_N\Delta a_N]+N\cdot \left\{k_g{a}_1{a}_{N-1}-k_N{a}_N\right\}\right]ds$$

$$=\sum _{j=1}^{N}j{D}_j\cdot {\int }_{\Gamma }\Delta a_{j}ds+{\int }_{\Gamma }(k_0+k_b{a}_N)u-k_d{a}_{1}ds$$

$$={\int }_{\Gamma }\left((k_0+k_b{a}_N)u-k_d{a}_1\right)ds$$

$$\frac{d}{dt}M={\int }_{\Omega }{\partial }_{t}udx+\sum _{j=1}^N\left\{j\cdot {\int }_{\Gamma }{\partial }_t{a}_{j}ds\right\}=0,$$

$$\hat{u}=\frac{u}{U},{\hat{a}}_j=\frac{a_j}{A_j}(j=1,\dots ,N),\hat{t}=\frac{t}{T},\text{and}\hat{x}=\frac{x}{L},$$

$$\nabla =\frac{1}{L}\hat{\nabla },{\nabla }^2=\frac{1}{L^2}{\hat{\nabla }}^2,\Delta =\frac{1}{L^2}\hat{\Delta },$$

$$\frac{U}{T}\frac{\partial \hat{u}}{\partial \hat{t}}$$

$$\frac{A_1}{T}\frac{\partial {\hat{a}}_1}{\partial \hat{t}}$$

$$\displaystyle\qquad-\hat{a}_{1}\left(\sum_{l=2}^{N-1}\frac{A_{j}k_{g}}{k_{d}}\hat{a}_{j}\right)+\sum_{j=3}^{N}\frac{k_{j}A_{j}}{k_{d}A_{1}}\hat{a}_{j}\Biggl{\}},\quad\hat{x}\in\hat{\Gamma},$$

$$\frac{A_2}{T}\frac{\partial {\hat{a}}_2}{\partial \hat{t}}$$

$$\frac{A_j}{T}\frac{\partial {\hat{a}}_j}{\partial \hat{t}}$$

$$\hat{x}\in \hat{\Gamma },j=3,\dots ,N$$

$$\frac{A_N}{T}\frac{\partial {\hat{a}}_N}{\partial \hat{t}}$$

$$-\frac{D_{u}U}{L}\left(\mathbf{n}\cdot \hat{\nabla }\hat{u}\right)=k_d{A}_1\left\{\left(\frac{k_{0}U}{k_d{A}_1}+\frac{k_b{A}_N}{k_d{A}_1}\widehat{a_N}\right)\hat{u}-\widehat{a_1}\right\}$$

$$U=\frac{M_0}{R\cdot |\Gamma |},A_j=\frac{M_0}{|\Gamma |}\text{for}j=1,2,\dots ,N,T=\frac{R^2}{D_1},L=R$$

$${\hat{k}}_0=\frac{k_{0}U}{k_d{A}_1},{\hat{k}}_b=\frac{k_b{A}_{N}U}{k_d{A}_1},{\hat{k}}_m=\frac{k_m{A}_1}{k_d},{\hat{k}}_j=\frac{k_j{A}_j}{k_d{A}_1}(j=2,\dots ,N),$$

$${\hat{k}}_g=\frac{k_g{A}_j}{k_d}(j=2,\dots ,N-1),,\gamma =\frac{k_d{R}^2}{D_1},\tilde{D}=\frac{D_u}{D_1},d_j=\frac{D_j}{D_1}(j=2,\dots ,N).$$

$$\frac{\partial \hat{u}}{\partial \hat{t}}=\tilde{D}{\hat{\nabla }}^2\hat{u},$$

$$-\tilde{D}\left(\mathbf{n}\cdot \hat{\nabla }\hat{u}\right)=\gamma \left\{\left[{\hat{k}}_0+{\hat{k}}_b{\hat{a}}_N\right]\hat{u}-{\hat{a}}_1\right\}.$$

$$\frac{\partial {\hat{a}}_1}{\partial \hat{t}}$$

$$-{\hat{k}}_g{\hat{a}}_1\left(\sum _{l=2}^{N-1}{\hat{a}}_j\right)+\sum _{j=3}^N{\hat{k}}_j{\hat{a}}_j\},$$

$$\frac{\partial {\hat{a}}_2}{\partial \hat{t}}$$

Full text

1Department of Mechanical and Aerospace Engineering, University of California, San Diego, La Jolla CA 92093-0411

2Scripps Institution of Oceanography, University of California, San Diego, La Jolla, CA 92093-0213

3Department of Mathematics, University of California, San Diego, La Jolla, CA 92093-0112

4Department of Physics, University of California, San Diego, La Jolla, CA 92093-0424

∗ Both these authors contributed equally

*∗∗*To whom correspondence should be addressed. e-mail: prangamani@ucsd.edu

1. 1 Introduction

2. 2 Model Development

3. 2.1 Assumptions

4. 2.2 Governing equations

5. 2.3 Mass conservation

6. 2.4 Non-dimensionalization

7. 2.5 System Reduction when $D_{u}\to\infty$

8. 3 Mathematical Analysis

9. 3.1 Homogeneous steady-states

10. 3.2 Linear Stability Analysis

11. 3.3 Special Case $N=2$: Necessary Conditions for Diffusion-Driven Instabilities

12. 3.4 Existence of Solutions to the General and Reduced Models

13. 4 Numerical Simulations

14. 4.1 Bistability under homogeneous perturbations

15. 4.2 Linear instability under non-homogeneous perturbations

16. 4.3 The emergence of the single-patch non-homogeneous steady-state

17. 4.4 Temporal evolution and pattern formation

18. 4.5 Change of the cell radius and single-patch area

19. 5 Discussion

20. 6 Acknowledgments

21. 7 Electronic Supplementary Material

1. Introduction

Biological membranes are unique two-dimensional structures that separate cellular contents from the extracellular environment and regulate the transport of material into and out of the cell [1, 2]. In addition to lipids and carbohydrates, these membranes contain a large proportion of proteins, the composition of which depends on the cell type [2, 3, 4, 5, 6]. One of the interesting features of membrane proteins is their ability to form clusters on the cell surface [7, 8, 9]. This clustering of proteins on the plasma membrane (PM) results in a spatial heterogeneity in the distribution of protein densities. Many factors can induce such a spatial heterogeneity, including lateral diffusion, physical barriers from the cytoskeleton [10], lipid raft affinity [11], and curvature differences along the membrane[8]. The formation of protein clusters is intimately related to various cellular phenomena such as polarization, membrane depolarization, receptor signaling, enzyme activity, and cytoskeletal regulation [12, 13, 14, 15, 16, 17].

A particular example of proteins forming clusters on the membrane is well-elucidated by amyloid-$\beta$ aggregation/fibrillation in the context of Alzheimer’s disease. It is thought that amyloid-$\beta$ can become cytotoxic is when it aggregates on the membrane at high levels [18]. Biophysical measurements show that amyloid-$\beta$ aggregates become more stable when oligomerized on the membrane surface [19, 20] and also can destabilize certain membrane compositions [21]. It is also thought that membrane components such as cholesterol may initiate aggregation of amyloid-$\beta$, which may then be bolstered by as a yet-unidentified secondary feedback mechanism [22].In addition to amyloid-$\beta$, surface receptors such as $\alpha$-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptor (AMPAR) [23, 24] and membrane-bound kinases such as Fyn [25] are known to cluster on the membrane; these clusters have been implicated in neuronal functioning in physiology [13, 11] and disease [21, 18].

One of the open questions in the field of protein aggregations is the role of the spatial organization of membrane proteins due to bulk-surface reactions and feedback mechanisms. Mathematical modeling has provided substantial insight into the geometric coupling of bulk-surface reaction-diffusion systems [26, 27, 28], including wave-pinning formulations [12, 29], spatial patterning [30, 31], and generalized stability analysis [32, 33, 34]. From a modeling perspective, several authors have proposed the classical Smoluchowski coagulation model [35, 36] as a suitable candidate for describing protein aggregation. The use of different aggregation kernels [37, 38] allowed a successful combination of experimental measurements with computational predictions. These models performed well in terms of comparisons to data and estimation of kinetic parameters such as the aggregation time and the asymptotic cluster distribution. However, by using the original Smoluchowski systems of ordinary differential equations (ODEs), these studies lack descriptions of the spatial protein organization, which can be crucial for the understanding of many cellular processes. To overcome this limitation, one can explicitly consider molecular diffusion and use systems of partial differential equations (PDEs) as has been done in the amyloid-$\beta$ aggregation models [39, 40, 41]. These studies have provided detailed theoretical estimates in terms of boundary conditions and homogenization tools. However, they have restricted the spatial scale to a small three-dimensional region of cerebral tissue and not described intracellular phenomena. Therefore, there is a need for mathematical models for protein aggregation in the PM with a proper spatial description to account for the numerous cellular processes that occur due to heterogeneous protein distribution.

In this work, our primary goal was to investigate the emergence of spatially heterogeneous steady-state profiles of membrane protein aggregates to identify how feedback between cytosolic and membrane components can drive pattern formation on the membrane. To this end, we merged the concept of bulk-surface reaction-diffusion systems with the Smoluchowski approach to introduce a new bulk-surface model for membrane protein clustering (Figure 1). The model equations describe a purely diffusive ligand in the cytosol which then undergoes membrane binding, without any cytosolic aggregation. The resultant membrane-bound protein can diffuse laterally and also form clusters with different oligomeric sizes. Finally, the oligomers of maximum size can further recruit more cytosolic proteins, resulting in a positive feedback for the membrane protein aggregates and stabilization of the oligomers [22, 19]. Following the approach of Ratz and Roger [33, 34], we then analyzed the model for diffusion-driven instabilities using the classical Turing mechanisms. We found these interactions allow diffusion-driven instabilities and pattern formation in the absence of a sustained localized stimulus.

In what follows, we present the model assumptions and derivation in Section 2, the mathematical analysis including stability analysis in Section 3, and conclude with numerical simulations (Section 4) and a discussion (Section 5) about our findings in the context of amyloid-$\beta$ and cluster of other membrane proteins.

2. Model Development

Here we present our bulk-surface reaction-diffusion model for protein aggregation, including feedback. We describe our assumptions (Section 2.1) and the governing equations (Section 2.2) in detail. In Section 2.3, we prove that the total mass of the system is conserved over time, and in Section 2.4, we non-dimensionalize the model. Finally, in section Section 2.5, we perform the systems reduction when the cytosolic diffusion goes to infinity, following the mathematical approach of Ratz and Roger [33, 34].

2.1. Assumptions

In our system, we assume that U represents the volume component, which can freely diffuse in the cytoplasm. Upon binding to the plasma membrane, it forms a surface monomer component $\textbf{A}_{1}$. The $\textbf{A}_{1}$ molecules laterally diffuse in the membrane and form the oligomeric components $\textbf{A}_{j}$. Here, $j$ denotes the number of $\textbf{A}_{1}$ molecules in the oligomer, which is at most $N\in\mathbb{N}$. In terms of chemical reactions, \ch{\textbf{U}$<=>[$f$]$\textbf{A}_{1}} denotes the binding of the cytosolic component to the plasma membrane with a reaction flux $f$. The subsequent oligomerization at the membrane is described by

[TABLE]

We also assume that the flux term $f$ describes ligand binding/unbinding to the cell surface, where the binding term will be linearly proportional to the concentrations of $U$ in the cytosol and $A_{N}$ in the plasma membrane. The oligomerization process is modeled as a particular version of the reversible Smoluchowski model for aggregation dynamics [42]. We also assume that the oligomerization process occurs only by monomer attachment in the mass action regime. Moreover, to keep the analysis tractable, we do not consider any cooperativity term such as Hill’s function [43]: the rate at which the different oligomers are formed is independent of their size.

2.2. Governing equations

We represent the cellular domain as the bounded region $\Omega$ with smooth boundary $\Gamma=\partial\Omega$. We define the concentrations $u(x,t):\Omega\times(0,\mathcal{T}]\to\mathbb{R}$ for the volume component and $a_{j}(x,t):\Gamma\times(0,\mathcal{T}]\to\mathbb{R}$ for the membrane oligomeric components. The molecular mechanisms underlying membrane protein aggregation and stabilization are quite complex. Therefore, we propose a mathematically tractable feedback term to represent this complex mechanism. The flux term is thus defined as

[TABLE]

for $k_{0}$, $k_{b}$ and $k_{d}$ positive constants, where $k_{0}$ is the basal binding rate, $k_{b}$ is the rate of $A_{N}$-dependent binding rate, and $k_{d}$ is the unbinding rate from the membrane into the cytosol. Then the governing equations for the spatiotemporal evolution of the different components are given by

[TABLE]

Here, $\nabla^{2}$ and $\Delta$ represent the Laplace and Laplace-Beltrami operators, respectively. The parameter $k_{m}$ represents the rate at which monomers bind to form dimers. The rate $k_{g}$ at which the oligomers of size greater than two are formed is assumed to be the same for all oligomerization reactions. Finally, $k_{j}$ represent the rates at which the oligomeric components of size $j$ will release a single monomer. The boundary condition for $a$ is periodic since the domain is closed and the boundary condition for $u$ is given by

[TABLE]

as a balance of the diffusive flux and the reaction rate at the membrane. All parameters and variables are non-negative real numbers.

2.3. Mass conservation

Let $n_{X}$ denote the number of molecules of the component $X$. For a closed system, we know that the total number of single molecules must be given by

[TABLE]

since each $A_{j}$ oligomer must have exactly $j$ molecules of $A_{1}$. From this fact, we define the total mass of the system, which accounts for spatial compartments (bulk and surface) and the different molecular size distributions. This is the content of the following.

Proposition 2.1**.**

Let $u$, $a_{1}$, $a_{2}$,…, $a_{N}$ be solutions of (2.2)–(2.7). Then the quantity

[TABLE]

represents the total mass of the system and is conserved over time, i.e, $M(t)=M_{0}\quad\forall t\geq 0$. In this case, $M_{0}$ denotes the initial mass which is given by $M_{0}=\int_{\Omega}u(x,0)dx+\sum_{j=1}^{N}\left\{j\cdot\int_{\Gamma}a_{j}(x,0)ds\right\}$

Proof.

By taking the time derivative of $M(t)$, and assuming $u$ and $a_{j}$ are $\mathcal{C}^{2}$ solutions for (2.2)–(2.7), we have

[TABLE]

For the integral $\int_{\Omega}\partial_{t}udx$ , we apply the divergence theorem and substitute Eq. (2.7) to obtain

[TABLE]

For the summation of surface integrals $\sum_{j=1}^{N}\left\{j\cdot\int_{\Gamma}\partial_{t}a_{j}ds\right\}$, we substitute the governing equations to obtain

[TABLE]

where the last equality comes from the fact that $\int_{\Gamma}\Delta a_{j}ds=0$ as a consequence of the First Green’s Theorem [44]. We therefore have

[TABLE]

from which we conclude that $M(t)=M(0)=:M_{0}$ for all $t\geq 0$ ∎

The mass conservation property for bulk-surface reaction-diffusion models has been established in different contexts [29, 33]. However, to the best of our knowledge, it has never been identified in the context of oligomerization reactions.

2.4. Non-dimensionalization

We introduce a non-dimensional version of the system that allow convenient qualitative interpretation of the system independent of the actual system size, but instead through the ratio of kinetic parameters to the diffusion contributions. We follow the approach in [33, 34] and define $U,A_{1},A_{2},...,A_{N}$ be the dimensional concentration quantities where $[U]=\text{mol}/\mu m^{3}$ and $[A]=\text{mol}/\mu m^{2}$ for $j=1,\ldots,N$. We also define $L$ and $T$ as the spatial and temporal quantities, where $[L]=\mu m$ and $[T]=s$. We then introduce the non-dimensional variables

[TABLE]

which lead to the transformed domains $\hat{\Omega}:=\{\xi\in\mathbb{R}^{3}|\xi L\in\Omega\}$ and $\hat{\Gamma}=\partial\hat{\Omega}$. By denoting $\hat{\nabla},\hat{\nabla}^{2}$, and $\hat{\Delta}$ as the dimensionless gradient, Laplace, and Laplace-Beltrami operators, respectively, and using

[TABLE]

we can apply the chain rule and rewrite the system (2.2)–(2.6) in the form

[TABLE]

The boundary conditions in (2.7) can be rewritten as

[TABLE]

Since $R>0$, we can define the characteristic concentrations $U$ and $A_{j}$ by dividing the total mass of the system per total volume and surface area, respectively. We also define the characteristic time with respect to the diffusion $D_{1}$ of the monomeric component across the cellular surface. Formally, we define

[TABLE]

and the dimensionless parameters

[TABLE]

As a result, (2.9) can be written as

[TABLE]

for $\hat{x}\in\hat{\Omega}$ with boundary condition

[TABLE]

for $\hat{x}\in\hat{\Gamma}$. Finally, for the surface components, (2.10)–(2.13) can be written as

[TABLE]

2.5. System Reduction when $D_{u}\to\infty$

We further reduce our system by assuming that the cytosolic diffusion coefficient is much larger than the lateral diffusion on the membrane as has been experimentally observed [45]. The resulting system is uniquely defined on the membrane surface, and the bulk variable $u$ will be represented by an integral operator also called a non-local functional. Our approach closely follows the work of Ratz and Roger [33, 34], though our system can be $N$-dimensional in principle. Formally, if we assume $D_{u}\to\infty$ and if the initial concentration of $\hat{u}$ is constant over $\hat{\Omega}$, then $\hat{u}$ no longer depends on space and $u=u(t)$. Therefore, the mass conservation law given by (2.8) implies

[TABLE]

where $\mathcal{M}_{0}=\hat{u}(0)|\hat{\Omega}|+\sum_{j=1}^{N}\left\{j\cdot\int_{\hat{\Gamma}}\hat{a}_{j}(s,0)ds\right\}$ is the total mass of the dimensionless system. We then define the non-local functional

[TABLE]

as in [33, 34]. Finally, we drop all the hats to obtain the reduced system

[TABLE]

where

[TABLE]

In the next sections, we will provide analytical estimates and numerical simulations to analyze the stability properties of the reduced system (2.23)–(2.24).

3. Mathematical Analysis

In this section, we present the mathematical framework for investigating diffusion-driven-instabilities in the system (2.23)–(2.24). We establish conditions that guarantee the existence and uniqueness of homogeneous steady-states, or the conditions for having multiple steady-states. We also present a characterization for the Jacobian Matrix in the case of homogeneous perturbations. For the non-homogeneous case, the linearization of the non-local functional yields a different Jacobian matrix, and a family of ordinary differential equations is derived to analyze the stability in terms of the eigenfunctions of the Laplace-Beltrami operator. We then apply our framework in the case $N=2$, where we obtain a necessary condition for diffusion-driven instabilities. Finally, we make some remarks about basic questions concerning the solution theory for the reduced system (2.23)–(2.24), as well as the more general model (2.2)–(2.7).

We start with the characterization of the homogeneous steady-states.

3.1. Homogeneous steady-states

The homogeneous solutions of (2.23)–(2.24) satisfy the ODE system

[TABLE]

and the steady-states in this case are given by $\textbf{a}^{*}=(a^{*}_{1},a^{*}_{2},a^{*}_{3},\dots,a^{*}_{N})$ such that

[TABLE]

for all $j=1,\ldots,N$. From $\mathcal{F}_{N}(\textbf{a}^{*})=0$, we obtain $a^{*}_{N}=\frac{k_{g}a^{*}_{1}a^{*}_{N-1}}{k_{N}}$ and proceeding recursively, it is easy to show that

[TABLE]

Hence $a^{*}_{j}=C_{j}(a^{*}_{1})^{j}$ where $C_{1}=1$ and

[TABLE]

Thus from $\mathcal{F}_{1}(\textbf{a}^{*})=0$, we must have

[TABLE]

By multiplying both sides by $|\Omega|$ and rearranging the $(a^{*}_{1})^{j}$ terms, we can define the polynomial

[TABLE]

where the roots of $\mathcal{P}_{N}$ are the steady-state values $a^{*}_{1}$. We then observe that the coefficient of $\alpha^{N}$ is a non-negative number if and only if

[TABLE]

which in this case implies that $\mathcal{P}_{N}(\alpha)$ has a unique positive root and therefore that the system has a unique steady-state. This is the case when $k_{b}=0$, which means that the largest oligomers do not promote ligand binding in the plasma membrane. On the other hand, if $k_{0}|\Gamma|N-\mathcal{M}_{0}k_{b}<0$, then multiple steady-states could exist.

3.2. Linear Stability Analysis

Linear stability is a traditional concept from the theory of dynamical systems that treats the study of the local behavior near a steady-state solution. The term “linear” stands for the analysis of the linear approximation of a nonlinear system, which can be sufficient to determine if a steady-state is stable or unstable. In the case of a system of ODEs, the analysis is carried out by evaluating the eigenvalues of the so-called Jacobian matrix. A similar analysis can be carried in the context of reaction-diffusion systems of PDEs with the analysis of the eigenvectors of the Laplace operator. A major contribution in this field is due to Alan Turing in the classic paper “The Chemical Basis of Morphogenesis” [46]. Turing established the notion of diffusion-driven instabilities and was the first to connect this mathematical idea with the formation of spatially heterogeneous patterns. In what follows, we first analyze the homogeneous perturbations of the steady-states by describing the Jacobian matrix of the system. Then we apply the Turing framework and define the conditions for diffusion-driven instabilities in our system (2.23)–(2.24).

3.2.1. Homogeneous perturbations

In this section, we investigate the linear stability of the steady-states $\textbf{a}^{*}$ against spatially homogeneous perturbations, that is in the absence of diffusion. Our study is an N-dimensional version of the approach taken in [33, 34] for a GTPase cylcling model. We need to compute the eigenvalues $\lambda$ of the Jacobian matrix

[TABLE]

for $\mathcal{F}_{j}$ defined in (2.23) and (2.24). If all the eigenvalues of $\mathcal{J}[\textbf{a}^{*}]$ have negative real parts, then the steady-state is called linearly stable [47]. That means that local perturbations will converge to the steady-state.

On the other hand, if at least one of the eigenvalues has a positive real part, then it is called linearly unstable, which that local perturbations will lead the system away from the steady-state. The next proposition generally characterizes $\mathcal{J}[\textbf{a}^{*}]-\lambda\textbf{I}$.

Proposition 3.1**.**

The matrix $\mathcal{J}[\textbf{a}^{*}]-\lambda\textbf{I}$ can be written in the form

[TABLE]

where $w_{0}$ and $\lambda$ are real numbers, $\textbf{w}\in\mathbb{R}^{N-1}$ is a row vector, $\textbf{v}\in\mathbb{R}^{N-1}$ is a column vector, and $H$ is a $(N-1)\times(N-1)$ tridiagonal matrix.

Proof.

We will first calculate $\frac{\partial\mathcal{F}^{*}_{j}}{\partial a_{i}}$ for $i,j=1,2,\ldots,N$ . For $j=1$ we obtain

[TABLE]

Now for $j=2$, we have

[TABLE]

and for $j=3$ to $j=N-1$, we obtain

[TABLE]

and

[TABLE]

otherwise, and finally for $j=N$,

[TABLE]

We then define $\mathcal{J}^{*}_{ij}:=\gamma\ \frac{\partial\mathcal{F}^{*}_{j}}{\partial a_{i}}$, $w_{0}:=\mathcal{J}^{*}_{11}-\lambda$, the vectors $\textbf{v},\textbf{w}\in\mathbb{R}^{N-1}$ such that

[TABLE]

and

[TABLE]

which proves the proposition. ∎

3.2.2. Non-homogeneous perturbations

We now consider a perturbation of the form $\textbf{a}_{s}=(a_{s,1},a_{s,2},...,a_{s,N})$ for $s\in(-1,1)$ of the homogeneous steady-state $\textbf{a}^{*}$ in the direction of $\Phi=\left(\varphi_{1},\varphi_{2},...,\varphi_{N}\right)$, for non-homogeneous $\varphi_{j}:\Gamma\times(0,T)\to\mathbb{R}$. Thus for each component, we assume

[TABLE]

so we may write the linear approximation $a_{s,j}\approx a^{*}_{j}+s\varphi_{j}$ for $j=1,\ldots,N$

[TABLE]

In particular, the linearization of the non-local functional yields \mathcal{U}[\textbf{a}_{s}]\approx\mathcal{U}[\textbf{a}^{*}]+s(\frac{d}{ds}\big{|}_{s=0}\ \mathcal{U}[\textbf{a}_{s}]) where

[TABLE]

Since we assume that $\varphi_{j}\in L^{2}(\Gamma)$ are orthogonal to the constant perturbations, which were analyzed in the previous section, we now consider

[TABLE]

which leads to a linearized system with a constant input $\mathcal{U}[\textbf{a}_{s}](t)=\mathcal{U}[\textbf{a}^{*}]$. For the approximation of the component $a_{1}$, we thus have

[TABLE]

where

[TABLE]

[TABLE]

The other terms of the Jacobian matrix remain the same as in the case of the homogeneous perturbations, so we omit the explicit calculations. In vector notation, we can then write the linearized system in the form

[TABLE]

where D is a diagonal matrix such that $\textbf{D}_{jj}=d_{j}$ where $d_{1}=1$ and $\tilde{\mathcal{J}}(\textbf{a}^{*})$ is the modified Jacobian matrix. We then define $\mathbb{N}_{0}:=\mathbb{N}\cup\{0\}$ and consider $(\omega_{l})_{l\in\mathbb{N}_{0}}\subset L^{2}(\Gamma)$, an orthonormal basis of infinitely smooth eigenfunctions of the Laplace-Beltrami operator, i.e,

[TABLE]

Then for each $j=1,\ldots,N$ we can express each component $\varphi_{j}$ as a linear combination

[TABLE]

where $\alpha_{jl}=\alpha_{jl}(t)$ for $l\in\mathbb{N}_{0}$ . Using vector notation, we can define the quantity $\mathcal{A}_{l}=(\alpha_{1l},\alpha_{2l},\ldots,\alpha_{Nl})^{T}$ such that

[TABLE]

By substituting the above expansion in (3.6), we obtain the linear ODE system

[TABLE]

and diffusion-driven instabilities occur if the above system is unstable for some $l\in\mathbb{N}_{0}$. This is true when at least one eigenvalue $\lambda$ of the matrix $-\eta_{l}\textbf{D}+\tilde{\mathcal{J}}(\textbf{a}^{*})$ has a positive real part. Therefore our target quantity is the so-called dispersion relation

[TABLE]

where $\operatorname{Re}(z)$ denotes the real part of a complex number $z$. Finally, the characteristic polynomials $p_{l}(\lambda):=\det(\lambda\textbf{I}-\gamma\tilde{\mathcal{J}}(\textbf{a}^{*})+\eta_{l}\textbf{D})$ can be written in the form

[TABLE]

where $b_{l,0}=det(-\gamma\tilde{\mathcal{J}}(\textbf{a}^{*})+\eta_{l}\textbf{D})$. Therefore, if $b_{l,0}<0$ for some $l\in\mathbb{N}$, then $p_{l}$ has a positive root and therefore $h(l)>0$.

3.3. Special Case $N=2$: Necessary Conditions for Diffusion-Driven Instabilities

We now fix $N=2$ and analyze the conditions for diffusion-driven instabilities. The equations are given by

[TABLE]

We provide a necessary condition in a special case where the system admits a unique spatially homogeneous steady-state. We prove that the system does not exhibit diffusion-driven instabilities provided that $k_{b}$ is sufficiently small.

Theorem 3.1**.**

Suppose $k_{b}\geq 0$ is such that

[TABLE]

for all $i\in\mathbb{N}$. Then the system admits a unique steady-state and no diffusion-driven instability exists.

Proof.

Let $\textbf{a}^{*}=\left(a^{*}_{1},a^{*}_{2}\right)$ be the spatially-homogeneous steady $\textbf{a}^{*}=\left(a^{*}_{1},a^{*}_{2}\right)$, which is obtained when $a^{*}_{2}=\frac{k_{m}(a^{*}_{1})^{2}}{k_{2}}$ and $a^{*}_{1}$ is a solution of $\mathcal{P}_{2}(\alpha)=0$, where

[TABLE]

Now, since $k_{b}\leq\frac{2k_{0}|\Gamma|}{\mathcal{M}_{0}}$, we have $2k_{0}|\Gamma|-\mathcal{M}_{0}k_{b}\geq 0$, and therefore the intermediate value theorem ensures that the system admits a unique positive steady-state. The Jacobian matrix with respect to homogeneous perturbations is then given by

[TABLE]

with a second-order characteristic polynomial $p(\lambda)=det(\lambda\textbf{I}-\mathcal{J}[\textbf{a}^{*}])$ given by $p(\lambda)=\lambda^{2}+b\lambda+c$, where

[TABLE]

and

[TABLE]

is also positive because $2k_{0}|\Gamma|-\mathcal{M}_{0}k_{b}\geq 0$. From that we conclude that both eigenvalues

[TABLE]

must have real negative parts and therefore the steady-states are linearly stable. We then perform a similar argument for non-homogeneous perturbations. From (3.5), we obtain the modified Jacobian matrix $\tilde{\mathcal{J}}(\textbf{a}^{*})$, and for a given $l\in\mathbb{N}_{0}$, we have

[TABLE]

with characteristic polynomials $p_{l}(\lambda):=\det(\lambda\textbf{I}-\gamma\tilde{\mathcal{J}}(\textbf{a}^{*})+\eta_{l}\textbf{D})$ given by the quadratics $p_{l}(\lambda)=\lambda^{2}+b_{l}\lambda+c_{l}$, where

[TABLE]

and

[TABLE]

In the case $l\in\mathbb{N}$, the $c_{l}$ terms are also positive, since we assume

[TABLE]

and in this case the $p_{l}(\lambda)$ have no roots with positive real parts. In the case $l=0$, we know that $\eta_{0}=0$, but the modified matrix $\tilde{\mathcal{J}}(\textbf{a}^{*})$ yields a different linearized system. Thus we have to analyze the stability of the linear equation given in (3.7) in the case where $l=0$, i.e,

[TABLE]

with characteristic polynomial $p_{0}(\lambda)=\lambda^{2}+b_{0}\lambda+c_{0}$ where

[TABLE]

and

[TABLE]

We now verify that $c_{0}\geq 0$. In fact, from (3.2) when $N=2$, we obtain

[TABLE]

and therefore by substituting the above equation on (3.11) and using that $a^{*}_{2}=\frac{k_{m}}{k_{2}}(a^{*}_{1})^{2}$, we obtain

[TABLE]

Finally, the hypothesis gives us $\frac{k_{0}}{k_{b}}\geq\frac{\mathcal{M}_{0}}{2|\Gamma|}$ and by using $\mathcal{M}_{0}-2|\Gamma|a^{*}_{2}\geq 0$ (total mass of $a_{2}$ at steady-state does not exceeds the total mass of the system), we obtain $\frac{\mathcal{M}_{0}}{2|\Gamma|}\geq a^{*}_{2}$ and therefore

[TABLE]

from which we conclude that $c_{0}\geq 0$. Therefore, the steady-state is stable against non-constant perturbations.

∎

3.4. Existence of Solutions to the General and Reduced Models

In the sections above, we developed a stability analysis of the reduced model (2.23)–(2.24), motivated by Turing in [46], and essentially following the technical approach in [33, 34]. Understanding the stability of solutions (if any exist) with respect to parameters in the reduced model is part of a larger set of basic mathematical questions regarding the solutions to the reduced system (2.23)–(2.24), as well as to the more general system (2.2)–(2.7). Knowing when these types of models are well-posed (solutions exist, are unique, and depend continuously on problem data) is critical for drawing scientific conclusions from mathematical estimates and numerical simulations of solutions. When models of critical phenomena correctly capture solution instability and even multiplicity, knowing when and how many solutions are permitted by the model is also crucial. In the context of dynamical problems, a standard approach to developing a solution theory is to prove that solutions exist for small time (local existence), and then endeavor to extend such results to large time (global existence). In developing existence results one generally first considers small data (initial and boundary conditions being small perturbations of zero), with again the hope of extending such results to large data (essentially no restriction on initial or boundary data, other than reasonable smoothness assumptions). Establishing that solutions are unique typically involves exploiting some inherent mathematical structure in the underlying equations.

In [48], an analysis of a reaction-diffusion system similar to our general model (2.2)–(2.7) is developed; see also [49] for related results. The main results in [48] include a global-in-time, large data, existence result based on a priori estimates for a linearized model and compactness techniques, together with a fixed-point argument using on a variant of the Leray-Schauder Theorem (Schaefer’s Theorem). They also obtain a uniqueness result using a Gronwall-type inequality. Their results require growth conditions on the nonlinearities appearing in the system (a local type of Lipschitz property, together with pointwise control of individual components in the solution). These conditions appear to be too restrictive to cover the nonlinearities arising here in (2.2)–(2.7), which include quadratic and more rapidly growing terms. Nevertheless, we expect that results similar to those in [48] can be shown to hold for (2.2)–(2.7) using similar arguments.

4. Numerical Simulations

We perform numerical simulations to complete our mathematical analysis. We start by analyzing the parameter regions of bistability (Section 4.1). Then we investigate whether the stable steady-states become linearly unstable under non-homogeneous perturbations (Section 4.2). From the linear instability analysis, we obtain the single-patch non-homogeneous steady-state (Section 4.3). Finally, we study the temporal dynamics of pattern formation (Section 4.4) and the single-patch dependence on the cell radius (Section 4.5). The numerical simulations were implemented in Matlab R2018a and Comsol Multiphysics 5.4.

4.1. Bistability under homogeneous perturbations

We begin by computing the homogeneous steady-states $\textbf{a}^{*}$ and the corresponding eigenvalues of the Jacobian matrix $\mathcal{J}[\textbf{a}^{*}]$ under homogeneous perturbations (cf. Section 3.2.1). We then explore the parameter regions of bistability where the system admits three steady-states, two of them stable and one unstable. In the case $N=2$, we obtain regions of bistability by change the basal binding rate $k_{0}$ and the $A_{2}$-dependent binding rate $k_{b}$ (Figure 2). For $k_{0}=0.015$, three steady-state values for $a^{*}_{1}$ emerge depending on $k_{b}$ (Figure 2 (A)). When $k_{0}$ also changes, we obtain both a bistability region (dark-gray) and a single steady-state region (light-gray) (Figure 2 (B)). A colored vertical line represents the region from Figure 2 (A). Other parameter choices also lead to bistability regions (see Figure S1 (A) for $N=2$ and Figure S2 (A) for $N=3$ in the ESM).

4.2. Linear instability under non-homogeneous perturbations

In this section, we numerically investigate which parameter values promote linear instability under non-homogeneous perturbations. We fix an eigenmode index $l\geq 1$ to explore diffusion-driven instabilities, and let $\textbf{a}^{*}$ be a stable steady-state under homogeneous perturbations. We can thus compute the dispersion relation $h(l)$ (3.8) defined in Section 3.2.2. If $h(l)<0$, the steady-state remains stable in the direction of the chosen eigenmode. In this case, the analysis is inconclusive, since we would also need to determine the stability for the other eigenmodes. If $h(l)>0$, the steady-state becomes unstable for the chosen eigenmode, and this is sufficient to ensure a diffusion-driven instability [50]. The case $h(l)=0$ is usually requires higher-order analysis, so we will not consider it in the context of linear stability.

Given a fixed eigenmode index $l$, we can then divide the parameter space into four regions. We will call them Regions [math], $1$, $2$, and $3$, where the numbers reflect the exact number of unstable steady-states. More precisely, we define:

• •

Region 0: The single steady-state region where $h(l)<0$; There are no unstable steady-states.

• •

Region 1: The single steady-state region where $h(l)>0$; There is only one unstable steady-state.

• •

Region 2: The bistability region where $h(l)>0$ for only one stable steady-state; a total of two unstable steady-states.

• •

Region 3: The bistability region where $h(l)>0$ for both stable steady-state; a total of three unstable steady-states.

Once we find the bistability and single steady-state regions, we can divide the same parameter space into Regions 0, 1, 2, and 3, by computing the number of unstable steady-states (Figure 3(A) for $N=2$ and (B) for $N=3$). The stability analysis in Region 0 is more subtle and requires further analysis since the stability criterion needs to be fulfilled for all eigenmodes. However, at least for $N=2$, 3.1 ensures that the system remains stable for sufficiently small $k_{b}$, which appears to be consistent with the numerical predictions. For higher $k_{b}$ values, the instabilities emerge in the bistability region (Regions $2$ and $3$) and also in the single steady-state regions Regions ([math] and $1$). We obtain a similar result for $N=3$ (Figure 3 (B) ). However, it should be noticed that the $k_{b}$ values promoting linear instabilities are higher (see $y$-axis ranging from $4$ to $14$) compared to with the case $N=2$. Regions 0, 1, 2, and 3 can be found with other parameter choices (see Figure S1 (B) for $N=2$ and Figure S2 (B) for $N=3$ in the ESM).

The Region 1 where a stable steady-state becomes unstable are known as Turing-type instability regions [46, 51], where the system may converge to a spatially non-homogeneous steady-state. We analyze these region when we increase both $k_{0}$ and $k_{b}$ ranges for different values of the diffusion coefficient $d_{2}$ (Figure 4 (A)). As $d_{2}$ decreases, Region 1 increases, which illustrates how the system becomes unstable as the discrepancies between diffusion become higher. A similar phenomenon occurs as we increase the dimensionless parameter $\gamma$ (Figure S3). On the other hand, as the eigenmode index $l$ increases, Region 1 significantly decreases (Figure 4 (B)). We exhibit the results for $l=2$, $l=6$, and $l=8$. Such a decrease allows us to explore the instability of the system by considering $(k_{0},k_{b})$ on Region 1 only for the first eigenmode $(l=1)$.

4.3. The emergence of the single-patch non-homogeneous steady-state

In this section, we investigate the spatio-temporal behavior of our system by numerically integrating the dimensionless equations. We consider a spherical domain of radius $R=1$ and, as in the previous sections, we fix $N=2$ or $N=3$. We avoid solving the surface system (2.23)–(2.24) due to the numerical complexity of the non-local functional. Instead, we solve the dimensionless bulk-surface equations (2.16)–(2.21) (dropping all the hats) for an extremely high cytosolic diffusion ($\tilde{D}=10^{8}$) on (2.16). In this way, our resulting system can be seen as an approximation of the reduced system when $\tilde{D}\to\infty$. We randomly perturbed the homogeneous steady-states by considering a small number $\varepsilon>0$ as the perturbation magnitude and a family $\{\xi(x)\}_{x\in\Gamma}$ of independent random variables uniformly distributed between $-\varepsilon$ and $\varepsilon$. In the case where $N=2$, we define the surface initial conditions

[TABLE]

where the $\frac{1}{2}$ accounts for mass conservation (see (2.22)). For the volume component, we define $u(x,0):=u^{*}$, where $u^{*}=\frac{1}{|\Omega|}\left[\mathcal{M}_{0}-|\Gamma|(a^{*}_{1}+2a^{*}_{2})\right]$ also because of the mass conservation property. For $N=3$ we define $a^{*}_{1}$ as in (4.1), $a_{j}(x,0):=a^{*}_{j}-\frac{1}{5}\varepsilon\ \xi(x)$ for $j=2$ and $j=3$, and $u^{*}=\frac{1}{|\Omega|}\left[\mathcal{M}_{0}-|\Gamma|(a^{*}_{1}+2a^{*}_{2}+3a^{*}_{3})\right]$.

Remark 4.1**.**

The element $(a^{*}_{1},a^{*}_{2},...,a^{*}_{N})$ is a homogeneous steady-state of the system (2.23)–(2.24) if and only if $(a^{*}_{1},a^{*}_{2},...,a^{*}_{N},u^{*})$ is a homogeneous steady-state of the system (2.16)–(2.21) provided that

[TABLE]

By the remark above, we can obtain the steady-states given the parameter choice in the reduced system (2.23)–(2.24). Then we can numerically integrate the bulk-surface PDE system (2.16)–(2.21) using the perturbation scheme described above. In order to associate the parameter regions that lead to instabilities with the formation of spatial patterns, we select four $(k_{0},k_{b})$ values in the four Regions 0, 1, 2, and 3 (Figure 5). We fix $N=2$ and the eigenmode index $l=1$. For each choice of $(k_{0},k_{b})$ , we integrate the system (2.16)–(2.21) to its final state by perturbing a homogeneous steady-state. We then plot the result for the $a_{1}$ component and visually inspect the results. For $(k_{0},k_{b})$ in Regions 1, 2, and 3, (colored in orange, red or black, respectively), a single-patch spatially heterogeneous steady-state emerges. On the other hand, when $(k_{0},k_{b})$ belong to Region 0, in which the system is stable for the eigenmode index $l=1$, the system converges to its homogeneous steady-state. This result indicates that the single-patch pattern is consistent across parameter choices in Regions 1, 2, and 3, once it remains unchanged in its circular shape and gradient of concentrations. Figure S4 in the ESM shows a similar result in the case $N=3$.

4.4. Temporal evolution and pattern formation

In this section, we further investigate the temporal evolution of the system. We consider $N=2$ and $(k_{0},k_{b})=(0.025,2.5)$ that belongs to Region 1 (see Figure 5). We then observe the spatial distribution of $a_{1}$ for different times (Figure 6 (A)). At $t=0$, We apply a random perturbation of magnitude $\varepsilon=10^{-10}$ around the unique homogeneous steady-state that is unstable under non-homogeneous perturbations. The system then smooths due to diffusion and the small random peaks continuously coalesce and react, until a few large domains emerge at $t=0.099$. At $t=0.114$ and $t=0.119$, multiple patches of higher $a_{1}$ concentration emerge. The feedback term (2.1) then plays its role, once the higher $a_{2}$ concentration location promotes the recruitment of more cytosolic component. This leads to the formation of the single-patch profile at $t=0.159$. From that time until the final time ($t=1$), the spatial configuration only changes in terms of concentration gradients. File F1 in the ESM contains a movie of the simulation shown in Figure 6 (A) for both monomeric ($a_{1}$) and dimeric components ($a_{2}$). In Figure S5 in the ESM, we show a similar result for $N=3$.

In order to quantify the single-patch size, we quantify the surface area of the high-concentration locations in the spherical domain. For this purpose, we define the function

[TABLE]

where $\varepsilon$ is the perturbation magnitude, $j$ is the index of the oligomeric component, and $\langle a_{j}\rangle(t)=\int_{\Gamma}a_{j}ds$ is the average concentration of $a_{j}$ across the sphere $\Gamma$. In Figure Figure S6 in the ESM, we illustrate how the percentage of $\mathcal{S}^{\varepsilon}_{j}$ with respect to the total surface area does not change significantly as $\varepsilon$ changes. We then evaluate the evolution of $\mathcal{I}^{\varepsilon}_{j}(t)$ over time (Figure 6 (B)). We exhibit the results of a single simulation for $N=2$ and $N=3$, and $\varepsilon=10^{-10}$. At early times, when the concentrations $a_{j}$ are close to the steady-state $a^{*}_{j}$ across the domain, $\mathcal{I}^{\varepsilon}_{j}(t)$ remains close to 0. Then the combination of diffusion and the feedback term makes the concentration gradients increase in a large portion of the domain, as illustrated in Figure 6 (A) for $t=0.099$, $t=0.114$, and $t=0.119$ . Finally, the formation of the single-patch promotes the decrease of $\mathcal{I}^{\varepsilon}_{j}(t)$, since the area of high concentration tends to be small in comparison with the total surface area. Moreover, the concentration outside the patch tends to be small, which makes the average $\langle a_{j}\rangle(t)$ assume lower values. Therefore, in the final times, the locations in the sphere where the concentrations remain above the average can be associated with the single patch. For this reason, we define the single-patch area

[TABLE]

where $t_{f}$ is the final simulation time. In this work, we avoid an analytical treatment for the temporal behavior of $\mathcal{I}^{\varepsilon}_{j}(t_{f})$. Instead, we base our definition of the single-patch area on visual inspection of the temporal evolution of $\mathcal{I}^{\varepsilon}_{j}(t_{f})$. In fact, we obtain the same temporal dynamics from Figure 6 (B) whenever the system forms the single-patch pattern. We also observe that $\mathcal{I}^{\varepsilon}_{1}>\mathcal{I}^{\varepsilon}_{2}$ for all times in the case where $N=2$, and also $\mathcal{I}^{\varepsilon}_{1}>\mathcal{I}^{\varepsilon}_{2}>\mathcal{I}^{\varepsilon}_{3}$ in the case where $N=3$. We conclude that $\mathcal{S}^{\varepsilon}_{1}>\mathcal{S}^{\varepsilon}_{2}$ (for $N=2$) and $\mathcal{S}^{\varepsilon}_{1}>\mathcal{S}^{\varepsilon}_{2}>\mathcal{S}^{\varepsilon}_{3}$ (for $N=3$). In order to better visualize this area shrinking as the oligomer size increases, we plot the final normalized concentration profiles (Figure 6 (C)). Given the arc-length parametrization of a geodesic curve crossing the single-patch region, the concentration distributions become tighter for $a_{2}$ compared to $a_{1}$ in the case $N=2$. The inset plot shows the non-normalized concentrations, where we see that $a_{2}>a_{1}$ in the single-patch location. A similar phenomenon occurs for $N=3$: the distribution and maximum value of $a_{j}$ becomes tighter and larger as $j$ increases from 1 to 3.

4.5. Change of the cell radius and single-patch area

We investigate how the single-patch area of a spherical cell depends on its radius $R$. From the non-dimensionalization of the bulk-surface system (see Section 2.4), we defined the characteristic quantities (2.15). In order to move through a dimensional system, we define the dimensionless parameters depending to $R$ to preserve a constant volume concentration. Therefore, we assume a constant $U$ such that $\mathcal{M}_{0}\propto R^{3}$, making the dimensionless parameters as functions of $R$. For this reason, each choice of $R$ will lead to a different solution of the non-dimensional system (2.16)–(2.21). In particular, it will also change the non-dimensional single-patch area $\mathcal{S}^{\varepsilon}_{j}$. We show the results for $R$ ranging from $0.5$ to $5$ and two different parameters (Figure 7):

[TABLE]

For $N=2$ and $N=3$, we provide the same total volume concentration for the system. For clarity, in this section we will refer to the non-dimensional system with the hat (\hat{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}a}}) notation. We define the initial conditions as a linear ramp of slope $\varepsilon$ around the steady-state

[TABLE]

where $\hat{x}=(\hat{x}_{1},\hat{x}_{2},\hat{x}_{3})\in\hat{\Gamma}$ and $\hat{x}_{1}$ is the position with the sphere centered in the origin. In the case $N=3$, we assumed $a_{3}^{*}=0$. See Figure 7 caption for details on the parameter choices. We then observe the same phenomena: the dimensional area (fig. 7 red; open circles) increases approximately linearly with $R$. On the other hand, the area percentage (fig. 7 black; closed circles) decreases with $R$. We can then conclude that although the dimensional area of clusters increases, the additional spherical area changes at a much faster rate since the area percentage varies with $\approx\frac{1}{R}$. The effect of $N$ on single-patch area shows both dimensional area and area percentage are higher for $N=3$ (fig. 7(B)) in comparison with $N=2$ (fig. 7(A)).

5. Discussion

Protein heterogeneity in the PM is of critical importance to cellular functions. Many factors influence this heterogeneity, including membrane composition, protein-protein interaction, phase separation, lateral diffusion, and possible feedback, resulting in the formation of spatial patterns [7, 8, 9]. For this reason, understanding the interplay of aggregation kinetics, lateral diffusion, and feedback in the formation of spatial patterns is an essential step towards developing a complete description of the mechanisms behind protein clustering on the cell surface. In this work, we developed a bulk-surface model for protein aggregation with positive feedback that exhibits a spatial heterogeneous single-patch steady-state. To the best of our knowledge, this is the first modeling attempt that merges the reaction-diffusion version of classical Smoluchowski dynamics with the modern bulk-surface geometrical setup.

A major result from our model is the role played by the feedback term $k_{b}\ a_{N}$ in the boundary conditions (2.1). If $k_{b}$ is low enough, the steady-state distribution is spatially uniform, and no protein heterogeneity exists. For $N=2$, we formally proved such a result (3.1), and for $N=3$, we used numerical simulations to observe a similar phenomenon. In particular, in the total absence of feedback ($k_{b}=0$), we observed that spatial heterogeneity is not achievable when we only considered protein-protein interaction. On the other hand, if $k_{b}$ is sufficiently high, we observed the emergence of linear instability and therefore patterning on the cellular surface. Experimental observations have shown that membrane proteins do organize in a spatially heterogeneous fashion [23, 24, 25]. However, the molecular mechanisms are still being investigated experimentally. The feedback mechanism we proposed here can also be interpreted in purely biological terms. The largest oligomers recruit ligands from the cytosol, which form ligand-receptor monomers. If the rate of recruitment of monomers is low, diffusive effects dominate, and the configuration of the system is homogeneous in space. On the other hand, a higher rate promotes a significant input of new monomeric components. Then, continued oligomerization generates higher concentrations of the largest components, which closes the positive feedback loop and drives pattern formation. The largest oligomers can thus be interpreted as self-activators of pattern formation. For this reason, our mechanism of pattern formation can be related to the classical Turing framework where self-activation is required to generate spatial patterns [46, 52]. Another interesting aspect of our model is the absence of an explicit description of cooperative binding. For the wave-pinning model [53, 54, 29], cooperativity is included with a Hill function, which accounts for the positive feedback. In contrast, our oligomerization reactions assume only mass-action kinetics, which seems to be insufficient for pattern formation without the feedback term.

Bistable systems are well known to promote diffusion-driven instabilities in the context of cell polarization [55, 56]. For the wave-pinning model [53, 54, 29], the structure of the Hill function is responsible for bistability. Other studies followed a similar approach, using a particular choice of reaction flux that is naturally associated with a bistable regime [57, 58]. In our model, bistability emerges by the combination of two key ingredients: positive feedback and mass conservation. This observation becomes clear as we carefully inspect the steady-state analysis of the reduced system (cf. Section 3.1). First, the equilibrium of the oligomerization reactions (driven by mass action kinetics only) provides the distribution across the different surface components. Then, the input from the non-local functional comes into play, as a consequence of the boundary conditions and mass conservation. The non-local functional at steady-state provides an extra equation, which gives the equilibrium solutions for the monomeric component. The particular contribution of the feedback comes from the coefficient $C_{N}\left(k_{0}|\Gamma|N-\mathcal{M}_{0}k_{b}\right)\alpha^{N}$ of the polynomial $\mathcal{P}_{N}(\alpha)$. If the coefficient is negative, then the existence of three roots, and therefore three steady-states, is achievable. In this case, we can compute their stability under homogeneous perturbations and verify bistability.

Under non-homogeneous perturbations, one or two stable steady-states may become unstable, and the system undergoes a diffusion-driven instability. Even more impressive is the emergence of a linear instability parameter region, called Region 1 in this study, when the system admits a single steady-state that becomes unstable. We note that in Getz et al. [59], the authors were able to find a region of linear instability for the Wave-Pinning model that is comparable with our Region 1. While the authors briefly discussed the changes in that parameter region for different wave-numbers, here we explicitly showed that the leading eigenmode exhibits a region of instability that shrinks as the eigenmode index increases. Such instability in the lower modes, which are associated to the smallest positive eigenvalues of the Laplace-Beltrami operator, has been ofter related to a single-patch steady-state pattern [60, 61], which is confirmed for our system.

The single-patch steady-state consistently appears for parameter values corresponding to the different instability regions (called as Regions 0, 1,2 and 3). Goryachev et al.  [61] found a similar spatial profile for the Cdc42 GTPase cycle, where the income of new cytoplasmic components maintained the cluster steady-state and compensated for its lateral diffusion. A similar phenomenon seems to happen in our system. An allegory that explains the stable existence of such heterogeneous steady-states is the so-called “rich get richer” competition [62]. In this case, larger domains outcompete for the smaller until only one stable domain arises. Our hypothesis about the existence of the single-patch is based on the role of the positive feedback term. We assume that the presence of high concentrations of the largest oligomer promotes ligand binding onto the PM in a linear fashion, without any saturation mechanism nor steric effects. As in [61], this assumption seems to account for a resource competition that excludes the possibility of multiple patches.

Based on the insights from our model, we identify future research directions that will enhance studies such as ours. In the current formulation, we lack a formal explanation for the emergence and robustness of the single patch steady-state. The spatial aspects of the model render such analysis hard, but it may be possible to obtain a formal proof by considering a one-dimensional version of our system as in [53]. Another interesting quantity to be computed in future studies is the so-called amplitude of the pattern, for which a formal calculation was recently developed [63]. Additionally, the mathematical challenge for a theoretical stability result lies in the increasing complexity of the system as $N$ increases. In this case, we have relied on numerical simulations for $N=3$ to identify the threshold phenomenon for diffusion-driven instabilities. However, future efforts in this direction could open up new mathematical avenues for stability analysis of increasingly complex systems. Finally, including the role of curvature and cytosolic diffusion in the formation of membrane protein aggregates would bring us closer to analyses of biological and biophysical systems. These are the focus on ongoing studies in our group.

6. Acknowledgments

This work was supported by Air Force Office of Scientific Research (AFOSR) Multidisciplinary University Research Initiative (MURI) grant FA9550-18-1-0051 to P. Rangamani. M. Holst was supported in part by NSF Awards DMS 1620366 and DMS 1345013.

7. Electronic Supplementary Material