\frac{\alpha_{(m-1)/2}}{1-x_{(m-1)/2}} + y_{(m-1)/2} + \frac{g^e}{2}(x_{[(m-3)/2] \bmod\zeta} + x_{[(m+1)/2] \bmod\zeta} - 2x_{(m-1)/2}), & x_{(m-1)/2} \leq 0 \vspace{0.5cm} \\
\alpha_{(m-1)/2} + y_{(m-1)/2} + \frac{g^e}{2}(x_{[(m-3)/2] \bmod\zeta} + x_{[(m+1)/2] \bmod\zeta} - 2x_{(m-1)/2}), & 0 < x_{(m-1)/2} < \alpha_{(m-1)/2} \\
& + y_{(m-1)/2} + \mathfrak{C}_{(m-1)/2} \vspace{0.5cm} \\
-1, & x_{(m-1)/2} \geq \alpha_{(m-1)/2} \\
& + y_{(m-1)/2} + \mathfrak{C}_{(m-1)/2} \\
\end{cases}
\end{split}
\label{eq:piecewise-ring-function}
\end{equation}

Because there are three distinct functions within this piecewise form of $F\e{m}$, we will consider each of the pieces in turn.

First, let us investigate the case where $x_{(m-1)/2} \leq 0$. Taking a look at this piece of the function, we see that the only variables present are $y_{(m-1)/2}$, $x_{(m-1)/2}$, $x_{[(m-3)/2]\bmod\zeta}$, and $x_{[(m+1)/2]\bmod\zeta}$. Then, the only non-zero derivatives of this piece of $F\e{m}$ are

\begin{align}
\frac{\partial F\e{m}}{\partial y_{(m-1)/2}} &= 1 \label{eq:piece-1-first-deriv} \\
\frac{\partial F\e{m}}{\partial x_{(m-1)/2}} &= \frac{\alpha_{(m-1)/2}}{(1-x_{(m-1)/2})^2} - g^e \label{eq:piece-1-sec-deriv} \\
\frac{\partial F\e{m}}{\partial x_{[(m-3)/2]\bmod\zeta}} &= \frac{g^e}{2} \label{eq:piece-1-third-deriv} \\
\frac{\partial F\e{m}}{\partial x_{[(m+1)/2]\bmod\zeta}} &= \frac{g^e}{2} \label{eq:piece-1-fourth-deriv}
\end{align}

We will now once again systematically determine the values of $j$ that give us these non-zero derivatives. For this case where $m$ is odd, checking the state vector in Equation <ref> indicates that $F\e{m}$ is the iteration function for $x_{(m-1)/2}$:

\begin{equation}
x_{(m-1)/2,\,k+1} = F\e{m}(x_{0,\,k},\,y_{0,\,k},\,x_{1,\,k},\,y_{1,\,k},\,\hdots\,,\,x_{\zeta-1,\,k},\,y_{\zeta-1,\,k})
\end{equation}

meaning $X\e{m} = x_{(m-1)/2}$. For the first non-zero derivative in Equation <ref>, we want to find $j$ such that $X\e{j} = y_{(m-1)/2}$, which is one dimension down from its associated fast variable $x_{(m-1)/2}$. In other words, $X\e{m+1} = y_{(m-1)/2}$, so $\frac{\partial F\e{m}}{\partial X\e{j}}=1$ if $j=m+1$ (given $m$ is odd). For the second non-zero derivative in Equation <ref>, we are differentiating with respect to $X\e{j} = x_{(m-1)/2}$, meaning $\frac{\partial F\e{m}}{\partial X\e{j}}=\alpha_{(m-1)/2}(1-x_{(m-1)/2})^{-2} - g^e$ if $j=m$. The slightly more challenging $j$ to find comes with the third non-zero derivative in Equation <ref>, where $X\e{j} = x_{[(m-3)/2]\bmod\zeta}$, or the “clockwise” neighbor. Similar to before, let us first consider the special case where $X\e{m} = x_{(m-1)/2}$ is the fast variable of the neuron with index $i=0$, which happens when $m=1$. In this case, the “clockwise” neighbor to $\mathbf{x}_0$ is $\mathbf{x}_{\zeta-1}$, which is associated with the fast variable $X\e{2\zeta-1} = x_{[(m-3)/2]\bmod\zeta} = x_{\zeta-1}$ for $m=1$. For the “normal” cases, going from $X\e{m} = x_{(m-1)/2}$ to its clockwise neighbor $X\e{j} = x_{(m-3)/2}$ requires moving up two dimensions: $X\e{m-2} = x_{(m-3)/2}$. Therefore, $\frac{\partial F\e{m}}{\partial X\e{j}}=g^e/2$ if $j=2\zeta-1$ and $m=1$, or $j=m-2$ and $m\neq 1$. The fourth non-zero derivative, concerning the “counterclockwise” neighbor $X\e{j} = x_{[(m+1)/2]\bmod\zeta}$, follows naturally from this, yielding that $\frac{\partial F\e{m}}{\partial X\e{j}}=g^e/2$ if $j=1$ and $m=2\zeta-1$, or $j=m+2$ and $m\neq 2\zeta-1$. We can now summarize the results from this paragraph into one piecewise function for the derivative $\frac{\partial F\e{m}}{\partial X\e{j}}$ that holds for $x_{(m-1)/2}\leq 0$ when $m\bmod 2=1$:

\begin{equation}
\frac{\partial F\e{m}}{\partial X\e{j}} =
\begin{cases}
1, & \hspace{0.325cm}\text{if }j=m+1 \\
\frac{\alpha_{(m-1)/2}}{(1-x_{(m-1)/2})^2} - g^e, & \hspace{0.325cm}\text{if }j=m\\
g^e/2, & \begin{cases}
\text{if }j=m-2\text{ and }m\neq 1 \\
\text{or }j=2\zeta-1\text{ and }m=1 \\
\text{or }j=m+2\text{ and }m\neq 2\zeta - 1 \\
\text{or }j=1\text{ and }m=2\zeta - 1
\end{cases} \\
0, & \hspace{0.325cm}\text{otherwise}
\end{cases}
\label{eq:ring-jacobian-part-2}
\end{equation}

Luckily, this is the peak of the difficulty mountain for this derivation of the Jacobian of our ring system, and it only goes downhill from here in complexity.

Continuing to dismantle the complex piecewise function in Equation <ref>, let us next consider the case where $0 < x_{(m-1)/2} < \alpha_{(m-1)/2} + y_{(m-1)/2} + \mathfrak{C}_{(m-1)/2}$. Taking a look at this piece of Equation <ref>, we can see that it contains the same variables as the piece where $x_{(m-1)/2}\leq 0$, namely, $y_{(m-1)/2}$, $x_{(m-1)/2}$, $x_{[(m-3)/2]\bmod\zeta}$, and $x_{[(m+1)/2]\bmod\zeta}$. The function itself, however, is evidently different, so the non-zero derivatives are also different:

\begin{align}
\frac{\partial F\e{m}}{\partial y_{(m-1)/2}} &= 1 \label{eq:piece-2-first-deriv} \\
\frac{\partial F\e{m}}{\partial x_{(m-1)/2}} &= - g^e \label{eq:piece-2-sec-deriv}\\
\frac{\partial F\e{m}}{\partial x_{[(m-3)/2]\bmod\zeta}} &= \frac{g^e}{2} \label{eq:piece-2-third-deriv}\\
\frac{\partial F\e{m}}{\partial x_{[(m+1)/2]\bmod\zeta}} &= \frac{g^e}{2} \label{eq:piece-2-fourth-deriv}
\end{align}

We might think that we now have to find the values of $j$ that give us these non-zero partial derivatives once again, but thankfully, we've lucked out in this case because we are dealing with the same state $X\e{m}=x_{(m-1)/2}$ and the same variables as the first piece of $F\e{m}$: $y_{(m-1)/2}$, $x_{(m-1)/2}$, $x_{[(m-3)/2]\bmod\zeta}$, and $x_{[(m+1)/2]\bmod\zeta}$. All of the relationships between these variables in state space are the same as before because the piece of $F\e{m}$ we are working with has no influence on the variables we are connecting. Therefore, all we have to do is plug the non-zero derivatives from Equations <ref>, <ref>, <ref>, and <ref> into the appropriate locations of the general form of the derivative in displayed in Equation <ref>. Making these substitutions yields the piecewise function for the derivative $\frac{\partial F\e{m}}{\partial X\e{j}}$ that holds for $0 < x_{(m-1)/2} < \alpha_{(m-1)/2} + y_{(m-1)/2} + \mathfrak{C}_{(m-1)/2}$ when $m\bmod 2=1$:

\begin{equation}
\frac{\partial F\e{m}}{\partial X\e{j}} =
\begin{cases}
1, & \hspace{0.325cm}\text{if }j=m+1 \\
- g^e, & \hspace{0.325cm}\text{if }j=m\\
g^e/2, & \begin{cases}
\text{if }j=m-2\text{ and }m\neq 1 \\
\text{or }j=2\zeta-1\text{ and }m=1 \\
\text{or }j=m+2\text{ and }m\neq 2\zeta - 1 \\
\text{or }j=1\text{ and }m=2\zeta - 1
\end{cases} \\
0, & \hspace{0.325cm}\text{otherwise}
\end{cases}
\label{eq:ring-jacobian-part-3}
\end{equation}

Finally, for the last part of the piecewise function $F\e{m}$ in Equation <ref>, where $x_{(m-1)/2} \geq \alpha_{(m-1)/2} + y_{(m-1)/2} + \mathfrak{C}_{(m-1)/2}$, we have been blessed with an extremely easy function to differentiate: $F\e{m}=-1$. Obviously, the derivative of $-1$ with respect to any variable is 0. Therefore, for the case where $x_{(m-1)/2} \geq \alpha_{(m-1)/2} + y_{(m-1)/2} + \mathfrak{C}_{(m-1)/2}$ for $m\bmod 2=1$, we can write the derivative $\frac{\partial F\e{m}}{\partial X\e{j}}$ very simply as

\begin{equation}
\frac{\partial F\e{m}}{\partial X\e{j}} = 0
\label{eq:ring-jacobian-part-4}
\end{equation}

which concludes all the possible cases of the derivative $\frac{\partial F\e{m}}{\partial X\e{j}}$.

So now, for the grand finale, we combine all of these calculations into one partial derivative to rule them all: $\frac{\partial F\e{m}}{\partial X\e{j}} = J_{mj}(\mathbf{X})$. This derivative represents the $m$th row and the $j$th column of the Jacobian of any state $\mathbf{X}$ of this symmetrically coupled ring lattice of Rulkov 1 neurons for any $m$ or $j$ between $1$ and $2\zeta$. To do this, we will combine all of the cases for the different values of $m$ and $j$ we have discussed: for even values of $m$ and all values of $j$, Equation <ref>; for odd values of $m$, values $x_{(m-1)/2}\leq 0$, and all values of $j$, Equation <ref>; for odd values of $m$, values $0 < x_{(m-1)/2} < \alpha_{(m-1)/2} + y_{(m-1)/2} + \mathfrak{C}_{(m-1)/2}$, and all values of $j$, Equation <ref>; and for odd values of $m$, values $x_{(m-1)/2} \geq \alpha_{(m-1)/2} + y_{(m-1)/2} + \mathfrak{C}_{(m-1)/2}$, and all values of $j$, Equation <ref>. Putting all of this together yields the generalized entry of the Jacobian $J_{mj}(\mathbf{X})$ shown on the next page. It is worth noting that when using this equation in practice, it should be read from right to left. And now, without further ado, here it is:

\begin{equation}
J_{mj}(\mathbf{X}) =
\begin{cases}
\begin{cases}

\begin{cases}
1, & \hspace{0.325cm}\text{if }j=m+1, \\
\frac{\alpha_{(m-1)/2}}{(1-x_{(m-1)/2})^2} - g^e, & \hspace{0.325cm}\text{if }j=m,\\
g^e/2, & \begin{cases}
\text{if }j=m-2,\\
\phantom{if }\text{ and }m\neq 1, \\
\text{or }j=2\zeta-1,\\
\phantom{if }\text{ and }m=1, \\
\text{or }j=m+2,\\
\phantom{if }\text{ and }m\neq 2\zeta - 1, \\
\text{or }j=1,\\
\phantom{if }\text{ and }m=2\zeta - 1,
\end{cases} \\
0, & \hspace{0.325cm}\text{otherwise,}
\end{cases} \text{for } x_{(m-1)/2} \leq 0, \\
\\
\\
\begin{cases}
1, & \hspace{2.225cm}\text{if }j=m+1, \\
- g^e, & \hspace{2.225cm}\text{if }j=m,\\
g^e/2, & \hspace{1.9cm}\begin{cases}
\text{if }j=m-2,\\
\phantom{if }\text{ and }m\neq 1, \\
\text{or }j=2\zeta-1,\\
\phantom{if }\text{ and }m=1, \\
\text{or }j=m+2,\\
\phantom{if }\text{ and }m\neq 2\zeta - 1, \\
\text{or }j=1,\\
\phantom{if }\text{ and }m=2\zeta - 1,
\end{cases} \stackunder{\text{for $0 < x_{(m-1)/2} < \alpha_{(m-1)/2}$}}{\text{$+ y_{(m-1)/2} + \mathfrak{C}_{(m-1)/2},$}} \\
0, & \hspace{2.225cm}\text{otherwise,}
\end{cases}\\
\\
\\
0, \hspace{6.685cm} \text{for} \stackunder{x_{(m-1)/2} \geq \alpha_{(m-1)/2}}{+ y_{(m-1)/2} + \mathfrak{C}_{(m-1)/2},}\\

\end{cases} & \text{when } m\bmod 2=1 \\
\\
\\
\begin{cases}
1, & \hspace{0.67cm}\text{if }j=m \\
-\eta(1+g^e), & \hspace{0.67cm}\text{if }j=m-1 \\
\eta g^e/2, & \hspace{0.345cm}
\begin{cases}
\text{if }j=m-3\\
\phantom{if }\text{ and }m\neq 2 \\
\text{or }j=2\zeta-1\\
\phantom{if }\text{ and }m=2 \\
\text{or }j=m+1\\
\phantom{if }\text{ and }m\neq 2\zeta \\
\text{or }j=1\\
\phantom{if }\text{ and }m=2\zeta
\end{cases} \\
0, & \hspace{0.67cm}\text{otherwise,}
\end{cases} & \text{when } m\bmod 2=0
\end{cases}
\label{eq:THE-jacobian-entry}
\end{equation}

Now that we have completed this Herculean task, there comes the relatively straightforward process of coding Equation <ref> into Python, which we do in Appendix <ref>. Along with a modification of the Lyapunov QR factorization code to account for the potentially much larger Jacobian matrices, this code calculates the Lyapunov spectrum of a ring lattice system orbit $O(\mathbf{X}_0)$ by using Equation <ref> to determine the Jacobian matrices $J(\mathbf{X})$ for all $\mathbf{X}\in O(\mathbf{X}_0)$. Now that we have established all the tools we need, we will apply them to exploring the system displayed in Figure <ref>: a ring of $\zeta=30$ Rulkov 1 neurons. Because each of these neurons has one slow and one fast variable (see Equation <ref>), we are dealing with a 60-dimensional state space, which is far bigger than any space we have dealt with so far in this paper. Within this ring system, we will explore the dynamics that emerge from three specific initial states: one where different neurons have different initial fast variable values $x_{i,\,0}$ but all have the same $y_{i,\,0}$, $\sigma_i$, and $\alpha_i$ values, one where different neurons have different $x_{i,\,0}$ and $\sigma_i$ values but the same $y_{i,\,0}$ and $\alpha_i$ values, and one where different neurons have different $x_{i,\,0}$, $\sigma_i$, and $\alpha_i$ values but the same $y_{i,\,0}$ values. We do not consider the case where different neurons start with different $y_{i,\,0}$ values because different evolutions of the slow variable are accounted for by different values of $\sigma_i$ (since we know from Section <ref> $\sigma_i$ determines the value of $x_i$ that will keep $y_i$ steady).

Before we start, let us quickly take care of some brief notational business. Since, in this section, we deal with many different distinct values of the parameters $\sigma_i$ and $\alpha_i$, let us collect all of the values of $\sigma_i$ and $\alpha_i$ into $\zeta$-dimensional vectors. We will denote the vector of all $\sigma_i$ values as

\begin{equation}
\boldsymbol{\sigma} = \begin{pmatrix}
\sigma_0 \\
\sigma_1 \\
\vdots \\
\sigma_{\zeta-1}
\end{pmatrix}
\end{equation}

and the vector of all $\alpha_i$ values as

\begin{equation}
\boldsymbol{\alpha} = \begin{pmatrix}
\alpha_0 \\
\alpha_1 \\
\vdots \\
\alpha_{\zeta-1}
\end{pmatrix}
\end{equation}

An interesting realization to make is that the parameter space of this ring lattice system, which is the space of all $\sigma_i$ and $\alpha_i$ values (not including $\eta$), is isomorphic to the state space of the system because they both can be represented in $\mathbb{R}^{2\zeta}$. For convenience, we will also introduce the one vector, $\mathbf{1}$, as a vector of all ones

\begin{equation}
\mathbf{1} = \begin{pmatrix}
1 \\
1 \\
\vdots \\
\end{pmatrix} \in \mathbb{R}^n
\end{equation}

for some arbitrary number of dimensions $n$ that will be clear from context.

$g^e=0$, $\lambda_1\approx -0.0938$

$g^e=0.05$, $\lambda_1\approx 0.0491$

$g^e=0.25$, $\lambda_1\approx 0.0595$

$g^e=1$, $\lambda_1\approx 0.1694$

Graphs of $x_{i,\,k}$ for eight neurons in a ring of $\zeta=30$ electrically coupled Rulkov 1 neurons, with $x_{i,\,0}\in (-1,\,1)$, $y_{i,\,0}=-3.25$, $\boldsymbol{\sigma} = -0.5\cdot\mathbf{1}$, and $\boldsymbol{\alpha} = 4.5\cdot\mathbf{1}$, visualized and $\lambda_1$ values calculated using the code in Appendix <ref>

Graph of the maximal Lyapunov exponent $\lambda_1$ against the electrical coupling strength $g^e$ for a ring of $\zeta=30$ electrically coupled Rulkov 1 neurons, with $x_{i,\,0}\in (-1,\,1)$, $y_{i,\,0}=-3.25$, $\boldsymbol{\sigma} = -0.5\cdot\mathbf{1}$, and $\boldsymbol{\alpha} = 4.5\cdot\mathbf{1}$, visualized using the code in Appendix <ref>

Now that all of this is established, let us begin exploring our first system. This system has $\zeta=30$ neurons with parameters $\boldsymbol{\sigma} = -0.5\cdot\mathbf{1}$ and $\boldsymbol{\alpha} = 4.5\cdot\mathbf{1}$. In other words, all the $\sigma_i$ values are $-0.5$ and all the $\alpha_i$ values are $4.5$. Additionally, we set the initial slow variable values for all of the neurons to be $y_{i,\,0} = -3.25$. However, setting the initial fast variable values for all of the neurons to be equal would be pointless because then the neurons would all be the same, resulting in no current flow between them. Instead, we will choose $x_{i,\,0}$ variables randomly from the interval $(-1,\,1)$. For the sake of reproducibility, we will use the initial state

\begin{equation}
\begin{split}
\mathbf{X}_0 &= \langle 0.68921784,\, -3.25,\, -0.94561073,\, -3.25,\, \\
&\mathrel{\phantom{=}} -0.95674631,\, -3.25,\,  0.91870134,\, -3.25,\, \\
&\mathrel{\phantom{=}} -0.32012381,\, -3.25,\, -0.23746836,\, -3.25,\, \\
&\mathrel{\phantom{=}} -0.43906743,\, -3.25,\, -0.48671017,\, -3.25,\, \\
&\mathrel{\phantom{=}} -0.37578533,\, -3.25,\, -0.00613823,\, -3.25,\, \\
&\mathrel{\phantom{=}} 0.25990663,\, -3.25,\, -0.54103868,\, -3.25,\, \\
&\mathrel{\phantom{=}} 0.12110471,\, -3.25,\,  0.71202085,\, -3.25,\, \\
&\mathrel{\phantom{=}} 0.689336,\, -3.25,\,   -0.03260047,\, -3.25,\, \\
&\mathrel{\phantom{=}} -0.90907325,\, -3.25,\,  0.93270227,\, -3.25,\, \\
&\mathrel{\phantom{=}} 0.51953315,\, -3.25,\, -0.46783677,\, -3.25,\, \\
&\mathrel{\phantom{=}} -0.96738424,\, -3.25,\, -0.50828432,\, -3.25,\, \\
&\mathrel{\phantom{=}} -0.60388469,\, -3.25,\, -0.56644705,\, -3.25,\, \\
&\mathrel{\phantom{=}} -0.42772621,\, -3.25,\,  0.7716625,\, -3.25,\, \\
&\mathrel{\phantom{=}} -0.60336517,\, -3.25,\,  0.88158364,\, -3.25,\, \\
&\mathrel{\phantom{=}} 0.0269842,\, -3.25,\,   0.42512831,\, -3.25 \rangle
\end{split}
\label{eq:big-initial-state}
\end{equation}

for all our analysis of this particular system. Using the code in Appendix <ref>, we can generate an orbit of this initial state $O(\mathbf{X}_0)$ and approximate the system's Lyapunov spectrum using this orbit. In Figure <ref>, we graph the first thousand iterations of the fast variable orbits of the first eight Rulkov 1 neurons in the ring. We start with uncoupled neurons $g^e=0$ in Figure <ref>, where we can see these neurons with identical parameters are all out of phase in the non-chaotic spiking domain. As expected, because there is no current flow and all of the individual Rulkov neurons are spiking regularly, the maximal Lyapunov exponent $\lambda_1$ is negative. When we up the electrical coupling strength to $g^e=0.05$, in Figure <ref>, the neurons still spike relatively periodically, although there are some irregularities when one voltage happens to catch onto another. This small amount of coupling conductance is enough to bump the system into chaos with a positive maximal Lyapunov exponent. Next, we bring the coupling strength up significantly to $g^e=0.25$, where we can see in Figure <ref> that the ring system now has synchronized chaotic bursts. This is reminiscent of the positive coupling of two symmetrically coupled chaotic bursting Rulkov 1 neurons in Section <ref> (Figure <ref>), where the bursts happen (mostly) in sync with each other but the spikes within the bursts are still chaotic and unsynchronized. Finally, we take the coupling strength to the extreme, with $g^e=1$ in Figure <ref>, where complete chaos ensues due to each Rulkov neuron having a tremendous influence over its nearest neighbors. Here, the maximal Lyapunov exponent jumps up to $\lambda_1=0.1694$, the highest Lyapunov exponent we have seen in this paper thus far.

A natural question to ask is how the maximal Lyapunov exponent varies as we vary $g^e$. This can be answered using the code in Appendix <ref>, which creates a graph of the maximal Lyapunov exponent of a ring lattice system against its electrical coupling strength. For this first system, the results of this code are displayed in Figure <ref>. Our first observation is that the maximal Lyapunov exponents are rather erratic for higher values of $g^e$, covering a wide range of values for similar values of $g^e$. However, there are some trends we can notice. Because the individual neurons of this system are non-chaotic, $\lambda_1$ values initially start below zero. As the neurons start to affect each other, we reach the range of chaotic spiking, where the $\lambda_1$ values quickly become positive and reach a maximum. Then, as we enter the synchronized bursting regime, the $\lambda_1$ values become much more varied and start a generally downward trend, which can be attributed to the silence between chaotic bursts of spikes being non-chaotic and lowering the Lyapunov exponents. As we reach the extreme values of $g^e$ towards the right side of the graph, the $\lambda_1$ values shoot up to extremely high and chaotic values.

$g^e=0$, $\lambda_1\approx 0.0644$

$g^e=0.05$, $\lambda_1\approx 0.0686$

$g^e=0.25$, $\lambda_1\approx 0.0663$

$g^e=1$, $\lambda_1\approx 0.2003$

Graphs of $x_{i,\,k}$ for eight neurons in a ring of $\zeta=30$ electrically coupled Rulkov 1 neurons, with $x_{i,\,0}\in (-1,\,1)$, $y_{i,\,0}=-3.25$, $\sigma_i\in(-1.5,\,-0.5)$, and $\boldsymbol{\alpha} = 4.5\cdot\mathbf{1}$, visualized and $\lambda_1$ values calculated using the code in Appendix <ref>

$g^e=0$, $\lambda_1\approx 0.0469$

$g^e=0.05$, $\lambda_1\approx 0.0563$

$g^e=0.25$, $\lambda_1\approx 0.0633$

$g^e=1$, $\lambda_1\approx 0.2053$

Graphs of $x_{i,\,k}$ for eight neurons in a ring of $\zeta=30$ electrically coupled Rulkov 1 neurons, with $x_{i,\,0}\in (-1,\,1)$, $y_{i,\,0}=-3.25$, $\sigma_i\in(-1.5,\,-0.5)$, and $\alpha_i\in (4.25,\,4.75)$, visualized and $\lambda_1$ values calculated using the code in Appendix <ref>

We will now examine our second and third Rulkov 1 neuron ring lattice systems, where we have different parameters between different neurons in the ring. The second system we will examine keeps the same randomly distributed $x_{i,\,0}$ values (Equation <ref>), the same $y_{i,\,0}=-3.25$ values, and the same $\boldsymbol{\alpha} = 4.5\cdot\mathbf{1}$ vector, but we randomly choose $\sigma_i$ values from the interval $(-1.5,\,-0.5)$. Referring back to the bifurcation diagram of Rulkov map 1 in Figure <ref> from Section <ref>, the set $\sigma\in(-1.5,\,-0.5)\:\cap\: \alpha=4.5$ gives us individual neurons in silence, bursting, and spiking domains, which can be seen in the uncoupled neuron system (Figure <ref>). The specific random $\boldsymbol{\sigma}$ vector that we use in this paper is

\begin{equation}
\begin{split}
\boldsymbol{\sigma} &= \langle -0.63903048,\, -0.87244087,\, -1.16110093,\, \\
&\mathrel{\phantom{=}} -0.63908737,\, -0.73103576,\, -1.23516699,\, \\
&\mathrel{\phantom{=}} -1.09564519,\, -0.57564289,\, -0.75055299,\, \\
&\mathrel{\phantom{=}} -1.01278976,\, -0.61265545,\, -0.75514189,\, \\
&\mathrel{\phantom{=}} -0.89922568,\, -1.24012127,\, -0.87605023,\, \\
&\mathrel{\phantom{=}} -0.94846269,\, -0.78963971,\, -0.94874874,\, \\
&\mathrel{\phantom{=}} -1.31858036,\, -1.34727902,\, -0.7076453,\, \\
&\mathrel{\phantom{=}} -1.10631486,\, -1.33635792,\, -1.48435264,\, \\
&\mathrel{\phantom{=}} -0.76176103,\, -1.17618267,\, -1.10236959,\, \\
&\mathrel{\phantom{=}} -0.66159308,\, -1.27849639,\, -0.9145025 \rangle
\end{split}
\label{eq:big-sigma-vector}
\end{equation}

Finally, the third system we simulate is one where we keep the randomly distributed $x_{i,\,0}$ and $\sigma_i$ values, keep $y_{i,\,0}=-3.25$, but randomly choose $\alpha_i$ values from the interval $(4.25,\,4.75)$. This even further varies the distribution of different possible behaviors between different neurons in the system. This can be seen in the uncoupled neuron system (Figure <ref>), with some neurons exhibiting very rapid spiking. The specific random $\boldsymbol{\alpha}$ vector we use in this paper is

\begin{equation}
\begin{split}
\boldsymbol{\alpha} &= \langle 4.31338267,\, 4.3882788,\,  4.6578449,\, \\
&\mathrel{\phantom{=}} 4.67308374,\, 4.28873181,\, 4.26278301,\, \\
&\mathrel{\phantom{=}} 4.73065817,\, 4.29330435,\, 4.44416548,\, \\
&\mathrel{\phantom{=}} 4.66625973,\, 4.26243104,\, 4.65881579,\, \\
&\mathrel{\phantom{=}} 4.68086764,\, 4.44092086,\, 4.49639124,\, \\
&\mathrel{\phantom{=}} 4.55500032,\, 4.33389054,\, 4.38869161,\, \\
&\mathrel{\phantom{=}} 4.57278526,\, 4.62717616,\, 4.62025928,\, \\
&\mathrel{\phantom{=}} 4.49780551,\, 4.46750298,\, 4.49561326,\, \\
&\mathrel{\phantom{=}} 4.66902393,\, 4.60858869,\, 4.6027906,\, \\
&\mathrel{\phantom{=}} 4.40563641,\, 4.54198743,\, 4.49388045 \rangle
\end{split}
\label{eq:big-alpha-vector}
\end{equation}

In Figures <ref> and <ref>, we graph the first thousand iterations of the fast variable orbits of the first eight neurons in the ring for the same electrical coupling strength values as the first system: $g^e=0$, 0.05, 0.25, and 1. Comparing both of these systems to the first system, we can see similar patterns emerging. For $g^e=0.05$, the adjacent neurons start to have some effect on each other, but the overall picture of the dynamics stays the same. Raising the electrical coupling strength up to $g^e=0.25$, the neurons all undergo synchronized chaotic bursting, and going to the extreme $g^e=1$, complete chaos ensues. One of the interesting things that is even clearer in these visualizations is the direct effect that neurons have on their adjacent partners. Specifically, we can see that spiking in one neuron is reflected in adjacent neurons through smaller spikes during a time of silence.

$x_{i,\,0}\in (-1,\,1)$, $y_{i,\,0}=-3.25$, $\sigma_i\in(-1.5,\,-0.5)$, $\boldsymbol{\alpha} = 4.5\cdot\mathbf{1}$

$x_{i,\,0}\in (-1,\,1)$, $y_{i,\,0}=-3.25$, $\sigma_i\in(-1.5,\,-0.5)$, $\alpha_i\in(4.25,\,4.75)$

Graphs of the maximal Lyapunov exponent $\lambda_1$ against the electrical coupling strength $g^e$ for our second and third ring lattice systems, visualized using the code in Appendix <ref>

In Figure <ref>, we generalize and visualize the maximal Lyapunov exponents of these two systems for many values of $g^e$. An immediate difference we notice when comparing these graphs to the graph in Figure <ref> is that all of $\lambda_1$ values remain above 0 since some of the individual neurons in these varied parameter systems are chaotic. However, the graphs of maximal Lyapunov exponents of all three of our systems have similar shapes, the major differences being in the beginning where the neurons are weakly coupled and operating under their own parameters. Past this weak coupling domain, all three graphs in Figures <ref> and <ref> follow the increase up to chaotic spiking, the swoop down as synchronized chaotic bursts occurring, and the shoot up as we approach the extreme values of $g^e$. This shows that despite making individual neurons more and more different from their neighbors, coupling all of them together makes the systems exhibit similar behaviors and chaotic dynamics.

Although these systems do exhibit multistability in the sense that some $\mathbf{X}_0$ states get attracted to a non-chaotic attractor and some get attracted to a chaotic attractor (like what we saw in the asymmetrically coupled system in the previous section), the multistability exhibited here isn't interesting because, practically, the only states that get attracted to a non-chaotic attractor are ones where all the neurons in the system have the same parameters and the same (or very close to the same) initial states. This is because, with $\zeta=30$ neurons, it is nearly impossible for all the neurons to sync up by chance like it does when the system has only two neurons. For this reason, in Section <ref>, we will not be examining complex geometries that come out of multistability, namely, basins and fractal basin boundaries; instead, we will be examining the fractal structure of these systems' chaotic attractors.

§ GEOMETRICAL ANALYSIS OF RULKOV NEURON SYSTEMS

In this section, we combine all of the research outlined in this paper, applying our theoretical and computational methods for quantifying the chaos emerging from the complex geometrical structure of attractors and basins (Section <ref>) to the chaotic neuron systems from our discussion of the Rulkov maps (Sections <ref>, <ref>, <ref>). We will begin by analyzing the fractal geometry of chaotic spiking and chaotic bursting attractors of a single Rulkov 2 neuron. Next, we will perform an in-depth analysis of the attractors, basins, and basin boundaries of a system composed of two asymmetrically electrically coupled Rulkov 1 neurons, which we discover exhibits some very interesting and unexpected geometrical properties. In doing so, we define the concept of a chaotic pseudo-attractor, develop visualizations of higher-dimensional basins, and analyze lower-dimensional slices of basins and basin boundaries. Last, we will conclude by calculating and analyzing the Lyapunov dimensions of the chaotic attractors generated by three systems of $\zeta=30$ electrically coupled Rulkov 1 neurons in a ring lattice. By the Kaplan-Yorke conjecture, these Lyapunov dimensions give a rough approximation of the fractal dimensions of these attractors that live in $60$-dimensional space.

$\alpha=4.5$, $\sigma=0$

$\alpha=4.1$, $\sigma=-0.5$

Spiking and bursting attractors of systems of one Rulkov 2 neuron in two-dimensional state space $\langle y,\,x\rangle$, visualized using the code in Appendix <ref>

§.§ One Rulkov 2 Neuron

Recall from Section <ref> that the attractor of a single Rulkov 1 neuron is not fractal because the resetting mechanism of the map brings states that deviate from each other back to the same value. However, the chaotic Rulkov map 2 doesn't have a built-in resetting mechanism; it isn't piecewise. Therefore, we suspect that the chaotic attractors of Rulkov map 2 are indeed fractal. In this section, we will examine a Rulkov 2 neuron that exhibits chaotic spiking ($\alpha=4.5$ and $\sigma=0$, shown in Figure <ref>) and a Rulkov 2 neuron that exhibits chaotic bursting ($\alpha=4.1$ and $\sigma=-0.5$, shown in Figure <ref>). The attractors generated by these systems in two-dimensional state space $\langle y,\,x\rangle$ are visualized in Figure <ref>, where they certainly look complex and fractal-like. Before we analyze the attractors though, one interesting thing to note about the bursting attractor in Figure <ref> is its similarity to the bursting attractor of Rulkov map 1. Comparing the chaotic bursting attractor in Figure <ref> to the diagram of the state space of Rulkov map 1 in Figure <ref>, we can see the high $x$ value spiking in both, chaotic in Rulkov map 2 and non-chaotic in Rulkov map 1, followed by a sudden attraction to the stable branch $B_{\text{stable}}$, a slow movement up $B_{\text{stable}}$, then a jump back up to the spiking branch $B_{\text{spikes}}$.

Now, to examine the possible fractal geometry of these Rulkov map 2 spiking and bursting attractors, we will first consider the Kaplan-Yorke conjecture. In Section <ref>, we calculated that the Lyapunov spectrum of the spiking attractor generated by the system with parameters $\alpha=4.5$ and $\sigma=0$ is $\lambda \approx \{0.5449,\, -2.070\times 10^{-4}\}$, and the Lyapunov spectrum of the bursting attractor generated by the system with parameters $\alpha=4.1$ and $\sigma=-0.5$ is $\lambda \approx \{0.5025,\, -0.03376\}$. Recall from the end of Section <ref> that the Kaplan-Yorke conjecture relies on finding an index $\kappa$ such that

\begin{equation}
\sum_{i=1}^{\kappa+1} \lambda_i < 0
\end{equation}

However, for both the spiking and bursting attractor, summing both Lyapunov exponents will give us a value greater than 0. Therefore, the Kaplan-Yorke conjecture will be of no use to us here; to calculate the dimension of these attractors, we must return to the method of box-counting.

2cSpiking Attractor

$\epsilon$     $N(\epsilon)$

$1/100$     1696

$1/200$     5577

$1/400$     17720

2cBursting Attractor

$\epsilon$     $N(\epsilon)$

$1/50$     1375

$1/100$     4982

$1/200$     17968

Some $N(\epsilon)$ values for a spiking attractor of Rulkov map 2 generated by the system with parameters $\alpha=4.5$ and $\sigma=0$ and a bursting attractor of Rulkov map 2 generated by the system with parameters $\alpha=4.1$ and $\sigma=-0.5$, both calculated approximately on the ranges shown in Figure <ref> using the code in Appendices <ref> and <ref>

In Section <ref>, we used the code in Appendix <ref> to count boxes on the Hénon attractor. However, this box-counting code works for box-counting on any attractor, so we can also use Appendix <ref> to determine $N(\epsilon)$ values on the spiking and bursting attractors of Rulkov map 2. The results from box-counting on the Rulkov 2 attractors in Figure <ref> are displayed in Table <ref>. Recalling from Section <ref> that we expect $N(\epsilon)$ to be proportional to $1/\epsilon^d$, we can take logarithms of the data in the table $N(\epsilon)$ vs. $1/\epsilon$ and perform a linear regression to extract the fractal dimension $d$, which produces stunningly good results. Considering the spiking attractor first, we get that

\begin{equation}
\log_2 N(\epsilon) = 1.693\log_2\epsilon - 0.509
\end{equation}

with an $R^2$ value of 0.9999. This indicates that the dimension of the spiking attractor is $d\approx 1.693$, which is not an integer, indicating that the spiking attractor shown in Figure <ref> is indeed a fractal. It is certainly believable that the attractor shown in Figure <ref> is closer to being two-dimensional than one-dimensional, as the points that make up the attractor seem to fill up a lot of two-dimensional state space. Performing a similar linear regression on the bursting attractor values in Table <ref> yields

\begin{equation}
\log_2 N(\epsilon) = 1.854\log_2\epsilon - 0.037
\end{equation}

with an $R^2$ value of 0.999999. This indicates that the bursting attractor is also fractal with a dimension of $d\approx 1.854$. This higher fractal dimension indicates that the bursting attractor in Figure <ref> is more “rough” than the spiking attractor in Figure <ref>, which is not immediately obvious. However, the main result from this section is that an isolated Rulkov neuron can indeed form a true strange attractor in two-dimensional state space, both for chaotic spiking and chaotic bursting systems.

§.§ Asymmetrical Electrical Coupling of Two Rulkov 1 Neurons

As we discovered in Section <ref>, multistability and chaotic dynamics appear to exist in a system of two asymmetrically electrically coupled Rulkov 1 neurons. In this section, we will analyze the many complex and interesting geometries of state space that appear in the system from Section <ref>, the asymmetrically coupled neuron system with $\beta^c_1=\sigma^c_1=\beta^c_2=\sigma^c_2=1$, $\sigma_1=\sigma_2=-0.5$, $\alpha_1=\alpha_2=4.5$, $g^e_1=0.05$, and $g^e_2=0.25$.

§.§.§ Attractor Analysis

In Figure <ref>, we can see that this system appears to have two distinct attractors: a non-chaotic spiking attractor and a chaotic attractor. However, upon further analysis of this system, we discover that this is not completely accurate. To see this, let us examine two seemingly chaotic orbits of this system. First, we will look at the system from Figure <ref>, which has initial states $\mathbf{x}_{1,\,0}=\langle -0.56,\, -3.25\rangle$ and $\mathbf{x}_{2,\,0}=\langle -1,\, -3.25\rangle$. Second, we will look at the system with initial states $\mathbf{x}_{1,\,0} = \langle 0.35,\, -3.25 \rangle$ and $\mathbf{x}_{2,\,0} = \langle -1.23,\, -3.25 \rangle$. Using the code in Appendix <ref>, we graph the fast variable orbits of these systems up to a large value of $k$ ($O(\mathbf{X_0})=\{\mathbf{X}_0,\,\mathbf{X}_1,\,\hdots\,,\,\mathbf{X}_{150000}\}$) in Figures <ref> and <ref>. As we can see, while the voltage orbits begin chaotic, they eventually fall into the non-chaotic spiking attractor. This is made especially clear in Figures <ref> and <ref>, where we zoom into the region where the orbits transition from chaotic to non-chaotic. We conjecture that these dynamics occur because while the neurons are exhibiting chaotic bursting, their voltages happen to line up and latch onto each other's spiking by chance, which propels the system's orbit to the non-chaotic spiking attractor. If this is true, then a true chaotic attractor of this system doesn't exist because everything will eventually, by chance, get attracted to the non-chaotic spiking attractor. For this reason, true multistability as we have defined it in Section <ref> doesn't exist here. However, we are still interested in analyzing this system because it exhibits sensitive dependence on initial conditions; although all orbits end up in the spiking attractor after enough time, initial conditions still have an effect on the qualitative dynamics of the system in the short term. For this reason, we will treat the short-term chaotic behavior exhibited by this system as being on a chaotic “pseudo-attractor,” which is where some orbits end up in before eventually being attracted to the non-chaotic spiking attractor.

$\mathbf{x}_{1,\,0} = \langle -0.56,\, -3.25 \rangle$, $\mathbf{x}_{2,\,0} = \langle -1,\, -3.25 \rangle$

$\mathbf{x}_{1,\,0} = \langle 0.35,\, -3.25 \rangle$, $\mathbf{x}_{2,\,0} = \langle -1.23,\, -3.25 \rangle$

(Zoom) $\mathbf{x}_{1,\,0} = \langle -0.56,\, -3.25 \rangle$, $\mathbf{x}_{2,\,0} = \langle -1,\, -3.25 \rangle$

(Zoom) $\mathbf{x}_{1,\,0} = \langle 0.35,\, -3.25 \rangle$, $\mathbf{x}_{2,\,0} = \langle -1.23,\, -3.25 \rangle$

Graphs showing the eventual attraction of two orbits of the system of two asymmetrically electrically coupled Rulkov 1 neurons to the system's non-chaotic spiking attractor, graphed using the code in Appendix <ref>

Four-dimensional state space $\langle x_1,\, y_1,\, x_2,\, y_2
\rangle$ projected onto two-dimensions $\langle x_1,\, y_1
\rangle$

Four-dimensional state space $\langle x_1,\, y_1,\, x_2,\, y_2
\rangle$ projected onto two-dimensions $\langle x_2,\, y_2
\rangle$

Projections of the four-dimensional non-chaotic attractor generated by the system of two asymmetrically coupled Rulkov 1 neurons with initial states $\mathbf{x}_{1,\,0}=\langle -0.54,\, -3.25 \rangle$ and $\mathbf{x}_{2,\,0}=\langle -1,\, -3.25 \rangle$ and Lyapunov spectrum $\lambda\approx\{-0.0057,\,-0.0125,\,-0.0126,\,-\infty\}$, graphed with the code in Appendix <ref>

Four-dimensional state space $\langle x_1,\, y_1,\, x_2,\, y_2
\rangle$ projected onto two-dimensions $\langle x_1,\, y_1
\rangle$

Four-dimensional state space $\langle x_1,\, y_1,\, x_2,\, y_2
\rangle$ projected onto two-dimensions $\langle x_2,\, y_2
\rangle$

Projections of the four-dimensional strange pseudo-attractor generated by two asymmetrically coupled Rulkov 1 neurons with parameters $\sigma_1=\sigma_2=-0.5$, $\alpha_1=\alpha_2=4.5$, coupling strengths $g^e_1=0.05$, $g^e_2=0.25$, initial states $\mathbf{x}_{1,\,0}=\langle -0.56,\, -3.25 \rangle$, $\mathbf{x}_{2,\,0}=\langle -1,\, -3.25 \rangle$, and Lyapunov spectrum $\lambda\approx\{0.0165,\,-0.0105,\,-0.0810,\,-\infty\}$, graphed with the code in Appendix <ref>

We are now interested in visualizing the non-chaotic spiking attractor and chaotic pseudo-attractor. Immediately, we run into the problem of the state space of this system being four-dimensional, which we cannot make a graph of. Despite this, as an attempt to visualize the attractors, we will project four-dimensional state space $\langle x_1,\, y_1,\, x_2,\, y_2\rangle$ onto two different two-dimensional surfaces, $\langle x_1,\, y_1\rangle$ and $\langle x_2,\, y_2\rangle$, by plotting only the $\mathbf{x}_1$ or $\mathbf{x}_2$ components of each state vector $\mathbf{X}$, respectively.[It is important to note that these projections do not capture the true geometry of these attractors in four-dimensional state space, as one of these projections is dimensionally analogous to compacting an entire region of three-dimensional state space onto a line segment. However, because we don't live in four spatial dimensions (at least, macroscopically), this is one of the best ways we have to visualize these attractors.] This is accomplished in Figures <ref> and <ref> for the initial conditions $\mathbf{X}_{0}=\langle -0.54,\, -3.25,\, -1,\, -3.25 \rangle$ (associated with Figure <ref>) and $\mathbf{X}_{0}=\langle -0.56,\, -3.25,\, -1,\, -3.25 \rangle$ (associated with Figure <ref>), respectively. It is immediately clear that the attractor in Figure <ref> is non-chaotic and the pseudo-attractor in Figure <ref> is chaotic.

However, the question now is whether or not the chaotic pseudo-attractor is strange. We know from Section <ref> that there is a possibility that it isn't because of the resetting mechanism of Rulkov map 1, which is what makes the attractor of a chaotic uncoupled Rulkov 1 neuron not fractal. We could try to use the Kaplan-Yorke conjecture again, but we should be cautious in this case because we are dealing with a pseudo-attractor, not a normal attractor. As a result, the Lyapunov exponents are only approximations inside this pseudo-attractor;[As shown in the code in Appendix <ref>, we calculate these Lyapunov exponents up to 65,000 iterations only, which is approximately when the orbit goes to the spiking attractor (see Figure <ref>).] calculating the spectrum for a larger number of iterations will contain information about the orbit on the non-chaotic spiking attractor. For this reason, it would be good to return to our good old-fashioned box-counting method. Using the code in Appendix <ref> for four-dimensional boxes of size $\epsilon=1/20,\,1/30,\,1/40$,[We calculate $N(\epsilon)$ for a small number of $\epsilon$ values because, since we are dealing with box-counting in four-dimensional state space, the code in Appendix <ref> has time complexity $O(\epsilon^4)$.] we get that the Minkowski–Bouligand dimension of the chaotic pseudo-attractor is $d\approx 1.84$ with an $R^2$ value of 0.9999. This is not an integer, so the chaotic pseudo-attractor is indeed fractal and strange.

Even though we suspect that the Kaplan-Yorke conjecture will not hold in this case, we can still calculate it for fun since it should still tell us that the pseudo-attractor is fractal even if it doesn't give us the correct dimension. The Lyapunov spectrum of the four-dimensional strange pseudo-attractor that $\mathbf{X}_{0}=\langle -0.56,\, -3.25,\, -1,\, -3.25 \rangle$ goes to is $\lambda\approx\{0.0165,\,-0.0105,\,-0.0810,\,-\infty\}$, which is calculated using the code in Appendix <ref>. Interestingly, in this case, the $\kappa$ index relevant to calculating the Lyapunov dimension $d_l$ is 2 since

\begin{equation}
\sum_{i=1}^2 \lambda_i \approx 0.006 > 0
\end{equation}

Therefore, the Lyapunov dimension of the strange pseudo-attractor is

\begin{equation}
\begin{split}
d_l &= \kappa + \frac{1}{|\lambda_{\kappa+1}|}\sum_{i=1}^{\kappa}\lambda_i \\
&= 2 + \frac{1}{|\lambda_{3}|}\sum_{i=1}^{2}\lambda_i \\
&\approx 2.07
\end{split}
\end{equation}

Although this is different from the true fractal dimension of the pseudo-attractor, the Kaplan-Yorke conjecture still indicates that the pseudo-attractor is fractal even though it indicated that the chaotic attractor of an isolated Rulkov 1 neuron wasn't (see Section <ref>).

§.§.§ Basin Analysis

Due to the short-term multistability of this asymmetrically coupled system, we are now interested in analyzing the basins of the spiking attractor and strange pseudo-attractor for interesting geometrical properties. Because the state space of our system is four-dimensional, these basins are four-dimensional sets. However, because of the difficulty in visualizing four-dimensional objects, as well as the fact that our coupled Rulkov 1 neurons tend to settle down in a small interval of $y$ values (see Figures <ref> and <ref>), we will first consider a two-dimensional slice of the basins, namely, the intersection of the four-dimensional basins with the set of state space

\begin{multline}
S_2'= \{\mathbf{X} = \langle x_1,\,y_1,\,x_2,\,y_2\rangle: \\-2<x_1<2,\,y_1=-3.25, \\ -2<x_2<2,\,y_2=-3.25\}
\label{eq:s2'}
\end{multline}

To determine which attractor a given state in $S_2'$ is initially attracted to, we will need to define a specific way of making the distinction between immediately going to the spiking attractor and spending some time in the chaotic pseudo-attractor first. In Figure <ref>, we see that states that are attracted to the chaotic pseudo-attractor spend on the order of $10^4$ iterations in the pseudo-attractor before going to the spiking attractor. Testing other initial states using the code in Appendix <ref> gives similar results. We will therefore say that an initial state that has a positive maximal Lyapunov exponent after 5,000 iterations was attracted to the chaotic pseudo-attractor, while an initial state that has a negative $\lambda_1$ after the same amount of time was attracted to the spiking attractor. It is worthwhile to emphasize that this way of detecting which attractor an initial state goes to is very different from the method we used before in Section <ref>, which was directly testing whether a state was in the vicinity of a given attractor's points. Instead, we are indirectly testing what attractor an initial state goes to using the power of Lyapunov exponents.

Maximal Lyapunov exponent color map, with red points indicating chaotic
dynamics and blue points indicating non-chaotic dynamics

Basin slices of the chaotic (white) and non-chaotic spiking (black) attractors

A slice of the basins of attraction of two asymmetrically electrically coupled Rulkov 1 neurons with $-2\leq x_{1,\,0}\leq 2$ on the horizontal axis and $-2\leq x_{2,\,0}\leq 2$ on the vertical axis, parameters $\sigma_1=\sigma_2=-0.5$, $\alpha_1=\alpha_2=4.5$, coupling strengths $g^e_1=0.05$, $g^e_2=0.25$, and fixed initial $y$ values $y_{1,\,0}=y_{2,\,0}=-3.25$, visualized using the code in Appendix <ref>

Using the code in Appendix <ref>, we test a large number of initial states in $S_2'$ using this method and present the results of these Lyapunov exponent calculations using a color map in Figure <ref>. An immediate observation is the beautiful complexity of the different maximal Lyapunov exponents in state space. This is further emphasized when we color the basins properly in Figure <ref>, with the basin of the chaotic pseudo-attractor (positive $\lambda_1$) shown in white and the basin of the non-chaotic attractor (negative $\lambda_1$) shown in black. Here, we can see that most of $S_2'$ is taken up by the basin of the white chaotic pseudo-attractor, with some black regions, curves, and scattered points of stability scattered throughout. The most obvious of these features is the black line that goes across the diagonal of Figure <ref> from the bottom left to the top right. This line and the areas surrounding it are places where the initial states of the two neurons are the same or very close to being the same. In this case, it makes sense that they will immediately be synchronized, spiking together in the non-chaotic attractor.

Now that we have visualized the basins in $S_2'$, we will now classify (using the method detailed in Section <ref>) the basin slices in the infinite two-dimensional plane containing $S_2'$. In other words, we are classifying the intersection of the four-dimensional basins with the infinite plane that lies parallel to and contains the square $S_2'$, which we will denote as $S_2$:

\begin{multline}
S_2 = \{\mathbf{X} = \langle x_1,\,y_1,\,x_2,\,y_2\rangle: \\x_1\in\mathbb{R},\,y_1=-3.25, \\ x_2\in\mathbb{R},\,y_2=-3.25\}
\end{multline}

To account for $S_2$ being two-dimensional, we will also need to specify what is meant by the mean and standard deviation of the four-dimensional attractors within the two-dimensional plane $S_2$. Let us denote a generic four-dimensional attractor as $A_4$ and its associated basin as $\hat{A}_4$. In this case, $A_4$ can represent either the spiking attractor or the chaotic pseudo-attractor. From Equation <ref>, if

\begin{equation}
a_1\e{1} \\[4pt]
a_1\e{2} \\[4pt]
a_1\e{3} \\[4pt]
\end{pmatrix},\,
\begin{pmatrix}
a_2\e{1} \\[4pt]
a_2\e{2} \\[4pt]
a_2\e{3} \\[4pt]
\end{pmatrix},\,\hdots
\right\}
\end{equation}

then the true mean of $A_4$, which we will denote as $\langle A_4\rangle$, is

\begin{equation}
\langle A_4\rangle = \lim_{j\to\infty}\frac{1}{j}\sum_{i=1}^j\mathbf{a}_i = \begin{pmatrix}
\langle x_1 \rangle \\[4pt]
\langle y_1 \rangle \\[4pt]
\langle x_2 \rangle \\[4pt]
\langle y_2 \rangle
\end{pmatrix}
\end{equation}

We introduce $\langle A_2\rangle$ to be the “effective mean” of the four-dimensional attractor $A_4$ that we will use when classifying the two-dimensional slice of a basin, or $\hat{A}_4\cap S_2$. We define it simply as a two-dimensional vector composed of the first and third entries of $\langle A_4\rangle$:

\begin{equation}
\langle A_2\rangle = \lim_{j\to\infty}\frac{1}{j}\sum_{i=1}^j\begin{pmatrix}
a_i\e{1} \\[4pt]
\end{pmatrix} = \begin{pmatrix}
\langle x_1 \rangle \\[4pt]
\langle x_2 \rangle
\end{pmatrix}
\end{equation}

Similarly, the true standard deviation of $A_4$ in four-dimensional state space is, from Equation <ref>,

\begin{equation}
\sigma_{A4} = \sqrt{\lim_{j\to\infty}\frac{1}{j}\sum_{i=1}^j|\mathbf{a}_i-\langle A_4\rangle|^2}
\end{equation}

It follows that the “effective standard deviation” $\sigma_{A2}$ of $A_4$ that we will use when classifying $\hat{A}_4\cap S_2$ is

\begin{equation}
\begin{split}
\sigma_{A2} &= \sqrt{\lim_{j\to\infty}\frac{1}{j}\sum_{i=1}^j\left|\begin{pmatrix}
a_i\e{1} \\[4pt]
\end{pmatrix}-\langle A_2\rangle\right|^2} \\
&= \sqrt{\lim_{j\to\infty}\frac{1}{j}\sum_{i=1}^j\left|\begin{pmatrix}
a_i\e{1}-\langle x_1 \rangle \\[4pt]
a_i\e{3}-\langle x_2 \rangle
\end{pmatrix}\right|^2}
\end{split}
\end{equation}

For our classification of the basin slices $\hat{A}_4\cap S_2$, we will consider the normalized two-dimensional distance $\xi_2$ from $\langle A_2\rangle$ of some

\begin{equation}
\mathbf{X} = \begin{pmatrix}
x_1 \\
-3.25 \\
x_2 \\
\end{pmatrix} \in S_2
\end{equation}

to be

\begin{equation}
\xi_2 = \frac{1}{\sigma_{A2}}\left|\begin{pmatrix}
x_1 \\ x_2
\end{pmatrix} - \langle A_2\rangle\right|
\end{equation}

Recall from Section <ref> that the goal of our basin classification method is to find a function $P(\xi)$ that represents the probability that a randomly selected initial state from an $n$-dimensional ball with radius $\xi$ centered at $\langle A\rangle$ is in a given basin. We know from Equation <ref> that in the limit $\xi\to\infty$, this function follows the power law

\begin{equation}
P(\xi) = \frac{P_0}{\xi^{\gamma}}
\label{eq:probability-function-power-law-2}
\end{equation}

where $P_0$ and $\gamma$ are parameters that we use to classify a basin. Our first goal for basin classification is then to find two functions $P_w(\xi_2)$ and $P_b(\xi_2)$, which are the probability functions associated with the white (chaotic) and black (non-chaotic) basins, respectively. In Appendix <ref>, we implement the Monte Carlo algorithm outlined in Section <ref> for classifying basins by altering the code in Appendix <ref>, which we used to classify the basin of the Hénon map. The first alteration we made was changing the way we detect which basin a given state goes to by using Lyapunov exponents rather than a direct location check. Another alteration we made was that because we are testing two basins, we added a boolean to be able to choose which basin we want to be testing for. One important numerical note is that, given we are testing values of $\xi_2=2^k$, we must choose the maximum value of $k$ carefully. This is because it must be large enough that we capture the nature of the basins far away from the attractors since we are interested in the limit $\xi\to\infty$, but it also must be small enough that states get attracted to the attractor fast enough for an accurate Lyapunov exponent calculation on the attractor. By experiment, we determine that a good maximum value of $k$ is 11.

$\xi_2$     $P_w(\xi_2)$     $P_b(\xi_2)$

$2^6=64$     0.8365     0.1530

$2^7=128$     0.8736     0.1237

$2^8=256$     0.8919     0.1010

$2^9=512$     0.9148     0.0839

$2^{10}=1024$     0.9419     0.0675

$2^{11}=2048$     0.9360     0.0544

Some approximate $P(\xi_2)$ values of the white and black basins of the asymmetrically coupled Rulkov 1 neuron system, both directly calculated using the code in Appendix <ref>

In Table <ref>, we show the results from running the code in Appendix <ref>, where we neglect values of $k$ from 0 to 5 as we are interested in large values of $\xi_2$. Our first observation is that $P_w(\xi_2)$ is behaving a bit oddly, as Equation <ref> indicates to us that $P(\xi)$ functions should either decrease or stay the same. For this reason, we turn our attention to $P_b(\xi_2)$, which is behaving as we expect. Running a linear regression calculation for $\log_2 P_b(\xi_2)$ against $\log_2\xi_2$, we get

\begin{equation}
\log_2 P_b(\xi_2) = -0.2957\log_2\xi_2 - 0.9360
\label{eq:asym-basin-2-lin-reg}
\end{equation}

Since $\gamma\approx0.2957$ is between 0 and $n=2$, we can conclude from Section <ref> that the basin slice of the non-chaotic spiking attractor is Class 3, meaning it extends to infinity in some directions but takes up an increasingly small fraction of state space. This makes sense because the farther away we start from the steady spiking attractor, the less chance there will be for the neurons to immediately synchronize with each other. In the limit $\xi_2\to\infty$, we expect that the only initial conditions that get attracted to the spiking attractor immediately are the ones where $x_{1,\,0}=x_{2,\,0}$, which is the diagonal shown in Figure <ref>. For this reason, we can conclude that the black basin slice in $S_2$ has finite measure since an infinitely long line contributes nothing to the measure of the two-dimensional black basin slice.[See Property 4 of measures (Equation <ref>) in Section <ref>.] To classify the white basin, we must make the realization that

\begin{equation}
P_w(\xi_2) + P_b(\xi_2) = 1
\label{eq:sum-the-basins}
\end{equation}

which is the case because of the way we defined the spiking attractor and chaotic pseudo-attractor: if an initial state does not immediately get attracted to a spiking orbit with a negative $\lambda_1$, then it goes to a chaotic orbit with a positive $\lambda_1$. Therefore, there is no other basin that would cause the probability of being in either the white or black basin to not be 100%. This fact can be confirmed by adding the Monte Carlo data in Table <ref>, which gives us values very close to 1. Therefore, even though $P_w(2^{11})$ goes down from $P_w(2^{10})$ due to the aforementioned difficulty of calculating Lyapunov exponents for orbits starting from initial states far away from the attractors, the values of $P_w(\xi_2)$ are indeed approaching 1 because from exponentiating both sides of Equation <ref>,

\begin{equation}
P_b(\xi_2) \approx \frac{0.5227}{\xi_2^{0.2957}}
\end{equation}

which goes to 0 as $\xi_2\to\infty$. Therefore, because we are interested in this very limit, the basin slice of the chaotic pseudo-attractor must be Class 1. Although it does not take up all of $S_2$, it does take up all of $S_2$ barring a set of finite measure, which is the black basin. If we were to choose a random initial state from the entirety of $S_2$, there would indeed be a 100% probability that the state gets attracted to the chaotic pseudo-attractor because the black basin has an infinite measure in $S_2$, whereas the white basin has only a finite measure.

Now that we have classified the basin slices that lie in $S_2$, we are now interested in the entire basins living in all of four-dimensional space $S_4$:

\begin{equation}
S_4 = \{\mathbf{X}=\langle x_1,\,y_1,\,x_2,\,y_2\rangle:\:\mathbf{X}\in\mathbb{R}^4\}
\end{equation}

It is evident that we should define normalized distance in this set as

\begin{equation}
\xi_4 = \frac{|\mathbf{X}-\langle A_4\rangle|}{\sigma_{A4}}
\end{equation}

The tricky part now is how to choose a random state in a given four-dimensional ball with radius $\xi_4$ centered at $\langle A_4\rangle$. In two dimensions, this is easy; we simply pick a random number between $0$ and $\xi$ that we say is the distance $r$ away from $\langle A\rangle$ and pick another random number between $0$ and $2\pi$ that we say is the angle $\phi$ from the positive $x_1$ axis. Then, our random initial state $\mathbf{x}_0$ from the two-dimensional disk $|\mathbf{x}-\langle A\rangle|<\xi$ is

\begin{equation}
\mathbf{x}_0 = \langle A_2\rangle +
\begin{pmatrix}
r\cos\phi \\
\end{pmatrix}
\end{equation}

This is the method we implement in Appendices <ref> and <ref> to choose a random state from a given two-dimensional ball. However, the way to do this in four dimensions is not at all obvious, so this challenge in computing $P(\xi_4)$ is worth discussion. If we temporarily take things down a dimension, the analogous way to pick a random state $\mathbf{x}_0$ in three dimensions is to use spherical coordinates. Using physicists' notation, any point in a three-dimensional ball with radius $\xi$ centered at some point $\langle A\rangle\in\mathbb{R}^3$ can be described in spherical coordinates using three numbers

\begin{equation}
\begin{cases}
r\in[0,\,\xi) \\
\theta\in[0,\,\pi) \\
\phi\in[0,\,2\pi)
\end{cases}
\end{equation}

Then, selecting random values of $r$, $\phi$, and $\theta$ in the intervals gives us a random initial state $\mathbf{x}_0$ that can be calculated with

\begin{equation}
\mathbf{x}_0 = \langle A\rangle + \begin{pmatrix}
r\sin\theta\cos\phi \\
r\sin\theta\sin\phi \\
\end{pmatrix}
\end{equation}

We now have a way of solving our original problem of choosing a random state from a four-dimensional ball with radius $\xi_4$ centered at $\langle A\rangle$ if we have a four-dimensional analogy to spherical coordinates. In a paper by blumenson, the $n$-dimensional coordinate system analogous to spherical coordinates is derived using linear algebra, so we can utilize these results for our four-dimensional case. The four coordinates used are

\begin{equation}
\begin{cases}
r\in[0,\,\xi_4) \\
\theta_1\in[0,\,\pi) \\
\theta_2\in[0,\,\pi) \\
\phi\in[0,\,2\pi)
\end{cases}
\end{equation}

Then, choosing these numbers randomly from the intervals, our random initial state $\mathbf{X}_0\in S_4$ is given by

\begin{equation}
\mathbf{X}_0 = \langle A_4\rangle + \begin{pmatrix}
r\sin\theta_1\sin\theta_2\cos\phi \\
r\sin\theta_1\sin\theta_2\sin\phi \\
r\sin\theta_1\cos\theta_2 \\
\end{pmatrix}
\end{equation}

Now that we know how to pick random states in $S_4$ to test, our Monte Carlo algorithm for classifying the full basins in four-dimensional space is a natural extension from our two-dimensional version, and we display the code used to do it in Appendix <ref>. However, we wish to make a few numerical notes before we detail our analysis. First, because $y$ is a slow variable, starting far away from where it eventually ends up (around $y=-3.25$) will inevitably result in the state taking even longer to reach its attractor. For this reason, we up the number of times we iterate test points for Lyapunov exponent calculation from 5000 to 20000, which raises computation time significantly. Additionally, recall from Section <ref> that we use a shell method for calculating higher values of $P(\xi)$. As a quick reminder, our shell method says that given we know $P(2^k)$ and $\Delta P(2^k)$, which is the probability that an initial state chosen from the $n$-dimensional shell with inner radius $2^k$ and outer radius $2^{k+1}$ centered at $\langle A\rangle$ is in the basin, we can calculate $P(2^{k+1})$ using Equation <ref>, which says that

\begin{equation}
P(2^{k+1}) = \frac{P(2^k)}{2^n} + \left(1-\frac{1}{2^n}\right)\Delta P(2^k)
\end{equation}

This depends on the dimension of the state space $n$ that the basin lives in because of the way objects scale differently based on their dimension.[See Section <ref> for a reminder.] Therefore, we must alter our shell method calculations because we are dealing with four-dimensional state space as opposed to a two-dimensional one.

In Appendix <ref>, we implement our spherical coordinate method of choosing random states in four-dimensional balls and the discussed numerical considerations to calculate $P(\xi_4)$ values for the white and black basins in $S_4$. The results of this code for the white basin $P_w(\xi_4)$ are shown in Table <ref>. In the table, $P_b(\xi_4)$ values are indirectly calculated by using the fact that $P_w(\xi_4)+P_b(\xi_4)=1$, which follows directly from Equation <ref>. Observing the values in Table <ref>, it appears that $P_w(\xi_4)$ and $P_b(\xi_4)$ stay relatively constant as we vary $\xi_4$. Given our analysis of the two-dimensional basin slices, this is a rather unexpected result, so we will proceed carefully with our analysis.

$\xi_4$     $P_w(\xi_4)$     $P_b(\xi_4) = 1-P_w(\xi_4)$

$2^0=1$     0.8480     0.1520

$2^1=2$     0.8649     0.1351

$2^2=4$     0.8491     0.1509

$2^3=8$     0.8537     0.1463

$2^4=16$     0.8690     0.1310

$2^5=32$     0.8549     0.1451

$2^6=64$     0.8766     0.1234

$2^7=128$     0.8760     0.1240

$2^8=256$     0.8629     0.1371

Some approximate $P(\xi_4)$ values of the white and black basins of the asymmetrically coupled Rulkov 1 neuron system, $P_w(\xi_4)$ values directly calculated using the code in Appendix <ref>, $P_b(\xi_4)$ indirectly calculated using the $P_w(\xi_4)$ values

Running linear regressions for the log-log plots of both $P_w(\xi_4)$ and $P_b(\xi_4)$, we get that

\begin{equation}
P_w(\xi_4) = 0.0042\xi_4 - 0.2315
\label{eq:asym-white-basin-4-lin-reg}
\end{equation}

with $R^2=0.392$, and

\begin{equation}
P_b(\xi_4) = -0.0260\xi_4 - 2.7540
\label{eq:asym-black-basin-4-lin-reg}
\end{equation}

with $R^2=0.386$. These low $R^2$ values confirm our suspicion that there is no upward or downward trend as we change $\xi_4$. This means that within the numerical variability of the Monte Carlo algorithm, the value of the basin classification parameter $\gamma$ is effectively 0. The value of the other basin classification parameter $P_0$ can be found by exponentiating both sides of Equations <ref> and <ref>, which yields that for the white basin, $P_0=0.8517$, and for the black basin, $P_0=0.1482$. These are both between 0 and 1, so by the classification method established in Section <ref>, both the white and black basins are Class 2, meaning they both occupy a fixed fraction of four-dimensional state space. This makes sense considering the data in Table <ref>, but again, it is unexpected. Our explanation for why the white and black basins are Class 2 is that since Rulkov map 1 is a slow-fast system, $x_1$ and $x_2$ adjust much faster than $y_1$ and $y_2$. Therefore considering some initial condition $\mathbf{X}_0=\langle x_{1,\,0},\,y_{1,\,0},\,x_{2,\,0},\,y_{2,\,0} \rangle$ some distance $\xi_4$ away from either attractor $A_4$, $y_1$ and $y_2$ will slowly drift towards $\langle y_1\rangle$ and $\langle y_2\rangle$ while $x_1$ and $x_2$ easily keep up with their slow attraction. By the time $y_1$ and $y_2$ find their way to $\langle y_1\rangle$ and $\langle y_2\rangle$, $x_1$ and $x_2$ will also be near $\langle x_1\rangle$ and $\langle x_2\rangle$ regardless of how much farther away they started compared to $y_1$ and $y_2$ since they evolve so much more quickly. Thus, once $y_1$ and $y_2$ get close to the attractors, we are effectively starting in $S_2'$ (or a square of similar size close and parallel to it). Indeed, comparing the fractions of state space the white ($P_0=0.8517$) and black ($P_0=0.1482$) basins take up to the $P_w(\xi_2)$ and $P_b(\xi_2)$ fractions for the small values of $\xi_2$ in Table <ref>, we find that these align quite well.

Maximal Lyapunov exponent color map, with red points indicating chaotic
dynamics and blue points indicating non-chaotic dynamics

Basin slices of the chaotic (white) and non-chaotic spiking (black) attractors

A slice of the basins of attraction of two asymmetrically electrically coupled Rulkov 1 neurons with $-256\leq y_{1,\,0}\leq 256$ on the horizontal axis and $-256\leq y_{2,\,0}\leq 256$ on the vertical axis, parameters $\sigma_1=\sigma_2=-0.5$, $\alpha_1=\alpha_2=4.5$, coupling strengths $g^e_1=0.05$, $g^e_2=0.25$, and fixed initial $x$ values $x_{1,\,0}=-1$ and $x_{2,\,0}=1$, visualized using the code in Appendix <ref>

Although we cannot truly visualize these four-dimensional basins, we can get a grasp of this Class 2 basin behavior by graphing a different two-dimensional slice of the basins, namely, one where we vary $y_1$ and $y_2$ but keep $x_1$ and $x_2$ fixed. In Figure <ref>, we graph a large range of initial $y$ values $-256\leq y_{1,\,0}\leq 256$ and $-256\leq y_{2,\,0}\leq 256$ while keeping the initial $x$ values fixed at $x_{1,\,0}=-1$ and $x_{2,\,0}=1$. In Figure <ref>, we can see a dispersed, seemingly random distribution of white and black basin points with no definite pattern in their arrangement. Similarly in Figure <ref>, we can see that the maximal Lyapunov exponent also appears random, perhaps with the exception of a region of darker red points. What is important to notice from these graphs is that the distribution of white and black points in Figure <ref> appears to remain the same no matter where we are on the graph even though we are spanning a huge range of different $y$ values. This is characteristic of Class 2 basins, as the fraction of state space they take up doesn't depend on the distance from their attractors. Something else that we notice in Figure <ref> is a lack of the black diagonal that appears in Figure <ref>. This line of stability doesn't appear in these graphs because, as previously mentioned, it is defined by the equation

\begin{equation}
x_{1,\,0} = y_{1,\,0} = x_{2,\,0} = y_{2,\,0}
\end{equation}

which the slice in Figure <ref> doesn't intersect since the initial states shown in it have different $x$ values.

§.§.§ Basin Boundary Analysis

From Figures <ref> and <ref>, it is clear that the basins of the non-chaotic spiking attractor and strange pseudo-attractor are not divided by a clear boundary. This naturally leads us to suspect that the basin boundary $\Sigma$ between the white and black basins might be fractal, which we know from Section <ref> leads to geometrical sensitivity to initial conditions with an uncertainty exponent $\mathfrak{u}<1$. Similar to how we classified the intersection of the system's basins with $S_2$ and $S_4$, we are interested in examining the uncertainty exponents and the fractalization of both the intersection of $\Sigma$ with $S_2'$ and the intersection of $\Sigma$ with $S_4'$. We have already defined $S_2'$ as a specific square-shaped subset of $S_2$ (Equation <ref>), so we will now define $S_4'$ as the subset of $S_4$

\begin{multline}
S_4'= \{\mathbf{X} = \langle x_1,\,y_1,\,x_2,\,y_2\rangle: \\-2<x_1<2,\,-1<y_1<-5, \\ -2<x_2<2,\,-1<y_2<-5\}
\label{eq:s4'}
\end{multline}

which is a four-dimensional cube, or tesseract.

First, let us examine the set $\Sigma\cap S_2'$, displayed in Figure <ref> as the boundary between the white and black basins. We will denote the uncertainty exponent of this basin boundary set as $\mathfrak{u}_2$ and its associated probability function as $\varrho_2(\epsilon)$. As a reminder, $\varrho_2(\epsilon)$ is the probability that a randomly selected state in $S_2'$ is uncertain, and we expect it to be proportional to $\epsilon^{\mathfrak{u}_2}$. Calculating $\varrho_2(\epsilon)$ for our coupled neuron system mirrors the Monte Carlo method we used for the Hénon map in Section <ref>: we pick a random initial state $\mathbf{X}_0\in S_2'$ and test four specific perturbed states

\begin{equation}
\begin{pmatrix}
x_{1,\,0}+\epsilon \\
-3.25 \\
x_{2,\,0} \\
\end{pmatrix},\,
\begin{pmatrix}
x_{1,\,0}-\epsilon \\
-3.25 \\
x_{2,\,0} \\
\end{pmatrix},\,
\begin{pmatrix}
x_{1,\,0} \\
-3.25 \\
x_{2,\,0}+\epsilon \\
\end{pmatrix},\,
\begin{pmatrix}
x_{1,\,0} \\
-3.25 \\
x_{2,\,0}-\epsilon \\
\end{pmatrix}
\end{equation}

to see if any of them end up in a different basin than $\mathbf{X}_0$ does. The only major difference in calculating $\mathfrak{u}_2$ for this system as compared to the Hénon map is that we detect which basin an initial state belongs to using maximal Lyapunov exponents, a distinction that we have already discussed at length.

Although we cannot visualize it, calculating the uncertainty exponent of the basin boundary living in four-dimensional state space $\Sigma\cap S_4'$ is a natural extension from calculating $\mathfrak{u}_2$. Similar to before, let us denote the uncertainty exponent of this four-dimensional basin boundary set as $\mathfrak{u}_4$ and its associated probability function as $\varrho_4(\epsilon)$. The step of picking a random initial state $\mathbf{X}_0\in S_4'$ is much easier than picking a random initial state from a four-dimensional ball (which we did for basin classification) since we chose our region of analysis $S_4'$ to be a simple tesseract. Instead of testing four perturbed initial states, however, we need eight of them because we are now considering a four-dimensional region. Explicitly, these perturbed states are

\begin{equation}
\begin{gathered}
\begin{pmatrix}
x_{1,\,0}+\epsilon \\
y_{1,\,0} \\
x_{2,\,0} \\
\end{pmatrix},\,
\begin{pmatrix}
x_{1,\,0}-\epsilon \\
y_{1,\,0} \\
x_{2,\,0} \\
\end{pmatrix},\,
\begin{pmatrix}
x_{1,\,0} \\
y_{1,\,0}+\epsilon \\
x_{2,\,0} \\
\end{pmatrix},\,
\begin{pmatrix}
x_{1,\,0} \\
y_{1,\,0}-\epsilon \\
x_{2,\,0} \\
\end{pmatrix}, \\
\begin{pmatrix}
x_{1,\,0} \\
y_{1,\,0} \\
x_{2,\,0}+\epsilon \\
\end{pmatrix},\,
\begin{pmatrix}
x_{1,\,0} \\
y_{1,\,0} \\
x_{2,\,0}-\epsilon \\
\end{pmatrix},\,
\begin{pmatrix}
x_{1,\,0} \\
y_{1,\,0} \\
x_{2,\,0} \\
\end{pmatrix},\,
\begin{pmatrix}
x_{1,\,0} \\
y_{1,\,0} \\
x_{2,\,0} \\
\end{pmatrix}\phantom{,}
\end{gathered}
\end{equation}

$\epsilon$     $\varrho_2(\epsilon)$     $\varrho_4(\epsilon)$

$2^0=1$     0.710     0.404

$2^{-1}=1/2$     0.714     0.416

$2^{-2}=1/4$     0.658     0.394

$2^{-3}=1/8$     0.571     0.392

$2^{-4}=1/16$     0.516     0.339

$2^{-5}=1/32$     0.409     0.349

$2^{-6}=1/64$     0.340     0.341

$2^{-7}=1/128$     0.255     0.328

$2^{-8}=1/256$     0.219     0.323

$2^{-9}=1/512$     0.172     0.308

$2^{-10}=1/1024$     0.126     0.306

$2^{-11}=1/2048$     0.092     0.302

Some approximate $\varrho(\epsilon)$ values for the basin boundary sets $\Sigma\cap S_2'$ and $\Sigma\cap S_4'$ of the asymmetrically coupled Rulkov 1 neuron system, calculated using the code in Appendix <ref>

Now that we have discussed the theory behind calculating $\mathfrak{u}_2$ and $\mathfrak{u}_4$, let us begin our analysis. Because $S_2'$ and $S_4'$ are relatively small and close to the attractors, we can safely determine which basin an initial state in either set belongs to by calculating its Lyapunov spectrum using a 5000 iteration long orbit. In Appendix <ref>, we implement the theory and numerical considerations we discussed to calculate the values of the probability functions $\varrho_2(\epsilon)$ and $\varrho_4(\epsilon)$, the results of which are displayed in Table <ref>. Taking a linear regression of the log-log plot of $\varrho_2(\epsilon)$ first (neglecting the first two points because we are interested in the limit $\epsilon\to 0$), we get that

\begin{equation}
\log_2\varrho_2(\epsilon) = 0.314\log_2\epsilon + 0.196
\end{equation}

with an $R^2$ value of 0.986. This indicates that $\mathfrak{u}_2\approx 0.314$, and since this is less than one, some amount of geometrical sensitivity to initial conditions does indeed exist in $S_2'$. Recalling from Equation <ref> that an uncertainty exponent is related to the fractal dimension of its associated basin boundary, we can say that the fractal dimension $d$ of $\Sigma\cap S_2'$ is

\begin{equation}
\begin{split}
d &= n - \mathfrak{u}_2 \\
&\approx 2 - 0.314 = 1.686
\end{split}
\end{equation}

Observing the visualization of $S_2'$ in Figure <ref>, we can see that there are some black regions that have boundaries that appear to be smooth, meaning that these particular subsets of the basin boundary have fractal dimensions close to 1. This accounts for the fractal dimension of $\Sigma\cap S_2'$ not being exceptionally close to 2: these boundaries that appear smooth make the overall basin boundary in this slice less “rough.”

Now, let us take a look at the function $\varrho_4(\epsilon)$. Running a linear regression on the values $\log_2\varrho_4(\epsilon)$ vs. $\log_2\epsilon$ from Table <ref> gives us

\begin{equation}
\log_2\varrho_4(\epsilon) = 0.037\log_2\epsilon - 1.341
\end{equation}

with an $R^2$ value of 0.967. This indicates $\mathfrak{u}_4\approx 0.037$, which is exceedingly small and the most extreme example of unpredictability emerging from basin geometry that we have encountered so far in this paper. To put this in perspective, the fact that $\varrho_4(\epsilon)$ is proportional to $\epsilon^{\mathfrak{u}_4}$ means that to reduce the uncertainty in which attractor an initial state in $S_4'$ will end up in by a factor of 10, we will need to reduce our initial uncertainty $\epsilon$ by a factor on the order of $10^{27}$. Using Equation <ref> again, we can calculate that the fractal dimension of the basin boundary $\Sigma\cap S_4'$ is

\begin{equation}
\begin{split}
d &= n - \mathfrak{u}_2 \\
&\approx 4 - 0.037 = 3.963
\end{split}
\end{equation}

which means that this basin boundary is extremely “rough,” behaving similarly to a four-dimensional object even though it divides four-dimensional space. This amount of roughness is analogous to the surface of a human lung ($d\approx 2.97$), just raised one dimension lung. We can see this high fractal dimension justified in Figure <ref>. Even though Figure <ref> doesn't show the four-dimensional set $S_4'$, it is clear that the boundary dividing the seemingly random distribution of white and black points must be close to the dimension of the state space the basins live in. One numerical note is that the $R^2$ value for the linear regression associated with $\varrho_4(\epsilon)$, while certainly telling of a correlation, is less than the $R^2$ value for the regression associated with $\varrho_2(\epsilon)$. We conjecture that this is the case because of the significantly higher uncertainty associated with $\varrho_4(\epsilon)$.

We will conclude this section with a short discussion on how the uncertainty exponents $\mathfrak{u}_2$ and $\mathfrak{u}_4$ change as we expand out from $S_2'$ and $S_4'$. Although the basin classification method we used to classify the white and black basins of this asymmetrically coupled Rulkov neuron system doesn't take into account basin boundaries, we conjecture that $\mathfrak{u}_2$ increases, approaching 1, as we expand the bounds of $x_1$ and $x_2$ away from the set $S_2'$ because the white basin dominates the farther away we go from the attractors. Additionally, we conjecture that $\mathfrak{u}_4$ stays relatively constant as we expand the bounds of $S_4'$ because the white and black basins are both Class 2, so we suspect that they stay similarly overlapped with each other. These conjectures are supported by running the code in Appendix <ref> for different ranges of $x$ and $y$ values.

§.§ Ring Lattice of TEXT Rulkov 1 Neurons

In Section <ref>, we examined three systems of $\zeta=30$ electrically coupled Rulkov 1 neurons in a ring lattice. As a reminder, they are as follows:

* Random fast variable initial states $x_{i,\,0}\in (-1,\,1)$, identical slow variable initial states $y_{i,\,0} = -3.25$, identical $\sigma_i$ values $\boldsymbol{\sigma} = -0.5\cdot\mathbf{1}$, and identical $\alpha_i$ values $\boldsymbol{\alpha} = 4.5\cdot\mathbf{1}$

* Random fast variable initial states $x_{i,\,0}\in (-1,\,1)$, identical slow variable initial states $y_{i,\,0} = -3.25$, random $\sigma_i$ values $\sigma\in(-1.5,\,-0.5)$, and identical $\alpha_i$ values $\boldsymbol{\alpha} = 4.5\cdot\mathbf{1}$

* Random fast variable initial states $x_{i,\,0}\in (-1,\,1)$, identical slow variable initial states $y_{i,\,0} = -3.25$, random $\sigma_i$ values $\sigma\in(-1.5,\,-0.5)$, and identical $\alpha_i$ values $\alpha_i\in(4.25,\,4.75)$

For the specific initial state $\mathbf{X}_0$, $\boldsymbol{\sigma}$ vector, and $\boldsymbol{\alpha}$ vector we use in this paper, see Equations <ref>, <ref>, and <ref>, respectively. As we discovered in Section <ref>, these systems nearly always exhibit chaotic behavior with positive Lyapunov exponents. Therefore, we can conclude that these systems get attracted to some chaotic attractor in 60-dimensional state space. A 60-dimensional space is so far beyond human comprehension that we will not even attempt to present visualizations of these attractors by projecting them onto two-dimensional surfaces. It would certainly be interesting to analyze the geometry of these 60-dimensional attractors by calculating their fractal dimension $d$, but unfortunately, we do not have the computing power to calculate $d$ in 60-dimensional space for a large number of $g^e$ values. Instead, in this section, we will analyze the geometry of these attractors by approximating their Lyapunov dimensions $d_l$, which we know from the Kaplan-Yorke conjecture at the end of Section <ref> gives a ballpark approximation (since we are calculating the Lyapunov spectrum for only a thousand iterations) for the true fractal dimensions of these attractors.

Since we are dealing with 60 Lyapunov exponents, it will be more convenient to calculate Lyapunov dimensions by implementing our method of calculating them from Section <ref> (Equations <ref> and <ref>) into code, which we do in Appendix <ref>.[The function that accomplishes this in the appendix looks a little strange compared to Equations <ref> and <ref> because of the fact that Python indexes from 0 and we index our Lyapunov exponents from 1.] Using this code, we can make graphs similar to the ones in Figures <ref> and <ref>, plotting the values of the Lyapunov dimension $d_l$ for many different values of $g^e$, which we do in Figure <ref>.

First ring lattice system

Second ring lattice system

Third ring lattice system

Graphs of the Lyapunov dimension $d_l$ against the electrical coupling strength $g^e$ for a ring of $\zeta=30$ electrically coupled Rulkov 1 neurons, visualized using the code in Appendix <ref>

There are many interesting observations we can make about the graphs shown in Figure <ref>. First, we can assume that all the chaotic attractors of these three systems for different values of $g^e$ are fractal since our approximations of the Lyapunov dimensions are spread out, not sticking to any defined integers. The only true integer dimensions in these graphs are at the very left of Figure <ref>, where we can see a small number of attractors have dimension 0. This is associated with the non-chaotic attractors we observed in Figure <ref>, which consist of a finite number of zero-dimensional points. Another notable observation is the sizable amount of state space that these attractors take up. Because the state space of this system is so large, we might expect the attractors to take up only a small fraction of state space, but instead, we see that the strange attractors take up a sizeable portion of it for many values of $g^e$, with some of the “roughest” of these attractors taking up close to 45 of the 60 total dimensions.

Comparing Figure <ref> to the graphs of $\lambda_1$ vs. $g^e$ in Figures <ref> and <ref>, we can see that the Lyapunov dimension $d_l$ follows a similar pattern of increasing as we enter the chaotic spiking domain, decreasing as the neurons start to burst in sync with each other, then increasing again as we enter complete chaos. This is to be expected because the Lyapunov dimension is calculated directly from the Lyapunov spectrum $\lambda$. We can also see a similarity in how the $d_l$ and $\lambda_1$ values are distributed across the different systems. Specifically, the $\lambda_1$ values are more erratic and spread out in the first system than they are in the second and third systems, which is also reflected in the $d_l$ values to some degree, observing that the values of $d_l$ in Figure <ref> are more spread out vertically in the synchronized bursting region. However, there are some very clear differences between the trends of the maximal Lyapunov exponent $\lambda_1$ and the Lyapunov dimension $d_l$ as we vary $g^e$. The most obvious difference is in comparing the peaks of the $\lambda_1$ vs. $g^e$ graphs and the $d_l$ vs. $g^e$ graphs, with both peaks in both graphs being associated with chaotic spiking around $g^e=0.1$ and complete chaos around $g^e=1$. In the $\lambda_1$ vs. $g^e$ graphs, the peak in the region of complete chaos is always higher than the peak in the region of chaotic spiking, a fact that is especially obvious in Figure <ref> (the second and third ring systems), where the peaks on the right dwarf the peaks on the left. However, in the graphs of $d_l$ vs. $g^e$, the peaks are similar in height, and in Figure <ref> (the first ring system), the left peak is actually higher than the right peak. This means that, for this system, the chaotic spiking attractor for the systems with a relatively small electrical coupling strength takes up more of state space and is even “rougher” than the attractor for the system of complete chaos $g^e=1$, which is a somewhat surprising result. The main thing to notice from making this comparison, however, is that the maximal Lyapunov exponent $\lambda_1$, the number we use to quantify how chaotic a system is, does not directly correlate to how “big” or fractal an attractor is. For that, as the Kaplan-Yorke conjecture indicates, we need the entire Lyapunov spectrum.

§ DISCUSSION AND CONCLUSION

We begin by reviewing the original aims of our research and our major results. Then, we discuss our suggestions for future research and the possible implications of our research in different fields.

§.§ Summary of Results

The first aim of our research was to analyze and quantify the chaotic dynamics of uncoupled Rulkov neurons. In Rulkov map 1, we accomplished this by expanding on the bifurcation diagram showing the distribution of silence, spiking, and bursting behavior originally presented by rulkov with an approximation of the bifurcation curve $C_{\text{bs}}$ (see Figure <ref>). We extended this bifurcation diagram to give information about the chaotic spiking-bursting that occurs for some parameter values of Rulkov map 1 and to show a quantification of how chaotic uncoupled Rulkov 1 neurons are in the region of parameter space $-2\leq\sigma\leq 0$ and $2\leq\alpha\leq 8$ by utilizing the maximal Lyapunov exponent $\lambda_1$ (see Figure <ref>). We also established a general algorithm for calculating the Lyapunov exponents for any Rulkov 2 neuron and used this to demonstrate the chaotic nature of the spiking and bursting orbits of Rulkov map 2.

The second aim of our research was to model the behavior and analyze the dynamics of complex Rulkov neuron systems. Expanding on the original work by rulkov involving the symmetrical coupling of two Rulkov 1 neurons, we explored a natural extension of this system, namely, a system of two identical spiking Rulkov 1 neurons with an asymmetrical coupling. The complex neuron systems we chose to perform an in-depth analysis of were based around a structure composed of $\zeta$ Rulkov 1 neurons arranged in a ring lattice with current flow in each neighboring connection. We established a method to calculate the $2\zeta$ Lyapunov exponents in a ring lattice system's Lyapunov spectrum through a detailed derivation of the system's $2\zeta\times 2\zeta$ Jacobian matrix, then we applied this method to analyze the dynamics of three different systems composed of $\zeta=30$ Rulkov 1 neurons arranged in a ring. As we varied the electrical coupling strength between the neurons, we discovered patterns in the dynamics of all three of these systems, including the gradual change from chaotic spiking to synchronized bursting to complete chaos.

The third aim of our research was to explore the possibility of the existence of multistability and fractal geometry in the Rulkov maps. This was a success, as we found at least one of these geometrical properties exhibited in three distinct Rulkov neuron systems. The first system where this arose was an uncoupled Rulkov 2 neuron, which we suspected had a chaotic attractor that exhibited fractal geometry because it lacked a resetting mechanism. The second system we found exhibiting these geometrical properties was the system of two asymmetrically electrically coupled Rulkov 1 neurons. From our analysis, we discovered this system could exhibit both non-chaotic and chaotic dynamics depending on the initial conditions of the two neurons. We also suspected that the system's four-dimensional chaotic attractors might be fractal despite the resetting mechanism in both neurons. Finally, the third distinct system that we suspected exhibited fractal geometry was the ring lattice of $\zeta=30$ Rulkov 1 neurons, whose chaotic dynamics we believed likely resulted in complex fractal attractors in 60-dimensional state space. With the three ring lattice systems we analyzed, we were interested in connecting our work with “temporal” chaotic dynamics to the “geometrical,” fractal form of chaos by means of the Kaplan-Yorke conjecture.

The fourth and final aim of our research was to detect, classify, and quantify the geometrical properties associated with chaos and unpredictability in Rulkov map systems. Using the detection methods we established, we first confirmed that multistability, fractal geometry, or both were indeed present in all three distinct systems that we suspected had these properties.

For our first system, the uncoupled Rulkov 2 neuron, we used the computer model from our previous analysis of the neuron's dynamics and the code we used to count boxes on the Hénon attractor to calculate the Minkowski-Bouligand dimensions of attractors generated by a single Rulkov 2 neuron. Our analysis yielded that both chaotic spiking and bursting attractors of Rulkov map 2 are fractal and between one-dimensional and two-dimensional, with the chaotic bursting attractor we analyzed having a higher dimension than the chaotic spiking attractor.

For our second system, the two asymmetrically electrically coupled Rulkov 1 neurons, we found that true multistability didn't exist in the system because after enough time, the neurons eventually synchronized, but a type of quasi-multistability associated with a non-chaotic spiking attractor and chaotic pseudo-attractor did exist. Analyzing the four-dimensional chaotic pseudo-attractor generated by the initial state $\mathbf{X}_0 = \langle-0.56,\,-3.25,\,-1,\,-3.25\rangle$, we found that the pseudo-attractor was indeed fractal with a Minkowski-Bouligand dimension of $d\approx 1.84$. Then, we classified the two basins of the system, looking first at a specific two-dimensional slice of the basins $S_2$, then all of four-dimensional state space $S_4$. We discovered that the two-dimensional basin slice associated with the non-chaotic spiking attractor was Class 3, while the two-dimensional basin slice associated with the chaotic pseudo-attractor was Class 1. Additionally, we discovered that both full four-dimensional basins were Class 2. Finally, we analyzed the basin boundary of the system, considering the intersection of the basin boundary $\Sigma$ with a subset of the two-dimensional slice $S_2'$ and the intersection of $\Sigma$ with a subset of all four-dimensional space $S_4'$. We discovered that the fractal dimension of the basin boundary in $S_2'$ was $d\approx 1.686$ and the fractal dimension of the basin boundary in $S_4'$ was $d\approx 3.963$. This fractal dimension is extremely high, taking up nearly all of four-dimensional state space and indicating high unpredictability regarding which attractor an initial state will end up in. Finally, we made some computer-experiment-based conjectures on the nature of the uncertainty exponents of the system's basin boundary for different subsets of state space.

For our final distinct system, the ring lattice of $\zeta=30$ neurons, we used our Lyapunov spectrum calculation method to analyze the Lyapunov dimensions of the chaotic attractors of our three ring lattice systems. According to the Kaplan-Yorke conjecture, these Lyapunov dimensions should approximate the attractors' true fractal dimensions. We found that all the chaotic attractors of the three systems were fractal and that for some electrical coupling strength values, the attractors took up significant portions of 60-dimensional state space. Comparing the Lyapunov dimensions of the ring lattice systems to their maximal Lyapunov exponents, we also found that while the Lyapunov dimensions followed a similar pattern of increasing and decreasing as we varied the electrical coupling strength, the two quantities were not directly associated with each other. This indicated that the “temporal” chaos measurement of the maximal Lyapunov exponent did not directly relate to the geometrical structure of the attractors.

§.§ Future Research and Possible Implications

The research we detail in this paper lends itself to a significant amount of future research that can improve upon and extend our analysis. To start, relating to Section <ref>, a more rigorous analysis can be done on the distribution of non-chaotic and chaotic dynamics in the parameter space of Rulkov map 1. Specifically, we suggest further analysis be done in dividing parameter space into more than non-chaotic and chaotic regions, perhaps finding a rigorous way to differentiate between the different regimes of behavior: silence, non-chaotic spiking, non-chaotic bursting, chaotic spiking, and chaotic bursting. We also suggest our research in the dynamics of Rulkov map 1 be applied to Rulkov map 2, especially concerning the distribution of non-chaotic and chaotic dynamics in parameter space.

Our research in ring lattice systems of Rulkov 1 neurons sets a precedent for how the chaotic dynamics of different Rulkov neuron lattices may be analyzed and quantified. Specifically, our detailed calculation of the Jacobian matrix of a ring lattice system can be naturally extended to more complex lattices of Rulkov 1 neurons, such as a mesh, torus, or sphere. With more current connections in these two-dimensional lattices, we suspect that more interesting dynamics may appear. A further extension to this is a system where every neuron has an electrical coupling with every other neuron in the system. Although all-to-all couplings have been studied in the context of a mean field of Rulkov 2 neurons,[See the review by ibarz.] we believe that this has never been done with the more experimentally-applicable electrical coupling of Rulkov neurons. This amount of coupling connections will likely produce even more complex dynamics and hyperchaotic behavior. In summary, our research regarding ring lattice systems detailed in Sections <ref> and <ref> provides another step towards modeling a biological neural network.

Additionally, our research in Rulkov ring lattice systems has possible implications for physics, most evidently in statistical mechanics and condensed matter theory. The chaotic dynamics exhibited in our Rulkov lattice models may be applicable to spin models in statistical mechanics, such as the Ising or XY models, as well as crystal lattice structures in condensed matter theory. Chaos in these neuron lattice models may provide a new perspective on the emergence of chaos in lattice models from these branches of physics.

There is also a significant amount of research within dynamical systems theory that could build on our geometrical analysis of Rulkov map systems. For example, research can be done to support or disprove our conjectures about the uncertainty exponents of our asymmetrically electrically coupled Rulkov 1 neuron system in larger regions of space. This may provide more information on how the basin classification method detailed in Section <ref> may or may not be correlated with the nature of the basin boundaries of a system. Additionally, since we used Lyapunov exponents to calculate the Lyapunov dimensions of our ring lattice systems' attractors, researchers with access to the computing power to box count in 60-dimensional space may be interested in confirming or rejecting our use of the Kaplan-Yorke conjecture in these systems. Finally, we recommend future research be done on the possibility of applying the somewhat recent development by basin-entropy of using basin entropy to measure geometrical sensitivity to initial conditions in these Rulkov systems. To our knowledge, neither basin entropy nor the uncertainty exponents we used in this paper have been used to explore multistability and geometrical sensitivity in these neuron systems before, so we recommend future research be done on how these compare with each other.

Overall, our research provides new insight into how fractal geometry, multistability, and geometrical sensitivity to initial conditions can appear in chaotic neuron systems. This research has large implications for both theoretical and experimental research in the behavior of neuronal models and biological neurons by demonstrating that neuron systems can be sensitive to initial conditions in both the temporal and geometrical sense. We recommend experiments be done with real biological neurons to test our theoretical results of slightly different initial neuron states leading to vastly different eventual behaviors.

§ ACKNOWLEDGEMENTS

We would first like to thank our lab director, Mr. Mark Hannum, for providing beneficial discussions and guidance throughout the research process. We would also like to acknowledge Dr. Paul So, who gave us the idea for this research project and provided useful information that kickstarted our research. As for our fellow students, we are extremely grateful to Arthur Prudius for his extensive review and practical suggestions in all aspects of our research. We would also like to thank Henry Stievater, Isaac Park, and Kohlen Farah for their peer review, as well as Leonard Schrag for his coding assistance. Finally, we would like to thank you, the reader, for taking the time to read our paper.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

journaltitleInternational Journal of Bifurcation and Chaos
titleRiddled basins

titleChaos: An Introduction to Dynamical Systems

Cambridge University Press

titleAn Introduction to Dynamical Systems

University of Groningen

titleThe Hénon map
typeBachelor's Thesis

journaltitleThe American Mathematical Monthly
titleA derivation of n-dimensional spherical coordinates

University of Padua

titleLyapunov theory for discrete time systems

Courier Corporation

titlePrinciples of Statistics

journaltitleReviews of Modern Physics
titleIntroduction to bifurcation theory

journaltitleScientific Reports
titleBasin entropy: a new tool to analyze uncertainty in dynamical systems

journaltitleScientific Reports
titleTesting for basins of Wada

journaltitleReviews of Modern Physics
titleErgodic theory of chaos and strange attractors

journaltitleJournal of Statistical Physics
titleQuantitative universality for a class of nonlinear transformations

World Scientific Publishing Company

titleDifferential Geometry Applied to Dynamical Systems

journaltitlePhysica D: Nonlinear Phenomena
titleMeasuring the strangeness of strange attractors

journaltitlePhysics Letters A
titleFinal state sensitivity: an obstruction to predictability

journaltitlePhysica D: Nonlinear Phenomena
titleStrange attractors that are not chaotic

journaltitleCommunications in Mathematical Physics
titleA two-dimensional mapping with a strange attractor

journaltitlePhysics Reports
titleMap-based models in neuronal dynamics

MIT Press

titleDynamical Systems in Neuroscience: The Geometry of Excitability and Bursting

journaltitleInternational Journal of Bifurcation and Chaos
titleNeural excitability, spiking and bursting

journaltitleIEEE Transactions on Neural Networks
titleSimple model of spiking neurons

journaltitleFunctional Differential Equations and Approximation of Fixed Points
titleChaotic behavior of multidimensional difference equations

journaltitlePhysica D: Nonlinear Phenomena
titleBasins of Wada

titleNeimark-Sacker bifurcation

Pearson Education

titleLinear Algebra and Its Applications

titleAn Introduction to Dynamical Systems and Chaos

Naval Research Laboratory

titleDescribing chaotic systems

University of Washington Press

titleThe Essence of Chaos

journaltitleProceedings of the Royal Society A
titleFractal geometry: what is it, and what does it do?

Times Books

titleThe Fractal Geometry of Nature

titleIntroduction to Discrete Dynamical Systems and Chaos

titleSimple mathematical models with very complicated dynamics

journaltitlePhysica D: Nonlinear Phenomena
titleFractal basin boundaries

titleDynamical systems

journaltitleCommunications in Mathematical Physics
titleOn the concept of attractor

journaltitlePhysica D: Nonlinear Phenomena
titleWada basin boundaries and basin cells

journaltitleThe European Physical Journal Special Topics
titleSynchronization of slow-fast systems

titleAttractor dimensions

titleBasin of attraction

Cambridge University Press

titleChaos in Dynamical Systems

journaltitlePhysica D: Nonlinear Phenomena
titleThe transition to chaotic attractors with riddled basins

titleMultistability in Physical and Living Systems

journaltitleNotices of the American Mathematical Society
titleWhat is... a strange attractor?

journaltitlePhysical Review E
titleModeling of spiking-bursting neural behavior using two-dimensional map

journaltitlePhysical Review Letters
titleRegularization of synchronized chaotic bursts

journaltitleJournal of Difference Equations and Applications
titleOn invariant surfaces and bifurcation of periodic solutions of ordinary differential equations

journaltitleThe Mathematica Journal
titleNumerical calculation of Lyapunov exponents

titleUniversalités et fractales

journaltitleInternational Journal of Bifurcation and Chaos
titleOrigin of chaos in a two-dimensional map modeling spiking-bursting neural activity

titleClassifying and quantifying basins of attraction

Cengage Learning

titleMultivariable Calculus

CRC Press

titleNonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, And Engineering

American Mathematical Society

titleAn Introduction to Measure Theory

journaltitlePhysica D: Nonlinear Phenomena
titleDetermining Lyapunov exponents from a time series

journaltitlePhysics Letters A
titleUnpredictability of the Wada property in the parameter plane

journaltitleScience China Technological Sciences
titleDelayed Hopf bifurcation in time-delayed slow-fast systems

§ NOTATION GUIDE

This appendix contains some of the notation used in this paper for reference. Notation is presented approximately in order of appearance.

2cBasic General Notation

$t$     time

$n$     a number of dimensions

$m$     a dimension

or dimension index

$i,\,j,\,k,\,\kappa,\,p$     indices

$r,\,a,\,b,\,\alpha$     basic parameters

$\sigma$     a scaling factor

$d$     the box-counting dimension

of a geometric object

2cDynamical Systems Fundamentals

$\mathbf{x}$     a generic state vector

$x$     a state of a one-dimensional system

or the component of $\mathbf{x}$

in the 1st dimension

$y$     the component of $\mathbf{x}$

in the 2nd dimension

$\mathbf{x}_0,\,x_0$     the initial state of a system

$\mathbf{x}_t,\,x_t$    the state of a system at time $t$

$\mathbf{x}_s,\,x_s$     a stationary state/fixed point

$\mathbf{x}_p,\,x_p$     a periodic point

$q$     the period of an orbit

$\mathbf{g},\,g$     a differential equation function

$\mathbf{f},\,f$     an iteration function

$\mathbf{f}^t,\,f^t$     $t$ iterations of an iteration function

$O(\mathbf{x})$     an orbit of $\mathbf{x}$

$O^q(\mathbf{x}_p)$     the smallest periodic orbit of $\mathbf{x}_p$

$O^+(\mathbf{x}_0)$     the forward orbit of $\mathbf{x}_0$

$x\e{m}$     the component of $\mathbf{x}$

in the $m$th dimension

$f\e{m}$     the iteration function in

the $m$th dimension of $\mathbf{f}$

$T,\,T'$     transformations

$x',\,y'$     the result of transforming

$x$ and $y$ with $T'$

$\delta x_0$     a perturbation of $x_0$

$\delta x_t$     the evolution of $\delta x_0$ after $t$ steps

2cDynamical Systems Fundamentals (continued)

$\delta \mathbf{x}_0$     a perturbation of $\mathbf{x}_0$

$\delta \mathbf{x}\e{m}_0$     the component of $\delta\mathbf{x}_0$

in the $m$th dimension

$\delta \mathbf{x}_t$     the evolution of $\delta \mathbf{x}_0$ after $t$ steps

$d\mathbf{x}_0$     an infinitesimal perturbation of $\mathbf{x}_0$

$d\mathbf{x}_t$     the evolution of $d\mathbf{x}_0$ after $t$ steps

$\mathbf{u}_0$     a unit vector in the direction of $d\mathbf{x}_0$

2cLyapunov Exponents and Related

$\lambda$     the Lyapunov exponent

in a one-dimensional system

or the Lyapunov spectrum

in a multi-dimensional system

$\lambda_i$     the $i$th Lyapunov exponent

$J(\mathbf{x})$     the Jacobian matrix of a system

$J(\mathbf{x}_t)$     the Jacobian matrix of a system at $\mathbf{x}_t$

$J_{mi}(\mathbf{x}_k)$     the $m$th row and $i$th

column of the matrix $J(\mathbf{x}_k)$

$J^t$     the matrix product

$J(\mathbf{x}_{t-1})J(\mathbf{x}_{t-2})\hdots J(\mathbf{x}_0)$

$p_i(t)$     the $i$th principal axis

of an $n$-dimensional ellipsoid

of perturbations at time $t$

$\mu_i$     the $i$th eigenvalue of the matrix $J^{t\intercal}J^t$

$\mathbf{w}_i$     the eigenvector associated with $\mu_i$

$\nu_i$     the $i$th eigenvalue of the matrix $J(\mathbf{x}_s)$

$Q^{(1)}R^{(1)}$     the QR factorization matrices of $J(\mathbf{x}_0)$

$J^*_k$     the matrix product $J(\mathbf{x}_{k-1})Q^{k-1}$,

defined recursively from $Q^{(1)}$

$Q^{(k)}R^{(k)}$     the QR factorization matrices of $J^*_k$

$\Psi(t)$     the matrix product

$R^{(t)}R^{(t-1)}\hdots R^{(1)}$

$\psi_{ij}(t)$     the entry in the $i$th row

and $j$th column of $\Psi(t)$

$r^{(k)}_{ij}$     the entry in the $i$th row

and $j$th column of $R^{(k)}$

$d_l$     the Lyapunov dimension

2cAttractors and Basins

$A,\,C$     attractors

$\mathcal{A},\,\mathcal{C}$     Milnor attractors

$\hat{A}$     the basin of attraction of $A$

2cAttractors and Basins (continued)

$\hat{\mathcal{A}}$     the basin of attraction of $\mathcal{A}$,

which could be riddled

$A'$     a proper subset of $A$

$\mathcal{A}'$     a proper subset of $\mathcal{A}$

$U$     an open set of states that $A$ attracts

$\sigma_A$     the standard deviation of $A$

$\xi$     normalized distance from an attractor

$S(\xi)$     the set of all $\mathbf{x}$

that lie in an $n$-dimensional ball

of radius $\xi$ centered at $\langle A\rangle$

$\hat{A}(\xi)$     the intersection of $\hat{A}$ and $S(\xi)$

$P(\xi)$     the probability that an initial state

$\mathbf{x}_0\in S(\xi)$ is also in $\hat{A}$

$\Delta S(2^k)$     the set of all $\mathbf{x}$ that lie in an

$n$-dimensional shell with inner radius

$\xi=2^k$ and outer radius $\xi=2^{k+1}$

centered at $\langle A\rangle$

$\Delta \hat{A}(2^k)$     the intersection of $\hat{A}$ and $\Delta S(2^k)$

$\Delta P(2^k)$     the probability that an initial state

$\mathbf{x}_0\in\Delta S(2^k)$ is also in $\hat{A}$

$P_0,\,\gamma$     parameters for

basin classification

$\xi_0$     the linear, normalized size

of a Class 4 basin

$\Sigma$     a basin boundary

$\varrho(\epsilon)$     the fraction of states

in a given region of state space

that lie at most a distance $\epsilon$

from a basin boundary

$\mathfrak{u}$     the uncertainty exponent

$w$     the number of attractors

in a system with more

than 3 attractors

$G$     a grid, or set of boxes

$s$     the number of boxes

on one side of a grid

$\Box_i$     a box, or an element of a grid $G$

$C(\mathbf{x})$     the basin $\mathbf{x}$ belongs to

$C(\Box_i)$     the basin that the

center point of $\Box_i$ belongs to

$b(\Box_i)$     the set of boxes consisting of $\Box_i$

and all the boxes sharing at least

one boundary point with $\Box_i$

$K(\Box_i)$     the number of distinct

colors $C(\Box)$ in $b(\Box_i)$

2cAttractors and Basins (continued)

$G_k^p$     the set of all $\Box_i$

for which $K(\Box_i)=k$

immediately after step $p$

$W_k$     the fraction of boxes on or

sufficiently close to a boundary

shared by $k$ basins

$S_2$     a specific, infinite,

two-dimensional slice

of four-dimensional state space

$S_2'$     a specific subset of $S_2$

$S_4$     all of four-dimensional state space

$S_4'$     a specific subset of $S_4$

$A_4$     a four-dimensional attractor

$\hat{A}_4$     the basin of $A_4$

$\langle A_4\rangle$     the mean of $A_4$

$\langle A_2\rangle$     the effective mean of $A_4$

used when classifying the

intersection of $\hat{A}_4$ and $S_2$

$\langle x_1\rangle,\,\langle y_1\rangle,$     slow and fast variable

$\langle x_2\rangle,\,\langle y_2\rangle$     mean values on $A_4$

$\sigma_{A4}$     the standard deviation of $A_4$

$\sigma_{A2}$     the effective standard

deviation of $A_4$

used when classifying the

intersection of $\hat{A}_4$ and $S_2$

$\xi_4$     normalized four-dimensional

distance from $\langle A_4\rangle$

$\xi_2$     normalized two-dimensional

distance from $\langle A_2\rangle$

$P_w(\xi)$     white basin classification

probability function

$P_b(\xi)$     black basin classification

probability function

$\mathfrak{u}_2$     the uncertainty exponent of the

intersection of $\Sigma$ and $S_2'$

$\varrho_2(\epsilon)$     the probability function

associated with $\mathfrak{u}_2$

$\mathfrak{u}_4$     the uncertainty exponent of the

intersection of $\Sigma$ and $S_4'$

$\varrho_4(\epsilon)$     the probability function

associated with $\mathfrak{u}_4$

2cSlow-Fast Systems and Neuron Basics

$\eta$     a very small parameter

$\chi,\,\omega$     continuous slow-fast functions

$f,\,g$     discrete slow-fast functions

$E_{\text{ion}}$     the equilibrium

potential of an ion

$I_{\text{ion}}$     the current of an ion

$g_{\text{ion}}$     the conductance of an ion

$\overline{g}$     the maximum

conductance of an ion

$C$     the capacitance of an axon

$V,\,n,\,m,\,h$     variables of the

Hodgkin-Huxley model

$\alpha_n,\,\beta_n,\,\alpha_m,$     functions of the

$\beta_m,\,\alpha_h,\,\beta_h$     Hodgkin-Huxley model

$v$     fast voltage variable

of the Izhikevich model

$u$     slow recovery variable

of the Izhikevich model

$a,\,b,\,c,\,d$     parameters of the

Izhikevich model

2cRulkov Basics and Current Injection

$x$     the fast or voltage variable

of the Rulkov maps

$y$     the slow variable

of the Rulkov maps

$\alpha,\,\sigma,\,\eta$     parameters of the Rulkov maps

$f_1$     the fast variable iteration

function of Rulkov map 1

$\mathbf{f}_1$     the iteration function

of Rulkov map 1

$f_2$     the fast variable iteration

function of Rulkov map 2

$\mathbf{f}_2$     the iteration function

of Rulkov map 2

$x_{s,\,\text{slow}}$     the value of $x$

that leaves $y$ fixed

$x_{s,\,\text{fast}}$     the fixed points of the

fast map of Rulkov map 1

$x_{s,\,\text{fast, stable}}$     the stable fixed point of the

fast map of Rulkov map 1

$x_{s,\,\text{fast, unstable}}$     the unstable fixed point of the

fast map of Rulkov map 1

2cRulkov Basics and Current Injection (continued)

$B_{\text{stable}}$     the stable branch

of the Rulkov maps

$B_{\text{unstable}}$     the unstable branch

of the Rulkov maps

$B_{\text{spikes}}$     the spiking branch

of the Rulkov maps

$\sigma_{\mathrm{th}}$     the threshold of excitation

of Rulkov map 1

$C_{\text{bs}}$     the bifurcation curve between

spiking and bursting behavior

in Rulkov map 1

$\sigma_{\text{n-s bif}}$     the Neimark-Sacker

bifurcation curve

of Rulkov map 1

$x_{s,\,1},\,x_{s,\,2},\,x_{s,\,3}$     the fixed points of the

fast map of Rulkov map 2

$I_k$     time-varying injected

direct current

$\beta_k,\,\sigma_k$     parameters for modeling

a time-varying current

$\beta^c,\,\sigma^c$     coefficients for DC response

2cNeuron Couplings

$\mathbf{x}_i$     the $i$th coupled neuron

$\mathbf{x}_{i,\,k}$     the state of the $i$th

coupled neuron at time $t=k$

$x_i$     the fast variable of

the $ith$ coupled neuron

$y_i$     the slow variable of

the $ith$ coupled neuron

$x_{i,\,k}$     the fast variable state of the

$i$th couple neuron at time $t=k$

$y_{i,\,k}$     the slow variable state of the

$i$th coupled neuron at time $t=k$

$\mathbf{X}$     the state vector of all the

coupled neurons in a system

$\mathbf{X}_k$     the state of a coupled

neuron system at time $t=k$

$\mathbf{F}$     the iteration function of

a coupled neuron system

$X\e{j}$     the $j$th dimension of $\mathbf{X}$

$F\e{m}$     the $m$th dimension of $\mathbf{F}$

$\mathfrak{C}_{i,\,x},\,\mathfrak{C}_{i,\,y}$     the coupling parameters of

the $i$th coupled neuron

2cNeuron Couplings (continued)

$\alpha_i,\,\sigma_i$     the Rulkov map parameters of

the $ith$ coupled neuron

$\beta^c_i,\,\sigma^c_i$     the coefficients for

DC response of

the $i$th coupled neuron

$g^e$     electrical coupling strength

or coupling conductance

$g^e_{i}$     the electrical coupling strength

of neuron $\mathbf{x}_i$

$g^e_{ji}$     the electrical coupling strength

from neuron $\mathbf{x}_j$ to neuron $\mathbf{x}_i$

$\mathcal{N}_i$     the set of neurons

adjacent to $\mathbf{x}_i$

$\zeta$     the number of neurons

in a coupled neuron lattice system

$J(\mathbf{X})$     the Jacobian matrix of a

coupled neuron system

$J_{mj}(\mathbf{X})$     the $m$th row and $j$th column

of the matrix $J(\mathbf{X})$

$J_{\text{dg},\,a},\,J_{\text{odg},\,b},$     partition matrices of the

$J_{\text{odg},\,c},\,J_{\text{dg},\,d}$     Jacobian matrix for two

coupled Rulkov 1 neurons

$\boldsymbol{\sigma}$     the vector of all $\sigma_i$ values

$\boldsymbol{\alpha}$     the vector of all $\alpha_i$ values

2cBasic Mathematical Notation

$|x|$     the absolute value of $x$

$z$     a complex number

$|z|,\,r$     the modulus of $z$

$\Arg(z),\,\varphi$     the argument of $z$

$z^*$     the complex conjugate of $z$

$|\mathbf{x}|$     the magnitude of $\mathbf{x}$

$\phi,\,\theta,\,\theta_i$     angles describing

spherical coordinates in $\mathbb{R}^n$

$M$     a matrix

$M\transpose$     the transpose of $M$

$I$     the identity matrix

$Q$     an orthogonal matrix

$R$     an upper triangular matrix

$\det M,\,\Delta$     the determinant of $M$

$\tr M,\,\tau$     the trace of $M$

$\frac{df}{dx}\big|_{x=x_t}\,,\,f'(x_t)$     the derivative of $f$ with

respect to $x$ evaluated at $x_t$

2cBasic Mathematical Notation (continued)

$\dot{\mathbf{x}}$     the derivative of $\mathbf{x}$

with respect to $t$

$\frac{\partial f\e{m}}{\partial x\e{i}}$     the partial derivative

of $f\e{m}$ with respect to $x\e{i}$

$\sum_{i=a}^b$     the sum from $i=a$ to $i=b$

$\prod_{i=a}^b$     the product from $i=a$ to $i=b$

$\sum_{i\in S}$     the sum over all indices $i$ in $S$

$\sim$     is proportional to

$\mathbf{1}$     a vector of all ones

2cSet Notation

$S$     a set

$\{x_0,\,x_1,\,\hdots\,,\,x_t\}$     the set with elements $x_i$

for $0\leq i\leq t$

$(a,\,b)$     the open interval from $a$ to $b$

or an ordered pair

representing a point

$\emptyset$     the empty set

$\mathbb{R}$     the set of real numbers

$\mathbb{N}$     the set of natural numbers

$\mathbb{C}$     the set of complex numbers

$\mathbb{R}^n$     the set of all vectors of real

numbers that have length $n$

$\in$     is an element of

$\subseteq$     is a subset of

$\subset$     is a proper subset of

$:$     such that

$S_1\cup S_2$     the union of $S_1$ and $S_2$

$\bigcup_{i=a}^b S_i$     the union set

$S_a\cup S_{a+1}\cup\hdots\cup S_b$

$S_1\cap S_2$     the intersection of $S_1$ and $S_2$

$S_1\times S_2$     the Cartesian

product of $S_1$ and $S_2$

$\mu(S)$     the measure of $S$

$\mu^n(S)$     the measure of $S\subset\mathbb{R}^n$

$\langle S\rangle$     the average of

all the elements of $S$

$|S|$     the cardinality of $S$

$N(\epsilon)$     the number of $n$-dimensional

boxes of side length $\epsilon$

needed to cover a set $S\subset\mathbb{R}^n$

§ DERIVATIONS AND PROOFS

This appendix contains derivations, proofs, and details of some equations, criteria, and methods utilized in this paper.

§.§ The One-Dimensional Lyapunov Exponent

From Equation <ref>,

\begin{equation}
\begin{split}
\lambda &= \lim_{t\to\infty}\lim_{\delta x_0\to 0}\frac{1}{t}\ln\left|\frac{f^t(x_0+\delta x_0)-f^t(x_0)}{\delta x_0}\right| \\
&= \lim_{t\to\infty}\frac{1}{t}\ln\left|\frac{df^t(x)}{dx}\bigg|_{x=x_0}\right| \\
\end{split}
\end{equation}

by the limit definition of a derivative. Using the chain rule, this can be rewritten as

\begin{equation}
\begin{split}
\lambda &= \lim_{t\to\infty}\frac{1}{t}\ln\Bigg|\frac{df^t(x)}{df^{t-1}(x)}\bigg|_{f^{t-1}(x)=x_{t-1}} \\
&\mathrel{\phantom{=}} \frac{df^{t-1}(x)}{df^{t-2}(x)}\bigg|_{f^{t-2}(x)=x_{t-2}}\hdots \\
&\mathrel{\phantom{=}} \frac{df^2(x)}{df(x)}\bigg|_{f(x)=x_1}\frac{df(x)}{dx}\bigg|_{x=x_0}\Bigg|
\end{split}
\label{eq:chain_rule_step_1dlyap_deriv}
\end{equation}

Since $f$ is iterative, we can say

\begin{equation}
\begin{split}
\lambda &= \lim_{t\to\infty}\frac{1}{t}\ln\Bigg|\frac{df(x)}{dx}\bigg|_{x=x_{t-1}}\frac{df(x)}{dx}\bigg|_{x=x_{t-2}} \\
\end{split}
\end{equation}

Condensing and using the product property of logarithms,

\begin{equation}
\begin{split}
\lambda &= \lim_{t\to\infty}\frac{1}{t}\ln\left|\,\prod_{i=0}^{t-1}\frac{df(x)}{dx}\bigg|_{x=x_i}\right| \\
&= \lim_{t\to\infty}\frac{1}{t}\ln\left(\prod_{i=0}^{t-1}\left|\frac{df(x)}{dx}\bigg|_{x=x_i}\right|\right) \\
&= \lim_{t\to\infty}\frac{1}{t}\sum_{i=0}^{t-1}\ln\left|\frac{df(x)}{dx}\bigg|_{x=x_i}\right| \\
&= \lim_{t\to\infty}\frac{1}{t}\sum_{i=0}^{t-1}\ln|f'(x_i)|
\end{split}
\end{equation}

This completes the derivation of Equation <ref>.

§.§ The Lyapunov Spectrum

Rewriting Equation <ref> for an infinitesimal perturbation at a general time step $k$,

\begin{equation}
d\mathbf{x}_k = \begin{pmatrix}
dx\e{1}_k \\[4pt]
dx\e{2}_k \\[1pt]
\vdots \\[1pt]
\end{pmatrix}
\end{equation}

For a general dimension $m$, Equation <ref> reads

\begin{equation}
x\e{m}_{k+1} = f\e{m}(x\e{1}_k,\, x\e{2}_k,\,\hdots\,,\, x\e{n}_k)
\end{equation}

From calculus [p. 944]multi, we know that

\begin{equation}
dx\e{m}_{k+1} = \frac{\partial f\e{m}}{\partial x\e{1}_k}dx\e{1}_k + \frac{\partial f\e{m}}{\partial x\e{2}_k}dx\e{2}_k + \hdots + \frac{\partial f\e{m}}{\partial x\e{n}_k}dx\e{n}_k
\end{equation}

\begin{equation}
dx\e{m}_{k+1} = \sum_{i=1}^n\frac{\partial f\e{m}}{\partial x\e{i}_k}dx\e{i}_k
\label{eq:appendixmulti}
\end{equation}

We can rewrite Equation <ref> for $J(\mathbf{x}_k)$ concisely as

\begin{equation}
J_{mi}(\mathbf{x}_k) = \frac{\partial f\e{m}}{\partial x\e{i}_k}
\end{equation}

Substituting into Equation <ref>,

\begin{equation}
dx\e{m}_{k+1} = \sum_{i=1}^nJ_{mi}(\mathbf{x}_k)dx\e{i}_k
\end{equation}

Expanding into the full vector $d\mathbf{x}_{k+1}$,

\begin{align}
\begin{split}
\begin{pmatrix}
dx\e{1}_{k+1} \\[4pt]
dx\e{2}_{k+1} \\[1pt]
\vdots \\[1pt]
\end{pmatrix}
&= \begin{pmatrix}
\sum_{i=1}^nJ_{1i}(\mathbf{x}_k)dx\e{i}_k \\[4pt]
\sum_{i=1}^nJ_{2i}(\mathbf{x}_k)dx\e{i}_k \\[1pt]
\vdots \\[1pt]
\sum_{i=1}^nJ_{ni}(\mathbf{x}_k)dx\e{i}_k
\end{pmatrix} \\
&= J(\mathbf{x}_k)\begin{pmatrix}
dx\e{1}_{k} \\[4pt]
dx\e{2}_{k} \\[1pt]
\vdots \\[1pt]
\end{pmatrix}
\end{split}
\end{align}

by the nature of matrix multiplication. Concisely, this equation yields Equation <ref>.

Following Equation <ref>, consider a generic vector $\mathbf{v}$. It is easy to see that

\begin{align}
\begin{split}
\mathbf{v}\transpose\mathbf{v} &=
\begin{pmatrix}
v_1 & v_2 & ... & v_n
\end{pmatrix}
\begin{pmatrix}
v_1 \\
v_2 \\
\vdots \\
\end{pmatrix} \\
&= v_1^2 + v_2^2 + \hdots + v_n^2 = |\mathbf{v}|^2
\end{split}
\label{eq:vtransposev}
\end{align}

where $M\transpose$ is the transpose of $M$. Using this,

\begin{align}
\begin{split}
|J^t\mathbf{u}_0|^2 &= (J^t\mathbf{u}_0)\transpose J^t\mathbf{u}_0 \\
&= \mathbf{u}_0\transpose (J^{t\intercal}J^t)\mathbf{u}_0
\end{split}
\end{align}

because $(MN)\transpose = N\transpose M\transpose$ [p. 101]linear. We can now use this to rewrite Equation <ref>:

\begin{align}
\begin{split}
\lambda &= \lim_{t\to\infty}\frac{1}{t}\ln\left|J^t\mathbf{u}_0\right| \\
&= \lim_{t\to\infty}\frac{1}{t}\ln\left[\mathbf{u}_0\transpose (J^{t\intercal}J^t)\mathbf{u}_0\right]^{1/2} \\
&= \lim_{t\to\infty}\frac{1}{2t}\ln\left[\mathbf{u}_0\transpose (J^{t\intercal}J^t)\mathbf{u}_0\right]
\end{split}
\label{eq:appendix-derivation-almost-done}
\end{align}

The matrix $J^{t\intercal}J^t$ is clearly symmetric:[This is why we took the time to transform Equation <ref>. In general, $J^t$ is neither symmetric nor diagonalizable, so it is much more convenient to work with $J^{t\intercal}J^t$, which, in addition to being symmetric and diagonalizable, has a set of $n$ orthonormal eigenvectors corresponding to real and positive eigenvalues.]

\begin{equation}
(J^{t\intercal}J^t)\transpose = J^{t\intercal}(J^{t\intercal})\transpose = J^{t\intercal}J^t
\end{equation}

Therefore, from linear algebra, we know that $J^{t\intercal}J^t$ has $n$ real eigenvalues (accounting for multiplicities) [p. 397]linear. Labelling these eigenvalues $\mu_i$ and their associated normalized eigenvectors $\mathbf{w}_i$ with $i = 1,\,2,\,\hdots\,,\,n$, we can take $\mathbf{u}_0$ to be a normalized eigenvector $\mathbf{w}_i$. Then, it follows from Equation <ref> that

\begin{align}
\begin{split}
\lambda_i &= \lim_{t\to\infty}\frac{1}{2t}\ln\left[\mathbf{w}_i\transpose (J^{t\intercal}J^t)\mathbf{w}_i\right] \\
&= \lim_{t\to\infty}\frac{1}{2t}\ln\left[\mathbf{w}_i\transpose \mu_i\mathbf{w}_i\right] \\
&= \lim_{t\to\infty}\frac{1}{2t}\ln\left[\mu_i\mathbf{w}_i\transpose\mathbf{w}_i\right] \\
&= \lim_{t\to\infty}\frac{1}{2t}\ln\mu_i
\end{split}
\end{align}

by the definition of eigenvalues and eigenvectors. By labeling the eigenvalues $\mu_i$ and the normalized eigenvectors $\mathbf{w}_i$ with the appropriate indices to satisfy $\lambda_1\geq\lambda_2\geq\hdots\geq\lambda_n$, this completes the derivation of Equation <ref>.

§.§ The Criteria for the Attractiveness of Fixed Points

For a possible fixed point attractor $A=\{\mathbf{x}_s\}$, Property 1 of attractors (see Section <ref>) is automatically satisfied by Equation <ref>. Furthermore, if Properties 1 and 2 are satisfied, Property 3 is satisfied because there is no proper subset of $A$ besides the empty set. Therefore, to determine whether or not a fixed point is an attractor, Property 2 is the only property we need to consider.

First, let us consider a very small neighborhood of $A$ as our open set of initial conditions $U$ that $A$ attracts. For some $\mathbf{x}_0\in U$, Property 2 of attractors says that if $A$ is an attractor, the distance from $\mathbf{f}^t(\mathbf{x}_0)=\mathbf{x}_t$ to $\mathbf{x}_s$ goes to 0 as $t$ goes to infinity, or $\lim_{t\to\infty}|\mathbf{x}_t-\mathbf{x}_s|=0$. Since $U$ is very small, the perturbation $\mathbf{x}_t-\mathbf{x}_s$ follows Equation <ref>:

\begin{equation}
\mathbf{x}_t-\mathbf{x}_s = J^td\mathbf{x}_s
\label{eq:jacobian-alltheway-fixed}
\end{equation}

where $d\mathbf{x}_s = \mathbf{x}_0-\mathbf{x}_s$. By Equation <ref>,

\begin{equation}
J^t = J\left(\mathbf{f}^{t-1}(\mathbf{x}_{s})\right)\,J\left(\mathbf{f}^{t-2}(\mathbf{x}_{s})\right)\,\hdots\, J\left(\mathbf{x}_s\right)
\end{equation}

By Equation <ref>, this collapses to

\begin{equation}
J^t = [J(\mathbf{x}_s)]^t
\label{eq:jtfixed}
\end{equation}

Now, we consider just the first iteration of $d\mathbf{x}_s$:

\begin{equation}
\mathbf{x}_1-\mathbf{x}_s = J(\mathbf{x}_s)d\mathbf{x}_s
\end{equation}

If we choose $d\mathbf{x}_s$ to be in the direction of an eigenvector of $J(\mathbf{x}_s)$ and $\nu_i$ to be its associated eigenvalue,[Because the Jacobian matrix $J$ doesn't have the nice properties of $J^{t\intercal}J^t$, which is symmetric and diagonalizable, the eigenvalues of $J(\mathbf{x}_s)$ might be complex. Therefore, in order to choose $d\mathbf{x}_s$ to be in the direction of an eigenvector of $J(\mathbf{x}_s)$, the entries of $d\mathbf{x}_s$ may need to be complex. This is a departure from our established intuition of the state space of dynamical systems, but it is still mathematically sound. For complex entries of $d\mathbf{x}_s$ and complex values of $\nu_i$, the absolute value signs should be interpreted as modulus functions.] this reads

\begin{equation}
\mathbf{x}_1-\mathbf{x}_s = \nu_i\,d\mathbf{x}_s
\end{equation}

By Equation <ref>, we can now rewrite Equation <ref> as

\begin{equation}
\mathbf{x}_t-\mathbf{x}_s = \nu_i^t\,d\mathbf{x}_s
\end{equation}

Taking the limit of the magnitude on both sides,

\begin{equation}
\lim_{t\to\infty}|\mathbf{x}_t-\mathbf{x}_s| = \lim_{t\to\infty}|\nu_i^t\,d\mathbf{x}_s| = |d\mathbf{x}_s|\lim_{t\to\infty}|\nu_i|^t
\end{equation}

For $|\nu_i|<1$,

\begin{equation}
\lim_{t\to\infty}|\mathbf{x}_t-\mathbf{x}_s| = |d\mathbf{x}_s|\lim_{t\to\infty}|\nu_i|^t = 0
\label{eq:eigenvalue-attractor}
\end{equation}

For $|\nu_i|>1$,

\begin{equation}
\lim_{t\to\infty}|\mathbf{x}_t-\mathbf{x}_s| = |d\mathbf{x}_s|\lim_{t\to\infty}|\nu_i|^t\to\infty
\label{eq:eigenvalue-repeller}
\end{equation}

Now, let us consider the matrix $J^{t\intercal}J^t$ at $\mathbf{x}_s$. Then,

\begin{equation}
J^{t\intercal}J^t = [J(\mathbf{x}_s)\transpose]^t[J(\mathbf{x}_s)]^t
\end{equation}

The eigenvalues $\mu_i$ of $J^{t\intercal}J^t$ are therefore

\begin{equation}
\mu_i = \nu_i^{2t}
\end{equation}

since $J$ and $J\transpose$ share the same characteristic polynomial. Equation <ref> then reads

\begin{equation}
\begin{split}
\lambda_i &= \lim_{t\to\infty}\frac{1}{2t}\ln\mu_i \\
&= \lim_{t\to\infty}\frac{1}{2t}\ln\nu_i^{2t} \\
&= \lim_{t\to\infty}\ln\nu_i
\end{split}
\label{eq:lyapunov-jacobian-eigenvalues}
\end{equation}

Thus, the eigenvalues of the Jacobian at $\mathbf{x}_s$ are directly related to the system's Lyapunov exponents. We know from Section <ref> that Lyapunov exponents determine the maximum and minimum growth of an infinitesimal $(n-1)$-dimensional sphere of perturbations from an initial condition, so $\nu_1$ and $\nu_n$ provide an upper and lower bound on the growth of $d\mathbf{x}_s$ in any direction.

Therefore, if $|\nu_i|<1$ for all $i = 1,\,2,\,\hdots\,,\,n$, Equation <ref> says that $\lim_{t\to\infty}|\mathbf{x}_t-\mathbf{x}_s|=0$ for all $d\mathbf{x}_s$, meaning Property 2 is satisfied and $\mathbf{x}_s$ is an attractor. Similarly, if every $|\nu_i|>1$, then Equation <ref> says that $\lim_{t\to\infty}|\mathbf{x}_t-\mathbf{x}_s|\to\infty$ for all $d\mathbf{x}_s$, meaning $\mathbf{x}_s$ is a repeller. If at least one $|\nu_i|<1$ and at least one $|\nu_i|>1$, then $\lim_{t\to\infty}|\mathbf{x}_t-\mathbf{x}_s|=0$ for some $d\mathbf{x}_s$ and $\lim_{t\to\infty}|\mathbf{x}_t-\mathbf{x}_s|\to\infty$ for others, meaning $\mathbf{x}_s$ is a saddle point. This completes the proof of the criteria for the attractiveness of fixed points. It is worthwhile to note that these criteria imply that if $\mathbf{x}_s$ is an attractor, all of the Lyapunov exponents in the Lyapunov spectrum of the forward orbit $O^+(\mathbf{x}_0)$ are negative by Equation <ref>. Similarly, if $\mathbf{x}_s$ is a repeller, all of the Lyapunov exponents of $O^+(\mathbf{x}_0)$ are positive. Finally, if $\mathbf{x}_s$ is a saddle points, some $\lambda_i$ are positive and some are negative.

For a one-dimensional system $x_{k+1} = f(x_k)$ with a fixed point $x_s$, the Jacobian matrix $J(x_s)$ is

\begin{equation}
J(x_s) = \left(\frac{df}{dx}\bigg|_{x=x_s}\right) = \big(f'(x_s)\big)
\end{equation}

The eigenvalue $\nu_1$ of $\big(f'(x_s)\big)$ is obviously $f'(x_s)$, so the criteria for the attractiveness of fixed points say that if $|f'(x_s)|<1$, $x_s$ is an attractor, and if $|f'(x_s)|>1$, $x_s$ is a repeller. This completes the proof of the criteria for the attractiveness of fixed points in one dimension.

§.§ A QR Factorization Method of Lyapunov Spectrum Calculation

QR factorization is a method of decomposing some matrix $M$ into the of an orthogonal matrix $Q$ and an upper triangular matrix $R$: $M=QR$ [p. 359]linear.[An orthogonal matrix is a square matrix whose column vectors are orthonormal (orthogonal and normalized). An upper triangular matrix is a matrix whose entries are zero below its main diagonal .] In this appendix, we present the algorithm of eckmann for calculating Lyapunov spectrums by taking advantage of QR factorization.[eckmann were the first to propose this method, but we loosely model our notation after the way it is presented in the paper by sandri.]

Recall from Section <ref> that when calculating the Lyapunov spectrum, we are concerned with the matrix

\begin{equation}
J^t = J(\mathbf{x}_{t-1})J(\mathbf{x}_{t-2})\hdots J(\mathbf{x}_0)
\label{eq:j^t_expansion}
\end{equation}

Let us first decompose the matrix $J(\mathbf{x}_0)$ using the QR factorization method, defining[In this case, the superscripts of $Q$ and $R$ should be thought of as an index. This is so that we can reference a given entry of one of these matrices by writing its row and column number in the subscript.]

\begin{equation}
J(\mathbf{x}_0) = Q^{(1)} R^{(1)}
\end{equation}

Now, for $k=2,\,3,\,\hdots\,,\,t$, we can recursively define

\begin{equation}
J(\mathbf{x}_{k-1})Q^{(k-1)} = J_k^*
\label{eq:recursive_jacobian}
\end{equation}

then decompose this matrix:

\begin{equation}
J_k^* = Q^{(k)} R^{(k)}
\label{eq:recursive_jacobian_qr}
\end{equation}

Since an orthogonal matrix $Q$ has the property that $Q\transpose Q = QQ\transpose = I$ [p. 345]linear,[$I$ is the identity matrix.] we can equate Equations <ref> and <ref> and rearrange:

\begin{align}
J(\mathbf{x}_{k-1})Q^{(k-1)} = Q^{(k)} R^{(k)} \\
J(\mathbf{x}_{k-1}) = Q^{(k)} R^{(k)} \left(Q^{(k-1)}\right)\transpose
\end{align}

Substituting into Equation <ref>,

\begin{equation}
\begin{split}
J^t &= Q^{(t)} R^{(t)} \left(Q^{(t-1)}\right)\transpose \\
&\mathrel{\phantom{=}} Q^{(t-1)} R^{(t-1)} \left(Q^{(t-2)}\right)\transpose \\
&\mathrel{\phantom{=}} \hdots Q^{(1)} R^{(1)}
\end{split}
\end{equation}

The adjacent $Q$ matrices cancel:

\begin{equation}
J^t = Q^{(t)} R^{(t)} R^{(t-1)} \hdots R^{(1)}
\label{eq:j^t_qr_factorized}
\end{equation}

Since the product of upper triangular matrices is an upper triangular matrix,[This can be seen explicitly in Equation <ref>.] we can define

\begin{equation}
\Psi(t) = \prod_{k=t}^1 R^{(k)} = \begin{pmatrix}
\psi_{11}(t) & \psi_{12}(t) & \hdots & \psi_{1n}(t) \\
0 & \psi_{22}(t) & \hdots & \psi_{2n}(t) \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \hdots & \psi_{nn}(t)
\end{pmatrix}
\label{eq:big_upsilon}
\end{equation}

where $\psi_{ij}$ is the entry in the $i$th row and the $j$th column of $\Psi(t)$. Substituting this into Equation <ref>,

\begin{equation}
J^t = Q^{(t)}\,\Psi(t)
\end{equation}

Substituting this representation of $J^t$ into Equation <ref>,

\begin{equation}
\lambda = \lim_{t\to\infty}\frac{1}{t}\ln|Q^{(t)}\,\Psi(t)\mathbf{u}_0|
\end{equation}

Since the eigenvalues of a triangular matrix are on its diagonal [p. 271]linear, taking $\mathbf{u}_0$ to be a normalized eigenvector of $\Psi(t)$ yields

\begin{equation}
\lambda_i = \lim_{t\to\infty}\frac{1}{t}\ln|Q^{(t)}\,\psi_{ii}(t)\mathbf{u}_0|
\end{equation}

for $i=1,\,2,\,\hdots\,,\,n$.[To match the indices, we arrange the column vectors of the $Q$ matrices so that the diagonal entries in their associated $R$ matrices are ordered from greatest to least.] Since $\psi_{ii}(t)$ is a constant within the limit, we can pull it out from between $Q^{(t)}$ and $\mathbf{u}_0$:

\begin{equation}
\lambda_i = \lim_{t\to\infty}\frac{1}{t}\ln\left(|\psi_{ii}(t)||Q^{(t)}\mathbf{u}_0|\right)
\label{eq:qr_factorization_step}
\end{equation}

Now, if we consider just the $|Q^{(t)}\mathbf{u}_0|$ factor, using Equation <ref> yields

\begin{equation}
\begin{split}
|Q^{(t)}\mathbf{u}_0|^2 &= \left(Q^{(t)}\mathbf{u}_0\right)\transpose Q^{(t)}\mathbf{u}_0 \\
&= \mathbf{u}_0\transpose \left(Q^{(t)}\right)\transpose Q^{(t)}\mathbf{u}_0 \\
&= \mathbf{u}_0\transpose\mathbf{u}_0 \\
&= 1
\end{split}
\end{equation}

because $Q^{(t)}$ is an orthogonal matrix and $\mathbf{u}_0$ is a unit vector. This shows that $|Q^{(t)}\mathbf{u}_0|=1$, so we can rewrite Equation <ref> as

\begin{equation}
\lambda_i = \lim_{t\to\infty}\frac{1}{t}\ln |\psi_{ii}(t)|
\label{eq:qr-factorization-method-of-lyap-spec-calc}
\end{equation}

This completes the derivation of the QR factorization method of calculating a Lyapunov spectrum.

In practice, when approximating a Lyapunov spectrum computationally using a large value of $t$, the entries of $\Psi(t)$ usually overflow because we are multiplying many $R$ matrices. To counter this, let us consider a simple case, namely, stopping at $t=2$. Then, Equation <ref> reads

\begin{equation}