\section{Introduction}
\label{sec:introduction}

The concentration of free calcium ions plays a significant role in a large variety of cell types, controlling cellular functions such as secretion, fertilisation, muscle contraction and cell death, amongst many others.
Consequently, the dynamics of calcium has been studied by many groups via a combination of theoretical and experimental approaches over many years \cite{keener2010mathematical,dupont2016models}.

Models for calcium dynamics often contain multiple timescales, and there are several underlying physiological reasons for this. Some calcium release channels, such as inositol 1,4,5-trisphosphate (${\rm IP_3}${}) receptors or ryanodine receptors, have two inherent timescales depending on calcium concentration: rapid activation is followed by slow inhibition, providing a fundamental mechanism for the appearance of complex calcium dynamics in many cell types.
Multiple timescales in calcium dynamics can also arise from the spatial organisation of the calcium channels \cite{coombes2004receptors}, or simply because some calcium transport pathways are slower than others. For example, the exchange of calcium with the extracellular medium is usually considerably slower than calcium transport into and out of the endoplasmic reticulum \cite{sneyd2017dynamical}.

In previous work, Harvey et al.~\cite{harvey2011multiple} used geometric singular perturbation theory (GSPT) in its standard form to analyse the dynamics of several different calcium models. Their analysis relied on the assumption that there exists a globally valid way to identify some variables as uniformly `fast' or `slow' throughout phase space, which is a standard assumption in GSPT \cite{hek2010geometric}. However, in calcium dynamics the speed of a reaction often depends on the concentration of calcium. Thus, in a corresponding model the speed of that reaction may vary through phase space, meaning that one or more variables could be slow in some regions of phase space and fast in others. Such a feature is not easy to capture within the standard GSPT framework.
In recent years, numerous authors have shown the applicability of GSPT to systems with \wm{\em non-standard} time-scale separations of this kind; see, for example, \cite{Lizarraga2020b,Wechselberger2015,jelbart2019two,Kristiansen2019d,Huber2005,Kruff2019,goeke2014constructive,Lizarraga2020c} and in particular the work by Kosiuk and Szmolyan in \cite{Gucwa2009,kosiuk2011scaling,kosiuk2016geometric}.

In this article, we present a systematic analysis of a calcium dynamics model by adopting \wm{a coordinate-independent} framework recently developed in \cite{wechselberger2018geometric}, which applies to the study of general \wm{multiple time-scale} systems
\[
	z' = H(z,\epsilon) , \qquad z \in \mathbb R^n , \qquad 0 < \epsilon \ll 1 ,
\]
where $n \geq 2$
and the set $S = \{z \in \mathbb R^n:H(z,0)=0\}$ defines a smooth $(k<n)-$dimensional submanifold of $\mathbb R^n$ known as the \textit{critical manifold}. Importantly for our purposes, there is no assumed global separation within the vector of variables, $z$, between slow and fast variables.
The model we consider is two-dimensional, and derives from a higher-dimensional model recently proposed in \cite{sneyd2017dynamical}, that contains features and dynamics common to a number of intracellular calcium dynamics models.
The reduced model is shown to exhibit relaxation oscillations with a distinctive three-timescale structure, with a geometry similar to the autocatolytic oscillations studied in \cite{Gucwa2009}.

The main analytical obstacles stem from a \wm{non-standard} timescale separation and the large number of small parameters in the system. These complications are overcome by a \textit{process-based} identification of small parameters.
We identify a total of five small parameters, which we then link to a single small parameter based on suitable order of magnitude estimates. We prove the existence
of the observed relaxation oscillations by the application of {\em successive} blow-up transformations \wm{in the coordinate-independent GSPT framework}.

\wm{We then consider possible mechanisms for the onset of the observed relaxation oscillations under important model parameter variations. One possibility involving the rapid onset of relaxation oscillations in a so-called \textit{complete or incomplete canard explosion} is briefly discussed. A second possibility involving variation of a system parameter over several orders of magnitude, with \textit{no canard explosion}, is also considered.

We believe that the process-based philosophy and methodology presented here is well suited to the analysis of a large class of important biological systems, and that the methods and formalism is readily transferable to the analysis of these problems. Our analysis also serves as a first step towards understanding the higher-dimensional model from which the reduced model under consideration here was derived.

\

The manuscript is structured as follows. In Section \ref{sec:motivations_for_the_closed_cell_model} we present the model from \cite{sneyd2017dynamical} and its reduction to the form we will analyze. In Section \ref{sub:non_dimensionalisation} we non-dimensionalise the model, \wm{which then allows us in Section \ref{sub:construction_of_the_singular_limit} to} identify candidate small parameters for a singular perturbation analysis and derive a suitable singular perturbation problem for analysis. Section \ref{sec:slow-fast_analysis_and_statement_of_the_main_result} is devoted to a GSPT analysis of this system, including the statement of our main results. The detailed blow-up analysis required to prove existence
of the relaxation oscillations is deferred to the Appendix. \SJ{In Section \ref{sec:effect_of_tau_} we consider possibilities for the onset of the relaxation oscillations under parameter variation.
In Section \ref{sec:discussion_conclusion} we summarise and conclude.}

\section{Physiological background and the model}
\label{sec:motivations_for_the_closed_cell_model}

\begin{figure}
	\centering
	\includegraphics[width=0.52\linewidth]{Images/Calcium_fluxes_diagram_new}
	\caption{Typical calcium fluxes in a simplified model of intracellular calcium dynamics. A solid line denotes a physical flux, while a dotted line denotes an interaction of some kind. The binding of external agonists, such as hormones or neurotransmitters, to receptors on the cell membrane results in a cascade of biochemical reactions that leads, ultimately, to the production of ${\rm IP_3}${}. This ${\rm IP_3}${} diffuses through the cytoplasm and binds to ${\rm IP_3}${} receptors (IPR) that, in turn, release calcium from the endoplasmic reticulum (ER) to the cytosol; the resulting calcium flux is denoted by $J_{\rm \text{IPR}}$. The calcium flux through the IPR is controlled by the cytoplasmic calcium concentration, in a complex system of sequential positive and negative feedbacks, here all denoted by a single dotted line. Calcium is transported from the cytosol to the ER by calcium ATPase pumps, with the flux denoted by $J_{\text{SERCA}}$. The overall calcium concentration inside the cell is controlled by exchange of calcium with the exterior of the cell, with the resulting flux denoted by $J_{\rm pm}$.  The line denoting $J_{\rm pm}$ is shown in grey, as in this paper this flux is mostly assumed to be zero.}
	\label{fig:Calcium_fluxes_diagram}
\end{figure}

Figure \ref{fig:Calcium_fluxes_diagram} shows, in a highly simplified form, some of the major calcium fluxes involved in many models of intracellular calcium dynamics. To trigger calcium oscillations, agonists bind to cell membrane receptors, thus initiating a chain of reactions that leads (among other things) to the formation of ${\rm IP_3}${}. ${\rm IP_3}${} then diffuses in the cell and binds to ${\rm IP_3}${} receptors at the surface of the endoplasmic reticulum (ER), which release calcium from the ER into the cytosol. The concentration of free calcium in the cytosol modulates the release of further calcium through the ${\rm IP_3}${} receptors, by a process of fast positive feedback followed by a slower negative feedback. Calcium is transported from the cytosol back to the ER by calcium ATPase pumps known as SERCA pumps. The fast positive feedback and slow negative feedback of calcium on the ${\rm IP_3}${} receptor can result in oscillations of calcium concentration without the need for variation in ${\rm IP_3}${} concentration. Models in which calcium oscillations can occur for a constant ${\rm IP_3}${} concentration are called Class I models \cite{de1992single,dupont1997simulations,sneyd2017dynamical}.

There are other possible mechanisms that can cause calcium oscillations. For example, the cytoplasmic calcium concentration also modulates the production and degradation of ${\rm IP_3}${}. In this case, variations in ${\rm IP_3}${} concentration change the rate of calcium release which in turn influences ${\rm IP_3}${} concentration. Such feedback, which can be either positive or negative or both, can lead to oscillations. Models in which calcium oscillations depend on the effects of calcium concentration on the production and/or degradation of ${\rm IP_3}${} are called Class II models \cite{politi2006models,gaspers2014hormone,sneyd2017modeling}. In some cells both Class I and Class II mechanisms are present, leading to a hybrid oscillatory mechanism \cite{harvey2010understanding}.

It was shown in \cite{sneyd2017dynamical} that a Class I model could reproduce many qualitative features of the oscillations seen in a variety of cell types.
Accordingly, and for simplicity, we will ignore Class II mechanisms and focus on a Class I model. In this case, the ${\rm IP_3}${} concentration can be viewed as a substitute for agonist concentration. In many analyses of Class I models, ${\rm IP_3}${} concentration (as a proxy for agonist concentration) is treated as a bifurcation parameter. However, our focus here is not on bifurcations of the system but on the existence of oscillations
for certain intervals of ${\rm IP_3}${} concentration and of other system parameters, and we thus fix ${\rm IP_3}${} concentration.

Let $c$ denote the calcium concentration in the cytosol, $c_e$ the calcium concentration in the ER, and $h$ a gating variable which represents the capacity of the ${\rm IP_3}${} receptors to allow the passage of calcium. We consider the following model developed in \cite{sneyd2017dynamical}:
\begin{equation}
	\begin{split}
		\displaystyle\tau_h(c)\frac{\text{d}h}{\text{d}t}&=\displaystyleh_{\infty} (c) - h,\\
		\displaystyle\frac{\text{d}c}{\text{d}t}&=\displaystyleJ_{\text{IPR}}-J_{\text{SERCA}} + J_{\rm pm},\\
		\displaystyle\frac{\text{d}c_e}{\text{d}t}&=\displaystyle\gamma(J_{\text{SERCA}}-J_{\text{IPR}}).

		\label{eq:Full_model_open}
	\end{split}
\end{equation}
The functional forms for the fluxes $J_{\text{IPR}}$, $J_{\text{SERCA}}$ and $J_{\rm pm}$ will be specified later, as will the forms of $h_{\infty}$ and $\tau_h$.  The parameter $\gamma$ is the ratio between the volume of the cytosol and the volume of the ER.  Default parameter values used in the model are specified in Table \ref{tab:params}.

It is useful to define $c_t := c + c_e / \gamma$, which represents the total moles of free calcium in the ER and the cytosol, divided by the cytoplasmic volume. This variable can be used as an alternative to $c_e$ in the model, yielding:
\begin{equation}
	\begin{split}
		\displaystyle\tau_h(c)\frac{\text{d}h}{\text{d}t}&=\displaystyle h_{\infty} (c) - h,\\
		\displaystyle\frac{\text{d}c}{\text{d}t}&=\displaystyleJ_{\text{IPR}}-J_{\text{SERCA}} + J_{\rm pm},\\
		\displaystyle\frac{\text{d}c_t}{\text{d}t}&= J_{\rm pm}.
		\label{eq:Full_model_open_ct}
	\end{split}
\end{equation}

System \eqref{eq:Full_model_open_ct} includes the flux denoted $J_{\text{pm}}$ in Figure~\ref{fig:Calcium_fluxes_diagram}, which represents the exchange of calcium between the cytosol and the extracellular medium.
This calcium exchange is typically slow compared to the rate of other calcium fluxes in the cell and is not necessary on either physiological \cite{bruce2002phosphorylation} or mathematical \cite{pages2019model} grounds in order for there to be oscillations. For example, there are cell types which exhibit calcium oscillations in the absence of any calcium flux across the plasma membrane, while there are other cell types for which calcium oscillations are entirely dependent on calcium entry. Some cell types can exhibit both kinds of behaviour, depending on the exact agonist used, and the concentration of the agonist. Thus, as a general rule, it is important to understand the behaviour of models both with and without calcium transport across the plasma membrane. However, in cases where $J_{\text{pm}}$ is very small relative to the other fluxes, the study of the closed-cell model (where $J_{\text{pm}}=0$) gives enormous insight into the behaviour of the full model, and thus a closed-cell analysis is an important first step for the analysis of any calcium model.

Thus, we set $J_{\text{pm}} = 0$, in which case $c_t$ is constant and the model reduces to a planar system:
\begin{align}
	\begin{split}
		\displaystyle\tau_h(c)\frac{\text{d}h}{\text{d}t}&=\displaystyleh_{\infty} (c) - h,\\
		\displaystyle\frac{\text{d}c}{\text{d}t}&=\displaystyleJ_{\text{IPR}}-J_{\text{SERCA}}.
		\label{eq:Full_model}
	\end{split}
\end{align}

We now specify the remaining flux terms. $J_{\text{IPR}}$ represents the flux of calcium from the ER to the cytosol through the IPR, and is modelled by
\begin{equation}\label{def:PO}
J_{\text{IPR}} = k_{\text{IPR}} P_O (c_e -c) = \gamma k_{\text{IPR}} P_O\left(c_t - (1 + 1/\gamma) c \right)\,,
\end{equation}
where $k_{\text{IPR}}$ is the maximum flux through the IPR and $P_O$ is the open probability of the IPR. There are numerous models for $P_O$ in the literature \cite{dupont2016models,sneyd2005models,politi2006models,atri1993single},
but one of the simplest and most recent is that of \cite{sneyd2017dynamical}, i.e.,
\begin{align}
\label{eq:P_O}
	P_O = \frac{\beta}{\beta + k_\beta \left( \beta + \alpha\right)},
\end{align}
where
\begin{align*}
	\alpha = \phi_{pdown}(1 - m_\alpha h_{\infty}), \qquad
	\beta  = \phi_p m_\alpha h,
\end{align*}
with
\begin{align*}
		m_\alpha     = \frac{c^4}{K_c^4 + c^4},  \qquad
		h_{\infty}        = \frac{K_h^4}{K_h^4+ c^4}, \qquad
		\phi_{pdown} = \frac{K_p^2}{K_p^2+ p^2}, \qquad
		\phi_p       = \frac{p^2}{K_p^2+p^2}.
\end{align*}

$J_{\text{SERCA}}$ is the flux from the cytosol to the ER via the ATPase pumps. Our SERCA pump model is bidirectional, and given by
\begin{equation*}
	J_{\text{SERCA}} = V_s \frac{c^2-K c_e^2}{K_s^2+c^2}.
\end{equation*}

\begin{table}[t!]
	\center
	\begin{tabular}{|c|l|l||c|l|l|}
		\hline
		Parameter    & Value & Units    & Parameter & Value      & Units           \\
		\hline
		$k_\beta$    & 0.4   & -        & $p$       & 0.05       & $\mu$M          \\
		$K_c$        & 0.2   & $\mu$M   & $K_\tau$  & 0.1        & $\mu$M          \\
		$K_h$        & 0.08  & $\mu$M   & $V_s$     & 0.9        & $\mu$Ms$^{-1}$ \\
		$K_p$        & 0.2   & $\mu$M   & $K$   & 0.00001957 & -               \\
		$k_{\text{IPR}}$      & 10    & $s^{-1}$ & $K_s$     & 0.2        & $\mu$M          \\
		$\tau_{max}$ & 1000   & s        & $\gamma$  & 5.5        & -               \\
		$c_t$        & 2     & $\mu$M   & -         & -          & -               \\
		\hline
	\end{tabular}
	\caption{Default (dimensional) parameter values for the model, equations (\ref{eq:Full_model}).}
	\label{tab:params}
\end{table}

A significant feature of our model is that, depending on the calcium concentration $c$, the speed of evolution of $h$ varies by several orders of magnitude. This feature is captured by the functional form of $\tau_h(c)$, which we choose to be a Hill function
\begin{equation}\label{tau-h-orig}
	\tau_h (c) = \tau_{max} \frac{K_\tau^4}{K_\tau^4+ c^4}\,,
\end{equation}
which is a monotonically decreasing switch-like function; see Figure~\ref{fig:tau_h_fig}. The value of $\tau_h(c)$ varies between $0$s and $1000$s depending on $c$, leading to significant variation in the speed of evolution of $h$.

\begin{figure}
	\centering
		\includegraphics[scale=0.6]{Images/tau_h_new}
	\caption{
	The function $\tau_h(c)$ defined in \eqref{tau-h-orig} varies over several orders of magnitude with $c$, which causes the speed of evolution of $h$ to vary significantly through phase space.}
	\label{fig:tau_h_fig}
\end{figure}

\section{Non-dimensionalisation}
\label{sub:non_dimensionalisation}

In this section we non-dimensionalise the model \eqref{eq:Full_model}.

Since $h\in [0,1]$ is already dimensionless,
we only need to non-dimensionalise the variables $c$ and $t$. We define
\begin{equation}
	c = Q_c \tilde{c}, \qquad
	t = T \tilde{t} ,
\end{equation}
where $Q_c$ and $T$ denote a reference concentration and timescale respectively, to be specified below.

We obtain the following dimensionless model:
\begin{equation}
	\label{eq:Non_dim}
	\begin{split}
		\tilde{\tau}_h(\tilde{c}) \displaystyle\frac{dh}{d\tilde{t}} &= \displaystyle \tilde{h}_{\infty}(\tilde c) - h ,\\
		\displaystyle \frac{d\tilde{c}}{d\tilde{t}} &= \displaystyle \tilde{J}_{\text{IPR}}(h, \tilde c) - \tilde{J}_{\text{SERCA}}(\tilde c),\\
	\end{split}
\end{equation}
where the new (dimensionless) expressions $\tilde \tau_h(\tilde c), \tilde h_\infty(\tilde c), \tilde{J}_{\text{IPR}}(h,\tilde c)$ and $\tilde{J}_{\text{SERCA}}(\tilde c)$ are defined analogously to their dimensional counterparts, except in terms of the new (dimensionless) parameters
\begin{equation}
\label{eq:nondim_pars}
\tilde{K}_i = \frac{K_i}{Q_c}, \ \ \tilde \tau_{max} = \frac{\tau_{max}}{T} , \ \
\tilde{k}_\text{\text{IPR}} = T k_{\text{IPR}}, \ \ \tilde{V}_s = \frac{T}{Q_c}V_s, \ \
\tilde{c}_t = \frac{c_t}{Q_c}, \ \ \tilde p = \frac{p}{Q_c} ,
\end{equation}
where $i = c, h, \tau, s$ or $p$.
We choose the following reference scales for our model:
\begin{equation}
Q_c = c_t = 2\,,
\qquad T = \frac{1}{\gamma k_{\text{IPR}}} = \frac{1}{55}\,.
\end{equation}
Thus, \eqref{eq:Non_dim} shows the dynamics of $c$ relative to the total concentration $c_t$, on the IPR timescale set by the uniform scaling $\gamma k_{\text{IPR}}$; this is the fastest timescale in the system.
Numerical values for the dimensionless model parameters consistent with Table \ref{tab:params} are shown in
Table \ref{tab:params_non_dim}. In what follows we refer only to the dimensionless model \eqref{eq:Non_dim} and parameters \eqref{eq:nondim_pars}, \wm{dropping the tildes for notational simplicity}.

Figure \ref{fig:non_dim} shows
relaxation oscillations in system \eqref{eq:Non_dim} with the parameter values shown in Table \ref{tab:params_non_dim}. Both time series contain components evolving on different time-scales, which is indicative of a \wm{non-standard} timescale separation.

\begin{table}[t!]
	\center
	\begin{tabular}{|c|l||c|l|}
		\hline
		Parameter         & Value & Parameter            & Value      \\
		\hline
		$\tilde{k}_\beta$ & 0.4   & $\tilde{K}_\tau$     & 0.05       \\
		$\tilde{K}_c$     & 0.1 & $\tilde{V}_s$        & 0.0081     \\
		$\tilde{K}_s$     & 0.1 & $\tilde{K}$             & 0.000019 \\
		$\tilde{K}_p$     & 0.1   & $\tilde{K}_h$        & 0.04      \\
		$\tilde{k}_{\text{IPR}}$          & 0.18     & $\gamma$             & 5.5        \\
		$\tilde{p}$       & 0.025  & $\tilde{\tau}_{max}$ & 55000       \\
		$\tilde{c}_t$     & 1  & -                    & -          \\
		\hline
	\end{tabular}
	\caption{Dimensionless parameters in system \eqref{eq:Non_dim}. \wm{The tilde notation is dropped in Section~\ref{sub:construction_of_the_singular_limit}.}}
	\label{tab:params_non_dim}
\end{table}

\begin{figure}[t!]
	\centering
	\begin{subfigure}[b]{0.48\linewidth}
		\includegraphics[width=\linewidth]{Images/timeseries_nondim.png}
	\end{subfigure}
	\begin{subfigure}[b]{0.48\linewidth}
		\includegraphics[width=\linewidth]{Images/timeseries_nondim_zoom.png}
	\end{subfigure}
	\caption{Relaxation oscillations in system \eqref{eq:Non_dim} with parameters from Table \ref{tab:params_non_dim}. Left: time series for $\tilde c$ and $\tilde h$. Right: enlargement of one spike. Both $\tilde c$ and $\tilde h$ exhibit \wm{multiple time-scale} dynamics, indicative of a \wm{non-standard} timescale separation.}
	\label{fig:non_dim}
\end{figure}

\section{Small parameters and timescales}
\label{sub:construction_of_the_singular_limit}

In this section we identify candidate small parameters for use in a singular perturbation analysis of the model \eqref{eq:Non_dim}. As usual, \wm{the (dimensionless) pre-factors of any process term (listed in Table~\ref{tab:params_non_dim})} play an important role in determining the magnitude of a given contribution to the dynamics, and so \wm{will} be considered as candidates for small parameters.
\wm{Furthermore, since the relaxation oscillations in Figure \ref{fig:non_dim} suggest a time-scale separation that is \wm{non-standard}, we are motivated to identify processes including a small (dimensionless) parameter
$\nu$ that are smooth for all $0 < \nu \ll 1$ but not {\em uniformly continuous} in the limit $\nu \to 0$. Such non-uniform perturbations can lead to distinct singular limits depending on the region of phase space, and are frequently observed to result in \wm{non-standard} time-scale separations and dynamics. Typical examples in the context of biological models are Michaelis-Menten or Hill-type functions with `small' half-activation/deactivation values; see e.g.~\cite{Kosiuk2016,wechselberger2018geometric}.}

\subsection{$h$ dynamics}
\label{ssub:timescales_as_a_function_of_calcium_concentration}

The speed of evolution of $h$ varies with calcium concentration due to the presence of the Hill function $\tau_h(c)$ in system \eqref{eq:Non_dim}. The switch-like profile of $\tau_h(c)$ is a reflection of the physiological observation that the rate at which $h$ changes depends dramatically on the calcium concentration at the mouth of the channel. The parameter or \textit{half-value} $K_\tau$ defines the $c$ value at which $\tau_h(c)$ takes on its midpoint value; see Figure \ref{fig:tau_h_fig}. For each fixed value of $K_{\tau}$, we identify two \wm{non-overlapping} regimes:

\smallskip

\begin{enumerate}
	\item[(R1)] $c \in \mathcal I_1=[a,b]$ for some fixed $a$ and $b$ with $K_{\tau}<a<b$;

	\smallskip

	\item[(R2)] $c \in \mathcal I_2=[0,K_\tau \nu]$ for some fixed $\nu>0$.
\end{enumerate}

\medskip

\noindent For sufficiently small $K_\tau$, $c = O(1)$ in (R1) and $c = O(K_\tau)$ in (R2). Correspondingly, the speed of evolution of $h$ is different in regimes (R1) and (R2) if $K_\tau$ is sufficiently small.

Note from Table~\ref{tab:params_non_dim} that the default value for $K_\tau$ is numerically small. With an aim of being able to properly distinguish between regimes (R1) and (R2), we define our first small parameter,
\begin{equation}
\label{eq:h_scale_tau}
\epsilon_\tau:= K_\tau .
\end{equation}
Thus, the dynamics of $h$ evolves on different speeds within regimes (R1) and (R2) in the limit $\epsilon_\tau \to 0$.

A similar switch-like mechanism occurs due to the term $h_\infty(c)$, which also takes the form of a Hill function. The magnitude of this contribution also varies over different orders of magnitude as $c$ varies, as shown in Figure \ref{fig:processes_fig}. As with our discussion above, distinct regimes can be identified by considering the magnitude of $c$ relative to the half-value parameter $K_h$, which is again numerically small. In the following we set
\begin{align}
\label{eq:Kh_scaling}
\epsilon_h:= K_h.
\end{align}

Finally we see that the evolution of $h$ scales with the constant (i.e., not $c$-dependent) factor $\tau_{max} K_\tau^4 \approx 0.34$. Since $\tau_{max}$ is numerically large, we define a third small parameter
\begin{equation}
\label{eq:h_scale_tau_max}
\epsilon_{\tau_{max}} := \tau_{max}^{-1} .
\end{equation}

We have identified three small parameters $\epsilon_\tau, \epsilon_h$ and $\epsilon_{\tau_{max}}$, each relevant to the timescale on which the variable $h$ is expected to evolve. We return to the question of how these small parameters scale relative to one another after considering the $c$ dynamics.

\begin{figure}
	\centering
	\includegraphics[scale=0.1]{Images/switches_new}
	\caption{Normalised switch-like functions in system \eqref{eq:Non_dim}, which play a role in determining the dominant processes in \eqref{eq:Non_dim} for a given value of $ c$. Corresponding half-values $K_i$ for $i=h,\tau,c,s$ are also shown. The $J^-_{SERCA}$ profile is less than $0.5$ at the half-value $K_s$ due to extra $c_t-$dependence in the numerator in \eqref{eq:Jserca}.}
	\label{fig:processes_fig}
\end{figure}

\subsection{$c$ dynamics}
\label{ssub:scaling_by_process}

In order to identify small parameters relevant for understanding the evolution of $c$, we consider the independent IPR and SERCA processes contributing to the right-hand-side of the second equation in system \eqref{eq:Non_dim}. As before, we consider candidates for small parameters which constitute either (i) uniform scalings to IPR or SERCA contributions, or (ii) non-uniform perturbations which give rise to switching mechanisms.

First, consider the IPR term $J_{\text{IPR}}(h,c)$. Nonlinearity in this term arises from the open probability factor $P_O(h,c)$ in \eqref{eq:P_O}, which is primarily comprised of switch-like Hill functions.
In particular, one must consider the role of the functions $h_\infty(c)$ and $m_\alpha(c)$. The former we have already considered above, while the latter is a Hill function with a half-value $K_c$. Looking at Figure \ref{fig:processes_fig}, however, we see that $m_{\alpha}(c)$ cannot be reasonably approximated asymptotically by a constant function in either regime (R1) or (R2). Rather, the nonlinearity associated with this term is expected to play a significant role over the entire $c-$domain corresponding to the oscillations in Figure \ref{fig:non_dim}. For this reason, we do not consider $K_c$ further as a candidate small parameter; see also Remark \ref{rem:scaling} below.

Now consider the SERCA contribution $J_{\text{SERCA}}({c})$.
This term models two separate processes: one positive, resulting from calcium ions being pumped from the cytosol to the ER, and one negative, corresponding to a leak of calcium ions from the ER into the cytosol. These two processes can be distinguished by writing
\begin{equation}
J_{\text{SERCA}}(c) = J_{\text{SERCA}}^+(c) - J_{\text{SERCA}}^-(c),
\end{equation}
where
\begin{equation}
\label{eq:Jserca}
J_{\text{SERCA}}^+(c) = V_s \left( \frac{{c}^2}{{K}_s^2+c^2} \right),  \qquad
J_{\text{SERCA}}^-(c) = K \gamma^2 {V}_s \left( \frac{({c}_t - {c})^2}{{K}_s^2+{c}^2} \right).
\end{equation}
Both terms have a switch-like contribution with a common half-value $K_s$. Here again, the nonlinear switching mechanism is dynamically significant over both (R1) and (R2), i.e.,~$J_{\text{SERCA}}^{\pm}(c)$ are not well approximated by constant functions in either regime, and so we do not consider
$K_s$ as an additional small parameter.
However, the magnitudes of $J_{\text{SERCA}}^+({c})$ and $J_{\text{SERCA}}^-({c})$ are also determined by the constant, uniform factors, $V_s$ and $V_s {K} \gamma^2$, respectively. Since both of these values are numerically small (see Table \ref{tab:params_non_dim}), we identify two more candidates for small parameters:
\begin{equation}
\label{eq:c_scale}
\epsilon_{s} := V_s , \qquad  \epsilon_K :=  K \gamma^2 .
\end{equation}
Both of these are relevant for determining the expected timescale of the evolution of $c$.

\begin{remark}
	\label{rem:scaling}
	The decision to consider the half-values $K_{\tau,h}$ as `small' and the half-values $K_{s,c}$ as `$O(1)$' can also be motivated on modelling grounds.

	Our analysis in Section \ref{sec:slow-fast_analysis_and_statement_of_the_main_result} will show that the problem is tractable without the need to consider the limits $K_{s,c} \to 0$, which would correspond to an approximation of $m_\alpha(c)$ and $J^{\pm}_{\text{SERCA}}(c)$ by constant functions in regime (R1). By retaining $O(1)$ values for $K_{s,c}$ we preserve dynamics associated with the nonlinearity of $m_\alpha(c)$ and $J^{\pm}_{\textrm{SERCA}}(c)$ over the entire domain.
\end{remark}

\subsection{A common scaling for singular perturbation analysis}
\label{ssub:a_common_scaling_and_natural_timescales}

In Sections \ref{ssub:timescales_as_a_function_of_calcium_concentration}-\ref{ssub:scaling_by_process},
we introduced five independent candidates for small parameters:
\[
\epsilon_\tau,\ \epsilon_{\tau_{max}}, \  \epsilon_h,\  \epsilon_s,\  \epsilon_K .
\]
For the sake of tractability in the following singular perturbation analysis, we will relate each small parameter to a \textit{single small parameter} $\epsilon$, based on order of magnitude comparisons. Specifically, we define rescaled parameters $\hat \tau_{max}, \hat K_h, \hat V_s, \hat K$ as follows:
\begin{equation}
{\hat \tau_{max}}=K_\tau^4 \tau_{max}, \ \ {\hat K_h}=K_h K_\tau^{-1}, \ \ {\hat V_s}=V_s K_{\tau}^{-2}, \ \ {\hat K}=K K_{\tau}^{-2}
\label{eq:hatpars}
\end{equation}
With the choice $\sqrt{\epsilon} := K_{\tau}=0.05$, we propose the following relationship between the five small parameters and define the common small parameter $\epsilon$:
\begin{equation}
\label{eq:common_scaling}
\epsilon := \epsilon_\tau^2 = \epsilon_{\tau_{max}}^{1/2}{ \hat \tau_{max}}^{1/2} = \epsilon_h^{2} \hat K_h^{-2} = \epsilon_s \hat V_s^{-1} = \epsilon_K \hat K^{-1} \gamma^{-2}.
\end{equation}
The parameters $\hat \tau_{max}, \hat K_h, \hat V_s, \hat K \gamma^2$ are of numerical order 1, and
$\epsilon = 2.5 \times 10^{-3}$; see Table \ref{tab:params_post_scaling}. In particular, the scaling proposed in \eqref{eq:common_scaling} can be motivated and justified by order of magnitude comparisons.
After applying the common scaling \eqref{eq:common_scaling}, we finally arrive at a version of the model in a form that is amenable to analysis, i.e., \eqref{eq:main2} defined in the following Lemma.

\begin{lemma}
	System \eqref{eq:Non_dim} can be written as
	\begin{equation}\label{eq:main2}
	\begin{split}
	h' & = \hat \tau_{max}^{-1} \left(\mathfrak h_{\infty}({c},\epsilon) - h\right) \left({c}^4 + \epsilon^2 \right), \\
	{c}' &= \mathfrak J_{\text{IPR}}(h,{c},\epsilon) - \epsilon \mathfrak{J}^+_{\text{SERCA}}({c},\epsilon) + \epsilon^2 \mathfrak{J}^-_{\text{SERCA}}({c},\epsilon),
	\end{split}
	\end{equation}
	where
	\begin{equation}
	\label{eq:mathfrak_h}
	\mathfrak h_\infty(c,\epsilon) = \frac{\epsilon^2 \hat K_h^4}{\epsilon^2 \hat K_h^4 + c^4} ,
	\end{equation}
	$\mathfrak J_{\text{IPR}}(h,c,\epsilon)$ is given by \eqref{def:PO} with explicit $\epsilon-$dependence due to $K_\tau = \sqrt{\epsilon}$,
	and
	\begin{align}\label{eq:serca_terms}
	\mathfrak{J}_{\text{SERCA}}^+(c) = \hat V_s \left( \frac{c^2}{K_s^2+c^2} \right), \quad
	\mathfrak{J}_{\text{SERCA}}^-(c) = \hat K \gamma^2 \hat V_s\left( \frac{(c_t - c)^2}{ K_s^2+c^2} \right) .
	\end{align}
	%

	By considering $\epsilon \ll 1$, system \eqref{eq:main2} can be written as the series expansion
	\begin{align}
	\label{eq:main1}
	\begin{split}
	\begin{pmatrix}
	h' \\
	{c}'
	\end{pmatrix}
	=&
	\begin{pmatrix}
	- \hat \tau_{max}^{-1} \\
	\mathfrak{J}^{(0)}_{\text{\text{IPR}}}(h,{c})
	\end{pmatrix}
	{c}^4h + \epsilon
	\begin{pmatrix}
	0 \\
	- \mathfrak J^+_{\text{SERCA}}({c})
	\end{pmatrix}
	\\
	&+ \epsilon^2
	\begin{pmatrix}
	\hat \tau_{max}^{-1}(\hat K_h^4 - h) \\
	\mathfrak{J}^{(1)}_{\text{\text{IPR}}}(h,{c}){c}^8h + \mathfrak J^-_{\text{SERCA}}({c})
	\end{pmatrix}
	+ \epsilon^4
	\begin{pmatrix}
	R_{h_\infty}({c},\epsilon) \\
	R_{\text{\text{IPR}}}(h,{c},\epsilon)
	\end{pmatrix} ,
	\end{split}
	\end{align}
	with IPR terms
	\[
	\mathfrak{J}^{(k)}_{\textrm{\text{IPR}}}(h,c) := \left( \frac{1}{c^{4(k+1)} h}\right) k_{\text{IPR}} P_O^{(k)}(h,c,0) ( \gamma c_t - (1+ \gamma) c), \qquad k = 0, 1,
	\]
	where
	\begin{equation}
	\begin{split}
	P_O^{(0)}(h,c) &= \frac{p^2 c^4 h}{p^2 (1 + k_\beta) c^4 h + k_\beta K_p^2(K_c^4 + c^4)}, \\
	P_O^{(1)}(h,c) &= \frac{k_\beta K_p^2 p^2 c^4 h}{\left(c^4 h(k_\beta+1) p^2+k_\beta K_p^2 \left(c^4+K_c^4\right)\right)^2} .
	\end{split}
	\end{equation}
	The remainder terms $R_{h_\infty}({c},\epsilon)$ and $R_{\text{\text{IPR}}}(h,{c},\epsilon)$ are $O(1)$ as $\epsilon \to 0$.
\end{lemma}
\begin{proof}
	System \eqref{eq:main2} is obtained from system \eqref{eq:Non_dim} after making the substitutions in \eqref{eq:h_scale_tau}, \eqref{eq:Kh_scaling}, \eqref{eq:h_scale_tau_max} and \eqref{eq:c_scale}, applying the common scaling \eqref{eq:common_scaling}, and dividing both sides of the equation for $h'$ by $\tau_h(c)$.

	The expansion in \eqref{eq:main1} follows after significant but standard algebraic manipulations and Taylor expansion in $\epsilon$.
\end{proof}

Systems \eqref{eq:main2} and \eqref{eq:main1} will be considered in detail in Section \ref{sec:slow-fast_analysis_and_statement_of_the_main_result}, where we present a GSPT analysis of the dynamics for $0 < \epsilon \ll 1$. Since only these systems will be considered in the remainder of the manuscript, we will drop the hat notation on rescaled parameters, for the sake of readability. Final system parameters are given in Table \ref{tab:params_post_scaling}.

\begin{table}[t!]
	\center
	\begin{tabular}{|c|l||c|l|}
		\hline
		Parameter         & Value & Parameter            & Value      \\
		\hline
		$k_\beta$ & 0.4   & $\epsilon$     & 0.0025      \\
		$K_c$     & 0.1 & $\hat{V}_s$        & 3.24    \\
		$K_s$     & 0.1 & $\hat K$             & 0.0076 \\
		$K_p$     & 0.1   & $\hat K_h$        & 0.8      \\
		$k_{\text{IPR}}$          & 0.18     & $\gamma$             & 5.5        \\
		$p$       & 0.025  & $\hat{\tau}_{max}$ & 0.34       \\
		$c_t$     & 1  &  -    &   -    \\
		\hline
	\end{tabular}
	\caption{Numerical parameter values for  \eqref{eq:main2} and \eqref{eq:main1} obtained by setting $\epsilon=K_\tau^2=0.0025$ and applying \eqref{eq:hatpars}  to the values of $\tau_{max}, V_s, K$ in Table \ref{tab:params_non_dim}. All other parameter values are unchanged. The hat notation will be dropped in Section \ref{sec:slow-fast_analysis_and_statement_of_the_main_result}.}
	\label{tab:params_post_scaling}
\end{table}

\begin{rem}
	By setting $\epsilon = 2.5 \times 10^{-3}$ in system \eqref{eq:main2}, one obtains precisely system \eqref{eq:Non_dim}. However, systems \eqref{eq:main2} and \eqref{eq:Non_dim} are not equivalent for other choices of $\epsilon \ll 1$.
	Our primary
	assumption in the following is that the dynamics of \eqref{eq:Non_dim} is well captured (qualitatively at least) by considering the asymptotic limit $\epsilon \to 0$
	of \eqref{eq:main2}.
\end{rem}

\begin{comment}
\begin{rem}
	Systems \eqref{eq:main1} and \eqref{eq:main2} are not strictly equivalent to \eqref{eq:Non_dim}, because of assumption \eqref{eq:common_scaling}. It is possible to do have a similar hypothesis to \eqref{eq:common_scaling} and obtain equivalent systems to \eqref{eq:Non_dim}, for example by introducing
	\begin{equation}
		\label{eq:common_scaling_remark}
		\epsilon := \hat K_\tau \epsilon_\tau^2 = \hat K_h  \epsilon_h^2 = \hat V_s \epsilon_s = \hat K \epsilon_K\,,
	\end{equation}
	where
	\begin{equation}
		\hat K_\tau, \hat K_h, \hat V_s, \hat K  = O (1).
	\end{equation}
	However, that would introduce more notations for the following analysis and the overall behaviour of systems \eqref{eq:main1} and \eqref{eq:main2} is similar to \eqref{eq:Non_dim}.
\end{rem}
\end{comment}

\subsection{Timescales}
\label{sec:timescales}

The presence of higher order terms in system \eqref{eq:main1} is suggestive of -- though not sufficient for -- multiple time-scale dynamics when $\epsilon \ll 1$. However, this system is not in the standard form for \wm{multiple time-scale} systems since there is no globally-valid separation of
variables {according to their speed of evolution}. From a modelling point of view, this is a consequence of the underlying physiology, which leads to a time-scale separation that depends on the position in $(h,c)-$phase space. \wm{From a mathematical point of view, the reason the system cannot be recast in standard form is that the singular limit we identify, $\epsilon \to 0$, is non-uniform in phase space due to the (non-uniform) limits $\epsilon_\tau, \epsilon_h \to 0$ implicit in the limit $\epsilon \to 0$. For example, as $\epsilon \to 0$ the function $\mathfrak h_\infty(c,\epsilon)$ in \eqref{eq:mathfrak_h} tends to $0$}
{when $c>0$ but tends to $1$ when $c=0$.}

The factorisation by $c^4h$ in the leading order term in system \eqref{eq:main1} in particular is expected to lead to multiple timescale dynamics, with a slower evolution in regions close to $c=0$ and $h=0$. Here we identify three dynamically important timescales, by a heuristic analysis to be formalised in later sections:
\smallskip
	\begin{enumerate}
	\item[(i)] $t = O(1)$ fast dynamics, in region $h, c = O(1)$. Bounded away from $c=0$ and $h=0$, the leading order dynamics is governed by
	\begin{equation}
	\begin{pmatrix}
	h' \\
	c'
	\end{pmatrix}
	=
	\begin{pmatrix}
	-\tau_{max}^{-1}\\
	\mathfrak{J}^{(0)}_{\text{\text{IPR}}}(h,c)
	\end{pmatrix}
	c^4h .
	\end{equation}
	IPR processes and the negative feedback of the calcium release on the gating variable $h$ dominate the dynamics in this regime.

\smallskip

	\item[(ii)] $t = O(\epsilon)$ slow dynamics, in region $h \sim 0, \ c = O(1)$. Near $h = 0$ but bounded away from $c=0$, the leading order dynamics is governed by
	\begin{equation}
	\begin{pmatrix}
	h' \\
	c'
	\end{pmatrix}
	=
	\epsilon
	\begin{pmatrix}
	0 \\
	- \mathfrak J^+_{\text{SERCA}}(c)
	\end{pmatrix}.
	\end{equation}
	The clearance of calcium ions by the SERCA pump  dominates the calcium dynamics in this regime.

\smallskip

	\item[(iii)] $t = O(\epsilon^2)$ infra-slow dynamics, in region $c \sim 0$. Near $c = 0$ (for any $h$), the leading order dynamics is governed by
	\begin{comment}
	\begin{equation}
	\begin{pmatrix}
	h' \\
	c'
	\end{pmatrix}
	=
	\epsilon^2
	\begin{pmatrix}
	\hat \tau_{max}^{-1} (\hat K_h^4 - h) \\
	\mathfrak{J}^{(2)}_{\text{\text{IPR}}}(h,c)c^8h + \mathfrak J^-_{\text{SERCA}}(c)
	\end{pmatrix}.
	\end{equation}
	\end{comment}
	\begin{equation}
	\begin{pmatrix}
	h' \\
	c'
	\end{pmatrix}
	=
	\epsilon^2
	\begin{pmatrix}
	\tau_{max}^{-1} (K_h^4 - h) \\
	\mathfrak J^-_{\text{SERCA}}(0)
	\end{pmatrix},
	\end{equation}
	where $\mathfrak J^-_{\text{SERCA}}(0) = K \gamma^2 V_s c_t^2 / K_s^2$. In this regime the leak from the ER associated with the SERCA is expected to dominate the calcium dynamics.
	\end{enumerate}

	\smallskip

These observations are supported by numerical simulations, in which distinct timescales are observed for (i) $h, c =  O(1)$, (ii) $h \sim 0, c = O(1)$, and (iii) $c \sim 0$; see again Figure \ref{fig:non_dim}. Remaining terms in \eqref{eq:main1} of order $O(\epsilon^n)$, $n > 3$, constitute regular perturbations.

\section{\wm{Multiple time-scale} analysis and statement of the main result}
\label{sec:slow-fast_analysis_and_statement_of_the_main_result}

In this section we present a singular perturbation analysis of syste
~\eqref{eq:main1}.
System \eqref{eq:main1} can be written in the general form
\begin{equation}\label{eq:general_form}
\begin{pmatrix}
h' \\
c'
\end{pmatrix}
=g(h,c)+\sum_{j\geq 1} \epsilon^j W_j(z) = N(h,c)f(h,c)+\epsilon G(h,c,\epsilon), \qquad 0 < \epsilon \ll1,
\end{equation}
where $g(h,c) := N(h,c)f(h,c)$, with
\[
N(h,c) =
\begin{pmatrix}
- \tau_{max}^{-1} \\
\mathfrak{J}^{(0)}_{\text{IPR}}(h,c)
\end{pmatrix},
\qquad
f(h,c)=c^4h,
\]
and $G(h,c,\epsilon) := \sum_{j\geq 1} \epsilon^j W_j(z)$ is given by
\[
G(h,c,\epsilon) =
\begin{pmatrix}
0 \\
- \mathfrak J^+_{\text{SERCA}}(c)
\end{pmatrix}
+ \epsilon
\begin{pmatrix}
\tau_{max}^{-1}(K_h^4 - h) \\
\mathfrak{J}^{(1)}_{\text{IPR}}(h,c)c^8h + \mathfrak J^-_{\text{SERCA}}(c)
\end{pmatrix}
+ \epsilon^3
\begin{pmatrix}
R_{h_\infty}(c,\epsilon) \\
R_{\text{IPR}}(h,c,\epsilon)
\end{pmatrix}.
\]

\begin{rem} Singular perturbation problems in the general form \eqref{eq:general_form} frequently arise in applications, and have been studied using a combination of GSPT and blow-up techniques in, e.g.,~\cite{Gucwa2009,Huber2005,kosiuk2016geometric,KuehnSzmolyan2015,Lizarraga2020,Lizarraga2020b}. For a formal introduction to the use of GSPT in analysing problems in the general form \eqref{eq:general_form}, we refer the reader to \cite{jelbart2019two} for the planar case, and to \cite{goeke2014constructive,wechselberger2018geometric} for the general (dimension $n\geq 2$) case. See also \cite{Kruff2019,Lizarraga2020b,Lizarraga2020c} for application to systems with $m>2$ timescales.
\end{rem}

Due to the non-uniformity of the singular limit, separate analyses of the dynamics for $\epsilon \to 0$ are required in regimes (R1) and (R2). Note that following from \eqref{eq:common_scaling}, regime (R2) can defined in terms of $\epsilon$ as follows:

\smallskip

\begin{enumerate}
	\item[(R2)] $c \in \mathcal I_2=[0,\sqrt{\epsilon} \nu]$ for fixed $\nu>0$, so that $c = O(\sqrt{\epsilon})$ as $\epsilon \to 0$.
\end{enumerate}

\smallskip

\noindent
Singular limit analyses in both regimes (R1) and (R2) will provide sufficient information for the statement of our main result in Section \ref{ssec:statement_of_the_main_result}, i.e.~the existence of an attracting relaxation cycle in \eqref{eq:main1}. The observed dynamics and corresponding analysis is qualitatively similar, but not identical, to the autocatolator model considered in \cite{Gucwa2009}. \SJ{Similar systems have also been studied in (R2) regime only in \cite{Kuehn2014}.}

\subsection{\wm{Multiple time-scale} analysis in regime (R1)}
\label{ssec:singular_limit_analysis_R1}

We first consider the dynamics in regime (R1). Setting $\epsilon=0$ in \eqref{eq:general_form} gives the \emph{layer problem},
\begin{equation}\label{eq:layer_R1}
\begin{pmatrix}
h' \\
c'
\end{pmatrix}
=
N(h,c) f(h, c)
=
\begin{pmatrix}
- \tau_{max}^{-1} \\
\mathfrak{J}^{(0)}_{\text{IPR}}(h,c)
\end{pmatrix}c^4h,
\end{equation}
which has two lines of equilibria, or \emph{critical manifolds}
\begin{equation}
\label{eq:manifolds_R1}
S_c = \left\{(0, c) : c \geq 0 \right\}, \qquad S_h = \left\{(h, 0) : h \geq 0 \right\} .
\end{equation}

\begin{comment}
\begin{rem}
	The presence of critical manifolds $S_c$ and $S_h$ \textit{define} system \eqref{eq:manifolds_R1} as a singular (as opposed to regular) perturbation problem in the GSPT sense, originally defined in Fenichel's seminal work \cite{fenichel1979geometric}.
\end{rem}
\end{comment}

\begin{rem}
	In contrast to standard form slow-fast systems (see,  e.g.,~\cite{Kuehn2015}),
	there is no globally defined slow variable in system \eqref{eq:main1}. Hence the layer problem \eqref{eq:layer_R1} is \emph{not} a bifurcation problem in either $h$ or $c$.
\end{rem}

\begin{comment}
\begin{rem}
	\label{rem:auxilliary_system}
	Orbits of the layer problem \eqref{eq:layer_R1} are described by the auxiliary system $(h',c')^T = N(h,c)$, which is orbitally equivalent to system \eqref{eq:layer_R1} on $\{(h,c):h>0,c>0\}$ via a transformation of time $d\hat t = f(h,c) dt$.
\end{rem}
\end{comment}

Evaluating the Jacobian at any given $p\in S_c$ or $S_h$ yields a trivial eigenvalue $\lambda_0 \equiv 0$, whose corresponding eigenvector spans the corresponding tangent space at $p$. Direct calculations
\wm{provide the} associated non-trivial eigenvalue $\lambda_c(c)$ (resp.~$\lambda_h(h)$) \wm{of $S_c$ (resp.~$S_h$)}, given by the Lie derivative expressions
\begin{equation}
\label{eq:evs_R1}
\lambda_c(c) = \langle \nabla f,N\rangle \big |_{S_c} = - \tau_{max}^{-1} c^4 \leq 0 ,  \qquad
\lambda_h(h) = \langle \nabla f,N\rangle \big |_{S_h} \equiv0,
\end{equation}
see also \cite[eqn.~(3.7)]{jelbart2019two}.
\begin{comment}
\begin{lemma} (\cite{jelbart2019two,wechselberger2018geometric}).
	Given a planar singular perturbation problem in the general form \eqref{eq:general_form} with critical manifold $S = \{(x,y) : f(x, y) = 0 \}$. Then the \wm{non-uniform} eigenvalue at a point $p \in S$ is given by the inner product
	\begin{equation}\label{eq:nontrivial_ev}
	\lambda(p) = \left\langle \nabla f(p), N(p) \right\rangle .
	\end{equation}
	The corresponding eigenvector lies within $\mathcal N_p = span (N(p))$.
\end{lemma}

\begin{proof}
	Linearisation; see also \cite[eqn.~(3.7)]{jelbart2019two}.
\end{proof}
\end{comment}
\begin{comment}
\begin{rem}
	\label{rem_contact}
	In the case that $p\in S$, $\nabla f(p) \neq 0$ and $N(p) \neq 0$, equation \eqref{eq:nontrivial_ev} implies that a loss of normal hyperbolicity occurs at a tangency between the fast dynamics and the critical manifold $S$ (c.f. the loss of normal hyperbolicity at a fold point in standard form problems \eqref{eq:stnd_form}).
\end{rem}
\end{comment}
Hence, the manifold $S_c$ is normally hyperbolic and attracting for all $c > 0$, and degenerate at $(0,0)$. The manifold $S_h$ is degenerate.

Fast fibers constitute non-trivial heteroclinic connections between $S_h$ and $S_c$ in the layer problem \eqref{eq:layer_R1}. In particular, they obey the equation
\begin{equation}\label{eq:fibers_R1}
\frac{\text{d}c}{\text{d}h} = - \tau_{max} \mathfrak{J}^{(0)}_{\text{IPR}}(h,c)
\begin{cases}
< 0 \qquad & \text{if} \qquad c < \left(\frac{\gamma}{1+\gamma} \right)c_t, \\
= 0 \qquad & \text{if} \qquad c = \left(\frac{\gamma}{1+\gamma} \right)c_t, \\
> 0 \qquad & \text{if} \qquad c > \left(\frac{\gamma}{1+\gamma} \right)c_t,
\end{cases}
\end{equation}
see Figure \ref{SingLimFig}.

\begin{figure}[t!]
	\centering
		\includegraphics[scale=0.23]{Images/calcium_sing_lim}
		\includegraphics[scale=0.23]{Images/fig_calcium_autocatolator_sing_lim.png}
	\caption{Singular limit dynamics. (Left) The (nonlinear) fast connection from $S_h$ to $S_c$ is determined by the layer problem \eqref{eq:layer_R1} and the flow on $S_c$ is given by the reduced problem \eqref{eq:reduced_R1}. The manifold $S_h$ is degenerate. \wm{(Right) The singular orbit $\Gamma=\Gamma_h \cup \Gamma_c \cup \Gamma_l$. Note, the flow indicated on $S_h$ is not formally defined, and included for illustrative purposes only; it represents the rightward motion along $\mathcal S_{a,\delta}$ in Figure \ref{fig:manifold_R2}, and helps to illustrate the three-timescale structure of the oscillations.}}\label{SingLimFig}
\end{figure}

\begin{rem}
	\label{rem:return}
	The expressions in \eqref{eq:fibers_R1} are not on their own sufficient to prove that orbits of the layer problem connect points on $S_h$ and $S_c$. A proof of this property is deferred to Appendix \ref{app:proof_theorem} Lemma \ref{lem:map12}, and omitted here for expository reasons.
\end{rem}

\wm{Next, we consider \eqref{eq:general_form} on the slow timescale $\tau = \epsilon t$, i.e.,
\begin{equation}\label{eq:general_form_slow}
\epsilon
\begin{pmatrix}
\dot{h} \\
\dot{c}
\end{pmatrix}
= N(h,c)f(h,c)+\epsilon G(h,c,\epsilon), \qquad 0 < \epsilon \ll1.
\end{equation}
The following result (see \cite{jelbart2019two}) considers the limit $\epsilon \to 0$ of \eqref{eq:general_form_slow}  called the \textit{reduced problem} which describes the leading order slow flow on $S_c$:}

\begin{lemma}
	\label{lem:reduced}
	Given a planar singular perturbation problem in the general form \eqref{eq:general_form}, with critical manifold $S$, and normally hyperbolic submanifold $S_n \subset S$. Then the reduced problem on $S_n$ is given by
	\begin{equation}\label{eq:reduced_eqn}
	\begin{pmatrix}
	\dot h \\
	\dot c
	\end{pmatrix}
	=
	\left[\frac{\det (N|G)}{\left\langle \nabla f, N \right\rangle}
	\begin{pmatrix}
	-D_cf \\
	D_hf
	\end{pmatrix}\right]
	\bigg|_{S_n} ,
	\end{equation}
	where the $(\dot{\ })$ notation denotes differentiation with respect to slow time $\tau$, and the determinant $\det (N|G)$ is taken with respect to the matrix with columns $N(h,c)$ and $G(h,c,0)$.
\end{lemma}

\begin{proof}
	The expression in \eqref{eq:reduced_eqn} is obtained by projecting the leading order perturbation vector $G(h,c,0)$ at each point $p \in S_c \setminus \{(0,0)\}$ onto its component in the corresponding tangent space $T_pS_c$; see \cite[Proposition 3.4]{jelbart2019two}; \wm{see also \cite{fenichel1979geometric,goeke2014constructive,wechselberger2018geometric} for similar methods applicable in the} case that the dimension $n$ is  greater than two.
\end{proof}

\begin{comment}
\begin{rem}
	Lemma \ref{lem:reduced} applies only for normally hyperbolic submanifolds of the critical manifold $S$. There is no well-defined reduced problem on the degenerate line $S_h$, which by expression \eqref{eq:evs_R1} is everywhere nonhyperbolic.
\end{rem}
\end{comment}

\begin{comment}
\begin{equation}\label{eq:reduced_R11}
\begin{pmatrix}
\dot h \\
\dot c
\end{pmatrix}
=
\left[\frac{1}{\lambda_h(c)}\det (N|G)
\begin{pmatrix}
-D_cf \\
D_hf
\end{pmatrix}\right]
\bigg|_{h=\epsilon=0}
=
\begin{pmatrix}
0 \\
- \mathfrak J^+_{SERCA}(c)
\end{pmatrix},
\end{equation}
\
where the overdot notation denotes differentiation with respect to slow time $\tau$, and the determinant $\det (N|G)$ is taken with respect to the matrix with columns $N(h,c)$ and $G(h,c,0)$.
\end{comment}
We obtain an expression for the reduced problem on $S_c \setminus \{0,0\}$ directly via equation \eqref{eq:reduced_eqn}, which gives
\begin{equation}\label{eq:reduced_R1}
\begin{pmatrix}
\dot h \\
\dot c
\end{pmatrix}
=
\begin{pmatrix}
0 \\
- \mathfrak J^+_{\text{SERCA}}(c)
\end{pmatrix} .
\end{equation}
In particular, $- \mathfrak J^+_{\text{SERCA}}(c) < 0$ for all $c >0$ and $- \mathfrak J^+_{\text{SERCA}}(0) = 0$; the reduced flow on $S_c$ is toward the origin. The origin itself can be considered as a non-hyperbolic equilibrium for the reduced problem extended to all of $S_c$; see Figure \ref{SingLimFig}.

\begin{comment}
\begin{rem}
	In contrast to standard form problems \eqref{eq:stnd_form}, one does not in general obtain a differential-algebraic problem from the system \eqref{eq:general_form} by first rewriting it on the slow timescale $\tau = \epsilon t$ and setting $\epsilon=0$. The usual `total derivative approach' (see, e.g.~\cite{Kuehn2015}) for calculating reduced flows will not apply in this more general setting.
\end{rem}
\end{comment}

\

The preceding analysis for both layer and reduced problems implies the following result
for $0 < \epsilon \ll 1$, which follows by classical results due to Fenichel \cite{fenichel1979geometric}.

\begin{lemma}\label{lem:slow_manifolds_R1}
	There exists an $\epsilon_0>0$ such that for all $\epsilon\in(0,\epsilon_0]$, compact submanifolds of $S_c$
	perturb to $O(\epsilon)$-close locally invariant \wm{one-}dimensional slow manifold
	\[
	S_{c,\epsilon}=\left\{(\tilde{\nu}(c,\epsilon),c):c\in[c_-,c_+]\right\},
	\]
	where $\tilde \nu(c,\epsilon)=O(\epsilon)$, and $c_\pm$ are positive constants satisfying $c_+ > c_- > 0$. The leading order flow on each $S_{c,\epsilon}$ is governed by the reduced problem \eqref{eq:reduced_R1}.
\end{lemma}

It remains to understand the dynamics near the degenerate line $S_h$. For this we must look in regime (R2).

\subsection{\wm{Multiple time-scale} analysis in regime (R2)}
\label{ssec:singular_limit_analysis_R2}

We now consider the dynamics in regime (R2), where $c = O(\sqrt{\epsilon})$.
\wm{
Defining a rescaled variable $C=O(1)$ via
$$
c=\delta C\,,\qquad\mbox{where }\delta := \sqrt{\epsilon}
$$
and introducing the `intermediate-slow' timescale $t_1 = \epsilon^{3/2} t = \delta^3 t$ in system  \eqref{eq:general_form},}
we obtain
\begin{equation}
\label{eq:main_R2}
\begin{pmatrix}
h' \\
C'
\end{pmatrix}
=
\tilde N(h,C) \tilde f_0(h,C) + \delta \tilde G(h,C,\delta), \qquad 0 < \delta \ll 1 ,
\end{equation}
where
\begin{equation}
\tilde N(h,C) =
\begin{pmatrix}
0 \\
1
\end{pmatrix},
\qquad
\tilde f_0(h,C) = A_{\text{\text{IPR}}} \gamma c_t C^4 h - A_{\text{SERCA}} \left(C^2 - c_t^2 \right) ,
\end{equation}
and
\begin{equation}
\label{eq:K2_rhs}
\tilde G(h,c,\delta) =
\begin{pmatrix}
\tilde g(h,C) \\
\tilde f_{rem}(C,h,\delta)
\end{pmatrix}
=
\begin{pmatrix}
 - \tau_h(C)^{-1} \left(h - h_\infty (C) \right)  \\
 - \left(A_{\text{\text{IPR}}} (1 + \gamma) C^4 h +2 A_{\text{SERCA}} c_t \right) C +  O\left(\delta\right)
\end{pmatrix} .
\end{equation}
In the above we defined
\[
A_{\text{\text{IPR}}} = \left(\frac{k_{\text{IPR}} p^2}{k_\beta K_c^4 K_p^2}\right) , \qquad
A_{\text{SERCA}}=\frac{K \gamma^2V_s}{K_s^2},
\]
and permitted a slight abuse of notation by letting the prime ($'$) denote differentiation with respect to the new \wm{intermediate-slow  time $t_1$}. In writing \eqref{eq:K2_rhs}, we also appeal to the fact that $\mathfrak h_\infty(\delta C, \delta^2) = h_\infty(C)$. Finally, note that $\tilde f_0(h,C)$ is just the leading order term in an expansion in $\delta$, and that in this regime,
IPR and SERCA processes compete at leading order.

\begin{rem}
	System \eqref{eq:main_R2} is in the so-called standard form for slow-fast systems. We present it here in the general form \eqref{eq:general_form}, and

	proceed via the same (more general) approach

	adopted in Section \ref{ssec:singular_limit_analysis_R1}. Our reasons for doing so are three-fold: (i) consistency with earlier sections is maintained; (ii) there are no additional technical difficulties, and (iii) this approach illustrates the relationship between the theory developed for standard form problems and its non-standard form generalisation.
\end{rem}

Setting $\delta = 0$ in \eqref{eq:main_R2} yields the layer problem
\begin{equation}
\label{eq:layer_R2}
\begin{pmatrix}
h' \\
C'
\end{pmatrix}
=
\tilde N \tilde f_0(h,C)
=
\begin{pmatrix}
0 \\
1
\end{pmatrix}
\left(A_{\text{\text{IPR}}} \gamma c_t C^4 h - A_{\text{SERCA}} \left(C^2 -  c_t^2 \right)\right) ,
\end{equation}
which has
critical manifold
\begin{equation}
\label{eq:S_R2}
\mathcal S = \left\{\left(\zeta(C), C \right) : C \geq 0 \right\} , \qquad
\zeta(C) = \frac{A_{\text{SERCA}}}{A_{\text{IPR}} \gamma c_t C^4}(C^2-\gamma^2c_t^2 K) ,
\end{equation}
see Figure \ref{fig:manifold_R2}. The non-trivial eigenvalue along $\mathcal S$ is given by the Lie derivative expression
\begin{equation}
\label{EV_R2}
\lambda_{\mathcal S}(C) = \left\langle \nabla f, N \right\rangle \big|_{\mathcal S} =  D_C\tilde f_0(\zeta(C), C) = \frac{2 A_{\text{SERCA}}}{C}\left(C^2-2 K\gamma^2 c_t^2\right) .
\end{equation}
The manifold $\mathcal S$ has a fold point at
\begin{equation}
\label{eq:Fold}
F = (\zeta(C_F),C_F) = \left(\frac{A_{\text{SERCA}}}{4 K \gamma^3 c_t^3 A_{\text{IPR}}}, \sqrt{2 K}\gamma c_t  \right),
\end{equation}
where $\lambda_{\mathcal S}(C_F)=0$ and the following nondegeneracy conditions are satisfied:
\begin{equation}\label{eq:fold_nondegeneracy}
D_C^2\tilde f_0(\zeta(C_F), C_F) = 4 A_{\text{SERCA}} > 0, \ \ \
D_h\tilde f_0(\zeta(C_F), C_F) = 4 A_{\text{IPR}} K^2 \gamma^5 c_t^5 >0.
\end{equation}
The first inequality in \eqref{eq:fold_nondegeneracy} implies that $\mathcal S=\mathcal S_a\cup\{F\}\cup \mathcal S_r$, where
\[
\mathcal S_a=\left\{(\zeta(C),C):0<C<C_F\right\}, \qquad \mathcal S_r=\left\{(\zeta(C),C):C>C_F\right\},
\]
and $\mathcal S_a$ (resp.~$\mathcal S_r$) is normally hyperbolic and attracting (resp.~repelling), as sketched in Figure \ref{fig:manifold_R2}.

\begin{comment}
\begin{rem}
	Recall by Remark \ref{rem_contact} that normal hyperbolicity is lost at points of tangency between the layer flow and the critical manifold. The loss of normal hyperbolicity at a regular fold point $F \in \mathcal S$ in particular is due to a \textit{quadratic} tangency with between the layer flow and the critical manifold $\mathcal S$. For problems given in the general form \eqref{eq:general_form}, such points are more generally referred to as \textit{order one contact points}. These are identified as points satisfying $\lambda_{\mathcal S}(F) = 0$, together with nondegeneracy conditions which reduce to those in \eqref{eq:fold_nondegeneracy} in the case of standard form problems \eqref{eq:stnd_form} (i.e. in the case of a regular fold); see \cite{jelbart2019two} for the dimension $n=2$ case, and \cite{wechselberger2018geometric} for the general case.
\end{rem}
\end{comment}

\begin{figure}[t!]
	\centering
		\includegraphics[scale=0.35]{Images/chart_K2_b.png}
	\caption{Sketch of the dynamics in regime (R2). The critical manifold $\mathcal S=\mathcal S_a\cup\{F\}\cup \mathcal S_r$ is shown with $\mathcal S_a$ in blue and $\mathcal S_r$ in dashed red. The fold point $F$, equilibrium $q$ and $h-$nullcline (in grey) are also indicated. Perturbed slow manifolds $\mathcal S_{a,\delta}$ and $\mathcal S_{r,\delta}$ are shown in shaded blue and red respectively, with the attracting slow manifold $\mathcal S_{a,\delta}$ extending through the neighbourhood of the fold point $F$ as described by Lemma \ref{lem:slow_manifolds_R2}.}
	\label{fig:manifold_R2}
\end{figure}

Now consider the reduced problem on the manifold $\mathcal S$, \wm{which describes the leading order dynamics on the infra-slow timescale $\tau_1 = \delta t_1 = \epsilon^2 t$}. As in Section \ref{ssec:singular_limit_analysis_R1}, we use the formula \eqref{eq:reduced_eqn} to \wm{derive} the reduced problem,
\begin{equation}\label{eq:reduced_R2}
\begin{pmatrix}
\dot h \\
\dot C
\end{pmatrix}
=
\left[\frac{\det (\tilde N| \tilde G)}{\langle \nabla \tilde f_0, \tilde N \rangle}
\begin{pmatrix}
-D_c \tilde f_0 \\
D_h \tilde f_0
\end{pmatrix}\right]
\bigg|_{\mathcal S}
=
\begin{pmatrix}
1 \\
-\frac{D_h\tilde f_0(\zeta(C),C)}{\lambda_{\mathcal S}(C)}
\end{pmatrix}
\tilde g(\zeta(C),C),
\end{equation}
where by another (slight) abuse of notation, $(\dot{\ })$ refers to differentiation with respect to
the infra-slow timescale $\tau_1$. Explicitly, we obtain the system
\begin{equation}\label{eq:reduced_R2_explicit}
\begin{pmatrix}
\dot h \\
\dot C
\end{pmatrix}
=
\begin{pmatrix}
- \tau_h(C)^{-1} \\
\frac{A_{\text{IPR}} \gamma c_t  C^5}{2 A_{\text{SERCA}} \tau_h(C) (C^2 - 2 K \gamma^2 c_t^2)}
\end{pmatrix}
\left(\zeta(C) - h_\infty (C) \right) .
\end{equation}

\begin{rem}
	\label{rem:timescales}
	In total we have identified \emph{four} \wm{distinct} timescales in system \eqref{eq:general_form} that are involved in our analysis \wm{(fast, slow, intermediate-slow, infra-slow)}:
	\[
	t
	\qquad \tau = \epsilon t,
	\wm{
	\qquad t_1 = \delta^3 t = \epsilon^{3/2} t ,
	\qquad \tau_1 = \delta t_1 = \epsilon^2 t .
	}
	\]

\end{rem}

\begin{comment}
\begin{rem}
	Equation \eqref{eq:reduced_R2} shows that the `usual' expression for the reduced problem associated with standard form problems \eqref{eq:stnd_form} (see, e.g. \cite{Kuehn2015}) is obtained directly from the more general expression \eqref{eq:reduced_eqn}.
\end{rem}
\end{comment}

System \eqref{eq:reduced_R2_explicit} can have up to three equilibria $(\varphi(C_{\ast}), C_{\ast}) \in \mathcal S$ in the physiological domain $h, C \geq 0$. Their location can be determined by solving the equation
\[
\zeta(C) = h_\infty (C) ,
\]
which reduces to the problem of identifying positive, real-valued roots of the cubic equation
\begin{equation}\label{eqs:reduced_equilibria_R2b}
P(m) = \alpha_3 m^3 + \alpha_2 m^2 + \alpha_1 m + \alpha_0 = 0,
\end{equation}
obtained after setting $m=C^2$.
\wm{
The coefficients $\alpha_i$, $i=0,1,2,3$ are given by}
\begin{equation}
\begin{aligned}
\alpha_3 &= \frac{A_{\text{SERCA}}}{\gamma c_t A_{\text{\text{IPR}}}}, && \alpha_2 = - K_h^4 - \frac{A_{\text{SERCA}}}{ A_{\text{\text{IPR}}}} K \gamma c_t , \\
\alpha_1 &= \frac{K_h^4 A_{\text{SERCA}}}{\gamma c_t A_{\text{\text{IPR}}}} , && \alpha_0 = - \frac{K K_h^4 \gamma c_t A_{\text{SERCA}}}{A_{\text{\text{IPR}}}} .
\end{aligned}
\end{equation}
Note that the number of real-valued roots for equation \eqref{eqs:reduced_equilibria_R2b} depends on the system parameters $p$ and $c_t$; \wm{see Figure~\ref{fig:cusp}, which identifies a cusp structure in $(p,c_t)$-space.}

\begin{figure}[t!]
	\centering
	\includegraphics[scale=0.8]{Images/p_ct_diag

	\caption{\wm{Number of equilibria in $(p,c_t)$-parameter space. The region with three equilibria  is bounded by two codimension-one saddle-node branches emanating from a codimension-two cusp point at $(p,c_t)\sim (0.02,0.8)$; outside this region there exists a unique equilibrium. Parameter values are as in Table \ref{tab:params_post_scaling}}.}
	\label{fig:cusp}
\end{figure}

Of course, the location and stability properties of the corresponding equilibria for system \eqref{eq:reduced_R2} will also depend on these parameters. We are interested
{here in the case in which} system \eqref{eq:reduced_R2} has a unique, unstable equilibrium $q \in \mathcal S_r$.

\begin{assumption}
	\label{ass:1}
	The parameters $p$, $c_t$ are chosen so that the cubic polynomial $P(m)$ in \eqref{eqs:reduced_equilibria_R2b} has a unique positive real-valued root, such that system \eqref{eq:reduced_R2} has a unique equilibrium $q = (\varphi(C_{\ast}), C_{\ast}) \in \mathcal S$. Moreover, we require that
	\begin{equation}\label{eq:stability_bounds}
	\frac{2 A_{\text{IPR}} K_h^4 \gamma c_t}{A_{\text{SERCA}}} > \frac{1}{C_\ast^8}\left(K_h^4 + C_\ast^4 \right) \left(C_\ast^2 - C_F^2 \right) > 0 ,
	\end{equation}
	implying that $q$ is unstable on $\mathcal S_r$.
\end{assumption}

The rightmost inequality in equation \eqref{eq:stability_bounds} implies that $C_\ast > C_F$, so that $q \in \mathcal S_r$. The leftmost inequality ensures that $D_c \dot C|_{C=C_\ast} > 0$, so that $q$ is unstable as an equilibrium on $\mathcal S$.
\begin{comment}
Explicitly, a little algebra yields the requirement for instability
\begin{equation}
\label{eq:stability_assumption}
D_c \dot C \big|_{C=C_\ast}
= \frac{1}{\tau_h(C_\ast)} \left(\frac{2 A_{\text{IPR}} K_h^4 C_\ast^8 \gamma c_t}{ A_{SERCA} (K_h^4 + C_\ast^4) (C_\ast^2 - C_F^2)} - 1 \right) > 0 .
\end{equation}
\end{comment}
Note that Assumption \ref{ass:1} is satisfied for the parameter values in Table \ref{tab:params_post_scaling}, for which we have a unique equilibrium $q \approx (0.074,1.48) \in \mathcal S_r$ satisfying $C_\ast \approx 1.48 > C_F \approx 0.68$ and
\[
\frac{2 A_{\text{IPR}} K_h^4 \gamma c_t}{A_{\text{SERCA}}} \approx 17.01 > \frac{1}{C_\ast^8}\left(K_h^4 + C_\ast^4 \right) \left(C_\ast^2 - C_F^2 \right) \approx 0.391 > 0 .
\]
Assumption \ref{ass:1} also implies
\wm{that the regularity conditions
\begin{equation}\label{eq:fold_regularity_R2}
\tilde g(\zeta(C_F),C_F)\neq0, \qquad D_C \tilde g(\zeta(C_F),C_F)\neq0 ,
\end{equation}
at the fold point $F$ are fulfilled.

Note that} the reduced problem \eqref{eq:reduced_R2} is not defined at $F$, since the flow undergoes a finite time blow-up (solutions reach $F$ in finite time). The conditions in \eqref{eq:fold_regularity_R2} imply that the reduced flow is oriented locally toward $F$, classifying the regular fold as a \emph{regular jump point} \cite{Kuehn2015}, where there is a transition from slow to fast.

\

The dynamics for $0 < \delta \ll 1$ in compact subdomains in regime (R2) is well-described by existing theory, and summarised in the following result.

\begin{lemma}\label{lem:slow_manifolds_R2}
	Given system \wm{\eqref{eq:main_R2}}
	and Assumption \ref{ass:1}, there exists $\delta_0>0$ such that for all $\delta\in(0,\delta_0]$ compact submanifolds
	of $\mathcal S_a$ and $\mathcal S_r$ perturb to $O(\delta)$-close locally invariant slow manifolds of the form
	\[
	\mathcal S_{a,\delta} = \left\{(\zeta(C)+O(\delta),C):C\in[C_-,C_+]\right\}, \ \
	\mathcal S_{r,\delta} = \left\{(\zeta(C)+O(\delta),C):C\in[C_l,C_r]\right\},
	\]
	respectively, where $C_\pm$ and $C_{l,r}$ are any positive constants satisfying $0<C_-<C_+<C_F$ and $C_F<C_l<C_r<\infty$. The leading order flow on $\mathcal S_{a/r,\delta}$ is given by the reduced flow on $\mathcal S_{a/r}$, and the reduced flow equilibrium $q \in \mathcal S_r$ perturbs to a nearby equilibrium $q_\delta \in \mathcal S_{r,\delta}$ of unstable node type.

	The attracting slow manifolds $\mathcal S_{a,\delta}$ extend through a neighbourhood of the regular jump point, exiting along the fast flow which is $O(\delta^{2/3})$-close to the fast fiber $\{(\zeta(C_F),C):C>C_F\}$; see Figure \ref{fig:manifold_R2}.
\end{lemma}

\begin{proof}
	The dynamics bounded away from a neighbourhood of $F$, including the persistence of the unstable node $q_\delta$, is described by Fenichel theory \cite{fenichel1979geometric}. The dynamics near $F$ is described by \cite[Theorem 2.1]{Krupa2001a}.
\end{proof}

It follows that trajectories are either repelled to infinity or exponentially attracted to slow manifolds $\mathcal S_{a,\delta}$, after which they follow the slow flow up to the neighbourhood of the regular jump point $F$, before
leaving via the fast flow. Once on the fast flow, the global separation of slow and fast variables in regime (R2) prohibits the existence of a return mechanism on any compact domain.
Limit cycles for $0 < \delta \ll 1$ must traverse both regimes (R1) and (R2).

\subsection{Statement of the main results}
\label{ssec:statement_of_the_main_result}

We are now able to state our main result: existence of a limit cycle of relaxation type in system \wm{\eqref{eq:general_form}} in the $(p,c_t)-$parameter region specified by Assumption \ref{ass:1}.

We consider system \wm{\eqref{eq:general_form}} in regime (R1), for which the fold point $F$ identified in regime (R2) `collapses' onto the point $Q_f = (h_F,0) = (\zeta(C_F),0)\in S_h$.
Having identified this point, we can construct the singular orbit
\[
\Gamma=  \Gamma_c \cup \Gamma_h \cup \Gamma_l ,
\]
where
\[
\Gamma_h=\left\{(h,0):h\in[0,h_F] \right\}, \qquad \Gamma_c=\left\{(0,c):c\in[0,c_d] \right\},
\]
and $\Gamma_l$ is the (unique) heteroclinic orbit of the layer problem \eqref{eq:layer_R1} connecting $Q_f$ to the `drop point' $Q_d = (0, c_d) \in S_c$; see \wm{Figure~\ref{SingLimFig} (Right).}
It should be noted that because of curvature in the layer flow, existence of a drop point $Q_d \in S_c$ must be shown explicitly (recall Remark \ref{rem:return}).
\begin{comment}
\begin{figure}[t!]
	\centering
		\includegraphics[scale=0.25]{Images/fig_calcium_autocatolator_sing_lim.png}
	\caption{The singular orbit $\Gamma=\Gamma_h \cup \Gamma_c \cup \Gamma_l$. Note that the flow indicated on $S_h$ is not formally defined, and included for illustrative purposes only; it represents the rightward motion along $\mathcal S_{a,\delta}$ in Figure \ref{fig:manifold_R2}, and helps to illustrate the three-timescale structure of the oscillations. \wm{[show true relaxation oscillation, e.g., from Figure 3. Maybe merge with Figure 5?] }}
	\label{fig:singular_orbit_R1}
\end{figure}
\end{comment}

\begin{theorem}\label{thm:main}
	Consider system \wm{\eqref{eq:general_form}} with $(p,c_t)$ fixed within the region specified by Assumption \ref{ass:1}.
	There exists an $\epsilon_0>0$ such that for all $\epsilon\in(0,\epsilon_0)$, there exists a locally unique
	relaxation cycle $\Gamma_\epsilon$ which is $O(\epsilon^{1/3})-$close to the singular orbit $\Gamma$ in the Hausdorff distance as $\epsilon \to 0$. The relaxation cycle $\Gamma_\epsilon$ is exponentially attracting, with Floquet exponent bounded above by $-\kappa/\epsilon^2$ for some constant $\kappa>0$.
\end{theorem}
Theorem \ref{thm:main} states existence for the relaxation cycles observed in Figure \ref{fig:non_dim}. Note that the additional timescale leads to a relaxation cycle $\Gamma_\epsilon$ which is $O(\epsilon^{1/3})$ from $\Gamma$, in contrast to the usual $O(\epsilon^{2/3})$ separation associated with two-timescale relaxation oscillations; see, e.g., \cite{Krupa2001b,Kuehn2015}. This is a consequence of the $O(\delta^{2/3}) = O (\epsilon^{1/3})$ separation near the fold in the regime (R2), as described in Lemma \ref{lem:slow_manifolds_R2}. The bound $-\kappa/\epsilon^2$ on the Floquet exponent also differs from the usual $-\kappa/\epsilon$ bound. This is a consequence of the total time spent in the vicinity of the attracting slow manifold $\mathcal S_{a,\delta}$, which is $O(1)$ on the \wm{infra-slow timescale $\tau_1 = \epsilon^2 t$}.

The proof is deferred to Appendix \ref{app:proof_theorem}, and utilises the so-called \textit{blow-up method} \cite{Dumortier1996} in the formulation of \cite{Krupa2001a,Krupa2001b}. The main task is to resolve the degeneracy associated with the non-hyperbolic line $S_h$ and, in particular, the point $Q_o = (0,0)$. This can be achieved in a two-step process, by means of a cylindrical blow-up along $S_h$, and a second (successive) spherical blow-up necessary to resolve a persistent degeneracy stemming from the point $Q_o$. We refer the interested reader to Figures \ref{fig:cyl_bu} and \ref{fig:sphere_bu} in particular, which illustrate the main dynamical features after blow-up, and further illustrate the similarity with the autocatolator model presented in \cite{Gucwa2009}.

\

Our analysis also allows for a leading order approximation of the oscillation period $\mathcal T$ as a function of the model parameter $\tau_{max}$.

\begin{proposition}\label{prop:period}
	Consider system \wm{\eqref{eq:general_form}} with $(p,c_t)$ fixed within the region specified by Assumption \ref{ass:1}. Then there exists an $\epsilon_0 > 0$ such that for all $\epsilon \in (0,\epsilon_0)$, the leading order approximation for the period $\mathcal T$ of the relaxation cycle $\Gamma_\epsilon$ is linear in $\tau_{max}$. More precisely, written

	in terms of the \wm{fast timescale $t$ of system \eqref{eq:general_form}} we have
	\[
	\mathcal T \sim \wm{\epsilon^{-2}} \tau_{max} v(p, c_t) \qquad \text{as} \qquad \epsilon \to 0,
	\]
	where $v(p, c_t)$ is smooth, positive and bounded on the relevant domain.
\end{proposition}

\begin{proof}

	Assuming Theorem \ref{thm:main} and recalling the discussion in Section \ref{sec:timescales}, the leading order approximation for the period $\mathcal T$ as $\epsilon \to 0$ is determined by the total time spent in the vicinity of the attracting \wm{infra-slow} manifold $\mathcal S_a$,
	which can be approximated using the expression for the corresponding reduced flow given in \eqref{eq:reduced_R2} \wm{which evolves on the infra-slow timescale $\tau_1 = \epsilon^2 t$}.
	{Let $C(\tau_1)$ denote a solution for \eqref{eq:reduced_R2} such that $\zeta(C(0)) = 0$ and $C(\tau_1^0) = C_F$ for some $\tau_1^0 = \epsilon^2 t_0 > 0$. }

	\SJ{This yields $C(0) = \sqrt{K} \gamma c_t$ and $C(\tau_1^0) = \sqrt{2 K} \gamma c_t$, and hence}
	\begin{equation}
	\label{eq:period}
	\epsilon^{2} \mathcal T \sim \int_{0}^{\wm{\tau_1^0}} \wm{d\tau_1}
	= \int_{\sqrt{K} \gamma c_t}^{\sqrt{2 K} \gamma c_t} \frac{1}{\dot C} dC ,
	\end{equation}
	%
	where $\epsilon^{2}$ appears in the left-hand-side since we aim to approximate the period $\mathcal T$ on the \wm{fast timescale $t = \epsilon^{-2} \tau_1$ of system \eqref{eq:general_form}}.

	Using \eqref{eq:reduced_R2}, the integrand can be written as $1/\dot C = \tau_{max} V(C; p , c_t)$, where the function $V(C ; p , c_t)$ is smooth with respect to $C, c_t, p,$ and independent of $\tau_{max}$. By Assumption \ref{ass:1} we have that $\dot C > 0$, and hence $V(C ; p , c_t) > 0$ on $\mathcal S_a$. Thus,
	\[
	\mathcal T \sim \wm{\epsilon^{-2}} \tau_{max} \int_{\sqrt{K} \gamma c_t}^{\sqrt{2 K} \gamma c_t} V(C ; p ; c_t) dC = \wm{\epsilon^{-2}} \tau_{max} v(p, c_t),
	\]
	from a smooth, positive and bounded function $v(p,c_t)$ as $\epsilon \to 0$.

\end{proof}

\begin{remark}\label{rem:num-period}
Figure~\ref{fig:non_dim} shows relaxation oscillations with an approximate period $\mathcal T \approx 2\times 10^{4}$, which compares well to the order estimate $\epsilon^{-2} \tau_{max} \approx 5\times 10^{4}$ derived above and evaluated using the parameter values given in Table~\ref{tab:params_post_scaling}.
\end{remark}

\begin{comment}
\section{Effect of $\tau_{max}$ on the system}
\label{sec:effect_of_tau_}

The analysis we have done so far has been for a specific choice of the parameter $\tau_{max}$. However, this parameter controls the period of oscillations in the model \cite{sneyd2017dynamical} and the period of calcium oscillations can vary significantly depending on the type of cell.
A detailed understanding of the behaviour of the model for values of $\tau_{max}$ on different magnitudes is of physiologically importance. Calcium oscillations can have a period that varies significantly depending on the type of cell. For example, in parotid acinar cells calcium oscillations have periods that are around 5 to 10 seconds \cite{pages2018Anatomically}. Those oscillations are not spiking oscillations and hardly look like relaxation oscillations. By contrast, calcium oscillations in HSY cells consist of spikes occurring with a period of around 40 seconds \cite{sneyd2017dynamical,tanimura1996inositol} and oscillations in hepatocytes can occur with a period of more than a minute \cite{bartlett2014calcium}. To model each cell, increasing values of $\tau_{max}$ can be used. Of course, changes in the period of oscillation may not be uniquely linked to changes in a deterministic process such as that modelled using $\tau_{max}$, and very likely are influenced by the stochastic nature of calcium oscillations \cite{skupin2008does}. Nevertheless, we wish to investigate the effect on our analysis of different choices of $\tau_{max}$, looking at a values of $\tau_{max}$ on orders of magnitude either side of the previously studied parameter regime.
\end{comment}

\section{Onset of oscillations}
\label{sec:effect_of_tau_}

Here, we briefly address the basic mechanisms leading to the onset of oscillations under parameter variation. We focus on three important model parameters: total calcium concentration $c_t$, the IP$_3$ concentration $p$ and the timescale parameter $\tau_{max}$ of the $h$-dynamics.

\subsection{Singular Andronov-Hopf bifurcation and canard explosion under variation of $c_t$ or $p$}
\label{sub:can_explosion}

The first mechanism involves a singular Andronov-Hopf (AH) bifurcation and a corresponding \textit{canard explosion} \cite{Dumortier1996,Krupa2001b,Kuehn2015}
and may arise in system  \eqref{eq:general_form} under variation in either $p$ or $c_t$ (but not under variation of $\tau_{max}$).
The occurrence of this mechanism is correlated with the passage of the equilibrium $q$ through the fold point $F$ in Figure \ref{fig:manifold_R2} under parameter variation. Note, that this implies a violation of the regularity condition \eqref{eq:fold_regularity_R2}.
Results in \cite{Krupa2001b,Dumortier1996} then imply the existence of a nearby singular AH bifurcation for a locally unique parameter value $p = p_{ah} + O(\delta)$ or $c_t = c_{t,ah} + O(\delta)$. We highlight that periodic orbits arising from the singular AH bifurcation
are expected to have an oscillation period $\mathcal T = O(\epsilon^{-7/4})$, which is intermediate between the intermediate-slow and infra-slow timescales $t_1 = \epsilon^{3/2} t$ and $\tau_1 = \epsilon^2t$ respectively.

\begin{figure}[t!]
	\centering
	\begin{subfigure}[b]{0.48\linewidth}
		\includegraphics[width=\linewidth]{Images/Can_Hopf_again}
	\end{subfigure}
	\quad
	\begin{subfigure}[b]{0.485\linewidth}
		\includegraphics[width=\linewidth]{Images/expl_period_again}
	\end{subfigure}
	\caption{{Bifurcation diagram for system \eqref{eq:general_form} under variation of $c_t$, computed numerically in MatCont \cite{MATCONT}, for $p=0.015$  and other parameters as in Table \ref{tab:params_post_scaling}:
	(Left) The onset of oscillations happens via a subcritical singular Andronov-Hopf (AH) bifurcation at $c_t\approx 0.96$.
	The transition to stable relaxation oscillations happens via a canard explosion in an exponentially small parameter regime.
	(Right) Period of the branch of relaxation oscillations.

	}}
	\label{fig:can_exp}
\end{figure}

Figure~\ref{fig:can_exp} shows a bifurcation diagram for system \eqref{eq:general_form} under variation of $c_t$, with $p=0.015$ and other parameters as in Table \ref{tab:params_post_scaling}. As expected, the onset of oscillations happens via a singular AH bifurcation, at $c_t\approx 0.96$. The AH bifurcation is observed to be subcritical; the criticality has been confirmed numerically in MatCont \cite{MATCONT} by showing that the corresponding first Lyapunov coefficient of the AH bifurcation is positive.

The rapid onset of the (non-standard) relaxation oscillations described in Theorem \ref{thm:main} is expected to occur over an exponentially small interval in $c_t$ (or $p$); this is referred to as a \textit{canard explosion} \cite{Dumortier1996,Krupa2001b,Gucwa2009} because of the dramatic (or explosive) growth of amplitude. The rapid onset of oscillations is observed as an almost vertical segment in Figure~\ref{fig:can_exp}.

\begin{rem}
Due to the subcritical nature of the singular AH bifurcation observed in Figure~\ref{fig:can_exp}, there exists a small parameter regime in which there is bistability between a stable relaxation/canard cycle and  a stable equilibrium; on the scale of Figure~\ref{fig:can_exp}, this regime is too small to be seen. The existence of this regime of bistability implies the occurrence of a saddle-node bifurcation of canard cycles (not indicated in Figure~\ref{fig:can_exp}). We refer to, e.g.,~\cite{Krupa2001b,Wechselberger2015} for details.

\end{rem}

\begin{figure}[t!]
	\centering
	\includegraphics[scale=0.62]{Images/ct_diag_again}
	\caption{{Bifurcation diagram for system \eqref{eq:general_form}, computed numerically in MatCont \cite{MATCONT} with parameters as in Table \ref{tab:params_post_scaling}; compare with Figure~\ref{fig:can_exp}.

	The onset of oscillations happens via a subcritical singular AH bifurcation. The transition to stable relaxation oscillations happens via an `incomplete' canard explosion involving canard homoclinic cycles; see, e.g., \cite{Wechselberger2015} for details.}}

	\label{fig:can_explosion}
\end{figure}

Figure~\ref{fig:can_explosion} shows another bifurcation diagram for system \eqref{eq:general_form} under variation of $c_t$, this time with all parameters as in Table \ref{tab:params_post_scaling}. Note the S-shaped branch of equilibria which indicates that the choice of $p=0.025$ puts the system into the regime in which there can be three equilibria;
see Figure~\ref{fig:cusp}. There is a subcritical singular AH bifurcation at $c_t\approx 0.63$.  As before, the criticality has been confirmed numerically in MatCont \cite{MATCONT}.
Importantly, the singular AH bifurcation resides on the middle branch (close to the lower fold) of the S-shaped equilibrium branch.

The details of the canard explosion in system \eqref{eq:general_form} can differ quite significantly, depending on the number of equilibria.
In the case that there are three equilibria, there is an \textit{incomplete canard explosion} in which canard cycles may terminate prematurely in a homoclinic bifurcation. Figure~\ref{fig:can_explosion} indicates the small zone where such an incomplete canard explosion happens. In this zone, the large relaxation oscillations turn into large canard cycles which terminate in a  large amplitude homoclinic bifurcation. Similarly, the small unstable oscillation cycles born out of the singular AH bifurcation turn into small canard cycles that terminate in a small amplitude homoclinic bifurcation.

\begin{rem}
A rigorous treatment of this incomplete canard explosion goes beyond the scope of this article and is deferred for future work. We refer to, e.g., \cite{Wechselberger2015}, which studies that phenomenon in the context of (neural) excitation.
\end{rem}

\subsection{\SJ{Regular Andronov-Hopf bifurcation under $\tau_{max}$ variation}}

From a physiological point of view, $\tau_{max}$ is an important parameter, as it is one of the major determinants of oscillation period. For example, although we do not have a detailed understanding of oscillation frequency as a function of $\tau_{max}$, we do know that an increase in $\tau_{max}$ leads to a decrease in oscillation frequency \cite{sneyd2017dynamical}. \SJ{Numerical studies of the corresponding open-cell model \eqref{eq:Full_model_open_ct} in \cite{sneyd2017dynamical} show the termination of relaxation oscillations
in a supercritical AH bifurcation.
These findings have been reproduced for system \eqref{eq:general_form} in Figure \ref{fig:period}, which shows the period of oscillations in system \eqref{eq:general_form} as a function of the original dimensionless quantity $\tilde \tau_{max}$ from Section \ref{sub:non_dimensionalisation}.}

\begin{figure}[t!]
	\centering

		\begin{subfigure}[b]{0.49\linewidth}
	\includegraphics[width=\linewidth]{Images/period_new}
		\end{subfigure}
	\begin{subfigure}[b]{0.49\linewidth}
		\includegraphics[width=\linewidth]{Images/CaCycsTau_1}

	\end{subfigure}
	\caption{Bifurcation diagrams for \eqref{eq:general_form} as a function of $\tilde \tau_{max}$, calculated with MatCont \cite{MATCONT}  for parameter values as in Table \ref{tab:params_post_scaling}: (Left) Period $\mathcal T$ as a function of $\tilde \tau_{max}$; (Right) Maximum and minimum values of $c$ as a function of $\tilde \tau_{max}$. Note that we plot with respect to $c = \sqrt \epsilon C$ here.
The \SJ{green} disk indicates the onset of oscillations at a supercritical AH bifurcation at $\tilde \tau_{max} = \tilde \tau_{max,ah} \approx 112.5$ or equivalently, $\nu_{max,ah} \approx 5.83\times10^{-2}$. }

	\label{fig:period}
\end{figure}

\begin{comment}
\SJ{Since we } we return to the original dimensionless value for $\tilde \tau_{max}$ from Section \ref{sub:non_dimensionalisation}. It follows from \eqref{eq:h_scale_tau_max} and \eqref{eq:common_scaling} that Theorem \ref{thm:main} and Proposition \ref{prop:period} describe the relaxation oscillations in system \eqref{eq:main1} for $\tilde \tau_{max} = O(\epsilon^{-2})$.
We are motivated by numerical findings on the corresponding open-cell model \eqref{eq:Full_model_open_ct} in \cite[Figure S3]{sneyd2017dynamical},
where a supercritical Andronov-Hopf bifurcation is identified for a (dimensional) value of $\tau_{max}$ of numerical order $\approx10$s, to consider the effect of varying $\tilde \tau_{max}$ in system \eqref{eq:main2}. These findings have been reproduced in Figure \ref{fig:period}, which shows the period of oscillations in system \eqref{eq:main2} as a function of (dimensionless) $\tilde \tau_{max}$.
\end{comment}

\SJ{In order to identify the AH bifurcation analytically, we consider system \eqref{eq:main_R2} except with
\begin{equation}
\label{eq:scale_tau_new}
\tau_{max} = \delta \nu_{max} = \epsilon^{1/2} \nu_{max} ,
\end{equation}
}and vary $\nu_{max} \in [0,\xi]$, where $\xi>0$ is fixed, as a bifurcation parameter.
The alternative scaling \eqref{eq:scale_tau_new} amounts to restriction to the scaling regime $\tilde \tau_{max} = O(\epsilon^{-3/2})$, instead of the scaling regime $\tilde \tau_{max} = O(\epsilon^{-2})$ corresponding to relaxation oscillations. In this case, the dynamics within regime (R2) are governed by
\begin{equation}\label{sys:main_iii_rescaled}
\begin{split}
h' &= - \frac{1}{\tau_h(C)} \left(h - h_\infty(C)\right) , \\
C' &= \tilde f_0(h,C) + \delta \tilde f_{rem}(h, C, \delta) .
\end{split}
\end{equation}
System \eqref{sys:main_iii_rescaled} is similar to system \eqref{eq:main_R2}, except that there is no small parameter factoring the equation for $h'$. Setting $\delta = 0$ in \eqref{sys:main_iii_rescaled} yields
\begin{equation}\label{sys:main_iii_rescaled_lo}
\begin{split}
h' &= - \frac{1}{\tau_h(C)} \left(h - h_\infty(C)\right) , \\
C' &= \tilde f_0(h,C) ,
\end{split}
\end{equation}
which has equilibria for $(h_\ast, C_\ast)$ satisfying $h_\ast = h_\infty(C_\ast) = \zeta(C_\ast)$,  where $\zeta(C)$ is given by \eqref{eq:S_R2}. In fact, it follows from our observations in Section \ref{ssec:singular_limit_analysis_R2} that system \eqref{sys:main_iii_rescaled_lo} can have one, two, or three equilibria, depending on the location in $(p,c_t)-$parameter space. We note that system \eqref{sys:main_iii_rescaled} is a \emph{regular perturbation problem}, with leading order dynamics determined by the limiting system \eqref{sys:main_iii_rescaled_lo}.

\SJ{The observed AH bifurcation is described in the following result. Note that neither $h_\infty(C)$ nor $\zeta (C)$ depend on $\tilde \tau_{max}$, so the number and location of equilibria is also independent of $\tilde \tau_{max}$.}

\begin{comment}
	\centering
\begin{subfigure}[b]{0.48\linewidth}
	\includegraphics[width=\linewidth]{Images/timeseries_nondim.png}
\end{subfigure}
\begin{subfigure}[b]{0.48\linewidth}
	\includegraphics[width=\linewidth]{Images/timeseries_nondim_zoom.png}
\end{subfigure}
\end{comment}

\begin{theorem}
	\label{thm:hopf}
	Consider system \eqref{sys:main_iii_rescaled_lo}, with a unique equilibrium $q = (h_\ast, C_\ast)$ such that
	\begin{equation}
	\label{eq:det_cond}
	C_\ast > C_F, \qquad  \big| D_C h_\infty(C_\ast) \big| > \frac{D_C \tilde f_0(q_\delta)}{D_h \tilde f_0(q_\delta)} ,
	\end{equation}
	where the coordinate $C_\ast$ can be written explicitly in terms of the model parameters. Then there exists $\delta_0 >0$ such that for each $\delta \in [0, \delta_0)$, system \eqref{sys:main_iii_rescaled} has a unique equilibrium $q_\delta$ such that $q_\delta \to q$ as $\delta \to 0$.  Furthermore, $q_{\delta}$ undergoes an AH bifurcation \SJ{for}
	\begin{equation}
	\label{eq:hopf_value}
	\nu_{max,ah}(\delta) = \frac{\left(K_h^4 + C_\ast^4 \right) C_\ast }{2 A_{\text{SERCA}} K_h^4(C_\ast^2 - 2 K \gamma^2 c_t^2)} + O(\delta) ,

	\end{equation}
	or, equivalently,
	\[
	\tilde \tau_{max,ah}(\epsilon) = \epsilon^{-3/2} \nu_{max,ah}(\sqrt{\epsilon}) = \frac{\left(K_h^4 + C_\ast^4 \right) C_\ast }{2 A_{\text{SERCA}} K_h^4(C_\ast^2 - 2 K \gamma^2 c_t^2)} \epsilon^{-3/2} + O(\epsilon^{-1}) .
	\]
\end{theorem}

\begin{proof}
	The proof is based on a fixed point analysis, and deferred to Appendix \ref{app:hopf_proof}. The conditions in \eqref{eq:det_cond} ensure satisfaction of the positive determinant condition associated with the AH point. The fact that $C_\ast$ can be written explicitly in terms of the model parameters follows from the fact that $C_\ast^2$ is a solution to the cubic polynomial \eqref{eqs:reduced_equilibria_R2b}.
\end{proof}

Numerical observations indicate a supercritical bifurcation and hence, the existence of nearby stable oscillations; see Remark \ref{rem:hopf_numerics} below.
The qualitative shape of the profile in Figure \ref{fig:period} can also be explained by our findings. \SJ{Let $\mathcal T$ denote the oscillation period on the fastest timescale $t$, as in Proposition \ref{prop:period}. Then}

\smallskip

\begin{enumerate}
	\item[(i)] \SJ{The fixed point analysis used to prove Theorem \ref{thm:hopf} in Appendix \ref{app:hopf_proof}, in particular the form of the determinant \eqref{eq:det},} \SJ{implies an oscillation period \SJ{$\mathcal T \propto \epsilon^{-3/2} \sqrt \nu_{max} = \epsilon^{-1/2} \sqrt{\tilde \tau_{max}}$}
		for parameter values $\tilde \tau_{max} = O(\epsilon^{-3/2})$ close to \SJ{the AH bifurcation.}
	\item[(ii)] Proposition \ref{prop:period} implies and oscillation period $\mathcal T = O(\epsilon^{-2})$ and leading order dependence $\mathcal T \propto \tilde \tau_{max}$ for parameter values $\tilde \tau_{max} = O(\epsilon^{-2})$, in the relaxation oscillatory regime.}
\end{enumerate}

\begin{remark}
	\label{rem:hopf_numerics}
	The conditions for the applicability of Theorem \ref{thm:hopf} are satisfied for system \eqref{sys:main_iii_rescaled} with parameter values taken from Table \ref{tab:params_post_scaling} and $\delta = \sqrt{2.5 \times 10^{-3}} = 5 \times 10^{-2}$. The conditions in \eqref{eq:det_cond} are satisfied since
	\[
	\big| D_C h_\infty(C_\ast) \big| \approx 0.20 > \frac{D_C \tilde f_0(q_\delta)}{D_h \tilde f_0(q_\delta)} \approx 8.3 \times 10^{-2} ,
	\]
	and the AH bifurcation occurs at $q_\delta = (h_\ast, C_\ast) \approx (0.0792, 1.48 )$, so $C_\ast > C_F \approx 0.68$. We obtain a numerical estimate of $\nu_{max,ah} \approx 5.83\times 10^{-2}$,
	in close agreement with the leading order estimate \eqref{eq:hopf_value} as $\delta \to 0$:
	\[
	\nu_{max,ah}(0) = \frac{\left(K_h^4 + C_\ast^4 \right) C_\ast }{2 A_{\text{SERCA}} K_h^4 (C_\ast^2 - 2 K \gamma^2 c_t^2)}
	\approx 5.42\times 10^{-2} ,
	\]
	where we have used the
	numerical value for $C_\ast$. A first Lyapunov coefficient $l_1 \approx -3.13 < 0$ was calculated numerically in the software package MatCont \cite{MATCONT}, indicating a supercritical bifurcation.
\end{remark}

\begin{comment}
\begin{remark}
	\label{rem:connection} {\color{red} This comment seems out of character with the rest of the paper. Could it go in the Appendices?} \SJ{[I'm not sure where in the appendix it could go (none of the blow-up analysis relevant for the onset of oscillations is presented). Still it seems to me like an important observation... I'm not sure... I will leave it to NZ to decide.]}
One can view system \eqref{sys:main_iii_rescaled} as being obtained in the `family rescaling chart' of a cylindrical blow-up along the degenerate line $S_h$, which also appears in system \eqref{sys:main_iii_rescaled} when considered in regime (R1). In this context, the small amplitude oscillations born in a \textit{singular Andronov-Hopf bifurcation} giving rise to oscillations on the blow-up cylinder for $\delta = 0$. A rigorous treatment of the `connection problem' for the limit cycles would require an understanding of the transition of these oscillations to the relaxation type oscillations perturbing from $\Gamma$ (see Figure \ref{fig:cyl_bu}), which contain distinct cycle segments off of the cylinder. This problem is \wm{non-uniform} insofar as it requires (i) detailed information regarding the growth of the cycles on the cylinder, and (ii) additional parameter space blow-ups to deal with variation of the parameter $\tilde \tau_{max}$ over varying orders of magnitude.
\end{remark}
\end{comment}

\begin{comment}
\subsection{Larger values $\tilde \tau_{max} \gg \epsilon^{-2}$}
\label{sub:the_large_tau_max_case}

An analysis of the behaviour of system \eqref{eq:main1} for which $\epsilon_{\tau_{max}}$ scales in a precise way with $\epsilon$ when $\tilde \tau_{max} \gg \epsilon^{-2}$ is beyond the scope of this article, however some insight can be gained by considering $\epsilon_{\tau_{max}}$ as an \textit{independent} small parameter and considering the system

\begin{align}
\begin{split}
\frac{\text{d}h}{\text{d}t} &=  \frac{\epsilon_{\tau_{max}}}{\epsilon^2} g(h,c), \\
\frac{\text{d}\tilde c}{\text{d}t} &= f(h,c),
\end{split}
\qquad 0 < \frac{\epsilon_{\tau_{max}}}{\epsilon^2} \ll 1 ,
\label{eq:tau_max_large}
\end{align}
with
\begin{equation}
\begin{split}
g(h,c) &= \left(\frac{1}{\hat K_\tau^4 \hat \tau_{max}}\right) \left(h_{\infty}(\tilde{c},\epsilon) - h\right) \left(\tilde{c}^4 + \epsilon^2 \hat K_\tau^4\right), \\
f(h,c) &= J_{\text{IPR}}(h,\tilde{c},\epsilon) - \epsilon \mathfrak{J}^+_{\text{SERCA}}(\tilde{c},\epsilon) + \epsilon^2 \mathfrak{J}_{\text{SERCA}}(\tilde{c},\epsilon).
\end{split}
\end{equation}

\begin{figure}[t!]
	\centering
	\begin{subfigure}[b]{0.75\linewidth}
		\includegraphics[width=\linewidth]{Images/phase_tau_large.png}
	\end{subfigure}
	\caption{Phase portrait an oscillation in system \eqref{eq:tau_max_large}, for parameter values as in Table \ref{tab:params_post_scaling}, except for $\tilde \tau_{max} = 2.1~\times~10^{7}$ and $\epsilon = 2.5~\times~10^{-3}$. In the left panel the orange curve represents the critical manifold $S_\tau$.}
	\label{fig:tau_max_large}
\end{figure}

Considering $\epsilon > 0$ as fixed, system \eqref{eq:tau_max_large} is in the standard form for slow-fast systems \eqref{eq:stnd_form}, with a critical manifold defined by $S_\tau := \{(h,c) | f(h,c)=0\}$. This is shown alongside the numerically observed relaxation oscillation in Figure \ref{fig:tau_max_large}.

In the phase space, the periodic orbit looks like a relaxation oscillation from a GSPT problem in standard form \cite{desroches2012mixed}. The periodic orbit contains four strokes, two straight fast strokes and two slow ones that follow an S-shaped critical manifold, see also Figure \ref{fig:phase_tau_max_large}, and jumps near the fold points. However, the time series, see Figure \ref{fig:time_tau_max_large} corresponding to the same oscillation looks unlike a relaxation oscillations arising from standard GSPT problems. The time series shows that the solution spends significantly more time near the branch with lower values of $c$ than near the upper branch. This difference of time spent on the slow manifold is due to the variable $h$ changing its speed in the phase space because of $\tau_{h}$. The system then seems to have three timescales: a fast time scale outside the critical manifold $S_\tau$, an intermediate timescale on the high c-valued stable branch of $S_\tau$, and a slow timescale on the lower branch of $S_\tau$. Nevertheless, only two timescales are apparent from the formulation in equation \eqref{eq:tau_max_large}. So, studying this system using the standard form of GSPT seems ill-suited because of the capacity of standard GSPT to see only two timescales. On the other hand, using the non-standard methods presented in this article would require a new ad hoc
non-dimensionalisation, which we do not consider further.
\end{comment}

\section{Discussion and Conclusion}
\label{sec:discussion_conclusion}

\wm{Non-standard} timescale separation is ubiquitous in biological and physiological phenomena, but literature on the mathematical modelling of such systems is relatively sparse. This is perhaps due to the significant mathematical complications inherent in such systems, and the perceived absence of a well adapted mathematical framework for their analysis. At present, however, there is a growing awareness that classical methods based on GSPT can be successfully adapted and applied in such contexts. This has been demonstrated in the context of chemical and biological models in the work of Kosiuk and Szmolyan \cite{Gucwa2009,kosiuk2011scaling,kosiuk2016geometric}. Formal developments in subsequent work by Wechselberger \cite{wechselberger2018geometric}, Goeke \& Walcher \cite{goeke2014constructive}, Lizarraga, Marangell \& Wechselberger \cite{Lizarraga2020,Lizarraga2020c} and Kruff \& Walcher \cite{Kruff2019} have aided in providing a systematic framework for \wm{a coordinate-independent} analysis of such problems via GSPT.

In this work, we studied a closed-cell model for intracellular calcium dynamics which is obtained via a reduction from a canonical open-cell model developed in \cite{sneyd2017dynamical}. By adopting a process-oriented philosophy and methodology, we identified a total of five small parameters in the model. The identified small parameters were either (i) uniform perturbations setting the scale of a particular process (e.g.~$V_s$), or (ii) non-uniform perturbations deriving from the switching behaviour associated with particular biological processes.
Small parameters of the latter kind
are responsible for \wm{non-standard} timescale separations, since the corresponding singular limit is non-uniform in phase-space. In our particular model, this led to distinct scaling regimes (R1) and (R2), with distinct timescale separations. In the context of biological switching in general, an increasing number of studies show that distinct scaling regimes of this kind are to be expected \cite{Gucwa2009,kosiuk2016geometric,Kristiansen2019d}. These references, along with the current manuscript, also demonstrate the suitability of blow-up methods for combining dynamical information obtained in distinct scaling regimes. Following the identification of small parameters, a common scaling \wm{with a single small parameter $\epsilon \ll 1$} was proposed in \eqref{eq:common_scaling}. This led to the singular perturbation problem \eqref{eq:main2}, which is given in a general form suited to the GSPT framework in \cite{wechselberger2018geometric}.
We proved existence
of the observed three-timescale relaxation cycles in Theorem \ref{thm:main}.
\wm{We then identified a singular Andronov-Hopf bifurcation together with a complete/incomplete canard explosion as an onset mechanism for oscillations under the variation of the total calcium concentration or IP$_3$ concentration.} We also described a relationship between the system parameter $\tilde \tau_{max}$, which sets the speed of evolution of $h$, and the period of oscillations. In particular, Proposition \ref{prop:period} describes a linear relationship for values $\tilde \tau_{max} = O(\epsilon^{-2})$, and Theorem \ref{thm:hopf} shows that a supercritical Andronov-Hopf provides \SJ{one possible} mechanism for the onset of oscillations with $\tilde \tau_{max} = O(\epsilon^{-3/2})$.

There are many directions in which one can take this work. For example, {it could be interesting to give a detailed description of the canard and incomplete canard explosions described in Section \ref{sub:can_explosion}.}
\SJ{The} current manuscript also paves the way for rigorous analytical treatment of the higher dimensional open-cell model from which it is derived. In fact, one could view system \eqref{eq:Full_model} as a layer problem for the open-cell but `almost-closed' model \eqref{eq:Full_model_open_ct} with $|J_{pm}| \ll 1$. In this context, a \wm{one}-parameter bifurcation analysis of system \eqref{eq:main2} in $c_t$ amounts to an analysis of the layer problem to system \eqref{eq:Full_model_open_ct}; a three-dimensional system with (at least) four timescales, due to an additional slow evolution in $c_t$. The case $|J_{pm}| = O(1)$ corresponding to `more open' cells is also of physiological significance. In both cases, this work provides a solid foundation upon which these significant and highly non-trivial analyses become feasible.

\begin{comment}

In this article, we have studied a model of calcium dynamics that originated from a model in \cite{sneyd2017dynamical} but which we reduced to a closed-cell model of two dimensions. This model is of great interest physiologically as it has been designed to represent the calcium dynamics in a large class of cells. Thus understanding the underlying dynamics of this particular model is of broad interest for calcium modelling. As this model exhibits periodic orbits with multiple timescales, we wished to study it using GSPT.

Because of the underlying physiology, the model could not be put into the standard form introduced by Fenichel \cite{fenichel1979geometric}. In previous work by Harvey \cite{harvey2011analysing}, calcium models were studied in the standard form of GSPT by considering some variables fast in a first analysis then slow in a second one; then both analysis would be combined. However, based on the work of Szmolyan, Kosiuk and Weschelberger \cite{kosiuk2011scaling,kosiuk2016geometric,wechselberger2018geometric}, we were able to rewrite the model in a form that fit a non-standard framework of GSPT. This framework allows capturing in one formulation the slow-fast structure of the model. The advantages of this approach are that having variables separated by the speed of evolution is not necessary, there is no limitation on the number of explicit timescales in the model, and the timescales of the model can be related to biological processes. However, not having the model in standard form comes with drawbacks as we cannot consider the layer problem as a bifurcation problem and the method to identify the timescales is not straightforward.

To put the model under the non-standard form of GSPT, we used an approach to identify small parameters based on the biological processes of the model. This approach is model-based, and thus small parameters and timescales identified give a better insight about the underlying physiology. We think this method can be relevant in a large variety of problem that can be divided into processes and not only in physiological models like calcium modelling. For example, the autocatalator previously studied by Gucwa and Szmolyan \cite{gucwa2009geometric}, is a problem in the non-standard form derived from the laws of mass actions. The three timescales of the problem are based on three different processes: the fastest timescale corresponds to the autocatalatic reaction, the intermediate timescale corresponds to the transformation of the intermediate into the product and the slow timescales corresponds to the transformation of the precursor into the first intermediate.

Once we identified all the small parameters, we linked them to one single small parameter to facilitate the analysis. Studying the model by considering each of the small parameters independently is well beyond the scope of this article. A direct consequence of assuming that all the small parameters scaled with a single small parameter is that the GSPT analysis carried in Section \ref{sec:effect_of_tau_} when $\tau_{max}$ is small or large is entirely different from the analysis carried in Section \ref{sec:slow-fast_analysis_and_statement_of_the_main_result} because the small parameters are not linked together in the same way.

We then used GSPT to understand the behaviour of the system on the fast and intermediate timescales. Next, we rescaled the system to explain the dynamics near the non-hyperbolic part of the critical manifold. While using those two approaches gives a good understanding of the system in the different regimes, a blow-up analysis is needed to stitch together the dynamics in the two regimes (see Appendix \ref{sec:blow_up}).

We also proved that in the parameter regime in which we were studying the model, the period of the oscillations grows linearly with the parameter $\tau_{max}$; and in the case of small values of $\tau_{max}$, the growth of the period is proportional to the square root of $\tau_{max}$. Those proofs are consistent with \cite{sneyd2017dynamical}.

The next natural step is to look at the transitions and bifurcations of the system. In Section \ref{sec:effect_of_tau_}, we showed the existence of different regimes for periodic orbits. However, for $\tau_{max}$ small, the periodic orbit exists only in the regime $c = O(\epsilon)$, whereas the case detailed in Section \ref{sec:slow-fast_analysis_and_statement_of_the_main_result} has parts of the periodic orbit in both $c = O(1)$ and $c = O(\epsilon)$. The transition between those two types of periodic orbits is not yet understood.

In calcium modelling, the ${\rm IP_3}${} concentration ($p$) and the total calcium of the cell ($c_t$) are often used as bifurcation parameters. An analysis of the bifurcation of the system using those parameters might reveal some unusual dynamics. For example, in the case of the rescaled regime \eqref{eq:main_R2}, one would expect canards when the equilibrium of the system, $C_*$, crosses the fold point $F$. However, contrary to the standard case, nothing appears to be limiting the growth of the canard.

Another possible direction for future work would be to look at models of calcium dynamics of higher dimension, like open-cells models or more complex models including dynamics of ${\rm IP_3}${}. For example, studying the full model of \cite{sneyd2017dynamical} in term of non-standard GSPT. Such analysis might be challenging, as it would add other small parameters that would then need to be linked meaningfully to the exisiting ones.

\end{comment}

\bibliographystyle{siamplain}