text:
\section{Introduction}
Cell growth and division require long-distance communication and coordination between two ends of the cell. To polarize, cells require the ability to form distinct cellular domains with different molecular components.  Establishment of cell polarity relies on local accumulation of signaling molecules on the membrane and has attracted considerable attention from mathematicians and physicists \cite{jilkine2011comparison,rappel2017mechanisms,halatek2018self}. Positive and negative feedback loops can result in spontaneous symmetry breaking of an initially homogeneous distribution of polarity regulators to form one or multiple clusters, as well as dispersal and reformation of these clusters in an oscillatory manner \cite{johnson2011}. Thus, polarization involves spatial diffusion processes coupled with biochemical reactions occurring within localized signaling compartments.

The master regulator of cell polarity in a variety of organisms, from yeast to humans, is the protein Cdc42 \cite{etienne2002,martin2014}. Previous mathematical modeling of Cdc42 in cell polarization has focused primarily on symmetry breaking and establishment of active Cdc42 cortical zones in budding yeast, often via a Turing mechanism (reviewed in \cite{goryachev2017}). Detailed and phenomenological models of cell polarization investigated the effects of single or multiple positive feedback loops and mass conservation on the formation of a unique polarity zone \cite{otsuji2007,goryachev2008, mori2011, jilkine2011, freisinger2013,chiou2018}. In some circumstances, two or more active Cdc42 domains can coexist for some time, and the basis for the switch from a single to multiple polarity zones is not yet fully understood \cite{wu2013,chiou2018}. Here we focus on a model of spatio-temporal oscillations of Cdc42 from pole-to-pole that have been observed in fission yeast \cite{das2012}, a model organism for understanding how cells integrate polarity and spatial coordination of growth.

The mathematical modeling of  cell polarization typically involves a system of reaction-diffusion equations for the concentrations of membrane-bound active Cdc42 molecules and cytosolic inactive Cdc42 molecules with no-flux or periodic boundary conditions \cite{goryachev2008, mori2008, mori2011, lo2014}. This approach does not take into effect the spatial segregation of polarity molecules in the cell. For Cdc42 and many other polarity molecules, most biochemical reactions of interest happen at the interface between the cytosol and the membrane. Because the active form of the protein is found on the membrane, while the inactive form is in the cytosol \cite{etienne2002}, a more appropriate formulation would be a membrane bulk-cytosol model with diffusion in the interior of the cell, and boundary conditions modeling the exchange between the membrane and the cytosol. Pattern formation in membrane-bulk models with applications to cell polarity have previously been analyzed in \cite{ratz2012turing, ratz2014, madzvamuse2015} through linear stability analysis and PDE numerical simulations. Stationary front solutions corresponding to a polarized cell were demonstrated to exist in both the reaction-diffusion \cite{mori2011} and membrane-bulk versions \cite{cusseddu2018} of the wave-pinning model of cell polarization. When considered on a 1-D bulk domain, coupled membrane-bulk reaction-diffusion systems become coupled PDE--ODE systems. Recent studies of coupled PDE--ODE models \cite{gomez2007, gou2015_SIAM, gou2017,xu2017} showed the possible collective synchronization of localized active units (or membranes) coupled through a linear bulk diffusion field, even if each unit is at rest when isolated from the group.

Our work focuses on the spatio-temporal dynamics of Cdc42 during cell growth of fission yeast. Fission yeast (\textit{S. pombe}) cells are rod-shaped and grow in length by tip extension, while maintaining a constant diameter of 3 $\mu$m. During their cell cycle, the cells switch from growing at one end to growing at both ends. The active Cdc42 molecules localize at the growth site. At the cell ends, Cdc42 exhibits temporal oscillations with a period around 5 minutes \cite{das2012}. Earlier mathematical models for cell polarity in fission yeast emphasize the role of positive feedbacks on symmetry breaking \cite{csikasz2008,cerone2012}. These models assume that some ``tip factors'' (delivered by microtubules) are present at the two ends of the fission yeast that leads to Cdc42 preferentially binding at the two ends of the cell \cite{kokkoris2014,martin2014}. Recent models for spatio-temporal oscillations of Cdc42 involve both positive and negative feedback loops, as well as time delays \cite{das2012,xu2016}. Several potential positive feedbacks are presented in a recent review paper \cite{goryachev2017} while potential negative feedbacks may involve a Cdc42 GEF (Guanine nucleotide exchange factor) \cite{das2015}. In \cite{xu2018}, we used a 1-D coupled PDE--ODE model to include possible positive and negative feedbacks between Cdc42 and its GEF, with our analysis of the resulting dynamics focusing on the fast bulk diffusion limit, for which the full model is reduced to an ODE model. We then explored in \cite{xu2018BMB} the effect of intrinsic noise on Cdc42 oscillations and compared the dynamics of the reduced ODE model with its stochastic counterpart.

In this work, we relax the fast bulk diffusion assumption and compare the dynamics of the coupled PDE--ODE system to its well-mixed ODE limit using linear stability analysis, time-dependent simulations and numerical bifurcation analysis. We study the effect of decreasing the bulk diffusion level on the stability of the symmetric steady state, at which both ends of the cell have the same amount of Cdc42. We find conditions for a Hopf bifurcation, and analyze how the Cdc42 and GEF bulk diffusion rates affect its criticality. We also explore the role of one key biochemical reaction rate, the GEF dissociation rate, on the emergence of oscillatory dynamics.  We present numerical simulations showing two different types of pole-to-pole oscillations, either weakly nonlinear or of relaxation type.

In the second part of the paper, we extend the 1-D bulk-cytosol to a two-dimensional circular domain, with pure diffusion in the interior, as well as surface diffusion and nonlinear reactions on the membrane. We also consider a reduced one-dimensional system of nonlocal reaction-diffusion equations that governs the dynamics of membrane-bound species when assuming a well-mixed bulk. Linear stability analysis is performed for both the full membrane-bulk model and its nonlocal reduction. Our aim is to investigate the symmetry-breaking effects of varying the bulk diffusion coefficients, the membrane diffusion coefficients and the GEF dissociation rate. Finally, we present numerical results that confirm our analysis and exhibit a variety of spatio-temporal patterns on a unit disk, such as stationary Turing patterns, traveling waves and standing waves.

\section{Mathematical models}\label{sec:models}

\subsection{One-dimensional PDE--ODE model}
We first model the cell as a one-dimensional bulk domain separating two well-mixed compartments. This simplified geometry is based on the competition of two growth zones of active Cdc42, localized at the cell tips, for a common substrate, corresponding to inactive Cdc42, that diffuses in the cytosol and can become active at the two tips in the presence of a Cdc42 GEF. We also consider the spatial distribution of the Cdc42 GEF, found at each cell tip in its active form or in the cytosol in its inactive form. A schematic representation of this 1-D model is shown in Fig.~\ref{fig:2Dmodel}a.

\begin{figure}[t!]
\centering
\includegraphics[width=0.7\linewidth]{figures/Schematic}
\caption{(a) A one-dimensional PDE--ODE model with bulk diffusion in the interior given by equation \eqref{eq:pde} and binding kinetics at two well-mixed ends given by the ODE system \eqref{eq:ode}. (b) A two-dimensional model given by equation \eqref{eq:2Dmodel} and equation \eqref{eq:consv2dMAIN} with diffusion along the membrane. Filled circle: membrane-bound molecule. Open circle: cytosolic molecule. Blue: Cdc42. Green: GEF.}
\label{fig:2Dmodel}
\end{figure}

Let $C(x,t)$ and $G(x,t)$ denote the concentrations of inactive Cdc42 and GEF in a cytosolic (or bulk) domain of length $L$, satisfying
\begin{equation}\label{eq:pde}
\frac{\partial C}{\partial t} = D_c \frac{\partial^2 C}{\partial x^2}\,, \quad \frac{\partial G}{\partial t} = D_g \frac{\partial^2 G}{\partial x^2}\,, \quad 0 < x < L.
\end{equation}
where $D_c$ and $D_g$ are two constant bulk diffusion coefficients. Near each of the two tips $(i=1,\,2)$, the binding and unbinding processes are modeled by flux boundary conditions given by, in $x=0$,
\begin{equation}\label{eq:BC_1}
\left. D_c\frac{\partial C}{\partial x}\right|_{x=0} = k^+(c_1(t),g_1(t)) C(0,t) - k^- c_1(t) \,, \quad \left. D_g\frac{\partial G}{\partial x}\right|_{x=0} = k^{\rm on}(c_1(t))G(0,t) - k^{\rm off} g_1(t)\,,
\end{equation}
while in $x=L$ the boundary conditions satisfy,
\begin{equation}\label{eq:BC_2}
\left. -D_c\frac{\partial C}{\partial x}\right|_{x=L} = k^+(c_2(t),g_2(t)) C(L,t) - k^- c_2(t) \,, \quad \left. -D_g\frac{\partial G}{\partial x}\right|_{x=L} = k^{\rm on}(c_2(t))G(L,t) - k^{\rm off} g_2(t)\,,
\end{equation}
where $c_i(t)$ and $g_i(t)$ denote the concentrations of active Cdc42 and GEF. We assume that active GEF promotes the binding of Cdc42 on each tip, and this results in the following association rate $k^+$:
\begin{equation}
k^+(c_i,g_i)=(k_0+k_{\rm cat} c_i^2)g_i\,.
\end{equation}
Here we assume a form of positive feedback that was used in previous models of fission yeast \cite{csikasz2008} and \cite{cerone2012}. The nonlinearity $c_i^2$ is necessary for symmetry breaking, allowing an asymmetric steady state to emerge. Potential chemical reactions that can lead to such kinetics have been discussed in \cite{goryachev2017}. Some form of negative feedback effect is necessary to obtain the pole-to-pole oscillations observed in fission yeast \cite{das2012}. For negative feedback, we assume that active Cdc42 inhibits the binding of GEF and that the association rate $k^{\rm on}$ is given by
\begin{equation}
k^{\rm on}(c_i)=\frac{k^{\rm on}}{1+\kappa c_i^2}\,.
\end{equation}
We also note that this specific implementation has been discussed in our previous work~\cite{xu2018}, where it was found that negative feedback is more likely to act through inhibition of GEF association rather than upregulation of GEF dissociation. Finally, the dissociation rates $k^-$ and $k^{\rm off}$ are chosen to be linear.

For the rest of the paper, we will employ the following notation for the binding kinetics:
\begin{equation}\label{eq:kinetics}
\mathcal{F}\left(C,c,g\right) = (k_0 + k_{\rm cat} c^2) g C - k^- c\,, \quad \mathcal{G}\left(G,c,g\right) = \frac{k^{\rm on}}{1+\kappa c^2} G - k^{\rm off} g\,.
\end{equation}
Using this notation, we define the nonlinear ODEs regulating the active species dynamics as
\begin{equation}\label{eq:ode}
\begin{split}
\frac{d c_i(t)}{dt}&= \mathcal{F}\left(C|_{x=(i-1)L},c_i(t),g_i(t)\right)\,, \quad i=1\,,2\,, \\
\frac{d g_i(t)}{dt}&= \mathcal{G}\left(G|_{x=(i-1)L},c_i(t),g_i(t)\right)\,, \quad i=1\,,2\,.
\end{split}
\end{equation}
Here, we highlight that protein production and decay are neglected by our model formulation, and that only the exchange between the active and inactive forms is considered. This results in the total amount of Cdc42 and GEF being constant over the timescale on which cell polarization takes place. Hence, our model formulation implies the following conservation laws:
\begin{equation}\label{eq:mass_cons}
 c_1(t) + c_2(t) + \int_0^L C(x,t) dx = C_{\rm tot}\,, \quad g_1(t) + g_2(t) + \int_0^L G(x,t) dx = G_{\rm tot}\,,
\end{equation}
where $C_{\rm tot}$ and $G_{\rm tot}$ are the total amount of Cdc42 and GEF.

Finally, we mention the limit of fast bulk diffusion, corresponding to the well-mixed regime, for which the concentrations of inactive species $(C(x,t), G(x,t))$ become spatially homogeneous $(C(t),G(t))$. Using the conservation laws \eqref{eq:mass_cons}, we can reduce the full coupled PDE--ODE model to a system of four nonlinear ODEs given by
\begin{equation}\label{eq:well_mixed}
\begin{split}
\frac{d c_i(t)}{dt}&= \mathcal{F}\left(C(t),c_i(t),g_i(t)\right)\,, \quad i=1,2,\\
\frac{d g_i(t)}{dt}&=\mathcal{G}\left(G(t),c_i(t),g_i(t)\right)\,,\quad i=1,2.
\end{split}
\end{equation}
where $C(t)$ and $G(t)$ are  defined by
\begin{equation}\label{eq:leading_order}
 C(t) \equiv \frac{C_{\rm tot} - c_1(t) - c_2(t)}{L}\,, \quad G(t) \equiv \frac{G_{\rm tot} - g_1(t) - g_2(t)}{L}\,.
\end{equation}

The model presented here is dimensionless. A justification of the various parameters (and their values) can be found in Appendix~\ref{sec:appendix_params}.

\subsection{Two-dimensional membrane-bulk reaction-diffusion model}\label{subsec:2_D_membrane_bulk}
In this section we consider the model on a two-dimensional domain $\Omega\in\mathbb{R}^2$ with boundary $\partial\Omega$ corresponding to the cell membrane. The geometry of the model is shown in Fig.~\ref{fig:2Dmodel}b. We assume that inactive Cdc42 in the cytosol only undergoes diffusion, with all the nonlinear binding kinetics taking place near the boundary. We also consider the lateral diffusion of active Cdc42 on the boundary, with a diffusion coefficient $D_c^m$. Let $C({\mathbf x}, t)$ denote the concentration of Cdc42 in the bulk, and $c({\mathbf x},t)$ denote the concentration of active Cdc42 on the membrane. These concentrations evolve according to a system of reaction-diffusion equations defined by
\begin{subequations}\label{eq:2Dmodel}
\begin{align}
\frac{\partial C({\mathbf x},t)}{\partial t}&=D_c\nabla^2C({\mathbf x},t),\quad {\mathbf x}\in\Omega,\\
\frac{\partial c({\mathbf x},t)}{\partial t}&=D_c^m\nabla_s^2 c({\mathbf x},t)+\mathcal{F}[C({\mathbf x},t), c({\mathbf x}, t), g({\mathbf x},t)],\quad {\mathbf x}\in\partial\Omega,\\
-D_c({\mathbf n}\cdot\nabla C)&=\mathcal{F}[C({\mathbf x},t), c({\mathbf x}, t), g({\mathbf x},t)],\quad {\mathbf x}\in\partial\Omega\,,
\end{align}
\end{subequations}
where $\nabla_s^2$ is the Laplace-Beltrami operator and ${\mathbf n}$ represents the outward unit normal vector to $\partial \Omega$. Similar reaction-diffusion equations are defined for GEF, the second species; see below.
\begin{subequations}\label{eq:consv2dMAIN}
\begin{align}
\frac{\partial G({\mathbf x},t)}{\partial t}&=D_g\nabla^2G({\mathbf x},t),\quad {\mathbf x}\in\Omega,\\
\frac{\partial g({\mathbf x},t)}{\partial t}&=D_g^m\nabla_s^2 g({\mathbf x},t)+\mathcal{G}[G({\mathbf x},t),c({\mathbf x},t),g({\mathbf x},t)],\quad {\mathbf x}\in\partial\Omega,\\
-D_g({\mathbf n}\cdot\nabla G)&=\mathcal{G}[G({\mathbf x},t),c({\mathbf x},t),g({\mathbf x},t)],\quad {\mathbf x}\in\partial\Omega.
\end{align}
\end{subequations}
The nonlinear functions $\mathcal{F}$ and $\mathcal{G}$ are as defined in \eqref{eq:kinetics}. Mass conservation arises from the formulation of the model, and we have the following conservation laws:
\begin{equation}\label{eq:consv2d}
\int_{\Omega} C({\mathbf x},t)d{\mathbf x}+\int_{\partial\Omega}c({\mathbf x},t)d{\mathbf x} = C_{\rm tot}\,, \quad \int_{\Omega} G({\mathbf x},t)d{\mathbf x}+\int_{\partial\Omega}g({\mathbf x},t)d{\mathbf x}=G_{\rm tot}\,.
\end{equation}
In Appendix \ref{subsec:appendix_mass}, we demonstrate that $C_{\rm tot}$ and $G_{\rm tot}$ must remain constant and show that the numerical method employed preserves mass over the timescale of cell polarization.

\subsection{Nonlocal PDE model on a circular domain}

Models of polarization often consider the cytosolic forms of Rho GTPases to be well-mixed \cite{goryachev2008, lo2014}.
We conclude this section with the introduction of a reduced nonlocal reaction-diffusion model governing the dynamics of membrane-bound species in the limit of fast bulk diffusion  $D \to\infty$. For simplicity, we consider a two-dimensional bulk domain $\Omega\in\mathbb{R}^2$  to be  a disk of radius $R$ and obtain a nonlocal reaction-diffusion model for membrane-bound species on a circle. In the well-mixed bulk, $C(t)$ and $G(t)$ are the spatially homogeneous concentrations of inactive Cdc42 and GEF, which are defined  using the conservation laws as
\begin{subequations}\label{eq:consv}
\begin{align}
C(t)&\equiv\frac{C_{\rm tot}}{|\Omega|}-\frac{1}{|\Omega|}\int_{\partial\Omega}c({\mathbf x},t)d{\mathbf x}=C_{\text{avg}}-\frac{1}{\pi R}\int_0^{2\pi}c(\theta,t)d\theta\,,\\
G(t)&\equiv\frac{G_{\rm tot}}{|\Omega|}-\frac{1}{|\Omega|}\int_{\partial\Omega}g({\mathbf x},t)d{\mathbf x}=G_{\text{avg}}-\frac{1}{\pi R}\int_0^{2\pi}g(\theta,t)d\theta\,,
\end{align}
\end{subequations}
where $C_{\text{avg}} = C_{\rm tot}/(\pi R^2)$ and $G_{\text{avg}} = G_{\rm tot}/(\pi R^2)$ are the average concentrations. Upon using polar coordinates to parametrize the Laplace-Beltrami operator, we write $\nabla_s^2 = R^{-2}\partial_{\theta\theta}$, with $\theta$ the angular variable, and find that the active species $c(\theta,t)$ and $g(\theta,t)$ must satisfy
\begin{subequations}\label{eq:nonlocal}
\begin{align}
\frac{\partial c(\theta,t)}{\partial t}&=\frac{D_c^m}{R^2}\frac{\partial^2 c(\theta,t)}{\partial \theta^2}+\mathcal{F}(C(t),c(\theta,t),g(\theta,t))\,, \quad 0\leq \theta < 2\pi\,,\\
\frac{\partial g(\theta,t)}{\partial t}&=\frac{D_g^m}{R^2}\frac{\partial^2 g(\theta,t)}{\partial \theta^2}+\mathcal{G}(G(t),c(\theta,t),g(\theta,t))\,, \quad 0\leq \theta < 2\pi\,.
\end{align}
\end{subequations}

\section{Results}

\subsection{One-dimensional PDE--ODE model}\label{sec:PDE_ODE}

\subsubsection{Steady states}\label{subsec:steady_states}
We first consider the steady states of the coupled PDE--ODE system \eqref{eq:pde}--\eqref{eq:mass_cons}. In the bulk domain, the Cdc42 and GEF profiles at a steady state must satisfy the following boundary value problems:
\begin{equation}
C''(x) = 0\,, \quad G''(x) = 0\,, \quad C'(0) = C'(L) = G'(0) = G'(L) = 0\,, \quad 0 < x < L\,.
\end{equation}
Hence, $C(x) \equiv C_a$ and $G(x) \equiv G_a$ are spatially homogeneous. From the conservation laws \eqref{eq:mass_cons} we obtain
\begin{equation}
C_a = \frac{C_{\rm tot} - c_{1a} - c_{2a}}{L}\,, \quad G_a = \frac{G_{\rm tot} - g_{1a} - g_{2a}}{L}\,,
\end{equation}
with the subscript $a$ referring to potential asymmetry $c_{1a} \neq c_{2a}$ and $g_{1a} \neq g_{2a}$ between both tips. Next, the active species steady state solutions satisfy the following system of nonlinear algebraic equations:
\begin{equation}\label{eq:steady_states}
\mathcal{F}(C_a,c_{1a},g_{1a}) = 0\,, \quad \mathcal{F}(C_a,c_{2a},g_{2a}) = 0\,, \quad \mathcal{G}(G_a,c_{1a},g_{1a}) = 0\,, \quad \mathcal{G}(G_a,c_{2a},g_{2a}) = 0\,,
\end{equation}
from which $g_{1a}$ and $g_{2a}$ are expressed as,
\begin{equation}
g_{1a} = \frac{L k^- c_{1a}}{(k_0+k_{\rm cat} c_{1a}^2)(C_{\rm tot} - c_{1a} - c_{2a})}\,, \quad g_{2a} = \frac{L k^- c_{2a}}{(k_0+k_{\rm cat} c_{2a}^2)(C_{\rm tot} - c_{1a} - c_{2a})}\,,
\end{equation}
with the corresponding Cdc42 values satisfying a system of two polynomial equations given by
\begin{equation}
H(c_{1a},c_{2a}) = 0\,, \quad H(c_{2a},c_{1a}) = 0\,,
\end{equation}
where $H(c_1,c_2)$ is defined by
\begin{equation}
\begin{split}
H(c_1,c_2) &\equiv \frac{\konG_{\rm tot}}{L}(C_{\rm tot} - c_1 - c_2)(k_0 + k_{\rm cat} c_1^2)(k_0 + k_{\rm cat} c_2^2) - k^-k^{\rm on} c_1(k_0 + k_{\rm cat} c_2^2) - k^-k^{\rm on} c_2(k_0 + k_{\rm cat} c_1^2) \\ & - L k^{\rm off} k^- c_1(1+\kappa c_1^2) (k_0 + k_{\rm cat} c_2^2)\,.
\end{split}
\end{equation}
Assuming the same binding kinetics for both cell tips causes reflection symmetry of solutions with respect to the midpoint $L/2$. As a result, interchanging $(c_{1a},g_{1a})$ with $(c_{2a},g_{2a})$ also satisfies \eqref{eq:steady_states}.

The steady state is symmetric when both tips possess equal solution values $(c_s,g_s)$, with $g_s$ defined by
\begin{equation}
g_s = \frac{L k^- c_s}{(k_0+k_{\rm cat} c_s^2)(C_{\rm tot} - 2c_s)},
\end{equation}
while $c_s$ satisfies a third degree polynomial equation given by
\begin{equation}\label{eq:cubic}
\left( L^2k^-k^{\rm off}\kappa + 2k^{\rm on}\kcatG_{\rm tot} \right)c_s^3 - \konk_{\rm cat}\GtotC_{\rm tot} c_s^2 + \left( 2Lk^-k^{\rm on} + L^2 k^-k^{\rm off} + 2 k_0\konG_{\rm tot} \right)c_s - k^{\rm on} k_0\GtotC_{\rm tot} = 0.
\end{equation}
By Descartes' rule of signs, it can be shown that equation \eqref{eq:cubic} possesses either one or three positive real roots, and therefore there always exists a symmetric steady state. Finally, we point out that the existence of the steady states does not depend on the diffusivity level, and that the structure of solutions remains the same for the reduced ODE system \eqref{eq:well_mixed} in the well-mixed regime.

\subsubsection{Linear stability analysis}\label{subsec:linear_stability_analysis}

We now proceed to a linear stability analysis of the symmetric steady state, leaving the linear stability of a general asymmetric steady state to be addressed via numerical bifurcation analysis. Consider a perturbation of the symmetric steady state given by
\begin{equation}\label{eq:perturbation}
C(x,t) = C_s + \varphi(x)e^{\lambda t}\,, \quad G(x,t) = G_s + \psi(x)e^{\lambda t}\,, \quad
\begin{pmatrix} c_1(t) \\ g_1(t) \end{pmatrix} = \begin{pmatrix} c_s \\ g_s \end{pmatrix} + \begin{pmatrix} u_1 \\ v_1 \end{pmatrix} e^{\lambda t}\,, \quad
\begin{pmatrix} c_2(t) \\ g_2(t) \end{pmatrix} = \begin{pmatrix} c_s \\ g_s \end{pmatrix} + \begin{pmatrix} u_2 \\ v_2 \end{pmatrix} e^{\lambda t}\,,
\end{equation}
where $\phi(x)$ and $\psi(x)$ are two eigenfunctions, while $(u_i,v_i)^T$ for $i=1,\,2$ are the eigenvectors of the linearized dynamics of active species at the left and right tips. Applying then the conservation laws \eqref{eq:mass_cons} yields
\begin{equation}\label{eq:mass_cons_eig}
u_1 + u_2 + \int_0^L \varphi(x) dx = 0, \quad v_1 + v_2 + \int_0^L \psi(x) dx = 0\,,
\end{equation}
and thus any perturbations can only result in a redistribution of total mass between active and inactive species. Next, upon inserting \eqref{eq:perturbation} within the PDE--ODE system \eqref{eq:pde}--\eqref{eq:ode} and linearizing about the symmetric steady state, we obtain an eigenvalue problem
\begin{equation}\label{eq:bvp}
0 = D_c \varphi_{xx} - \lambda \varphi\,, \quad 0 = D_g \psi_{xx} - \lambda \psi\,, \quad 0 < x < L\,, \\
\end{equation}
subject to boundary conditions at $x=0$ given by
\begin{equation}\label{eq:BC_1_eig}
D_c \varphi_x(0) = \frac{\partial\mathcal{F}}{\partial C} \varphi(0) + \frac{\partial\mathcal{F}}{\partial c} u_1 + \frac{\partial\mathcal{F}}{\partial g} v_1\,, \quad D_g \psi_x(0) = \frac{\partial\mathcal{G}}{\partial G} \psi(0) + \frac{\partial\mathcal{G}}{\partial c} u_1 + \frac{\partial\mathcal{G}}{\partial g} v_1\,,
\end{equation}
while on $x=L$ the boundary conditions are
\begin{equation}\label{eq:BC_2_eig}
D_c \varphi_x(L) = - \frac{\partial\mathcal{F}}{\partial C} \varphi(L) - \frac{\partial\mathcal{F}}{\partial c} u_2 - \frac{\partial\mathcal{F}}{\partial g} v_2\,, \quad D_g \psi_x(L) = - \frac{\partial\mathcal{G}}{\partial G} \psi(L) - \frac{\partial\mathcal{G}}{\partial c} u_2 - \frac{\partial\mathcal{G}}{\partial g} v_2\,.
\end{equation}
Finally, from the linearized ODE dynamics we find that the eigenvectors $(u_i,v_i)^T$ satisfy
\begin{equation}\label{eq:ode_eig}
\begin{pmatrix} \lambda - \frac{\partial\mathcal{F}}{\partial c} & - \frac{\partial\mathcal{F}}{\partial g} \\ - \frac{\partial\mathcal{G}}{\partial c} & \lambda - \frac{\partial\mathcal{G}}{\partial g} \end{pmatrix} \begin{pmatrix} u_1 \\ v_1 \end{pmatrix} = \begin{pmatrix} \frac{\partial\mathcal{F}}{\partial C} \varphi(0) \\ \frac{\partial\mathcal{G}}{\partial G} \psi(0) \end{pmatrix}\,, \qquad
\begin{pmatrix} \lambda - \frac{\partial\mathcal{F}}{\partial c} & - \frac{\partial\mathcal{F}}{\partial g} \\ - \frac{\partial\mathcal{G}}{\partial c} & \lambda - \frac{\partial\mathcal{G}}{\partial g} \end{pmatrix} \begin{pmatrix} u_2 \\ v_2 \end{pmatrix} = \begin{pmatrix} \frac{\partial\mathcal{F}}{\partial C} \varphi(L) \\ \frac{\partial\mathcal{G}}{\partial G} \psi(L) \end{pmatrix}\,.
\end{equation}
All the partial derivatives in equations \eqref{eq:BC_1_eig}--\eqref{eq:ode_eig} are evaluated at the symmetric steady state.

\subsubsection{Hopf bifurcation}\label{subsubsec:cond_hopf}
We first suppose that $\lambda \neq 0$. Hence, the only possible bifurcating eigenvalues are purely imaginary and lead to a Hopf bifurcation. Then, upon solving equation \eqref{eq:bvp}, we find that the eigenfunctions may be written as the sum of an even $(+)$ and an odd $(-)$ parts,
\begin{equation}\label{eq:phi_psi}
\varphi(x) = \varphi_+^0 \frac{\cosh\left(\omega_c\left(\frac{L}{2}-x\right)\right)}{\cosh\left(\omega_c \frac{L}{2}\right)} + \varphi_-^0 \frac{\sinh\left(\omega_c\left(\frac{L}{2}-x\right)\right)}{\sinh\left(\omega_c \frac{L}{2}\right)}\,, \quad
\psi(x) = \psi_+^0 \frac{\cosh\left(\omega_g\left(\frac{L}{2}-x\right)\right)}{\cosh\left(\omega_g \frac{L}{2}\right)} + \psi_-^0 \frac{\sinh\left(\omega_g\left(\frac{L}{2}-x\right)\right)}{\sinh\left(\omega_g \frac{L}{2}\right)}\,,
\end{equation}
where $\omega_c = \sqrt{\lambda/D_c}$ and $\omega_g = \sqrt{\lambda/D_g}$. Next, the application of the boundary conditions \eqref{eq:BC_1_eig} and \eqref{eq:BC_2_eig} yields the following expressions for the even coefficients $\varphi_+^0$ and $\psi_+^0$:
\begin{equation}
\varphi_+^0 = -\frac{1}{2}\frac{\frac{\partial\mathcal{F}}{\partial c}(u_1+u_2) + \frac{\partial\mathcal{F}}{\partial g}(v_1+v_2)}{D_c\omega_c \tanh\left(\omega_c\frac{L}{2}\right) + \frac{\partial\mathcal{F}}{\partial C}}\,, \quad
\psi_+^0 = -\frac{1}{2}\frac{\frac{\partial\mathcal{G}}{\partial c}(u_1+u_2) + \frac{\partial\mathcal{G}}{\partial g}(v_1+v_2)}{D_g\omega_g \tanh\left(\omega_g\frac{L}{2}\right) + \frac{\partial\mathcal{G}}{\partial G}}\,,
\end{equation}
while the odd coefficients $\varphi_-^0$ and $\psi_-^0$ are given by
\begin{equation}
\varphi_-^0 = -\frac{1}{2}\frac{\frac{\partial\mathcal{F}}{\partial c}(u_1-u_2) + \frac{\partial\mathcal{F}}{\partial g}(v_1-v_2)}{D_c\omega_c \coth\left(\omega_c\frac{L}{2}\right) + \frac{\partial\mathcal{F}}{\partial C}}\,, \quad
\psi_-^0 = -\frac{1}{2}\frac{\frac{\partial\mathcal{G}}{\partial c}(u_1-u_2) + \frac{\partial\mathcal{G}}{\partial g}(v_1-v_2)}{D_g\omega_g \coth\left(\omega_g\frac{L}{2}\right) + \frac{\partial\mathcal{G}}{\partial G}}\,.
\end{equation}
Since any solutions of the linearized system can be decomposed into an odd and an even parts, we treat the two different modes separately.

\paragraph{Odd mode $(-)$.} For the odd mode, we have that $u_2 = - u_1$, $v_2 = - v_1$ and thus
\begin{equation}\label{eq:phi_psi_odd}
\varphi(x) = -\dfrac{\left(\frac{\partial\mathcal{F}}{\partial c} u_1 + \frac{\partial\mathcal{F}}{\partial g} v_1\right)\sinh\left(\omega_c\left(\frac{L}{2}-x\right)\right)}{\left(D_c\omega_c\coth\left(\omega_c\frac{L}{2}\right) + \frac{\partial\mathcal{F}}{\partial C}\right)\sinh\left(\omega_c\frac{L}{2}\right)}\,, \quad
\psi(x) = -\dfrac{\left(\frac{\partial\mathcal{G}}{\partial c} u_1 + \frac{\partial\mathcal{G}}{\partial g} v_1\right)\sinh\left(\omega_g\left(\frac{L}{2}-x\right)\right)}{\left(D_g\omega_g\coth\left(\omega_c\frac{L}{2}\right) + \frac{\partial\mathcal{G}}{\partial G}\right)\sinh\left(\omega_g\frac{L}{2}\right)}\,,
\end{equation}
which automatically satisfies the constraints \eqref{eq:mass_cons_eig}. We then substitute \eqref{eq:phi_psi_odd} within \eqref{eq:ode_eig} and find that $(u_1,\,v_1)^T$ must be a non-trivial solution of
\begin{equation}
\begin{pmatrix}
\lambda - \dfrac{\frac{\partial\mathcal{F}}{\partial c} D_c \omega_c\coth\left(\omega_c\frac{L}{2}\right)}{D_c\omega_c\coth\left(\omega_c\frac{L}{2}\right) + \frac{\partial\mathcal{F}}{\partial C}} &
- \dfrac{\frac{\partial\mathcal{F}}{\partial g} D_c\omega_c\coth\left(\omega_c\frac{L}{2}\right)}{D_c\omega_c\coth\left(\omega_c\frac{L}{2}\right) + \frac{\partial\mathcal{F}}{\partial C}} \\[5pt]
- \dfrac{\frac{\partial\mathcal{G}}{\partial c} D_g\omega_g\coth\left(\omega_g\frac{L}{2}\right)}{D_g\omega_g\coth\left(\omega_g\frac{L}{2}\right) + \frac{\partial\mathcal{G}}{\partial G}} &
\lambda - \dfrac{\frac{\partial\mathcal{G}}{\partial g} D_g\omega_g\coth\left(\omega_g\frac{L}{2}\right)}{D_g\omega_g\coth\left(\omega_c\frac{L}{2}\right) + \frac{\partial\mathcal{G}}{\partial G}}
\end{pmatrix}
\begin{pmatrix} u_1 \\ v_1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\,,
\end{equation}
which, upon dividing through $\lambda \neq 0$, can be rearranged into
\begin{equation}
\begin{pmatrix}
\lambda + \frac{\omega_c}{\coth\left(\omega_c \frac{L}{2}\right)}\frac{\partial\mathcal{F}}{\partial C} - \frac{\partial\mathcal{F}}{\partial c} & - \frac{\partial\mathcal{F}}{\partial g} \\
- \frac{\partial\mathcal{G}}{\partial c} & \lambda + \frac{\omega_g}{\coth\left(\omega_g \frac{L}{2}\right)}\frac{\partial\mathcal{G}}{\partial G} - \frac{\partial\mathcal{G}}{\partial g} \\
\end{pmatrix}
\begin{pmatrix} u_1 \\ v_1\end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\,.
\end{equation}
The eigenvalue parameter $\lambda$ is therefore a root of the following characteristic equation:
\begin{equation}\label{eq:trans_odd}
 F_-(\lambda) \equiv \left[ \lambda + \frac{\omega_c}{\coth\left(\omega_c \frac{L}{2}\right)}\frac{\partial\mathcal{F}}{\partial C} - \frac{\partial\mathcal{F}}{\partial c} \right]\left[\lambda + \frac{\omega_g}{\coth\left(\omega_g \frac{L}{2}\right)}\frac{\partial\mathcal{G}}{\partial G} - \frac{\partial\mathcal{G}}{\partial g}\right] - \frac{\partial\mathcal{F}}{\partial g}\frac{\partial\mathcal{G}}{\partial c} = 0\,.
\end{equation}
In Section \ref{sec:results1-D}, equation \eqref{eq:trans_odd} will be solved numerically for a Hopf bifurcation, i.e. for a pair of roots $\lambda = \pm i\lambda_I$ with $\lambda_I > 0$. In the well-mixed regime, corresponding to the infinite bulk diffusion limit $D_{c,\,g} \to \infty$, equation \eqref{eq:trans_odd} reduces to
\begin{equation}\label{eq:trans_odd_WM}
\lim_{D_{c,\,g} \to \infty} F_-(\lambda) = \left| \lambda I - J_- \right| = 0\,,
\end{equation}
where $I$ is the $2\times 2$ identity matrix and $J_-$ is a Jacobian matrix defined by
\begin{equation}\label{eq:jac_odd}
 J_- = \begin{pmatrix} \frac{\partial\mathcal{F}}{\partial c} & \frac{\partial\mathcal{F}}{\partial g} \\[5pt] \frac{\partial\mathcal{G}}{\partial c} & \frac{\partial\mathcal{G}}{\partial g} \end{pmatrix}\,.
\end{equation}
In computing the limit, we have employed the following approximation of the hyperbolic cotangent function for a small argument:
\begin{equation*}
\coth\left(\omega_{c,\,g}\frac{L}{2}\right) \approx \frac{2}{L\omega_{c,\,g}}\,, \qquad \text{with} \qquad \omega_{c,g} = \sqrt{\frac{\lambda}{D_{c,g}}} \ll 1\,.
\end{equation*}
In the well-mixed regime, necessary conditions for the symmetric steady state to undergo a Hopf bifurcation are given by
\begin{equation}
{\rm tr}(J_-) = \frac{\partial\mathcal{F}}{\partial c} + \frac{\partial\mathcal{G}}{\partial g} = 0\,, \quad \det(J_-) = \frac{\partial\mathcal{F}}{\partial c}\frac{\partial\mathcal{G}}{\partial g} - \frac{\partial\mathcal{G}}{\partial c}\frac{\partial\mathcal{F}}{\partial g} > 0\,.
\end{equation}

\paragraph{Even mode $(+)$.} For the even mode, we have that $u_1=u_2$, $v_1=v_2$ and thus
\begin{equation}\label{eq:phi_psi_even}
\varphi(x) = -\dfrac{\left(\frac{\partial\mathcal{F}}{\partial c} u_1 + \frac{\partial\mathcal{F}}{\partial g} v_1\right)\cosh\left(\omega_c\left(\frac{L}{2}-x\right)\right)}{\left(D_c\omega_c\tanh\left(\omega_c\frac{L}{2}\right) + \frac{\partial\mathcal{F}}{\partial C}\right)\cosh\left(\omega_c\frac{L}{2}\right)}\,, \quad
\psi(x) = -\dfrac{\left(\frac{\partial\mathcal{G}}{\partial c} u_1 + \frac{\partial\mathcal{G}}{\partial g} v_1\right)\cosh\left(\omega_g\left(\frac{L}{2}-x\right)\right)}{\left(D_g\omega_g\tanh\left(\omega_c\frac{L}{2}\right) + \frac{\partial\mathcal{G}}{\partial G}\right)\cosh\left(\omega_g\frac{L}{2}\right)}\,.
\end{equation}
Upon substituting \eqref{eq:phi_psi_even} within \eqref{eq:ode_eig}, we find that $(u_1,\,v_1)^T$ is a non-trivial solution of
\begin{equation}\label{eq:cond_even_lambda}
\begin{pmatrix}
\lambda - \dfrac{\frac{\partial\mathcal{F}}{\partial c} D_c\omega_c\tanh\left(\omega_c \frac{L}{2}\right)}{D_c\omega_c\tanh\left(\omega_c \frac{L}{2}\right) + \frac{\partial\mathcal{F}}{\partial C}} &
- \dfrac{\frac{\partial\mathcal{F}}{\partial g} D_c\omega_c\tanh\left(\omega_c \frac{L}{2}\right)}{D_c\omega_c\tanh\left(\omega_c \frac{L}{2}\right) + \frac{\partial\mathcal{F}}{\partial C}} \\[10pt]
- \dfrac{\frac{\partial\mathcal{G}}{\partial c} D_g\omega_g\tanh\left(\omega_g \frac{L}{2}\right)}{D_g\omega_g\tanh\left(\omega_g \frac{L}{2}\right) + \frac{\partial\mathcal{G}}{\partial G}} &
\lambda - \dfrac{\frac{\partial\mathcal{G}}{\partial g} D_g\omega_g\tanh\left(\omega_g \frac{L}{2}\right)}{D_g\omega_g\tanh\left(\omega_g \frac{L}{2}\right) + \frac{\partial\mathcal{G}}{\partial G}} \\
\end{pmatrix}
\begin{pmatrix} u_1 \\ v_1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\,.
\end{equation}
Similarly as for the odd mode, we can divide through $\lambda \neq 0$ in \eqref{eq:cond_even_lambda} to obtain
\begin{equation}
\begin{pmatrix}
\lambda + \frac{\omega_c}{\tanh\left(\omega_c \frac{L}{2}\right)}\frac{\partial\mathcal{F}}{\partial C} - \frac{\partial\mathcal{F}}{\partial c} & - \frac{\partial\mathcal{F}}{\partial g} \\[5pt]
- \frac{\partial\mathcal{G}}{\partial c} & \lambda + \frac{\omega_g}{\tanh\left(\omega_g \frac{L}{2}\right)}\frac{\partial\mathcal{G}}{\partial G} - \frac{\partial\mathcal{G}}{\partial g} \\
\end{pmatrix} \begin{pmatrix} u_1 \\ v_1\end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\,,
\end{equation}
which corresponds to the linear system obtained from applying the constraints \eqref{eq:mass_cons_eig} onto an even eigenfunction. The characteristic equation is therefore given by
\begin{equation}\label{eq:trans_even}
F_+(\lambda) = \left[ \lambda + \frac{\omega_c}{\tanh\left(\omega_c \frac{L}{2}\right)}\frac{\partial\mathcal{F}}{\partial C} - \frac{\partial\mathcal{F}}{\partial c} \right] \left[ \lambda + \frac{\omega_g}{\tanh\left(\omega_g \frac{L}{2}\right)}\frac{\partial\mathcal{G}}{\partial G} - \frac{\partial\mathcal{G}}{\partial g} \right] - \frac{\partial\mathcal{F}}{\partial g} \frac{\partial\mathcal{G}}{\partial c} = 0\,.
\end{equation}
Once again, in the limit of fast bulk diffusion, the conditions for a Hopf bifurcation associated with the even mode is easily derived. Upon taking the limit $D_{c,\,g} \to \infty$ in equation \eqref{eq:trans_even} and using the following approximation of the hyperbolic tangent function:
\begin{equation*}
\tanh\left(\omega_{c,\,g}\frac{L}{2}\right) \approx \omega_{c,\,g}\frac{L}{2}\,, \qquad \text{with} \quad \omega_{c,g}=\sqrt{\frac{\lambda}{D_{c,g}}} \ll 1\,,
\end{equation*}
we obtain a quadratic equation given by
\begin{equation}\label{eq:trans_even_WM}
\lim_{D_{c,\,g} \to \infty} F_+(\lambda) = \left| \lambda I - J_+ \right| = 0\,,
\end{equation}
where the Jacobian matrix $J_+$ is defined by
\begin{equation}\label{eq:jac_even}
J_+ = \begin{pmatrix}
-\frac{2}{L}\frac{\partial\mathcal{F}}{\partial C} + \frac{\partial\mathcal{F}}{\partial c} & \frac{\partial\mathcal{F}}{\partial g} \\
\frac{\partial\mathcal{G}}{\partial c} & - \frac{2}{L}\frac{\partial\mathcal{G}}{\partial G} + \frac{\partial\mathcal{G}}{\partial g}
\end{pmatrix} \,.
\end{equation}
In the well-mixed regime, necessary conditions for a Hopf bifurcation associated with the even mode are given by ${\rm tr}(J_+) = 0$ and $\det(J_+) > 0$. We will however see that no Hopf bifurcation of the even mode occurs for the nonlinear functions and the parameter values used in our study, and thus only anti-phase oscillations are observed in both the well-mixed and finite diffusion regimes. The absence of in-phase oscillations was  also noted in a PDE--DDE model with mass conservation \cite{xu2016}. Whether a mass-conserved PDE-ODE(DDE) model can have in-phase oscillations remains an open question.

\subsubsection{Zero-crossing eigenvalues}\label{subsubsec:zero_lambda}
We now investigate the occurrence of zero-crossing eigenvalues leading to either pitchfork or saddle-node bifurcations. We remark that a rigorous treatment of the eigenvalue problem defined in equations \eqref{eq:mass_cons_eig}--\eqref{eq:ode_eig} would require the introduction of a branch cut in the complex $\lambda\text{-plane}$ along $\infty < \Re(\lambda) \leq 0$ and $\Im(\lambda) = 0$. Fortunately, the limit of $\lambda \to0$, with $\lambda$ real, has the same effect as the infinite bulk diffusion limit: all the square roots from the characteristic equations \eqref{eq:trans_odd} and \eqref{eq:trans_even} cancel, and no continuous spectrum is introduced. Therefore, the solution for $\lambda = 0$ can be readily obtained by taking the well-defined limit $\lambda \to 0$.

\paragraph{Odd mode $(-)$.} After taking the limit $\lambda \to 0$ in equation \eqref{eq:trans_odd}, we obtain the following condition for the linearized system to have a zero-crossing eigenvalue associated with the odd mode:
\begin{equation}
\lim_{\lambda \to 0} F_-(\lambda) = \det(J_-) = \frac{\partial\mathcal{F}}{\partial c}\frac{\partial\mathcal{G}}{\partial g} - \frac{\partial\mathcal{G}}{\partial c}\frac{\partial\mathcal{F}}{\partial g} = 0.
\end{equation}
When this condition is satisfied, the system undergoes a pitchfork bifurcation with two branches of asymmetric steady states emerging from the symmetric steady state.

\paragraph{Even mode $(+)$.} We proceed similarly for the even mode. The symmetric steady state undergoes a saddle-node bifurcation when the following condition is satisfied:
\begin{equation}
\lim_{\lambda \to 0} F_+(\lambda) = \det(J_+) = \frac{4}{L^2}\frac{\partial\mathcal{F}}{\partial C}\frac{\partial\mathcal{G}}{\partial G} - \frac{2}{L}\left(\frac{\partial\mathcal{F}}{\partial c}\frac{\partial\mathcal{G}}{\partial G} + \frac{\partial\mathcal{G}}{\partial g}\frac{\partial\mathcal{F}}{\partial C}\right) + \left(\frac{\partial\mathcal{F}}{\partial c}\frac{\partial\mathcal{G}}{\partial g} - \frac{\partial\mathcal{G}}{\partial c}\frac{\partial\mathcal{F}}{\partial g}\right) = 0\,.
\end{equation}

We highlight that bulk diffusion has no effect on the existence of steady states and on bifurcations resulting from zero-crossing eigenvalues. However, as will be seen in the next section, diffusion does influence the stability of symmetric steady states and possible oscillatory dynamics.

\subsubsection{Numerical bifurcation analysis}\label{sec:results1-D}

Our previous work on the well-mixed version of the model pointed out the importance of membrane residence times of GEF proteins (set by $1/k^{\rm off}$) in determining whether Cdc42 proteins were at steady state or undergoing pole-to-pole oscillations \cite{xu2018}. Hence, we select $k^{\rm off}$ and the bulk diffusion coefficients as bifurcation parameters, and numerically investigate their effects on the stability of the symmetric steady state ($c_1=c_2, \, g_1=g_2$). Details of the numerical bifurcation analysis are given in Appendix~\ref{sec:numerical_bifurcation_analysis}.

We start with the dynamics of the ODE model (\ref{eq:ode}). The intricate succession of bifurcations that results when increasing $k^{\rm off}$ is presented in Fig.~\ref{fig:bif_diag_koff_WM}. The symmetric steady state is unstable for low $k^{\rm off}$ values, with linear stability gained through a supercritical Hopf bifurcation when $k^{\rm off} \approx 0.99$. Further beyond this threshold, two fold bifurcations in $k^{\rm off} \approx 6.23$ and $k^{\rm off} \approx 9.17$ cause hysteresis and coexistence of up to three symmetric steady states. We observe that bistability of symmetric steady states occurs in the range $(6.24, 7.86)$, with such thresholds corresponding to subcritical pitchfork bifurcations connecting branches of asymmetric and symmetric steady states. Finally, stable and unstable branches of asymmetric steady states are connected via fold bifurcations at $k^{\rm off} \approx 0.56$ and $k^{\rm off} \approx 11.22$.

\begin{figure}[H]
\centering
\begin{subfigure}{0.36\linewidth}
\includegraphics[width=\linewidth]{figures/coco_c12_WM_koff}
\caption{Active Cdc42 levels.}
\end{subfigure}
\begin{subfigure}{0.36\linewidth}
\includegraphics[width=\linewidth]{figures/coco_g12_WM_koff}
\caption{Active GEF levels.}
\end{subfigure}
\caption{\label{fig:bif_diag_koff_WM} Symmetric steady state (\textcolor{red}{red}), asymmetric steady state (\textcolor{black}{black}) and periodic solution branches (\textcolor{green}{green}) for the reduced ODE model in the well-mixed regime, as a function of $k^{\rm off}$. Other parameter values are given in Table~\ref{table:parm}. A full line indicates linear stability, while unstable branches are dashed.}
\end{figure}

For the parameter regime considered, stable anti-phase oscillations are observed in the limit of fast bulk diffusion when $k^{\rm off}$ is small ($k^{\rm off} < 1$). A closer view of this regime is shown in panel (a) of Fig.~\ref{fig:hopf_WM}, where we observe that oscillations emerge either via supercritical Hopf or homoclinic bifurcations. Numerical simulations exhibiting the role of the asymmetric equilibrium on the shape of the oscillations are shown in panel (c) of the same figure. From $k^{\rm off} = 0.9$ to $k^{\rm off} = 0.5$, we observe a transition from weakly to highly nonlinear oscillations. The intermediate range $k^{\rm off} \, \in \, (0.56, 0.79)$ is characterized by the absence of stable oscillations, with all the trajectories attracted to the asymmetric steady states.

\begin{figure}[H]
\centering
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/AUTO_g12_WM_edited}
\caption{Active GEF levels.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/AUTO_period_WM}
\caption{Oscillatory period.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/g1_WM}
\caption{$g_1(t)$ for $k^{\rm off} = 0.9,\, 0.7, \, 0.5$.}
\end{subfigure}
\caption{\label{fig:hopf_WM} Oscillatory dynamics for the reduced ODE model (\ref{eq:well_mixed}) in the well-mixed regime. (a) A zoomed-in view of the bifurcation diagram in Fig.~\ref{fig:bif_diag_koff_WM}(b). (b) Period of oscillations in the $k^{\rm off}$ regimes $(0, 0.56)$ and $(0.79, 0.99)$. For the first branch, the period has a minimum near $k^{\rm off} \approx 0.3$ and rapidly increases near both $k^{\rm off}=0$ and $k^{\rm off} \approx 0.56$, suggesting the occurrence of homoclinic bifurcations as limit cycles collide with asymmetric equilibria. For the branch of periodic solutions emerging from a supercritical Hopf bifurcation at $k^{\rm off} \approx 0.99$, the period increases rapidly at a suspected homoclinic bifurcation near $k^{\rm off} \approx 0.79$. (c) Numerical solutions for $k^{\rm off}=0.9, 0.7, 0.5$. For a high rate $k^{\rm off}=0.9$, weakly nonlinear oscillations are observed. For an intermediate rate $k^{\rm off}=0.7$, numerical solutions evolve to an asymmetric steady state. As $k^{\rm off}$ further decreases to $0.5$, highly nonlinear relaxation oscillations are observed. For each of the three simulations, the initial conditions correspond to the symmetric steady state (red horizontal line) plus a small perturbation stimulating the odd mode.}
\end{figure}

We consider next the effect of finite bulk diffusion on the emergence of oscillatory dynamics, focusing on the parameter regime where there is a unique symmetric steady state. We numerically solve the eigenvalue relation \eqref{eq:trans_odd} for a Hopf bifurcation associated with the anti-phase mode, i.e. for a pair of roots $\lambda = \pm i\lambda_I$ with $\lambda_I > 0$. The corresponding Hopf stability boundaries in the $D_c$ versus $D_g$ parameter plane are given in panel (a) of Fig.~\ref{fig:Dg_Dc}, where the region of linear stability of the symmetric steady state is located below the curves. We point out that fast cytosolic diffusion is sufficient to achieve oscillations and that decreasing $k^{\rm off}$ shrinks the linear stability region. Since $1/k^{\rm off}$ sets the membrane dwelling time for the GEF, we conclude that increasing this quantity pushes the system into oscillatory dynamics even for relatively modest diffusion coefficients. Numerical simulations exhibiting anti-phase oscillations with qualitatively different shapes are shown in panel (b) of the same figure.

\begin{figure}[H]
\centering
\begin{subfigure}{0.36\linewidth}
\includegraphics[width=\linewidth]{figures/stability_diagram_model2_Dg_Dc}
\caption{$(D_g,\, D_c)$ with $k^{\rm off} = 0.5,\, 0.9$.}
\end{subfigure}
\begin{subfigure}{0.36\linewidth}
\includegraphics[width=\linewidth]{figures/c_Dg1_Dc20}
\caption{$D_g = 1,\, D_c = 20$.}
\end{subfigure}
\caption{\label{fig:Dg_Dc} Panel (a): Hopf stability boundaries in the $D_c$ versus $D_g$ plane in the full PDE--ODE model for $k^{\rm off} = 0.9$ (\textcolor{black}{black}) and $k^{\rm off} = 0.5$ (\textcolor{red}{red}). Other parameter values are given in Table~\ref{table:parm}. Linear stability regions are located below each curve. Panel (b): for $D_c = 20$ and $D_g = 1$, increasing the membrane residence time results in a transition from weakly to highly nonlinear oscillations, a finding consistent with the dynamics in the well-mixed regime. The phase shift of half a period between $c_1(t)$ (\textcolor{blue}{blue}) and $c_2(t)$ (\textcolor{red}{red}) is a clear indication of anti-phase oscillations.}
\end{figure}

Shown next in Fig.~\ref{fig:stability_boundaries} are stability diagrams for the symmetric steady state that use $k^{\rm off}$ as a bifurcation parameter. We note that in panel (a), the bulk diffusion coefficients are set to be equal, while in panels (b) and (c) only $D_c$ or $D_g$ is allowed to vary. The branching behavior of the Hopf bifurcation, i.e. whether it is super- or subcritical, was resolved numerically with the software AUTO (cf.~\cite{doedel2007}) and a nonlocal formulation of the PDE--ODE model (see Appendix~\ref{sec:numerical_bifurcation_analysis}). Letting $D \equiv D_c = D_g$ or letting $D_c$ free with $D_g$ fixed yields qualitatively similar stability boundaries, with the second case seemingly converging to the first one as $D_g$ increases. In both cases, the bifurcation is supercritical when the diffusivity level is above some threshold and it is subcritical in the low diffusion regime. However, stability boundaries in the $D_g$ versus $k^{\rm off}$ parameter plane have a qualitatively different shape, with the critical diffusivity threshold being inversely proportional to $k^{\rm off}$. Here, increasing $D_c$ tenfold pushes the system further into instability and near the onset of oscillations. Hence, we conclude that $D_c$ has a stronger effect on the oscillatory dynamics than $D_g$.

\begin{figure}[H]
\centering
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/AUTO_stabCurve_koff_D}
\caption{$(k^{\rm off}, D)$ with $D \equiv D_c = D_g$.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/AUTO_stabCurve_koff_Dc}
\caption{$(k^{\rm off}, D_c)$ with $D_g=1,\, 10$.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/AUTO_stabCurve_koff_Dg}
\caption{$(k^{\rm off}, D_g)$ with $D_c=1,\, 10$.}
\end{subfigure}
\caption{\label{fig:stability_boundaries} Hopf stability boundaries for the symmetric steady state in the full PDE--ODE model in the $D$ versus $k^{\rm off}$ (panel (a)), $D_c$ versus $k^{\rm off}$ (panel (b)) and $D_g$ versus $k^{\rm off}$ (panel (c)) parameter planes. Full and dashed curves indicate respectively super- and subcritical Hopf bifurcations, while triangular points indicate a criticality change (generalized Hopf point). Other parameter values are given in Table~\ref{table:parm}.}
\end{figure}

Next, we compute global bifurcation diagrams along two horizontal slices from panel (a) of Fig.~\ref{fig:stability_boundaries}. The results are shown in Fig.~\ref{fig:hopf_finite_diffusion} below for $D=10$ (upper panels (a)-(c)) and $D=1$ (lower panels (d)-(f)). When $D=10$, the symmetric steady state loses stability through a supercritical Hopf bifurcation at $k^{\rm off} \approx 0.81$. Further below this threshold, the stable weakly nonlinear anti-phase oscillations are annihilated by the asymmetric steady state at $k^{\rm off} \approx 0.71$, with the oscillatory period undergoing a sharp increase to infinity. As in the well-mixed regime, relaxation oscillations exist in the range $k^{\rm off} \, \in \, (0,0.56)$.

As predicted by panel (a) of Fig.~\ref{fig:stability_boundaries}, we find that the symmetric steady state undergoes a subcritical Hopf bifurcation on the horizontal slice $D=1$, near $k^{\rm off} \approx 0.17$. Furthermore, the branch of limit cycles emerging from the Hopf bifurcation gains stability at a fold point near $k^{\rm off} \approx 0.27$. Finally, panel (f) of Fig.~\ref{fig:hopf_finite_diffusion} illustrates a typical hard loss of stability near a subcritical Hopf bifurcation. Notice the absence of weakly nonlinear oscillations and the direct transition to relaxation oscillations.

\begin{figure}[H]
\centering
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/AUTO_g12_D10}
\caption{Active GEF levels, $D=10$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/AUTO_period_D10}
\caption{Oscillatory period, $D=10$.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/g1_D10}
\caption{$D=10$ and $k^{\rm off} = 0.75,\, 0.6,\, 0.5$.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/AUTO_g12_D1}
\caption{Active GEF levels, $D=1$.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/AUTO_period_D1}
\caption{Oscillatory period, $D=1$.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/g1_D1}
\caption{$D=1$ and $k^{\rm off} = 0.3,\, 0.2,\, 0.1$.}
\end{subfigure}
\caption{\label{fig:hopf_finite_diffusion} Oscillatory dynamics in the full PDE--ODE model for $D_c = D_g = 10$ (panels (a)-(c)) and for $D_c = D_g = 1$ (panels (d)-(f)). Other parameter values are given in Table~\ref{table:parm}. The case $D=10$ yields qualitatively similar dynamics to that of the well-mixed regime, with the Hopf bifurcation being supercritical. However, when $D=1$ the Hopf bifurcation is subcritical and a hard loss of stability of the symmetric steady state is observed; see panel (f).}
\end{figure}

We conclude this section with Fig.~\ref{fig:pole_to_pole}, which gives some explanations behind the emergence of anti-phase relaxation oscillations in the limit of long GEF membrane residence time. Relaxation oscillations are characterized by long periods of rest followed by sudden sharp variations and are typical of slow-fast systems. Since $k^{\rm off} < k^-$, GEF dissociation is slow compared to Cdc42 dissociation, and active GEF can thus be seen as a slow variable while active Cdc42 as a fast variable. This results in a more rapid flux on (or off) each pole for Cdc42 (see panel (c)).

\begin{figure}[H]
\centering
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/C_bulk_Dc1_Dg1_koff1e-01}
\caption{Inactive Cdc42 $C(x,t)$.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/phase_plane_Dc1_Dg1_koff1e-01}
\caption{Phase plane trajectories.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/derivative_Dc1_Dg1_koff1e-01}
\caption{Flux on/off each pole.}
\end{subfigure}
\caption{\label{fig:pole_to_pole} Relaxation oscillations in the full PDE--ODE model for $D_c = D_g = 1$, $k^{\rm off} = 0.1$ and other parameter values as in Table~\ref{table:parm}. Panel (a): kymograph of inactive Cdc42 level, with time on the horizontal axis and space on the vertical axis (see Supplemental Movie S1). Panel (b): trajectories in the phase plane $(c_i(t),g_i(t))$. Panel (c): Cdc42 (blue) and GEF (red) fluxes on/off each pole as a function of time.}
\end{figure}

\subsection{Nonlocal PDE model on a circular domain}\label{sec:nonlocal}
In this section, we consider a two-dimensional membrane-bulk model with circular geometry and well-mixed bulk dynamics. This results in the reduced nonlocal PDE model given by system~\eqref{eq:nonlocal}. Our approach combines linear stability analysis and numerical simulations to predict and explore the formation of spatio-temporal patterns near stability thresholds and away from them, in the highly nonlinear regime.

We first define the spatially uniform (or trivial) steady state as $C(r,\theta) \equiv C^*,\, G(r,\theta) \equiv G^*,\, c(\theta) \equiv c^*$ and $g(\theta) \equiv g^*$. Similarly to the one-dimensional PDE--ODE model, $C^*,\, G^*$ and $g^*$ satisfy
\begin{equation*}
C^*=C_{\text{avg}}-\gamma c^*\,, \quad G^*=G_{\text{avg}}-\gamma g^*\,, \quad g^*=\frac{k^- c^*}{[k_0 + k_{\rm cat} (c^*)^2](C_{\text{avg}}-\gamma c^*)}\,,
\end{equation*}
where $G_{\text{avg}}=G_{\rm tot}/|\Omega|$, $C_{\text{avg}}=C_{\rm tot}/|\Omega|$ and $\gamma = |\partial\Omega|/|\Omega|$ are respectively the average masses and the ratio of the perimeter over the area. The steady state concentration of active Cdc42 is then found to satisfy the following cubic equation:
\begin{equation}\label{eq:ss2D}
[k^{\rm off}\kappa k^- +\gamma\konk_{\rm cat}G_{\text{avg}}](c^*)^3-\konk_{\rm cat}\GavgC_{\text{avg}}(c^*)^2 + [k^{\rm off} k^-+k^{\rm on}\gamma(k^-+k_0G_{\text{avg}})]c^*-k^{\rm on} k_0\GavgC_{\text{avg}}=0,
\end{equation}
where $\gamma=2/R$ for a circle of radius $R$. In the one-dimensional case, the average masses are defined as $C_{\text{avg}}=C_{\rm tot}/L$ and $G_{\text{avg}}=G_{\rm tot}/L$, and equation \eqref{eq:ss2D} reduces to equation \eqref{eq:cubic} derived for the symmetric steady state of the PDE--ODE model.

\subsubsection{Linear stability analysis}
We consider the following nonuniform, spatially periodic, perturbations of the trivial steady state:
\begin{equation}\label{eq:pertb}
 c(\theta, t)=c^* + u_n{\rm e}^{\lambda t + in\theta}\,, \quad g(\theta, t)=g^* + v_n{\rm e}^{\lambda t + in\theta}\,, \quad n\neq 0\,,
\end{equation}
where once again $\lambda \in \mathbb{C}$ is the eigenvalue parameter. After substitution of \eqref{eq:pertb} within the system of nonlocal PDEs \eqref{eq:nonlocal} and upon linearizing around the trivial steady state, we find that the eigenvector $(u_n,v_n)^T$ must satisfy
\begin{equation}\label{eq:eigen_nonlocal}
\left[ \lambda I - J_n\right] \begin{pmatrix} u_n \\ v_n \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix},
\end{equation}
where the Jacobian matrix $J_n$ is defined by
\begin{equation}\label{eq:J_n}
J_n = \begin{pmatrix} - D_c^m\frac{n^2}{R^2} + \frac{\partial\mathcal{F}}{\partial c} & \frac{\partial\mathcal{F}}{\partial g} \\ \frac{\partial\mathcal{G}}{\partial c} & - D_g^m\frac{n^2}{R^2} + \frac{\partial\mathcal{G}}{\partial g} \end{pmatrix} =
\begin{pmatrix} -D_c^mn^2/R^2+2k_{\rm cat} c^*g^*C^*-k^- &C^*[k_0+k_{\rm cat} (c^*)^2]\\[5pt] -{2\kappak^{\rm on} c^*G^*}/{[1+\kappa(c^*)^2]^2} &-D_g^mn^2/R^2-k^{\rm off} \end{pmatrix}.
\end{equation}
For each mode $n$, the corresponding pair of eigenvalues is given by
\begin{equation*}
\lambda_{\pm}(n)=\frac{{\rm tr} J_n\pm\sqrt{{\rm tr}^2J_n-4\det J_n}}{2},
\end{equation*}
where the trace and the determinant are
\begin{equation*}
{\rm tr} J_n=-n^2D_c^m/R^2-n^2D_g^m/R^2+2k_{\rm cat} c^*g^*C^*-(k^{-}+k^{\rm off})\,, \quad \det J_n=[n^2D_c^m/R^2-2k_{\rm cat} c^*g^*C^*+k^-][n^2D_g^m/R^2+k^{\rm off}]+\eta\,,
\end{equation*}
with $\eta$ defined as
\begin{equation}\label{eq:eta}
\eta = \frac{2\kappak^{\rm on} c^*[k_0+k_{\rm cat} (c^*)^2]C^*G^*}{ {[1+\kappa(c^*)^2]^2}}.
\end{equation}
As usual, linear stability requires a negative trace and a positive determinant of the Jacobian matrix for each $n$, therefore yielding the following inequality to be satisfied by $D_c^m$:
\begin{align}\label{cond:nonlocal}
&D_c^m > \max\left\{-D_g^m+\frac{R^2}{n^2}(2k_{\rm cat} c^*g^*C^*-k^{-}-k^{\rm off})\,, -\frac{\eta R^2/n^2}{D_g^mn^2/R^2+k^{\rm off}}+\frac{R^2}{n^2}(2k_{\rm cat} c^*g^*C^*-k^-)\right\}.
\end{align}
This inequality places a bound on the minimal diffusion coefficient of Cdc42 for the linear stability of the trivial steady state with respect to non-spatially homogeneous perturbations. This boundary is plotted in Fig.~\ref{fig:stb_nonlocal} and verified numerically.

Notice that because of the integral terms within the nonlocal PDE system \eqref{eq:nonlocal}, we cannot simply set $n=0$ in \eqref{eq:J_n} to recover the Jacobian matrix associated with spatially uniform perturbations. A similar observation was made by R{\"a}tz and R{\"o}ger in their stability analysis of a coupled membrane-bulk reaction-diffusion model for polarity of small GTPases \cite{ratz2012turing}. Furthermore, stability with respect to spatially uniform perturbations, or equivalently, in the absence of membrane diffusion, is required for Turing diffusion-driven instability (cf.~\cite{turing1952}). For this special case, the Jacobian matrix of the linearized system satisfies
\begin{equation}\label{eq:J_0}
J_0 = \begin{pmatrix} -\frac{\partial\mathcal{F}}{\partial C}\gamma + \frac{\partial\mathcal{F}}{\partial c} & \frac{\partial\mathcal{F}}{\partial g} \\ \frac{\partial\mathcal{G}}{\partial c} & -\frac{\partial\mathcal{G}}{\partial G}\gamma + \frac{\partial\mathcal{G}}{\partial g} \end{pmatrix}
\end{equation}
and necessary conditions for diffusion-driven instability are given by
\begin{equation}\label{eq:stability_condition}
{\rm tr} J_0 = \frac{\partial\mathcal{F}}{\partial c}+\frac{\partial\mathcal{G}}{\partial g} - \gamma\left( \frac{\partial\mathcal{F}}{\partial C} + \frac{\partial\mathcal{G}}{\partial G} \right)<0,\quad \det J_0 = \gamma^2 \frac{\partial\mathcal{F}}{\partial C} \frac{\partial\mathcal{G}}{\partial G} - \gamma\left( \frac{\partial\mathcal{F}}{\partial c} \frac{\partial\mathcal{G}}{\partial G} + \frac{\partial\mathcal{G}}{\partial g} \frac{\partial\mathcal{F}}{\partial C} \right) + \frac{\partial\mathcal{F}}{\partial c}\frac{\partial\mathcal{G}}{\partial g}-\frac{\partial\mathcal{F}}{\partial g}\frac{\partial\mathcal{G}}{\partial c}>0,
\end{equation}
which have been verified numerically for the parameter values employed in our study (see Appendix \ref{sec:appendix_nonlocal}). Restoring membrane diffusion, two different symmetry-breaking mechanisms may cause loss of stability of the trivial steady state: the stationary or the oscillatory Turing instabilities, which we describe below.
\begin{itemize}
\item The \textbf{stationary Turing instability}, which happens when the Jacobian matrix possesses an eigenvalue at the origin: $\det J_n = 0$, with $n \neq 0$. We note that because of circular symmetry (or $O(2)$ symmetry), an unstable mode always comes in a pair $\pm n$, and thus the stationary Turing instability corresponds to a pitchfork bifurcation.
\item The \textbf{oscillatory Turing instability}, which happens when the Jacobian matrix possesses a pair of critical eigenvalues on the imaginary axis: ${\rm tr} J_n = 0$ and $\det J_n > 0$, with $n \neq 0$. This instability corresponds to a Hopf bifurcation with $O(2)$ symmetry (cf.\cite{vangils1986}), and traveling waves as opposed to stationary patterns are expected to form.
\end{itemize}
More intricate bifurcations resulting from the interaction of two spatial modes may also be detected when performing linear stability analysis. One example is the Bogdanov--Takens bifurcation with $O(2)$ symmetry (cf.~\cite{dangelmayr1987}), which occurs when the stationary and oscillatory Turing instabilities coincide for the same (nonzero) spatial modes $\pm n$. At such a bifurcation point, we have ${\rm tr} J_n = 0$ and $\det J_n = 0$, and thus the linearized system possesses two zero-eigenvalues of multiplicity two.

\subsubsection{Spatio-temporal pattern formation in the nonlocal PDE model}\label{sec:nonlocal_stb}

Here, we compare our linear stability analysis against numerical simulations of the nonlocal system of PDEs \eqref{eq:nonlocal}. Simple finite differences are employed to spatially discretize diffusion terms while nonlocal integral terms are handled using the trapezoidal quadrature rule. As initial conditions, we take in this section
\begin{equation}
c(\theta,0) = c^* + u_n\cos\left(n\theta+\varphi\right)\,, \qquad g(\theta,0) = g^* + v_n\cos\left(n\theta+\varphi\right)\,,
\end{equation}
where the eigenvector $(u_n,\,v_n)^T$ is a non-trivial solution of \eqref{eq:eigen_nonlocal} and $\varphi$ is an arbitrary angular phase shift.

In Fig.~\ref{fig:stb_nonlocal}, we show three stability diagrams in the $(D_g^m, D_c^m)$ parameter plane for different $k^{\rm off}$ values. As expected, the trivial steady state is linearly stable when Cdc42 diffuses fast. Fig.~\ref{fig:stb_nonlocal} also confirms the destabilizing effect of long GEF membrane residence times, since the linear stability region shrinks as $k^{\rm off}$ decreases. In each diagram, the different spatio-temporal patterns observed via numerical simulations are labeled with a variety of symbols. For instance, the star and the filled circle respectively indicate a traveling wave and a stationary pattern, while the empty circle shows the absence of patterns. This non-exhaustive parameter exploration reveals that the spatio-temporal dynamics thresholds clearly correspond to the stability boundaries shown in Fig.~\ref{fig:stb_nonlocal}, suggesting that in the limit of fast bulk diffusion, all the identified bifurcations are supercritical.

\begin{figure}[H]
\centering
\includegraphics[width=0.9\linewidth]{figures/StabilityNonlocalPDEJun27}
\caption{\label{fig:stb_nonlocal} Stability boundaries for the nonlocal PDE model in the $(D_g^m, D_c^m)$ parameter plane with $k^{\rm off} = 1$ (panel (a)), $k^{\rm off}=0.5$ (panel (b)) and $k^{\rm off}=0.1$ (panel (c)) and other parameters the same as in Table \ref{table:parm}. For illustrative purposes, parameter values where we observe a stationary pattern with a single peak, a stationary pattern with two peaks, a traveling and a standing waves are labeled by a filled circle, a square, a star and a polygon. A stationary pattern that evolves into a traveling wave is indicated by a star superposed on a square. Some secondary stability boundaries associated with higher spatial modes have been omitted.}
\end{figure}

Some numerically computed spatio-temporal patterns obtained when $k^{\rm off} = 1$ are shown in Fig.~\ref{fig:WM_pattern}. For this case, the trivial steady state is expected to lose stability via pitchfork bifurcations associated with the modes $n=1,2$ and $3$ (indicated by the blue, red and golden curves). Nevertheless, our simulations revealed that some patterns that are initially stationary can evolve into traveling waves, and this phenomenon occurs in the absence of any oscillatory Turing instabilities (see panels (b), (d) and (e)). One potential explanation comes from the vicinity of multiple stability curves, which causes two or more spatial modes to interact and ultimately leads to a stationary pattern undergoing a secondary Hopf bifurcation. We also observed ``quasi-stationary'' patterns with multiple peaks. For instance in Fig.~\ref{fig:WM_pattern}(c), two peaks are seen to collapse and split again immediately after, resulting in a different pattern configuration. Finally, when the ratio $D_c^m/D_g^m$ is very small, a large number of modes become unstable and linear stability analysis cannot be employed to predict the final number of peaks of a stationary pattern. This phenomenon precisely corresponds to what happens in panel (f), where for $D_c^m/D_g^m = 0.1$, an initial pattern with three peaks rapidly evolves into a highly localized pattern with a single peak.

\begin{figure}[H]
\centering
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/WM_cs_TW_n1_Dcm2e-01_Dgm1_koff1_phi1e-02_tend1000}
\caption{$D_g^m = 1,\, D_c^m = 0.2$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/WM_cs_TW_n1_Dcm1p8e-01_Dgm1_koff1_phi1e-01_tend1000}
\caption{$D_g^m = 1,\, D_c^m = 0.18$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/WM_cs_TW_n2_Dcm1p1e-01_Dgm5e-01_koff1_phi1e-01_tend500}
\caption{$D_g^m = 0.5,\, D_c^m = 0.11$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/WM_cs_TW_n2_Dcm1e-01_Dgm5e-01_koff1_phi1e-01_tend500}
\caption{$D_g^m = 0.5,\, D_c^m = 0.1$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/WM_cs_TW_n3_Dcm4e-02_Dgm2e-01_koff1_phi1e-01_tend300}
\caption{$D_g^m = 0.2,\, D_c^m = 0.04$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/WM_cs_TW_n3_Dcm4e-02_Dgm4e-01_koff1_phi1e-01_tend300}
\caption{$D_g^m = 0.4,\, D_c^m = 0.04$}
\end{subfigure}
\caption{\label{fig:WM_pattern} A gallery of numerical simulations of the nonlocal model with $k^{\rm off} = 1$ and membrane diffusion coefficients taken from Fig.~\ref{fig:stb_nonlocal}(a). Shown above are kymographs of membrane-bound active Cdc42, with time on the horizontal axis and space on the vertical axis. Panel (a): weakly nonlinear stationary pattern near a mode $n=1$ pitchfork stability boundary (see Supplemental Movie S2). Panel (b): for a slightly lower $D_c^m$ value, the same pattern evolves into a traveling wave. Panel (c): quasi-stationary pattern with two peaks near a mode $n=2$ pitchfork stability boundary. Panel (d): transition to a traveling wave from an initial pattern with two peaks. Panel (e): transition to a traveling wave from an initial pattern with three peaks. Panel (f): transition to a localized stationary pattern with a single peak from an initial pattern with three peaks.}
\end{figure}

At smaller $k^{\rm off}$ values, as indicated by the black curve in panels (b-c) of Fig.\ref{fig:stb_nonlocal}, the loss of stability of the trivial steady state can arise via a Hopf bifurcation (see Appendix \ref{sec:appendix_nonlocal}). This is similar to what is seen in the one-dimensional PDE--ODE model from Section \ref{sec:PDE_ODE}, where oscillatory dynamics also require a small GEF dissociation rate. Near those oscillatory Turing instabilities, the formation of traveling waves is expected. We distinguish between rotating and standing waves, with the latter type often referred to as pole-to-pole oscillations. As shown in Fig.~\ref{fig:WM_wave}, our numerical experiments revealed that the two types coexist near the Hopf stability boundaries, while away from them in the highly nonlinear regime the standing waves are only transient.

\begin{figure}[H]
\centering
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/WM_cs_TW_n1_Dcm2p8e-01_Dgm2e-01_koff5e-01_phi1e-01_tend1000}
\caption{$D_g^m = 0.2,\, D_c^m=0.28$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/WM_cs_SW_n1_Dcm2p8e-01_Dgm2e-01_koff5e-01_phi1e-01_tend1000}
\caption{$D_g^m = 0.2,\, D_c^m=0.28$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/WM_cs_SW_n1_Dcm2e-01_Dgm2e-01_koff5e-01_phi1e-01_tend300}
\caption{$D_g^m = 0.2,\, D_c^m=0.2$}
\end{subfigure}
\caption{\label{fig:WM_wave} Numerically computed traveling waves with $k^{\rm off} = 0.5$ and membrane diffusion coefficients taken from Fig.~\ref{fig:stb_nonlocal}(b). In panels (a-b), there is coexistence of rotating and standing waves near the mode $n=1$ Hopf stability boundary. Further below this threshold, in $D_c^m = 0.2$, standing waves are only transient and evolve into rotating waves. Each solution is shown in Supplemental Movies S3, S4 and S5.}
\end{figure}

Finally, for both $k^{\rm off} = 0.5$ and $k^{\rm off} = 0.1$, we notice in Fig.~\ref{fig:stb_nonlocal} the presence of Bogdanov--Takens bifurcation points associated with the $n=1$ spatial mode and indicated by the branching of the black from the blue curves. Traveling waves are expected to interact with stationary patterns in the neighborhood of such a bifurcation point, and further details about its effect on the dynamics will be given in Section \ref{sec:2D_membrane_bulk} in the context of finite bulk diffusion.

\subsection{Two-dimensional membrane-bulk reaction-diffusion model}\label{sec:2D_membrane_bulk}

In this section, we consider the full two-dimensional membrane-bulk reaction-diffusion model with circular geometry. Our aim is to determine the effects of a finite bulk diffusion field on the spatio-temporal dynamics from Section \ref{sec:nonlocal}. Linear stability analysis is performed first, followed then by numerical simulations.

\subsubsection{Linear stability analysis}

We consider a perturbation of the spatially uniform steady state given by
\begin{equation}\label{eq:pert_disk}
C(r,\theta,t)=C^*+{\rm e}^{\lambda t}\Phi(r,\theta)\,, \quad G(r,\theta,t)=G^*+{\rm e}^{\lambda t}\Psi(r,\theta)\,, \quad c(\theta,t)=c^*+{\rm e}^{\lambda t}\phi(\theta)\,, \quad g(\theta,t)=g^*+{\rm e}^{\lambda t}\psi(\theta),
\end{equation}
where $\Phi(r,\theta),\,\Psi(r, \theta)$ and $\phi(\theta),\, \psi(\theta)$ are respectively the bulk and membrane eigenfunctions, while $(C^*,\, G^*,\, c^*,\, g^*)$ is such as introduced in Section \ref{sec:nonlocal}. Next, upon inserting \eqref{eq:pert_disk} within equations \eqref{eq:2Dmodel}--\eqref{eq:consv2dMAIN} and after linearization, we obtain the following eigenvalue problem in the bulk:
\begin{subequations}\label{eq:linearsys2d}
\begin{align}
& D_c\left[\Phi_{rr}(r,\theta) +\frac{1}{r}\Phi_r(r,\theta)+\frac{1}{r^2}\Phi_{\theta\theta}(r,\theta) \right] - \lambda\Phi(r,\theta) = 0,\quad 0<r<R,\quad 0<\theta<2\pi,\\
& D_g\left[\Psi_{rr}(r,\theta) +\frac{1}{r}\Psi_r(r,\theta)+\frac{1}{r^2}\Psi_{\theta\theta}(r,\theta) \right] - \lambda\Psi(r,\theta) = 0,\quad 0<r<R,\quad 0<\theta<2\pi,
\end{align}
and on the membrane:
\begin{align}\label{eq:linear_memb}
& \frac{D_c^m}{R^2}\phi_{\theta\theta}(\theta) + \frac{\partial\mathcal{F}}{\partial C}\Phi(R,\theta) + \left(\frac{\partial\mathcal{F}}{\partial c} - \lambda\right)\phi(\theta) + \frac{\partial\mathcal{F}}{\partial g}\psi(\theta) = 0\,, \quad 0<\theta<2\pi\,, \\
& \frac{D_g^m}{R^2}\psi_{\theta\theta}(\theta) + \frac{\partial\mathcal{G}}{\partial G}\Psi(R,\theta) + \frac{\partial\mathcal{G}}{\partial c}\phi(\theta) + \left(\frac{\partial\mathcal{G}}{\partial g} - \lambda\right)\psi(\theta) = 0\,, \quad 0<\theta<2\pi\,,
\end{align}
with boundary conditions given by
\begin{align}\label{eq:linear_BC2D}
&-D_c\Phi_r(R,\theta)=\frac{\partial\mathcal{F}}{\partial C}\Phi(R,\theta)+\frac{\partial\mathcal{F}}{\partial c}\phi(\theta)+\frac{\partial\mathcal{F}}{\partial g}\psi(\theta),\quad -D_g\Psi_r(R,\theta)=\frac{\partial\mathcal{G}}{\partial G}\Psi(R,\theta)+\frac{\partial\mathcal{G}}{\partial c}\phi(\theta)+\frac{\partial\mathcal{G}}{\partial g}\psi(\theta).
\end{align}
Here each of the partial derivatives is evaluated at the spatially uniform steady state. Also, the conservation laws \eqref{eq:consv2d} yield
\begin{equation}\label{eq:zero_mass}
\int_0^{2\pi}\int_0^R \Phi(r,\theta)rdrd\theta + \int_0^{2\pi} \phi(\theta)R d\theta = 0\,, \quad \int_0^{2\pi}\int_0^R \Psi(r,\theta)rdrd\theta + \int_0^{2\pi} \psi(\theta)R d\theta = 0.
\end{equation}
\end{subequations}

Our aim is to derive, for each spatial mode, the characteristic equations satisfied by the eigenvalue parameter. We consider first the general case where $\lambda$ is nonzero, followed then by the zero-crossing eigenvalue case.

\paragraph{Nonzero eigenvalue.} When $\lambda \neq 0$, separation of angular and radial variables yields the following expressions for the eigenfunctions:
\begin{equation}\label{eq:ansatz}
\Phi_n(r,\theta)=A_n \frac{I_n(\omega_c r)}{I_n(\omega_c R)} e^{in\theta}\,, \quad \Psi_n(r,\theta)=B_n \frac{I_n(\omega_g r)}{I_n(\omega_g R)} e^{in\theta}\,, \quad \phi_n(\theta) = u_n e^{in\theta}\,, \quad \psi_n(\theta) = v_n e^{in\theta}\,,
\end{equation}
where $\omega_c=\sqrt{\lambda/D_c}$, $\omega_g=\sqrt{\lambda/D_g}$ and $I_n(z)$ for $n \in \mathbb{Z}$ are the usual modified Bessel functions. Next, from the linearized boundary conditions \eqref{eq:linear_BC2D}, $A_n$ and $B_n$ are readily found to satisfy
\begin{equation*}
A_n = - \frac{\frac{\partial\mathcal{F}}{\partial c} u_n + \frac{\partial\mathcal{F}}{\partial g} v_n}{D_c \omega_c \frac{I_n'(\omega_c R)}{I_n(\omega_c R)} + \frac{\partial\mathcal{F}}{\partial C}}\,, \quad B_n = - \frac{\frac{\partial\mathcal{G}}{\partial c} u_n + \frac{\partial\mathcal{G}}{\partial g} v_n}{D_g \omega_g \frac{I_n'(\omega_g R)}{I_n(\omega_g R)} + \frac{\partial\mathcal{G}}{\partial G}}.
\end{equation*}
Then upon substituting \eqref{eq:ansatz} within the linearized membrane reaction-diffusion system, we derive the following linear system to be satisfied by the eigenvector $(u_n, v_n)^T$:
\begin{equation}\label{eq:nonzero}
\begin{pmatrix}
\lambda + \frac{n^2 D_c^m}{R^2} - p_n(\lambda) \frac{\partial\mathcal{F}}{\partial c} & - p_n(\lambda) \frac{\partial\mathcal{F}}{\partial g} \\
- q_n(\lambda) \frac{\partial\mathcal{G}}{\partial c} & \lambda + \frac{n^2 D_g^m}{R^2} - q_n(\lambda) \frac{\partial\mathcal{G}}{\partial g}
\end{pmatrix}
\begin{pmatrix} u_n \\ v_n \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\,,
\end{equation}
where $p_n(\lambda)$ and $q_n(\lambda)$ are defined by
\begin{equation*}
p_n(\lambda) = \frac{D_c\omega_cI_n'(\omega_cR)}{D_c\omega_cI_n'(\omega_cR) + \frac{\partial\mathcal{F}}{\partial C} I_n(\omega_cR)}\,, \quad q_n(\lambda) = \frac{D_g\omega_gI_n'(\omega_gR)}{D_g\omega_gI_n'(\omega_gR)+\frac{\partial\mathcal{G}}{\partial G} I_n(\omega_gR)}.
\end{equation*}
Finally, setting the determinant of the matrix in \eqref{eq:nonzero} to zero yields the following characteristic equation:
\begin{equation}\label{eq:lambdaFull}
F_n(\lambda) \equiv \Bigg[\lambda+\frac{n^2 D_c^m}{R^2}-p_n(\lambda)\frac{\partial\mathcal{F}}{\partial c}\Bigg] \Bigg[\lambda+\frac{n^2 D_g^m}{R^2}-q_n(\lambda)\frac{\partial\mathcal{G}}{\partial g}\Bigg] - p_n(\lambda)q_n(\lambda)\frac{\partial\mathcal{F}}{\partial g}\frac{\partial\mathcal{G}}{\partial c} = 0.
\end{equation}
For an oscillatory Turing instability, purely imaginary roots of $F_n(\lambda)$ for $n \neq 0$ are necessary. Fortunately when treating the spatially uniform mode, the integral constraints \eqref{eq:zero_mass} and the linear system \eqref{eq:nonzero} with $n=0$ lead to an equivalent solvability condition. A rigorous Turing stability analysis would require us to prove the absence of roots $\lambda$ with $\Re(\lambda) \geq 0$ for the characteristic equation $F_0(\lambda) = 0$, which cannot simply be verified numerically with the argument principle from complex analysis because of the singularity at the origin. This challenging task is therefore left for further study. Similar issues were encountered when treating the in-phase mode in the PDE--ODE model.

\paragraph{Zero-crossing eigenvalue.} For this special case, a suitable ansatz for the eigenfunctions is given by
\begin{equation}\label{eq:ansatz_zero}
\Phi_n(r,\theta) = A_n \left(\frac{r}{R}\right)^n e^{in\theta}\,, \quad \Psi_n(r,\theta) = B_n \left(\frac{r}{R}\right)^n e^{in\theta}\,, \quad \phi_n(\theta) = u_n e^{in\theta}\,, \quad \psi_n(\theta) = v_n e^{in\theta}\,,
\end{equation}
where from the linearized boundary conditions \eqref{eq:linear_BC2D}, $A_n$ and $B_n$ are found to satisfy
\begin{equation*}
A_n = -\frac{\frac{\partial\mathcal{F}}{\partial c} u_n + \frac{\partial\mathcal{F}}{\partial g} v_n}{D_c\frac{n}{R} + \frac{\partial\mathcal{F}}{\partial C}}\,, \quad B_n = -\frac{\frac{\partial\mathcal{G}}{\partial c} u_n + \frac{\partial\mathcal{G}}{\partial g} v_n}{D_g\frac{n}{R} + \frac{\partial\mathcal{G}}{\partial G}}\,.
\end{equation*}
The homogeneous linear system satisfied by the eigenvector $(u_n,v_n)^T$ is then given by
\begin{equation}\label{eq:zero}
\begin{pmatrix}
\frac{n^2 D_c^m}{R^2} - p_n(0)\frac{\partial\mathcal{F}}{\partial c} & - p_n(0)\frac{\partial\mathcal{F}}{\partial g} \\
- q_n(0) \frac{\partial\mathcal{G}}{\partial c} & \frac{n^2 D_g^m}{R^2} - q_n(0)\frac{\partial\mathcal{G}}{\partial g}
\end{pmatrix} \begin{pmatrix} u_n \\ v_n \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\,,
\end{equation}
where $p_n(0)$ and $q_n(0)$ are defined by
\begin{equation*}
p_n(0) = \frac{D_c n}{D_c n + R\frac{\partial\mathcal{F}}{\partial C}}\,, \quad q_n(0) = \frac{D_g n}{D_g n + R\frac{\partial\mathcal{G}}{\partial G}}\,.
\end{equation*}
After setting the determinant of the matrix in \eqref{eq:zero} to zero, we obtain the following necessary condition for a zero-crossing eigenvalue:
\begin{equation}\label{eq:charlambda0}
F_n(0) \equiv \left[\frac{n^2 D_c^m}{R^2} - p_n(0)\frac{\partial\mathcal{F}}{\partial c} \right] \left[\frac{n^2 D_c^m}{R^2} - q_n(0)\frac{\partial\mathcal{G}}{\partial g} \right] - p_n(0)q_n(0)\frac{\partial\mathcal{F}}{\partial g}\frac{\partial\mathcal{G}}{\partial c} = 0.
\end{equation}
When equation \eqref{eq:charlambda0} is satisfied for $n \neq 0$, the trivial steady state undergoes a pitchfork bifurcation, corresponding to the stationary Turing instability.

\subsubsection{Spatio-temporal pattern formation in the two-dimensional model}

We now explore the effect of diffusion and GEF membrane residence time on the formation of spatio-temporal patterns on a circular bulk domain when the inactive species diffuse at finite rates. As bifurcation parameters, we therefore select the bulk and membrane diffusion coefficients, as well as $k^{\rm off}$, the GEF dissociation rate. The spatial discretization in our numerical method combines a finite element approach in the bulk with a finite difference approach on the boundary (see Appendix~\ref{2dnumerics} for further details). The initial conditions stimulate a specific spatial mode $n$ by adding  the corresponding eigenfunction to the uniform steady state. Finally, the stability of the uniform steady state with respect to small perturbations is predicted from a numerical eigenvalue computation, directly by solving \eqref{eq:charlambda0} and \eqref{eq:lambdaFull}.

We present in panel (a) of Fig.~\ref{fig:wave} a stability diagram in the $(k^{\rm off},D)$ parameter plane, where we assume that both inactive Cdc42 and GEF diffuse at the same rate $D \equiv D_c = D_g$. We remark that the stability boundaries are similar to those arising in the one-dimensional PDE--ODE model as shown in panel (a) of Fig.~\ref{fig:stability_boundaries}. In both cases, the simultaneous increase of $D_c$ and $D_g$ has a destabilizing effect and results in a reduction of the region of linear stability of the trivial steady state. However, a key difference between the two diagrams is the curve of primary instabilities. In 2-D circular geometry, this curve is composed of a pitchfork (\textcolor{blue}{blue}) and a Hopf (\textcolor{black}{black}) bifurcation segments, both associated with the $n=1$ spatial mode, that connect in a Bogdanov--Takens bifurcation point in $(k^{\rm off}, D) \approx (0.06, 4.83)$. Hence, decreasing $k^{\rm off}$ when the bulk diffusion level is above this threshold leads to a loss of stability of the trivial steady state through a Hopf bifurcation and traveling (or rotating) waves are expected to form. Below the Bogdanov--Takens bifurcation point, our linear stability analysis predicts a loss of stability through a pitchfork bifurcation, leading to the formation of stationary patterns with accumulation of proteins in a single location.

As the Hopf curve approaches the pitchfork stability boundary, the pair of critical eigenvalues tends to zero and a slowing down of the phase velocity of the traveling wave is expected. As shown in panels (b-c) of Fig.~\ref{fig:wave}, this phenomena is observed in our numerical simulations. Eventually, the wave stops and becomes a stationary pattern (see Fig.~\ref{fig:pattern_n1}). Furthermore, in the vicinity of the Bogdanov--Takens bifurcation point, it is possible to derive a four-dimensional system of ODEs approximating the dynamics of the full bulk-membrane reaction-diffusion system. Although, such a dimensionality reduction process is beyond the scope of this study, we remark that the analysis of the resulting normal form could reveal the precise boundary delimiting the regions of existence of traveling waves and stationary patterns.

\begin{figure}[H]
\centering
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/stability_diagram_disk_koff_D_Dcm1e-01_Dgm1e-01}
\caption{$(k^{\rm off},D)$ with $D_c^m = D_g^m = 0.1$.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_TW_n1_Dc30_Dg30_Dcm1e-01_Dgm1e-01_koff6e-01_phi1e-01_tend300}
\caption{$k^{\rm off} = 0.6,\, D=30$.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_TW_n1_Dc6_Dg6_Dcm1e-01_Dgm1e-01_koff1e-01_phi1e-01_tend300}
\caption{$k^{\rm off} = 0.1,\, D=6$.}
\end{subfigure}
\caption{\label{fig:wave} Panel (a): stability boundaries in the $D$ versus $k^{\rm off}$ parameter plane. Panels (b-c): numerically computed traveling wave solutions of active Cdc42 along the circular membrane. The blue stars in panel (a) indicate where each simulation is performed. As the Bogdanov--Takens bifurcation point is approached, we remark a slowing down of the phase velocity of the traveling wave. Membrane diffusion coefficients are equal to $D_c^m = D_g^m = 0.1$, while other parameter values are given in Table~\ref{table:parm}.}
\end{figure}

For the stationary pattern in Fig.~\ref{fig:pattern_n1}, we notice a small GEF dissociation rate ($k^{\rm off} = 0.05$), which causes the proteins to mostly accumulate on the membrane. Also, in contrast with coupled bulk-membrane reaction-diffusion systems without mass conservation (cf.~\cite{levine2005membrane,gomez2019,paquin2018}), we note that high (low) concentrations on the membrane match with low (high) concentrations in the bulk.

\begin{figure}[H]
\centering
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_TW_n1_Dc4_Dg4_Dcm1e-01_Dgm1e-01_koff5e-02_phi1e-01_tend300}
\caption{Kymograph of active Cdc42 level.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/C_cs_TW_n1_Dc4_Dg4_Dcm1e-01_Dgm1e-01_koff5e-02_phi1e-01_tend300}
\caption{Active and inactive Cdc42 levels.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/G_gs_TW_n1_Dc4_Dg4_Dcm1e-01_Dgm1e-01_koff5e-02_phi1e-01_tend300}
\caption{Active and inactive GEF levels.}
\end{subfigure}
\caption{\label{fig:pattern_n1} Stable polarized state with $k^{\rm off} = 0.05$ and $D = 4$ (other parameters are the same as in the caption of Fig.~\ref{fig:wave}). In panels (b-c), the level of membrane-bound (active) proteins is given by the blue curve. The black dot in the same panels indicates the position on the disk with zero azimuthal coordinate. Because of circular symmetry, any rotations and reflections of this pattern is also a solution. See Supplemental Movie S6.}
\end{figure}

Our numerical simulations also revealed the possible subcritical branching of the Hopf bifurcation associated with the $n=1$ spatial mode. In Fig.~\ref{fig:bistability}, we show the result of two simulations performed near the Hopf stability boundary, where the trivial steady state is expected to be linearly stable, and observe that a traveling wave can form given a large enough perturbation of the trivial steady state. We highlight that such a bistable behavior was not observed in Section \ref{sec:nonlocal}, when considering the nonlocal PDE system governing the dynamics in the well-mixed regime. This would be consistent with our findings in the one-dimensional PDE--ODE model, that the Hopf bifurcation is subcritical for low level of bulk diffusion, while it is supercritical in the well-mixed regime. The precise determination, through a normal form computation, of the diffusion level at which the criticality transition occurs, is left for further study.

\begin{figure}[H]
\centering
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/stability_diagram_disk_koff_D_Dcm1e-01_Dgm1e-01_zoom}
\caption{$k^{\rm off} = 0.2$, $D = 6$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_TW_n1_Dc6_Dg6_Dcm1e-01_Dgm1e-01_koff2e-01_phi1e-01_tend1000}
\caption{Trivial steady state.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_TW_n1_Dc6_Dg6_Dcm1e-01_Dgm1e-01_koff2e-01_phi2e-01_tend1000}
\caption{Traveling waves.}
\end{subfigure}
\caption{\label{fig:bistability} Bistability between the trivial steady state and traveling wave solutions. The blue star just below the stability boundary in panel (a) indicates where the simulation is performed. In panel (b), we observe a small perturbation of the trivial steady state that slowly vanishes. However, as shown in panel (c), applying a bigger perturbation leads to a traveling wave. This suggests a loss of stability through a subcritical Hopf bifurcation with respect to the $n=1$ spatial mode. Other parameters are the same as in the caption of Fig.~\ref{fig:wave}.}
\end{figure}

We now address the existence of standing waves (or pole-to-pole oscillations) in the full coupled membrane-bulk reaction-diffusion model, and recall that such oscillatory behavior is essential during the process of cell division. Numerical simulations from Fig.~\ref{fig:p2p} suggest the coexistence of pole-to-pole oscillations and traveling waves near the Hopf stability boundary. The parameter values employed here are the same as for panel (c) of Fig.~\ref{fig:wave}, although we remark that similar simulation results were obtained for different sets of parameters near the stability threshold. As shown in Fig.~\ref{fig:p2p_to_wave}, clear transitions between standing and traveling waves are observed away from the Hopf stability boundary. Overall, those findings are consistent with the dynamics in the limit of fast bulk diffusion (see Section~\ref{sec:nonlocal}).

\begin{figure}[H]
\centering
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_TW_n1_Dc6_Dg6_Dcm1e-01_Dgm1e-01_koff1e-01_phi1e-01_tend1000}
\caption{Clockwise rotating waves.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_TW_n-1_Dc6_Dg6_Dcm1e-01_Dgm1e-01_koff1e-01_phi1e-01_tend1000}
\caption{Anticlockwise rotating waves.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_SW_n-1_Dc6_Dg6_Dcm1e-01_Dgm1e-01_koff1e-01_phi1e-01_tend1000}
\caption{Standing waves.}
\end{subfigure}
\caption{\label{fig:p2p} Coexistence of rotating and standing waves when $k^{\rm off} = 0.1$ and $D=6$ (other parameters are the same as in the caption of Fig.~\ref{fig:wave}). In the vicinity of a Hopf stability boundary, different initial conditions can either lead to clockwise rotating waves (panel (a)), anticlockwise rotating waves (panel (b)) or standing waves (pole-to-pole oscillations, see panel (c)). The traveling waves in each panel are shown in Supplemental Movies S7, S8 and S9.}
\end{figure}

\begin{figure}[H]
\centering
\begin{subfigure}{0.36\linewidth}
\includegraphics[width=\linewidth]{figures/stability_diagram_disk_koff_D_Dcm1e-01_Dgm1e-01_away}
\end{subfigure}
\begin{subfigure}{0.42\linewidth}
\includegraphics[width=\linewidth]{figures/cs_SW_n1_Dc30_Dg30_Dcm1e-01_Dgm1e-01_koff4e-01_phi1e-01_tend300}
\end{subfigure}
\caption{\label{fig:p2p_to_wave} Transition from standing to rotating waves in the highly nonlinear regime, far from any stability boundaries. The blue star in the left panel indicates where the simulation is performed (other parameters are the same as in the caption of Fig.~\ref{fig:wave}).}
\end{figure}

Next, we investigate the formation of patterns with multiple protein accumulation sites, which happens when higher spatial modes become unstable. For the purpose of avoiding any oscillatory instabilities, we set $k^{\rm off} = 1,\, D=10$, and select the membrane diffusion coefficients as bifurcation parameters. Shown in panel (a) of Fig.~\ref{fig:multipeak} are four stability curves in the $D_c^m$ versus $D_g^m$ parameter plane. Each of these curves is associated with a distinct mode and computed directly from the zero-crossing eigenvalue relation \eqref{eq:charlambda0}. The resulting diagram is qualitatively similar to its counterpart in the limit of fast bulk diffusion (see panel (a) of Fig.~\ref{fig:stb_nonlocal}). Here also, we note the stabilizing effect of simultaneously decreasing the bulk diffusion coefficients. Finally, in panels (b-c) of Fig.~\ref{fig:multipeak} we show Cdc42 stationary distributions in the bulk and on the membrane when multiple accumulation sites form.

\begin{figure}[H]
\centering
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/stability_diagram_disk_Dgm_Dcm_koff1_Dc10_Dg10}
\caption{$(D_g^m,D_c^m)$ with $k^{\rm off} = 1,\,D=10$.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/C_cs_TW_n2_Dc10_Dg10_Dcm1e-01_Dgm6e-01_koff1_phi1e-01_tend300}
\caption{Stationary pattern with two peaks.}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/C_cs_TW_n3_Dc10_Dg10_Dcm4e-02_Dgm2e-01_koff1_phi1e-01_tend300}
\caption{Stationary pattern with three peaks.}
\end{subfigure}
\caption{\label{fig:multipeak} Panel (a): pitchfork stability boundaries in the $D_c^m$ versus $D_g^m$ parameter plane with $k^{\rm off} = 1$, $D=10$ and other parameters taken from Table~\ref{table:parm}. Panel (b): the stationary pattern is computed near a mode $n=2$ stability boundary, as indicated by the square in panel (a) (see Supplemental Movie S10). Panel (c): the stationary pattern is computed near a mode $n=3$ stability boundary, as indicated by the triangle in panel (a) (see Supplemental Movie S11).}
\end{figure}

We conclude this section with simulations showing how the number of peaks of a stationary pattern can be reduced as the GEF membrane diffusion coefficient increases and as the mode $n=1$ pitchfork stability boundary is approached. More precisely, in panels (a-c) of Fig.~\ref{fig:pattern_n23}, we consider the horizontal slice $D_c^m=0.1$ and observe a transition from a two-peak to a single peak pattern. In panels (d-f) a transition between a three-peak and a single peak pattern is observed on the slice $D_c^m = 0.04$. We remark that the intermediate patterned states involve circular motions and peak switching that occur in the absence of any Hopf bifurcations. We believe this phenomena to result from nonlinear interactions between two or more spatial modes. This was also observed in the nonlocal version of this model, although here the relatively low level of bulk diffusion allows for smoother switching. Our model suggests that the number of polarized sites is dependent on membrane diffusion coefficients. As both $D_c^m$ and $D_g^m$ decrease, more polarity sites can be observed. Our findings for all three model variants are summarized in Table~\ref{table:compare}.

\begin{figure}[H]
\centering
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_TW_n2_Dc10_Dg10_Dcm1e-01_Dgm6e-01_koff1_phi1e-01_tend300}
\caption{$D_g^m = 0.6,\,D_c^m = 0.1$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_TW_n2_Dc10_Dg10_Dcm1e-01_Dgm8e-01_koff1_phi1e-01_tend300}
\caption{$D_g^m = 0.8,\,D_c^m = 0.1$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_TW_n2_Dc10_Dg10_Dcm1e-01_Dgm1_koff1_phi1e-01_tend300}
\caption{$D_g^m = 1,\,D_c^m = 0.1$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_TW_n3_Dc10_Dg10_Dcm4e-02_Dgm2e-01_koff1_phi1e-01_tend300}
\caption{$D_g^m = 0.2,\,D_c^m = 0.04$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_TW_n3_Dc10_Dg10_Dcm4e-02_Dgm4e-01_koff1_phi1e-01_tend300}
\caption{$D_g^m = 0.4,\,D_c^m = 0.04$}
\end{subfigure}
\begin{subfigure}{0.31\linewidth}
\includegraphics[width=\linewidth]{figures/cs_TW_n3_Dc10_Dg10_Dcm4e-02_Dgm6e-01_koff1_phi1e-01_tend300}
\caption{$D_g^m = 0.6,\,D_c^m = 0.04$}
\end{subfigure}
\caption{\label{fig:pattern_n23} Effect of reducing the membrane diffusivity ratio $D_c^m/D_g^m$ on pattern formation. In panels (a-c), we observe a transition from a pattern with two peaks to a localized pattern with a single peak. The same transition is observed in panel (d-f) for an initial weakly nonlinear pattern with three peaks. Other parameters are the same as in the caption of Fig.~\ref{fig:multipeak}.}
\end{figure}

\section{Conclusions and discussion}

\begin{table}[h!]
\caption{Comparison of three different model geometries}
\centering
\begin{tabular}{p{3.5cm}| p{4.3cm}| p{4.3cm}| p{4.3cm}}
\hline\hline
 & \footnotesize{\textbf{1-D PDE--ODE model}} & \footnotesize{\textbf{1-D nonlocal PDE model (circle)}} & \footnotesize{\textbf{2-D membrane-bulk model (disk)}} \\
\hline
\footnotesize{\textbf{Assumptions}} & \footnotesize{finite diffusion in the bulk, well-mixed in each end compartment} & \footnotesize{well-mixed in the bulk, finite diffusion on the membrane} & \footnotesize{finite diffusion in the bulk, finite diffusion on the membrane}\\
\hline
\footnotesize{\textbf{Patterns}} & \footnotesize{symmetric or asymmetric steady states, anti-phase (pole-to-pole) oscillations} & \footnotesize{stationary patterns with one or multiple peaks, traveling or standing waves} & \footnotesize{similar as for the nonlocal model} \\
\hline
\footnotesize{\textbf{Bifurcations}} & \footnotesize{ pitchfork,  subcritical or supercritical Hopf} & \footnotesize{stationary Turing (pitchfork),  oscillatory Turing (Hopf)}  & \footnotesize{similar as for the nonlocal model} \\
\hline
\footnotesize{\textbf{Slower bulk diffusion $(D\equiv D_c = D_g)$}}& \footnotesize{(1) stability region of symmetric steady state increases, (2) supercritical Hopf bifurcation becomes subcritical} & \footnotesize{NA} & \footnotesize{(1) stability region of homogeneous steady state increases, (2) traveling waves slow down near a Bogdanov--Takens bifurcation} \\
\hline
\footnotesize{\textbf{Slower Cdc42 membrane diffusion}} & \footnotesize{NA} & \footnotesize{(1) homogeneous steady state loses stability, (2) stationary patterns or traveling waves emerge} & \footnotesize{similar as for the  nonlocal model}  \\
\hline
\footnotesize{\textbf{Slower GEF membrane diffusion}} & \footnotesize{NA} & \footnotesize{(1) stationary patterns are less localized, (2) formation of traveling waves even without Hopf bifurcations} & \footnotesize{similar as for the nonlocal model, but stationary patterns with multiple peaks can form} \\
\hline
\footnotesize{\textbf{Low  GEF dissociation rate $k^{\rm off} \approx 0.1$}} & \footnotesize{anti-phase relaxation oscillations} & \footnotesize{(1) stationary patterns or traveling waves can form, (2) that interact near Bogdanov--Takens bifurcations} & \footnotesize{similar as for the nonlocal model}  \\
\hline
\footnotesize{\textbf{High  GEF dissociation rate $k^{\rm off} = 1$}} & \footnotesize{no oscillations} & \footnotesize{(1) no Hopf bifurcations, (2) stationary patterns with one or several peaks, that can evolve into traveling waves} & \footnotesize{similar as for the nonlocal model, but stationary patterns are less likely to evolve into traveling waves}  \\

\end{tabular}
\label{table:compare}
\end{table}

\subsection{Conclusions}

In conclusion, we have analyzed three spatial model variants of a cell polarization model in order to explore the effect of spatial dimensionality on pattern formation. We focused on a membrane-bulk model for the Cdc42-GEF system in fission yeast, where the same reaction kinetics were considered for three, consecutively more complicated, geometries: a one-dimensional PDE--ODE model, a nonlocal one-dimensional PDE model obtained when the diffusion coefficients of the cytosolic forms $D \to\infty$ (which is often a limit considered in models for polarization \cite{goryachev2008, lo2014}), and a full two-dimensional membrane-bulk model with finite diffusion in the cytosol and on the membrane. For these models, we analyze how bulk diffusion and lateral diffusion along the membrane affect the linear stability of a homogeneous steady state. Interestingly, fast membrane diffusion of Cdc42 stabilizes the homogeneous steady state, while fast bulk diffusion of Cdc42 has a destabilizing effect. Decreasing the GEF dissociation constant $k^{\rm off}$ shrinks the regions of linear stability.

For the 1-D PDE--ODE model, we observe both symmetric and asymmetric distributions of Cdc42 at the two ends, as well as pole-to-pole oscillations. Hopf bifurcations occur when GEF dissociation is slow compared to Cdc42 dissociation, with bulk diffusion affecting the amplitude and the shape of the oscillations near the stability thresholds. When bulk diffusion is slow, relaxation oscillations emerge from a subcritical Hopf bifurcation, while weakly nonlinear oscillations emerge from a supercritical Hopf bifurcations in faster bulk diffusion regime. For the 1-D nonlocal PDE model and the 2-D membrane-bulk model, stationary or oscillatory Turing instabilities, either corresponding to pitchfork or Hopf bifurcations, can occur when membrane-bound Cdc42 diffuses slowly. All the observed patterns have a symmetry-breaking effect, either consisting of stationary patterns or traveling waves. Stationary patterns with multiple peaks are observed when the diffusion rates of both membrane-bound species decrease. These stationary patterns are sensitive to perturbations and can undergo sudden rotations or evolve into traveling waves even in the absence of Hopf bifurcations. The main results are summarized in Table ~\ref{table:compare}.

Note that there are two different limits that can result in a one-dimensional model. If the reactive compartments at the two cell ends are assumed to be well-mixed, while finite diffusion is assumed in the bulk, we obtain a one-dimensional PDE--ODE model. If we assume that the bulk domain is well-mixed and consists of a disk, with lateral diffusion allowed, we obtain a one-dimensional nonlocal PDE model on a circle.  No temporal oscillations with a uniform spatial profile are observed in the nonlocal model. Instead, the dominant patterns are traveling waves and stationary patterns. In comparison to these 1-D model reductions, the 2-D model exhibits similar spatio-temporal patterns as the 1-D nonlocal PDE model. On a 2-D circular bulk domain, there are an infinite number of spatial modes associated to the linearized system, all corresponding to the circular harmonics, and this therefore increases the likelihood of multiple interacting unstable modes. The Bogdanov--Takens bifurcation with $O(2)$ symmetry is one notable example of this richer bifurcation structure that is absent from the 1-D PDE--ODE model. One last difference consists in the absence of diffusion-driven instability in the 1-D PDE--ODE model, partly due because of the well-mixed assumption of the two end compartments.

\subsection{Discussion}
The role that spatial segregation between components plays in giving rise to spatio-temporal patterns in cell biology is beginning to be appreciated. By taking into account spatial segregation between nucleus and cytoplasm, oscillations can emerge in the Hes1 and p53 pathways \cite{sturrock2011}. A rigorous analysis showing that diffusion can give rise to oscillations (i.e., treating the diffusion coefficient as a bifurcation parameter) was presented in \cite{chaplain2015hopf}. Pole-to-pole oscillations of Min proteins in \textit{E. coli} have typically been studied using standard reaction-diffusion equations (reviewed in \cite{kruse2007,loose2011}). However, bulk-surface models for Min oscillations have also been proposed \cite{huang2003,halatek2012} and are more realistic in terms of protein localization. An often overlooked problem in mathematical modeling of signaling networks is the effect of cellular geometry. Initial models of spatio-temporal phenomena are often studied in an idealized 1-D setting \cite{jilkine2011comparison}. If a certain model gains acceptance, then its extension to 2-D or 3-D domains may be done, often with some corrections of parameter values or reaction terms. However, a systematic analysis of the effect of dimensionality on spatio-temporal behavior of models of cell polarization as dimensionality increases from 1-D to 2-D is not usually performed. (One notable exception is the 1-D reaction-diffusion \cite{mori2011} and 2-D membrane-bulk versions \cite{cusseddu2018} of the wave-pinning model of cell polarization.)

Pattern formation resulting from coupling a diffusion process in a bounded 2-D domain to a nonlinear reaction-diffusion process on the domain boundary has been studied in \cite{levine2005membrane,ratz2012turing,ratz2015turing,paquin2018}. It was first shown in \cite{levine2005membrane} that the formation of membrane-bound stationary patterns is possible even when the two species diffuse in the bulk at the same rate (membrane diffusion was neglected). R{\"a}tz and R{\"o}gers \cite{ratz2012turing,ratz2014} analyzed a coupled membrane-bulk model and its nonlocal reduction for small GTPase. Their model includes a diffusion equation for inactive Cdc42 in the bulk and a pair of reaction-diffusion equations for active and inactive Cdc42 on the membrane. From a linear stability analysis, they showed that differences between bulk and lateral diffusion rates can trigger a Turing instability. Fast bulk diffusion and slow membrane diffusion of Cdc42 can also lead to pattern formation in our model, even in the case of equal membrane diffusivities and equal bulk diffusivities. The major difference lies in the modeling assumptions about the underlying role of GEF in the activation of Cdc42 proteins. In \cite{ratz2012turing,ratz2014}, GEF is assumed to be at quasi-steady state and is expressed in terms of the concentration of active Cdc42, whereas our model includes a separate variable for GEF. Furthermore, in \cite{ratz2012turing,ratz2014} it is assumed that GEF only activates Cdc42, while our model has both  positive and negative feedback loops that involve GEF. Thus, Hopf bifurcations and oscillatory patterns can occur in our model, but are not seen in the model of R{\"a}tz and R{\"o}gers \cite{ratz2012turing,ratz2014}. Our results can also be compared with patterning in regular reaction-diffusion systems, with neither membrane-bulk coupling nor mass conservation. For a disk domain with no-flux boundary conditions, stationary patterns are expected to form over the full domain, with their building blocks in the linear regime corresponding to the eigenfunctions of the Laplacian on the disk. In our model, because diffusion was the only process assumed in the bulk, we merely observed patterning in some layer near the boundary of the bulk. Another scenario is to consider a 1-D reaction-diffusion with periodic boundary conditions, with no nonlocal coupling. There, one can verify if stationary patterns with multiple peaks exist in a wider parameter regime, and if they are more stable than in their mass-conserved membrane-bulk counterpart.

Finally, we discuss some limitations and possible extensions of our work. From a biological perspective, our model predicts that spatio-temporal patterning is more likely to emerge when the GEF dissociation rate $k^{\rm off}$ is small compared to the Cdc42 dissociation rate. In the one-dimensional PDE--ODE model, Cdc42 and GEF oscillate out of phase, while the spatial model with diffusion of membrane-bound forms predicts that Cdc42 and GEF spatially segregate along the membrane. Our model also predicts that the number of polarized sites should increase as the membrane diffusion coefficients of both Cdc42 and GEF proteins decrease. While our model only involves one Cdc42 GEF, fission yeast cells have multiple GEFs. Experimental studies in fission yeast suggest that the GEF Scd1 oscillates from one pole to the other during cell growth, while no obvious oscillations of another GEF, Gef1, have been reported \cite{das2012,das2015}. Another study suggests that Gef1 is distributed in the cytosol \cite{tay2018}. Recent work suggests that positive feedback involves  Scd1, while negative feedback is  achieved through Cdc42 inhibiting Gef1 \cite{hercyk2019}. To take these findings into account, our model with a single GEF regulator would need to be replaced by a more complicated model with two GEFs, which would be harder to analyze. To limit active GEF localization to the cell tips in a full two-dimensional model, a space-dependent binding rate $k^{\rm on}$ should be assumed. It may also be necessary to include GAP dependent feedback to fully model the fission yeast system \cite{goryachev2017,tay2018}. While the 2-D model geometry may not be applicable to rod-shaped fission yeast, where Cdc42 seems to preferentially bind to the poles instead of all over the membrane, it is applicable to spherical budding yeast, where positive and negative feedbacks involving Cdc42 result in stable polarization and oscillatory behaviors. Note that traveling waves of Cdc42 have been observed in budding yeast under some experimental conditions \cite{ozbudak2005}.

For the convenience of mathematical analysis, we studied the 2-D model on a circular bulk domain, though fission yeast cells are rod-shaped. One follow-up question is to consider the effect of domain geometry and growth on the full 2-D membrane-bulk model. Another question to consider is how the patterning will change if model geometry is assumed to be three-dimensional. In a recent pattern formation study on spherical domains \cite{diegmiller2018}, Diegmiller et al.\ reduced a 3-D membrane-bulk model to a nonlocal reaction-diffusion equation on a 2-D surface by assuming that bulk diffusion is fast. Their model considers autocatalysis and predicts that a single polarity site would form. Following their framework, we can extend our 1-D nonlocal model to a 2-D surface and explore how the negative feedback may affect the number of polarity sites. Another potential extension is to take into account anomalous diffusion in the cytoplasm \cite{jeon2011} which raises challenges due to the Robin boundary conditions. Previous work \cite{henry2002} used fractional diffusion for modeling anomalous diffusion assuming an infinite domain. In a bounded domain, appropriate formulation of boundary conditions is an active area of research \cite{baeumer2018}.

From a dynamical systems perspective, a weakly nonlinear analysis could be performed to determine the criticality of the bifurcations in each of the three model variants, and to investigate how stationary patterns and traveling waves interact near Bogdanov--Takens bifurcation points. Our numerical results also revealed the formation of membrane-bound localized stationary patterns in the limit of small membrane diffusion coefficients ratio, similar to those studied in \cite{gomez2019}. There, another type of linear stability analysis can be performed, which consists of linearizing the membrane-bulk PDE system around a localized pattern constructed using matched asymptotic expansion, as opposed to linearizing around a uniform steady state. Finally, the analytical framework presented here can be used to study the spatio-temporal dynamics of signaling proteins in other cell biology contexts.

\section*{Acknowledgments}
FPL was supported by the NSERC Doctoral Award $\#6564$. BX is supported by the Robert and Sara Lumpkins Endowment for Postdoctoral Fellows in Applied and Computational Math and Statistics at the University of Notre Dame. This work was supported by the NSF (AEL by DMS-1815216, AJ by DMS-1615800 and KLD by DGE-1313583). We would also like to thank Michael Ward and Wayne Nagata for helpful discussions.