# Dynamics of the complex rational delay recursive sequence
$\displaystyle{z_{n+1}=\frac{\alpha+\beta z_{n-k}}{\gamma-z_{n}}}$

Sk. Sarif Hassan
Department of Mathematics
University of Petroleum and Energy Studies
Bidholi, Dehradun, India
Email<EMAIL_ADDRESS>

###### Abstract

Dynamics of the delay rational difference equation
$\displaystyle{z_{n+1}=\frac{\alpha+\beta z_{n-k}}{\gamma-z_{n}}}$ with
complex parameters $\alpha$, $\beta$, $\gamma$ and arbitrary complex initial
conditions is investigated. Existence of prime period two solutions and higher
order periods are ensured in the complex parameters unlike in the case of real
parameters of the same rational difference equation. In addition, a new
dynamical behavior, chaotic solutions of the difference equation are ensured
computationally.

Keywords: Rational difference equation, Local asymptotic stability, Chaotic
trajectory and Periodicity.

Mathematics Subject Classification: 39A10 & 39A11

##  1 Introduction and Background

Consider the rational difference equation,

$\displaystyle{z_{n+1}=\frac{\alpha+\beta z_{n-k}}{\gamma-
z_{n}}},n=0,1,2,\ldots$ (1)

where the parameters $\alpha$, $\beta$, $\gamma$ and the initial conditions
$z_{-k},z_{-k+1},\dots z_{-1},z_{0}$ are arbitrary complex numbers.

This (k+1)th order rational difference equation Eq.(1) and closely similar
equations are studied when the parameters $\alpha\geq 0$, $\beta>0$ and
$\gamma>0$ and the initial conditions are real numbers while $k$ is any
natural number [1], [2] and [3]. They have shown that the positive fixed point
of Eq.(1) is a global attractor with a basin that depends on certain
conditions posed on the coefficients $\alpha$, $\beta$ and $\gamma$ [1].

In this present article, an attempt has been made to study the dynamics of the
rational difference equation (1) under the assumption that the parameters and
the initial conditions are arbitrary complex numbers. Similar works can be
found in [4], [5] and [6]. Here some of the important basics are being
discussed in the following [7], [8] and [9]:

Definition 1: A difference equation of order $(k+1)$ is of the form

$\displaystyle{z_{n+1}=f(z_{n},z_{n-1},\dots,z_{n-k})},n=0,1,2,\ldots$ (2)

where $f$ is a continuous function which maps a subset $\mathbb{D}^{k+1}$ into
$\mathbb{D}$ and $\mathbb{D}\subset\mathbb{C}$. A fixed point $\bar{z}$ of the
difference equation Eq. (2) is a point that satisfy the condition
$\bar{z}=f(\bar{z},\bar{z},\dots,\bar{z})$.

Definition 2: Let $\bar{z}$ be a fixed point of the Eq.(2), then $\bar{z}$ is
_locally asymptotically stable_ if for every $\epsilon>0$, there exist a
$\delta(\epsilon)>0$ such that, if $z_{-k},\dots,z_{-1},z_{0}\in$ $D$ with
$\left\lvert z_{-k}-\bar{z}\right\rvert+\dots+\left\lvert
z_{-1}-\bar{z}\right\rvert+\left\lvert
z_{0}-\bar{z}\right\rvert<\delta(\epsilon)$, then $\left\lvert
z_{n}-\bar{z}\right\rvert<\epsilon$ for all $n\geq-k$.

Definition 3: A sequence ${z_{n}}_{n=−k}^{\infty}$ is said to be _periodic
with period $p$_ if $z_{n+p}=z_{n}$ for all $n\geq-k$. A sequence
${z_{n}}_{n=−k}^{\infty}$ is said to be _periodic_ with prime period _p_ if
$p$ is the smallest positive integer having this property.

Definition 4: An open ball $B(a,r)\in\mathbb{C}$ is called an _invariant open
ball_ of Eq.(2) if $z_{-k},\dots,z_{-1},z_{0}\in B(a,r)$ then $z_{n}\in
B(a,r)$ for all $n>0$. That is every solution of Eq.(2) with initial
conditions in $B(a,r)$ remains in $B(a,r)$.

Definition 5: The difference equation Eq.(2) is said to be _permanent and
bounded_ if there exist positive real numbers $M$ and $N$ with $0<M\leq
N<\infty$ such that for any initial conditions $z_{-k},\dots,z_{-1},z_{0}$
there exist a positive integer P which depends on the initial conditions such
that $M\leq\left\lvert z_{n}\right\rvert\leq N$ for all $n\geq P$.

The linearized equation associated with Eq.(2) about the equilibrium point
$\bar{z}$ is

$y_{n+1}=\sum_{i=0}^{k}\frac{\partial
f(\bar{z},\bar{z},\dots,\bar{z})}{\partial u_{i}}y_{n-i}$

Its characteristic equation is

$\lambda^{k+1}=\sum_{i=0}^{k}\frac{\partial
f(\bar{z},\bar{z},\dots,\bar{z})}{\partial u_{i}}\lambda^{n-i}$

where $n=0,1,2,\dots$.

######  Theorem 1.1.

Assume that $f$ is a $C^{1}$-function and let $\bar{z}$ a fixed point of
Eq.(2). Then the following statements are true:

  * •

If all the roots of the characteristic equation lie in the open unit disk
$\left\lvert\lambda\right\rvert<1$, then the fixed point $\bar{z}$ of Eq.(2)
is locally asymptotically stable.

  * •

If at least one root of the characteristic equation has the absolute value
greater than one, then the fixed point $\bar{z}$ of Eq.(2) is unstable.

  * •

If all the roots of the have the characteristic equation absolute value
greater than one, then the fixed point $\bar{z}$ of Eq.(2) is a source.

######  Theorem 1.2.

Assume that $p$, $q$ $\in\mathbb{C}$ and $k\in\mathbb{N}$. Then $\left\lvert
p\right\rvert+\left\lvert q\right\rvert<1$ is a sufficient condition for
asymptotically stability of the difference equation

$z_{n+1}-pz_{n}+qz_{n-k}=0,n=0,1,2,3,\dots$

Now we shall use these basic theorems to explore the local stability of the
fixed points of the Eq.(1).

##  2 Local Asymptotic Stability of the Fixed Points and Boundedness

The fixed points of Eq.(1) are the solutions of the quadratic equation

$\bar{z}=\frac{\alpha+\beta\bar{z}}{\gamma-\bar{z}}$

Eq.(1) has the two fixed points $\bar{z}_{1}$ and $\bar{z}_{2}$ when
$\left\lvert\alpha\right\rvert\neq\left\lvert\frac{({\gamma-\beta})^{2}}{4}\right\rvert$
…
$\frac{1}{2}\left(-\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$
and
$\frac{1}{2}\left(\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$
respectively.

It is noted that if $\alpha=\frac{({\gamma-\beta})^{2}}{4}$, then there is
only one fixed point, $\frac{(\gamma-\beta)}{2}$.

The linearized equation of the rational difference equation Eq.(1) with
respect to the fixed points $\bar{z}$ is

$\displaystyle{z_{n+1}-\frac{\alpha+\beta\bar{z}}{(\gamma-\bar{z})^{2}}z_{n}-\frac{\beta}{\gamma-\bar{z}}z_{n-k}=0,n=0,1,\ldots}$
(3)

with associated characteristic equation

$\lambda^{k+1}-\frac{\alpha+\beta\bar{z}}{(\gamma-\bar{z})^{2}}\lambda^{k}-\frac{\beta}{\gamma-\bar{z}}=0.\\\
$ (4)

The linearized equation for the fixed point $\frac{(\gamma-\beta)}{2}$ is

$\displaystyle{z_{n+1}-\frac{\gamma-\beta}{\gamma+\beta}z_{n}-\frac{2\beta}{\beta+\gamma}z_{n-k}=0,n=0,1,\ldots}$
(5)

with associated characteristic equation

$\lambda^{k+1}-\frac{\gamma-\beta}{\gamma+\beta}\lambda^{k}-\frac{2\beta}{\beta+\gamma}=0.$
(6)

It is found that there does not exist any $\alpha$, $\beta$ and $\gamma$ in
$\mathbb{C}$ such that $\alpha=\frac{({\gamma-\beta})^{2}}{4}$ in
$\mathbb{D}\in\mathbb{C}$ for which the condition

$\left\lvert\frac{\gamma-\beta}{\gamma+\beta}\right\rvert+\left\lvert\frac{2\beta}{\beta+\gamma}\right\rvert<1$

Hence the fixed point $\frac{(\gamma-\beta)}{2}$ is not _locally
asymptotically stable_.

The following result gives the local asymptotic stability of the fixed point …
$\frac{1}{2}\left(-\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$
and
$\frac{1}{2}\left(\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$
of the Eq.(1).

######  Theorem 2.1.

The fixed point
$\frac{1}{2}\left(-\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$
of Eq.(1) is locally asymptotically stable if

$\left\lvert\frac{2\beta}{\sqrt{(\beta-\gamma)^{2}-4\alpha}+\beta+\gamma}\right\rvert+\left\lvert\frac{\gamma\left(-\sqrt{(\beta-\gamma)^{2}-4\alpha}-\beta+\gamma\right)-2\alpha}{2(\alpha+\beta\gamma)}\right\rvert<1$

######  Theorem 2.2.

The fixed point
$\frac{1}{2}\left(\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$
of Eq.(1) is locally asymptotically stable if

$\left\lvert\frac{2\beta}{-\sqrt{(\beta-\gamma)^{2}-4\alpha}+\beta+\gamma}\right\rvert+\left\lvert\frac{\gamma\left(\sqrt{(\beta-\gamma)^{2}-4\alpha}-\beta+\gamma\right)-2\alpha}{2(\alpha+\beta\gamma)}\right\rvert<1$

Proof of these two theorems are straightforward from the result stated in the
Theorem 1.2. Here we go with few examples which illustrate the asymptotic
behavior of these two fixed points.

S.n. |  Parameters $\alpha$, $\beta$, $\gamma$ |  Remark
---|---|---
1 |  $\alpha\to-107.7-374.6i,\beta\to 56-91.2i,\gamma\to 147.4+210i$ |  Here modulus of the zeros of the characteristic equation are both less than one and so the fixed point is _locally asymptotically stable_. It is found that for all $k=1,2,3,\dots$, The trajectories are convergent to the fixed point $\frac{1}{2}\left(-\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$.
2 |  $\alpha\to-297.4-75.3i,\beta\to 78-71.6i,\gamma\to-220.6-128.5i$ |  Here modulus of the zeros of the characteristic equation are less than one and so the fixed point $\frac{1}{2}\left(\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$ is _locally asymptotically stable_. It is found that for all $k=1,2,3,\dots$, The trajectories are convergent to the fixed point.
Table 1: Parameters $\alpha$, $\beta$ and $\gamma$ for which the fixed points
$\frac{1}{2}\left(\pm\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$
of Eq.(1) are _locally asymptotically stable_ for arbitrary initial values.

---

Figure 1: Local Asymptotical Stability of the fixed points. In the top of the
figure, for different values of $k$, complex trajectories plot (left) and
their corresponding time series plot with real and imaginary part (right) are
given for the fixed point
$\frac{1}{2}\left(-\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$.
Same is given in the bottom for the other fixed point.

In the Table-1 as stated, for the parameters given in the serial number $1$
and $2$, the trajectories are convergent and converge to the fixed points
$\frac{1}{2}\left(-\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$
and
$\frac{1}{2}\left(\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$
of Eq.(1) respectively for all delay term $k$.

Further, a matrix of $20$ rows of parameters $\alpha$, $\beta$ and $\gamma$
for which the fixed point
$\frac{1}{2}\left(-\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$
of Eq.(1) is _locally asymptotically stable_ , is given here. In each row the
first, second and third complex number referred as $\alpha$, $\beta$ and
$\gamma$ respectively.

$\left(\begin{array}[]{ccc}2213.8\,-627.2i&-1158.7-591.8i&4355.\,-33.7i\\\
-1755.8-435.3i&1657.2\,-99.3i&3223.4\,-1018.8i\\\
-880.8+1780.7i&-1282.1-813.5i&-555.-3870.8i\\\
2860.1\,+1430i&-1262.7-860i&-643.8+3893.9i\\\
-2501.1+1648.8i&1508.5\,+1652.8i&1773.1\,-4323.1i\\\
-60.8-492i&7.9\,-90.7i&307.5\,+659i\\\
2750.5\,-1175.1i&-953.7+392.9i&3987.6\,+1450.2i\\\
338.9\,+77.1i&-203.4-84.8i&750.7\,+226.7i\\\
-635.5+187.8i&-224.9-840.4i&1347.8\,+1812.5i\\\
18\,+69.1i&1.2\,+45.3i&276.5\,-26.2i\\\
1055.8\,+3679.6i&501.5\,-369.8i&1400.3\,+1747.7i\\\
2337.3\,+4233.9i&-377.3-451.5i&120.7\,+1373.2i\\\
3231.5\,+700.7i&-968.6-1000.5i&3070.4\,-968.5i\\\
2021.8\,+1066.1i&-361.+739i&3616.2\,+1451.6i\\\
-1058.1+1749.3i&357.7\,+4.8i&788.9\,-3083.i\\\
-890.5+1361.9i&-1421.6-670.8i&116.4\,-3391.3i\\\
-3635.7-4440.8i&-1692.9+1093i&3719.2\,-3308.7i\\\
-76.2+2700.8i&-733.+1500i&3174.1\,+1390i\\\
-1883.8-4906.7i&-1490.2+742.2i&2899.5\,-4163.9i\\\
301\,+1028.6i&-112.9-1099.1i&2190.1\,-3783.9i\\\ \end{array}\right)$

Similarly, the fixed point
$\frac{1}{2}\left(\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$
of Eq.(1) is _locally asymptotically stable_ for the set of parameters which
is given as matrix as follows:

$\left(\begin{array}[]{ccc}2686.1\,-1250.3i&629.7\,-417i&-1594.+4310.2i\\\
1419.2\,-2938.1i&588.9\,-568.9i&-2005.2+421i\\\
763.8\,+3108i&-646.-533.3i&-1073.2-3104.1i\\\
2816.4\,-1018.7i&-503.9-401.1i&-2364.7-2712.i\\\
-106.2-138.8i&231.4\,+106.1i&-631.5+854.9i\\\
2608.3\,-3767.4i&-363.6-366.9i&-1909.8+299.2i\\\
3435.3\,+2539.7i&535.3\,+20.1i&-618.1+2127.1i\\\
-984.2+955.5i&666.8\,-431.4i&26.\,-2208.2i\\\
-2113.7+3096.8i&991.3\,+434.3i&-326.8-3814.8i\\\
-669.8-794.3i&-86.8-1545.9i&-3677.5+1627.6i\\\
2663.8\,-2211i&1366.2\,+1397.7i&-3862.-2238.9i\\\
575\,+221i&1515.2\,-8.2i&1087.6\,+4131.8i\\\
992.7\,-772i&180.9\,-1155.2i&-2175.5-885.1i\\\
3337.7\,+2555i&-1931.1+111.4i&-2914.1-4788i\\\
982.5\,+492.9i&807.6\,+666.8i&-3412.6-1495.6i\\\
1908.8\,+3639.9i&-527.4+1553.7i&-3264.5+645.1i\\\
1393.1\,-853.1i&261.1\,+1276.7i&-2089.8-4087.5i\\\
-1931.3-1221.9i&340.2\,-527.6i&-1151.6-2448i\\\
-1319.1+724.3i&394.1\,+798.2i&-4187.5-992.1i\\\
1582.8\,+350.7i&251.5\,+37.2i&-2431.8+892.6i\\\ \end{array}\right)$

The $100$ set of such parameters $\alpha$, $\beta$ and $\gamma$ are plotted in
the complex plane in the following Fig. 2.

---

Figure 2: Parameters $\alpha$, $\beta$ and $\gamma$ plot for which the fixed
points are locally asymptotically stable. Here the top figure stands for the
fixed point
$\frac{1}{2}\left(\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$
and bottom one stands for the fixed point
$\frac{1}{2}\left(\sqrt{-4\alpha+\beta^{2}-2\beta\gamma+\gamma^{2}}-\beta+\gamma\right)$.

###  2.1 Local Asymptotic Stability in the case $\alpha=0,\in\mathbb{C}$

When the parameter $\alpha=0$, then the the difference equation Eq.(1) would
reduce to

$\displaystyle{z_{n+1}=\frac{\beta z_{n-k}}{\gamma-
z_{n}},n=0,1,\ldots,k\in{1,2,\dots}}$ (7)

The fixed points $z_{\pm}$ of the equation Eq.(7) are $0$ and $\gamma-\beta$
respectively.

In the similar fashion we did earlier, the fixed point $z_{+}=0$ of the
difference equation Eq.(7) is _locally asymptotically stable_ , _unstable_ if
$\left\lvert\frac{\beta}{\gamma}\right\rvert<1$,
$\left\lvert\frac{\beta}{\gamma}\right\rvert>1$ respectively.

Similarly, the fixed point $z_{-}=\gamma-\beta$ is _locally asymptotically
stable_ if $\left\lvert\frac{\gamma}{\beta}-1\right\rvert<0$ i.e.
$\left\lvert\frac{\gamma}{\beta}\right\rvert<1$. If
$\left\lvert\frac{\gamma}{\beta}\right\rvert>1$ then the fixed point is
_unstable_.

###  2.2 Local Asymptotic Stability in the case $\alpha=\beta=\gamma$

When all the three parameters are same $(\beta=\gamma=\alpha)$, then the
difference equation Eq.(1) becomes

$\displaystyle{z_{n+1}=\frac{\alpha+\alpha z_{n-k}}{\alpha-
z_{n}}},n=0,1,2,\ldots$ (8)

The fixed points of the difference equation Eq.(8) are $\pm i\sqrt{\alpha}$.

The fixed points $\pm i\sqrt{\alpha}$ of Eq.(8) are _locally asymptotically
stable_ if

$\left\lvert\frac{2}{1+\mp\frac{i}{\sqrt{\alpha}}}\right\rvert<1$

Here we shall look for the subset $\mathbb{S}$ of $\mathbb{C}$ of the
parameter $\alpha$ for which the above stated condition does hold good. Before
seeking the set of our interest, let us try to envisage the region and we
found it to be the following in the Fig. 3 where the real part and imaginary
part of the parameter $\alpha$ are assumed to varied over the region
$(-1,1)\times(-1,1)$.

Now we can see what we are to find, i.e. dependence of the boundary of this
set, i.e. we should find a few functions yielding $y$ as a function of $x$ on
the boundary. where $\alpha=x+iy$. We have the following $x$ and $y$ which
satisfy the boundary
$\left\lvert\frac{2}{1+\mp\frac{i}{\sqrt{\alpha}}}\right\rvert=1$.

---

---

Figure 3: Region of the parameter $\alpha$ while its real and imaginary part
are varying along the (-1,1) and (-1,1) respectively.

To figure out why the solution looks slightly involved, we plot all the roots
of the underlying polynomial then the curves will form more symmetric pattern:

---

Figure 4: Boundary of the roots of the underlying polynomials. Different
colors denote four different polynomials.

Here we go with a set of examples of parameter $\alpha$ for which the fixed
points $\pm i\sqrt{\alpha}$ of the difference equation Eq.(8) are _locally
asymptotically stable_.

_A=_ $\\{\\{\alpha\to 0.294548\,+0.27933i\\},\\{\alpha\to
0.303387i\\},\\{\alpha\to-0.111111\\},\\{\alpha\to-0.111111+0.288786i\\},\\{\alpha\to-0.12042i\\}\\}$

_B=_
$\\{\\{\alpha\to-0.428571i\\},\\{\alpha\to-0.526596-0.0285714i\\},\\{\alpha\to-0.0725055\\},\\{\alpha\to
0.325133\,-0.0233645i\\},\\{\alpha\to 0.28391\,+0.0116959i\\}\\}$

For these two sets _A_ and _B_ of parameters the fixed points $\mp
i\sqrt{\alpha}$ is _locally asymptotically stable_ respectively.

###  2.3 Permanence and Boundedness of Solutions

######  Theorem 2.3.

Let
$\left\lvert\alpha\right\rvert<\left\lvert\frac{(\gamma-\beta)^{2}}{4}\right\rvert$
and the initial values $z_{-k},z_{-k+1},\dots,z_{-1},z_{0}\in
B(0,\left\lvert\bar{z}_{1}\right\rvert)$. If $\\{z_{n}\\}$ is any solution of
Eq.(1), then $z_{n}\in B(0,\left\lvert\bar{z}_{1}\right\rvert)$ for all $n\geq
1$, that is $\\{z_{n}\\}$ is permanent and bounded.

Proof: If the the initial values $z_{-k},z_{-k+1},\dots,z_{-1},z_{0}\in
B(0,\left\lvert\bar{z}_{1}\right\rvert)$, then $\left\lvert
z_{-k}\right\rvert<\left\lvert\bar{z}_{1}\right\rvert,\left\lvert
z_{-k+1}\right\rvert<\left\lvert\bar{z}_{1}\right\rvert,\dots,\left\lvert
z_{-1}\right\rvert<\left\lvert\bar{z}_{1}\right\rvert$ and $\left\lvert
z_{0}\right\rvert<\left\lvert\bar{z}_{1}\right\rvert$, we have

$0<\left\lvert z_{1}\right\rvert=\left\lvert\frac{\alpha+\beta z_{-k}}{\gamma-
z_{0}}\right\rvert<\left\lvert\frac{\alpha+\beta\bar{z}_{1}}{\gamma-\bar{z}_{1}}\right\rvert=\left\lvert\bar{z}_{1}\right\rvert$

Similarly,

$0<\left\lvert z_{2}\right\rvert=\left\lvert\frac{\alpha+\beta
z_{-k+1}}{\gamma-
z_{1}}\right\rvert<\left\lvert\frac{\alpha+\beta\bar{z}_{1}}{\gamma-\bar{z}_{1}}\right\rvert=\left\lvert\bar{z}_{1}\right\rvert$

Using _Mathematical Induction_ we have $0<\left\lvert
z_{n}\right\rvert<\left\lvert\bar{z}_{1}\right\rvert$ for all $n\geq 1$.
Therefore, the open ball $B(0,\left\lvert\bar{z}_{1}\right\rvert)$ is
invariant of Eq.(1) and the solution $\\{z_{n}\\}$ of Eq.(1) is permanent and
bounded.

######  Theorem 2.4.

Let
$\left\lvert\alpha\right\rvert=\left\lvert\frac{(\gamma-\beta)^{2}}{4}\right\rvert$
and the initial values $z_{-k},z_{-k+1},\dots,z_{-1},z_{0}\in
B(0,\left\lvert\frac{\gamma-\beta}{2}\right\rvert)$ where
$\left\lvert\beta\right\rvert<\left\lvert 3\gamma\right\rvert$. If
$\\{z_{n}\\}$ is any solution of Eq.(1), then $z_{n}\in
B(0,\left\lvert\frac{\gamma-\beta}{2}\right\rvert)$ for all $n\geq 1$, that is
$\\{z_{n}\\}$ is _permanent and bounded_.

Proof: Considering $\left\lvert\beta\right\rvert<\left\lvert
3\gamma\right\rvert$, we have
$0<\left\lvert\frac{\gamma-\beta}{2}\right\rvert<\left\lvert\gamma\right\rvert$
and the initial values $z_{-k},z_{-k+1},\dots,z_{-1},z_{0}\in
B(0,\left\lvert\frac{\gamma-\beta}{2}\right\rvert)$, then $\left\lvert
z_{-k}\right\rvert<\left\lvert\frac{\gamma-\beta}{2}\right\rvert,\left\lvert
z_{-k+1}\right\rvert<\left\lvert\frac{\gamma-\beta}{2}\right\rvert,\dots,\left\lvert
z_{-1}\right\rvert<\left\lvert\frac{\gamma-\beta}{2}\right\rvert$ and
$\left\lvert z_{0}\right\rvert<\left\lvert\frac{\gamma-\beta}{2}\right\rvert$,
we have

$0<\left\lvert z_{1}\right\rvert=\left\lvert\frac{\alpha+\beta z_{-k}}{\gamma-
z_{0}}\right\rvert<\left\lvert\frac{\alpha+\beta\frac{\gamma-\beta}{2}}{\gamma-\frac{\gamma-\beta}{2}}\right\rvert=\left\lvert\frac{\gamma-\beta}{2}\right\rvert$

Similarly,

$0<\left\lvert z_{2}\right\rvert=\left\lvert\frac{\alpha+\beta
z_{-k+1}}{\gamma-
z_{1}}\right\rvert<\left\lvert\frac{\alpha+\beta\frac{\gamma-\beta}{2}}{\gamma-\frac{\gamma-\beta}{2}}\right\rvert=\left\lvert\frac{\gamma-\beta}{2}\right\rvert$

Using _Mathematical Induction_ we have $0<\left\lvert
z_{n}\right\rvert<\left\lvert\frac{\gamma-\beta}{2}\right\rvert$ for all
$n\geq 1$. Therefore, the open ball
$B(0,\left\lvert\frac{\gamma-\beta}{2}\right\rvert)$ is invariant of Eq.(1)
and the solution $\\{z_{n}\\}$ of Eq.(1) is _permanent and bounded_.

An example is taken here as an evidence of the Theorem 2.4. Consider
$\gamma=5+6i$, $\beta=2+3i$ and so
$\alpha=\left\lvert\frac{(\gamma-\beta)^{2}}{4}\right\rvert=0+4.5i$. Here the
fixed point is $\frac{\gamma-\beta}{2}=1.5+1.5i$ with modulus $2.1213$. Also
we assume k=2 and the initial values $z_{-2}=i,z_{-1}=1+i$ and $z_{0}=1.5+i$
note that $\left\lvert z_{-2}\right\rvert,\left\lvert z_{-1}\right\rvert$ and
$\left\lvert z_{0}\right\rvert$ are all less than
$\left\lvert\frac{\gamma-\beta}{2}\right\rvert=2.1213$ and
$\left\lvert\beta\right\rvert<3\left\lvert\gamma\right\rvert$. So the Theorem
2.4 is applicable here. What we expect is that all the solution $\\{z_{n}\\}$
of Eq.(1) is _permanent and bounded_ , that is $0<\left\lvert
z_{n}\right\rvert<\left\lvert\frac{\gamma-\beta}{2}\right\rvert$ for all
$n\geq 1$. The complex trajectory plot including its time series plot are
given in the Fig. 5.

---

Figure 5: Complex trajectory plot including its time series plot.

In the Fig. $5$, it is seen that all the solutions $\\{z_{n}\\}$ of Eq.(1) are
lying in the disk $B(0,\left\lvert\frac{\gamma-\beta}{2}\right\rvert)$.

##  3 Periodic of Solutions

We shall first look for the _prime period two solutions_ of the three
difference equation Eq.(1).

Let $\ldots,\phi,~{}\psi,~{}\phi,~{}\psi,\ldots$, $\phi\neq\psi$ be a _prime
period two solution_ of the difference equation Eq.(1).
If $k$ is even, then $z_{n}=z_{n-k}$ and then
$\psi=\frac{\alpha+\beta\phi}{\gamma-\phi}$ and
$\phi=\frac{\alpha+\beta\psi}{\gamma-\psi}$. By solving these two equations we
get $(\gamma+\beta)(\psi-\phi)=0$. So $\gamma=-\beta$ will lead to _prime
period two solutions_ of the difference equation Eq.(1).
if $k$ is odd, then $z_{n+1}=z_{n-k}$ and then
$\psi=\frac{\alpha+\beta\psi}{\gamma-\phi}$ and
$\phi=\frac{\alpha+\beta\phi}{\gamma-\psi}$. By solving these two equations we
get $(\gamma-\beta)(\psi-\phi)=0$. Therefore $\gamma=\beta$ will lead to
_prime period two solutions_ of the difference equation Eq.(1).

It is noted that for $\gamma>\beta>0$ in the real line, there was no _prime
period two solution_ of the difference equation Eq.(1). Here we list a few
examples where we found the higher order periods viz. $9,10,13$ and etc.

Parameters: $\alpha$, $\beta$, $\gamma$ and Delay term: $k$ |  Periodic Solutions |  Remarks
---|---|---
$\alpha=0.2729+0.0372i$, $\beta=0.7690+0.3960i$, $\gamma=-\beta$ and $k=1$ |  $\psi=-0.0444+0.1575i$ $\phi=0.2333+1.6669i$, Period: 2 |
$\alpha=0.4024+0.9828i$, $\beta=0.0740+0.6841i$, $\gamma=\beta$ and $k=2$ |  $\psi=-0.5589-0.1317i$ $\phi=-0.1991+0.9915i$, Period: 2 |
$\alpha=-0.1111+0.2337i$, $\beta=\alpha$, $\gamma=\beta$ and $k=6$ |  Period: 9 |
$\alpha=-0.1111+0.2337i$, $\beta=\alpha$, $\gamma=\beta$ and $k=8$ |  Period: 11 |
$\alpha=-0.1111+0.2337i$, $\beta=\alpha$, $\gamma=\beta$ and $k=12$ |  Period: 15 |
$\alpha=-0.1111+0.2337i$, $\beta=\alpha$, $\gamma=\beta$ and $k=17$ |  Period: 20 |
Table 2: Higher Order Periodic Solutions where $\gamma=\pm\beta$, of the equation Eq.(1) for different initial values. In the right most column the corresponding periodic trajectory plots are given. Parameters: $\alpha$, $\beta$, $\gamma$ and Delay term: $k$ |  Periodic Solutions |  Remarks
---|---|---
$\alpha=-0.1111+0.2337i$, $\beta=\alpha$, $\gamma=\beta$ and $k=38$ |  Period: 41 |
$\alpha=-0.1111+0.2337i$, $\beta=\alpha$, $\gamma=\beta$ and $k=40$ |  Period: 43 |
$\alpha=-0.1111+0.2337i$, $\beta=\alpha$, $\gamma=\beta$ and $k=52$ |  Period: 55 |
$\alpha=-0.1111+0.2337i$, $\beta=\alpha$, $\gamma=\beta$ and $k=80$ |  Period: 83 |
$\alpha=-0.1111+0.2337i$, $\beta=\alpha$, $\gamma=\beta$ and $k=100$ |  Period: 103 |
$\alpha=-0.1111+0.2337i$, $\beta=\alpha$, $\gamma=\beta$ and $k=5100$ |  Period: 5103 |
Table 3: Higher Order Periodic Solutions where $\gamma=\pm\beta$, of the
equation Eq.(1) for different initial values. In the right most column the
corresponding periodic trajectory plots are given.

In the Table $2$ and Table $3$, a set of examples of higher order periods are
given for different values of $k$ and parameters. What is interesting to note
here is that the period is associated to the delay term $k$ as $p=k+3$ ($p$
denotes period) where all the parameters are taken as same
($\alpha=\beta=\gamma$). This inspired us to propose the following
conjectures.

######  Conjecture 3.1.

There exist as many as higher order periods of the difference Eq.(1) is
demanded.

######  Conjecture 3.2.

For $\alpha=\beta=\gamma$, period ($p$) is increased by $3$ with the delay
term $k$. i.e. $p=k+3$.

##  4 Chaotic Solutions

A new dynamical behavior of the rational difference equation Eq.(1) is
chaoticity which was not present in the real set of parameters. It is really
hard to determine the set of all parameters $\alpha$, $\beta$ and $\gamma$ for
which the solutions of the Eq.(1) are chaotic but computationally we have
encountered some chaotic solutions for some values of the parameters which are
given in the following Table. 4.
The Lyapunov characteristic exponents serves as a useful tool to quantify
chaos. Specifically Lyapunav exponents measure the rates of convergence or
divergence of nearby trajectories. Negative Lyapunov exponents indicate
convergence, while positive Lyapunov exponents demonstrate divergence and
chaos. The magnitude of the Lyapunov exponent is an indicator of the time
scale on which chaotic behavior can be predicted or transients decay for the
positive and negative exponent cases respectively. In this present study, the
largest Lyapunov exponent is calculated for a given solution of finite length
numerically [10].
From computational evidence, it is arguable that for complex parameters
$\alpha$, $\beta$ and $\gamma$ which are stated in the following table the
solutions are chaotic for every initial values.

Parameters $\alpha$, $\beta$, $\gamma$ and Delay term $k$ |  Lyapunav exponent
---|---
$\alpha=(0.8003,0.1419)$, $\beta=(0.1576+0.9706)$, $\gamma=0.9572+0.4854)$ and $k=1$ |  $0.0543$
$\alpha=(0.2217,0.1174)$, $\beta=(0.6028,0.7112)$, $\gamma=(0.2217,0.1174)$ and $k=2$ |  $1.062$
$\alpha=(0.2564,0.6135)$, $\beta=(0.6620,0.4162)$, $\gamma=(0.8419,0.8329)$ and $k=3$ |  $1.314$
$\alpha=(0.5349,0.7210)$, $\beta=(0.4795,0.6393)$, $\gamma=(0.557,0.6473)$ and $k=4$ |  $1.373$
Table 4: Chaotic solutions of the equation Eq.(1 for different choice of
parameters and initial values.

The chaotic trajectory plots including corresponding complex plots are given
the following Fig. 6.

---

Figure 6: Chaotic Trajectories of the equation Eq.(1 of four different cases
as stated in Table 4.

In the Fig. 6, for each of the four cases ten different initial values are
taken and plotted in the left and in the right corresponding complex plots are
given. From the Fig. 6, it is evident that for the four different cases the
basin of the chaotic attractor is neighbourhood of the centre $(0,0)$ of
complex plane. It is noted that the _Conjecture_ $3.1$ and $3.2$ together
suggest that the chaos is simply proportional to the delay term $k$ of the
rational dynamical system since the existence of all possible periods suggest
that the trajectory is eventually chaotic.

## References

  * [1] E.M.E. Zayed and M.A. El-Moneam, (2009), On the rational recursive sequence $x_{n+1}=\frac{\alpha+\beta x_{n-k}}{\gamma-x_{n}}$, _J Appl Math Comput_ , 31: 229 - 237.
  * [2] Xing-Xue Yan and Wan-Tong Li (2003), Global attractivity in the recursive sequence $x_{n+1}=\frac{\alpha-\beta x_{n}}{\gamma-x_{n-1}}$ _Applied Mathematics and Computation_ , 138, 415 - 423.
  * [3] Chenquan Gan, Xiaofan Yang, and Wanping Liu, (2013) Global Behavior of $x_{n+1}=\frac{\alpha+\beta x_{n-k}}{\gamma+x_{n}}$, _Discrete Dynamics in Nature and Society, Hindawi Publishing Corporation_ , Volume 2013, Article ID 963757.
  * [4] Sk. S. Hassan, E. Chatterjee, (2015), Dynamics of the equation $\displaystyle{z_{n+1}=\frac{\alpha+\beta z_{n}}{A+z_{n-1}}}$ in the Complex Plane, _Cogent Mathematics, Taylor & Francis_, 2 (1), 1-12.
  * [5] Sk. S. Hassan, Dynamics of Delay Logistic Difference Equation in the Complex Plane, _arXiv:1507.02964_ [math.DS].
  * [6] Sk. S. Hassan, Sk. S. Hassan, Complex Dynamics of $z_{n+1}=\frac{\alpha+\beta z_{n}+\gamma z_{n-1}}{A+Bz_{n}+Cz_{n-1}}$ _arXiv:1509.00850_ [math.DS].
  * [7] Saber N Elaydi, Henrique Oliveira, José Manuel Ferreira and João F Alves, (2007) _Discrete Dyanmics and Difference Equations, Proceedings of the Twelfth International Conference on Difference Equations and Applications_ , World Scientific Press.
  * [8] M.R.S. Kulenovi$\acute{c}$ and G. Ladas, (2001) _Dynamics of Second Order Rational Difference Equations; With Open Problems and Conjectures_ , Chapman & Hall/CRC Press.
  * [9] V.L. Kocic and G. Ladas, (1993) _Global Behaviour of Nonlinear Difference Equations of Higher Order with Applications_ , Kluwer Academic Publishers, Dordrecht, Holland.
  * [10] A. Wolf, J. B. Swift, H. L. Swinney and J. A. Vastano, (1985), Determining Lyapunov exponents from a time series Physica D, 126, 285-317.


# gACSON software for automated segmentation and morphology analyses of
myelinated axons in 3D electron microscopy

Andrea Behanova 111Shared first authorship. Ali Abdollahzadeh 222Shared first
authorship. Ilya Belevich  Eija Jokitalo  Alejandra Sierra333Shared senior
authorship<EMAIL_ADDRESS>Jussi Tohka 444Shared senior
authorship<EMAIL_ADDRESS>Department of Information Technology, Uppsala
University, Uppsala, Sweden.  Biomedical Imaging Unit, A.I.Virtanen Institute
for Molecular Sciences, University of Eastern Finland, Kuopio, Finland.
Electron Microscopy Unit, Institute of Biotechnology, University of Helsinki,
Helsinki, Finland.

(Matlab-based; the software; It; we estimated that manual instance
segmentation of axons in a white matter EM with 300 million voxels of the size
$15\times 15\times 50nm^{3}$ would require 2400 h; webKnossos [11]; VAST [16],
NeuroMorph [17]; open source; reading and writing a variety of image formats;
many image processing tools; excellent alignment features.; includes all MIB
image processing tools; user friendly GUI for training and inference of deep
models.; 3D visualization of large EM volumes, not limited to the RAM size;
visualizing large EM volumes; skeletonization by manually placing and
connecting nodes; manual processing of pre-segmented data.; fast browsing
speeds; efficient distributed EM annotation in petabytes-scale; provides
flight mode, a single-view egocentric reconstruction method to accelerate
annotation.; extracts morphometric information, e.g., axon diameter or myelin
g-ratio.; easy training procedure.; image stitching, registration (using scale
invariant feature transform), editing, and annotation features (floating text
labels).; multiple linked image stack display; many tools for collaborative
microcircuit reconstruction; based with Direct3D graphics for optimal
performance; petabyte-scale image data browsing; many general labeling tool;
annotating arbitrary structures, regions and locations, depending on the
user’s needs; proof-reading of segmentations.; in Blender; import and
visualize mesh models; cross-sectional analysis of axons and their synaptic
connections; measuring density of organelles; mesh objects should be generated
by other 3D image segmentation software.; object tracking; interactive
processing of image sub-volumes, followed by offline batch mode for complete
volume analysis.; quantification of myelin g-Ratio.; automated segmentation
and cross-sectional morphology analysis of myelinated axons in EM volumes;
instance segmentation of myelin; myelin g-ratio measurement; 3D volume
visualization, provides proof-reading tools.; Matlab-based; It; the software;
the new tool; All the parameters of gACSON are summarized in Supplementary
Table 1.; We suppress these complex noise patterns because the segmentation
algorithms in steps 2 and 3 are affected by it. Specifically, region growing
algorithms are sensitive to image noise as a voxel belonging to a region may
seem dissimilar from the average intensity of the region simply because of
noise. For denoising, we apply the non-local BM4D algorithm [20]. Unlike local
denoising filters operating in a small window around the target voxel, non-
local denoising considers all the voxels in the image and weights them by the
similarity of their neighbourhood to the neighbourhood of the target voxel.
This results in better preservation of the image texture and edges and reduces
the loss of details in the image compared with local denoising algorithms.; in
particular; The magnification factor of the zoomed-in images in the orange
boxes is 300 %.; The magnification factor of the zoomed-in images in the
orange boxes is 300 %.; 2D slices; we; Scribbles can be drawn on single or
multiple 2D slices of the 3D volume leading to a training set to train the
random forest using the above described features for each voxel. The trained
random forest then automatically extends training set to predict the myelin
class in the whole 3D EM volume. The user can scroll between the 2D images and
add new scribbles if needed. It is not necessary to adapt means or contrast of
the image, as we assume that contrast normalization has been applied as a pre-
processing step to 2D images of the 3D-EM volume.; noisier; noisier; (see
Figure 4a for an example of semantic segmentation of myelin); See Figure 4b
for an example of instance segmentation of the intra-axonal space); The BVG;
In more detail, we compute; and generate; this seeding method; the method; we
compute; we assign; The algorithm; we; the software; The morphometry
algorithm; it; Thereafter, it; is calculated; The manual segmentation panel
allows; e.g., a leakage at nodes of Ranvier or if staining artifacts were
mistakenly segmented as myelin.; The user can also discard a segmentation
label, for example, when a myelinated axon is partially outside the image
domain, and replace it with a null label representing the background;
evaluating; The morphometry panel is used; The morphometry algorithm; We used
three adult male Sprague-Dawley rats (10-weeks old, weight of
$320\text{\,}\mathrm{-}$$380\text{\,}\mathrm{g}$, Harlan Netherlands B.V.,
Horst, Netherlands). The rats were accommodated in separated cages in an
environment of $22\text{\,}\mathrm{\SIUnitSymbolCelsius}$ temperature, 60$\%$
humidity, 12 h of light per day, and unlimited water and food. Animal Care and
Use Committee of the Provincial Government of Southern Finland authorized all
the actions performed on the animals accordingly to the ground rules addressed
by the European Community Council Directive 86/609/EEC. TBI was induced by
lateral fluid percussion injury in two rats, as described in [30] using a
fluid-percussion device. The impact pressure was adjusted to a 3.2-3.4 atm to
induce a severe injury. The sham-operated rat underwent all the surgical
procedures except the impact. After five months, the rats were transcardially
perfused using 0.9$\%$ NaCl followed by 4$\%$ paraformaldehyde and the brains
were post-fixed in 4$\%$ paraformaldehyde 1$\%$ glutaraldehyde overnight at
$4\text{\,}\mathrm{\SIUnitSymbolCelsius}$. Then, the brains were sectioned
into 1-mm thick coronal sections and a section at approximately -3.80 mm from
bregma from each brain were further dissected into smaller samples containing
the areas of interest. The samples were stained using an enhanced staining
protocol described in [31]. After selecting the area of interest within the
samples, the blocks were further trimmed and imaged in a scanning electron
microscope (Quanta 250 Field Emission Gun; FEI Co., Hillsboro, OR, USA),
equipped with the 3View system (Gatan Inc., Pleasanton, CA, USA) using a
backscattered electron detector (Gatan Inc.). The face of the blocks was the
x-y plane, and the z-axis was the cutting direction. All blocks were imaged
using a beam voltage of $2.5\text{\,}\mathrm{kV}$ and a pressure of 0.15 Torr.
Microscopy Image Browser [8] (MIB; http://mib.helsinki.fi) was used to find
the bounding box of the collected 2D-EM images and run the voxel-size
calibration. The collected images were further aligned by measuring the
translation between the consecutive images using the cross-correlation cost
function (MIB, drift correction). Finally, the contrast of the collected
images was normalized such that the mean and standard deviation of the
histogram of each image matches the mean and standard deviation of the whole
image stack. Information about the 3D-EM data are shown in Table 2.; We
collected two samples from each brain: the ipsilateral and contralateral
hemispheres of the layer VI of the primary somatosensory cortex.; along the
length of the axon (L); along the length of the axon (L); The reported
p-values have not been corrected for multiple comparisons.; The reported
p-values have not been corrected for multiple comparisons.; It would be
possible to exclude axons also manually, but we preferred to set a length-
threshold to increase the objectivity and repeatability of the analysis.; We
provide here uncorrected p-values for clarity as 1) our focus on
quantification, not hypothesis testing, 2) there are several alternatives to
set the number of tests (i.e., should we correct for tests between Sham and
TBI and the tests between the hemispheres or should we consider difference
between the Sham and TBI and hemispheric differences as separate experiments),
and 3) knowing that p-values are uncorrected, the readers can adjust their
expectations [33]; Supplementary; Supplementary; although if any multiple
comparisons correction took place that effect would disappear; Supplementary;
Supplementary; Especially in TBI rats, the hemisphere effects were strong and
would survive multiple comparison correction.; The BVG algorithm is written in
“plain” Matlab-code optimized for speed, while BM4D
(https://webpages.tuni.fi/foi/GCF-BM3D), SLIC
(https://github.com/fk128/SLICSupervoxels), and multi-stencils fast marching
and Runge-Kutta method used in the skeletonization process
(https://github.com/cbm755/fast_marching_kroon) rely on compiled MEX files.
Note that to reduce the computation time advantages of compiled programming
languages such as C and C++ compared to interpreted programming languages such
as Matlab, the Matlab execution engine uses just-in-time compilation of all
Matlab code with a single execution pathway that considerably speeds-up the 3D
voxel-based processing.; The peak of the memory consumption of the pipeline
occurred after calculating SLIC supervoxels when the pipeline consumed 11 GB
of memory from which approximately 6.7 GB was consumed by SLIC and the rest
was due to maintaining multiple copies of data in the memory.; Matlab-based;
software; The software; Importantly, it; They; Also, despite the volumes
included in this study were selected in a specific area, the layer VI of the
primary somatosensory cortex, variations in the cellular content of these
small volumes may occur, which might affect the comparative results.; Certain
limitations and challenges need to be addressed in future software updates.
For example, the size of datasets cannot exceed the RAM size because the
segmentation pipeline loads the whole dataset into memory. For a very large EM
volume this requires tiling EM datasets into sub-volumes and then performing a
segmentation task. This task should be automated. The instance segmentation of
myelin may include inaccuracies when adjacent myelin sheaths have a
substantial difference in thickness as the current watershed-based
segmentation divides the myelin content equally between two adjacent axons. We
plan to use distance transform with a propagation speed relative to myelin
thickness to remedy such scenarios.; Supplementary Information)

###### Abstract

Background and Objective: Advances in electron microscopy (EM) now allow
three-dimensional (3D) imaging of hundreds of micrometers of tissue with
nanometer-scale resolution, providing new opportunities to study the
ultrastructure of the brain. In this work, we introduce a freely available
gACSON software for visualization, segmentation, assessment, and morphology
analysis of myelinated axons in 3D-EM volumes of brain tissue samples.
Methods: The software is equipped with a graphical user interface (GUI). It
automatically segments the intra-axonal space of myelinated axons and their
corresponding myelin sheaths and allows manual segmentation, proofreading, and
interactive correction of the segmented components. gACSON analyzes the
morphology of myelinated axons, such as axonal diameter, axonal eccentricity,
myelin thickness, or g-ratio. Results: We illustrate the use of by segmenting
and analyzing myelinated axons in six 3D-EM volumes of rat somatosensory
cortex after sham surgery or traumatic brain injury (TBI). Our results suggest
that the equivalent diameter of myelinated axons in somatosensory cortex was
decreased in TBI animals five months after the injury. Conclusions: Our
results indicate that gACSON is a valuable tool for visualization,
segmentation, assessment, and morphology analysis of myelinated axons in 3D-EM
volumes. is freely available at https://github.com/AndreaBehan/g-ACSON under
the MIT license.

###### keywords:

myelinated axons, electron microscopy, segmentation, brain, morphology

††journal: Computer Methods and Programs in Biomedicine

##  1 Introduction

Assessing the structure of the brain is critical to better understanding its
normal and abnormal functioning. Advances in electron microscopy (EM) now
allow three-dimensional (3D) imaging of hundreds of micrometers of tissue with
nanometer-scale resolution, providing new opportunities to study the
ultrastructure of the brain [1, 2]. Quantitative analysis of 3D-EM data, such
as morphological assessment of ultrastructure, spatial distribution or
connectivity of cells, requires the instance segmentation of individual
ultrastructural components [3, 4, 5]. Performing this segmentation manually is
tedious, if not impossible, due to the large size and enormous number of
components in typical 3D-EM data. For example, manual labeling of 215 neurites
in a volume of 500 million voxels required 1500 h [6], and [7]. Therefore, the
analyses of 3D-EM data of brain tissue require the development of advanced
software tools that allow neuroscientists to automatically visualize, segment,
and extract the geometric and topological features of brain ultrastructures.

There are several software tools for analyzing 3D-EM data, including open
source software packages such as Microscopy Image Browser (MIB) [8], DeepMIB
[9], Knossos [10], , AxonSeg [12], AxonDeepSeg [13], TrackEM2 [14], CATMAID
[15], , SegEM [6], Ilastik [18], or MyelTracer [19], and proprietary software
such as Amira (Thermofisher ScientificTM, MA, www.thermofisher.com/amira-
avizo) or Imaris (Oxford Instruments, Abingdon, UK, www.imaris.oxinst.com).
Table LABEL:tbl:software summarizes the main features of the software
packages. However, the segmentation tools in these packages have mainly been
designed for manual or semi-automatic segmentation, which makes their use very
labor intensive for the user. The supervised learning-based segmentation
methods available in these packages, as in DeepMIB [9] or AxonDeepSeg [13],
typically require manually annotated training data, which can be difficult and
tedious to annotate. Moreover, while some of these packages have excellent
visualization capabilities, they are not geared towards interactive
correction, proofreading, and validation of automatic segmentations.

Table 1: Examples of open-source software tools for 2D and 3D-EM image analysis and their main features. Software |  Main features
---|---
MIB [8] |  Matlab-based; manual and semi-automated segmentation of large EM volumes (graph cuts method, supervoxels);
DeepMIB [9] |  Matlab-based; supervised deep learning-based image segmentation of EM and light microscopy volumes;
Knossos [10] |  Python-based;
webKnossos [11] |  In-browser;
AxonSeg [12] |  Matlab-based; automated 2D segmentation of axons and myelin on histology images (the extended minima algorithm; active contours to detect the outer and the inner boundaries of myelin sheath);
AxonDeepSeg [13] |  Python-based; supervised deep learning-based software for semantic segmentation of axons and myelin from microscopy images (3-class semantic segmentation task);
TrackEM2 [14] |  ImageJ plugin; manual and semi-automated segmentation of large EM volumes (brushing, skeletonization);
CATMAID [15] |  Python-based; fast terabyte-scale image data browsing; user-interactive image annotation; WebGL neuronal morphology viewer;
VAST [16] |  C++
NeuroMorph [17] |  Operating
SegEM [6] |  Matlab-based; semi-automated segmentation of 3D-EM data using convolutional neural networks and manually provided neuronal skeletons.
Ilastik [18] |  Python-based; automated and semi-automated segmentation of EM volumes (graph cut methods for instance segmentation);
MyelTracer [19] |  Python-based; 2D semi-automated segmentation of axons and myelin (intensity thresholding and user’s drawn lines to generate contours for myelin segmentation);
gACSON |  Matlab-based;

In this context, we developed a Graphical user interface for AutomatiC
Segmentation of axONs, called gACSON, and axonal morphology analysis in large
3D-EM data of the brain. The gACSON builds on our previous pipelines ACSON [7]
and DeepACSON [3] for 3D-EM segmentation and morphometrics and provides an
easy-to-use graphical user interface (GUI) and additional segmentation
algorithms for myelin and the intra-axonal space of myelinated axons. The
software allows for automatic unsupervised segmentation that requires no
manually annotated training data. While ACSON and DeepACSON are automated
segmentation pipelines, gACSON is an interactive tool, providing user-friendly
tools to review and correct the segmentations. also allows segmentation of a
subset of all myelinated axons so that the user does not need to wait for the
completion of the segmentation of all myelinated axons. More importantly,
provides instance segmentation of myelin that maps myelin sheaths to the
appropriate intra-axonal space instance. This feature makes unique compared to
our existing ACSON and DeepACSON pipelines [7]. With the software, users can
visualize and review the segmented ultrastructure or correct the results of an
automated segmentation by providing a dynamic environment for interactive
merge or split operations on the segmented components. Like the ACSON and
DeepACSON pipelines, gACSON extracts morphological parameters of segmented
myelinated axons, such as axonal diameter or eccentricity, but extends the
morphological parameters to measure myelin thickness and g-ratio by enabling
instance segmentation of myelin. Moreover, visualization, segmentation
evaluation, or correction of automatic segmentation results is not limited to
EM but can be applied to other biomedical imaging modalities, such light
microscopy images.

Figure 1: gACSON segmentation pipeline. The pipeline takes 3D-EM volume as
input and 1) denoises the EM volume using the BM4D method, 2) segments the
myelin in the denoised volume using a user-defined method (thresholding,
seeded region growing (SRG) or an ML-based approach), 3) performs the instance
segmentation of the intra-axonal space of myelinated axons, and 4) performs
the instance segmentation of myelin. The numbers in parentheses refer to
subsections where the steps are detailed.

The structure of this paper is as follows: Sections 2.1 and 2.2 present the
gACSON segmentation and morphometry pipelines, respectively. Section 2.3
introduces the GUI including the parameter settings and the tools for
proofreading and correcting the automatic segmentation and evaluating the
accuracy of the segmentation results. In Section 3, we demonstrate the
application of the software to segment a new dataset of the six serial block-
face scanning electron microscopy (SBEM) volumes of rat somatosensory cortex.
We further make these volumes and segmentations publicly available.

##  2 Methods

###  2.1 gACSON segmentation pipeline

gACSON includes an automated segmentation pipeline that traces myelinated
axons in 3D-EM data of brain tissue. The segmentation pipeline, summarized in
Figure 1, consists of following steps: 1) denoising of EM volumes by the
block-matching and 4D filtering (BM4D) algorithm [20] (Section 2.1.1), 2)
semantic segmentation of myelin (Section 2.1.2), 3) instance segmentation that
provides a separate label to each myelinated axon (Section 2.1.3), and 4)
assignment of myelin to individual intra-axonal spaces, referred to as
instance segmentation of myelin (Section 2.1.4).

####  2.1.1 BM4D denoising

SBEM images are affected by primary beam noise, secondary emission noise, and
noise in the final detection system [21]. BM4D, , enhances a sparse
representation in the transform domain by grouping similar 3D image patches
(i.e., continuous 3D blocks of voxels) into 4D data arrays. BM4D automatically
estimates noise type and variance and we use BM4D’s default parameter values.
Figure 2 illustrates a representative SBEM image from the somatosensory cortex
of a sham-operated rat and its BM4D-denoised version.

Figure 2: a A representative SBEM image from the somatosensory cortex of a
sham-operated rat. (b) BM4D denoising of a.

####  2.1.2 Semantic segmentation of myelin

gACSON is equipped with four different methods for semantic segmentation of
myelin: 1) interactive segmentation of myelin by intensity thresholding
(Figure 3a); 2) use of seeded region growing initiating with an automatic seed
selection; 3) use of seeded region growing with user-defined seeds (manual
seed selection); and 4) use of weakly supervised machine learning (ML) by
training a random forest [22] with a scribble annotation (see [23]) drawn by
the user on of 3D-EM data. The user can perform semantic segmentation of
myelin using various algorithms that enable accurate, automated, or
computationally efficient segmentation. For example, a ML -based segmentation
produces a more accurate segmentation, but requires user input and can be more
computationally intensive than a thresholding approach.

In ML-based segmentation, use the following features to classify voxels: 1)
the image intensity, 2) the Gaussian-smoothed image intensity ($\sigma=2)$, 3)
the Laplacian of the Gaussian image ($\sigma=0.5$), 4) the Gaussian gradient
magnitude image ($\sigma=2$), and 5) the difference of Gaussians ($\sigma=5$).
$\sigma$ is the standard deviation of the isotropic Gaussian smoothing kernel
in voxels. Figure 3b-left panel shows an annotation example drawn by the user
to segment myelin against other components, and Figure 3b-right panel is the
resulting semantic segmentation.

Regarding seeded region growing, gACSON automatically selects the seeds for
seeded region growing by either randomly distributing the seeds over the image
region or placing seeds based on the intensity values in the images.

Figure 3: (a) Semantic segmentation of myelin by a user-defined threshold. (b)
Machine learning (ML)-based semantic segmentation of myelin. Left panel: user
specifies myelin (magenta) against other ultrastructures (green) by drawing a
scribble on an EM image. These scribbles provide data to train a random forest
model. Right panel: the semantic segmentation of myelin (magenta) against
other ultrastructures (green). The thresholding-based segmentation of myelin
incorrectly includes cell membranes and is therefore than the ML-based
segmentation. Both methods segment mitochondria as myelin.

####  2.1.3 Instance segmentation of the intra-axonal space of myelinated
axons

This step segments the complement of the myelin segment into instances of
myelinated axons using the Bounded Volume Growing algorithm (BVG) (.
algorithm, that we presented in [7], integrates seeded region growing [24] and
the Canny edge detector [25] to segment the intra-axonal space of myelinated
axons. BVG iteratively segments an EM volume into myelinated axons
$V_{1},\ldots,V_{N}$, where $V_{0}$ denotes the myelin segment. Each
myelinated axon $V_{i}$ is segmented based on region growing starting from a
seed voxel and adding more voxels from the neighborhood of the current $V_{i}$
to $V_{i}$ until no more voxels can be added. For a voxel to be added to
$V_{i}$, it must meet three criteria (the voxel must not be an edge voxel; the
voxel must not belong to any $V_{j}$ with $j<i$; the voxel must have an
intensity similar to the average intensity of $V_{i}$). BVG grows one seed at
a time and then moves on to the next seed. At the completion of an axon, all
remaining seeds within the segmented axon are discarded. For a detailed
technical description of the algorithm, see [7]. The key to successful
segmentation by BVG is the selection of seed voxels. gACSON allows the user to
choose between three seed generation techniques:

  1. 1.

Initial seeds by the regional maxima method uses the semantic segmentation of
the myelin to define the positions of the seeds as the myelin encloses the
axons. the Euclidean distance transform of the semantic segmentation of the
myelin for each 2D image slice along the $z$ axis the seeds by extracting the
regional maxima from the distance transform images.

  2. 2.

Initial seeds as the watershed centroids: the regional maxima method for
generating seeds may overestimate the number of initial seeds, which increases
the computation time. To save computation time, uses the H-maxima transform of
the distance transform of myelin segmentation as a marker for a marker-
controlled watershed segmentation [26] and defines the centroids of the
watershed segments as starting seeds for the BVG algorithm. The H-maxima
transform suppresses all local maxima whose height is less than $H$. The
marker-controlled watershed transform divides the image into regions one-to-
one based on the positions of the markers, such that the number of markers
equals the final number of watershed regions.

  3. 3.

Manual placement of start seeds: the user can manually select start seeds for
the bounded volume growing algorithm. Manual seed selection is preferred for
fast segmentation of a limited number of axons.

The accuracy, automation level, and computational efficiency of a segmentation
are determined by the choice of seeding method. For example, defining seed
positions using watershed centroids is computationally more efficient than
defining them using regional maxima because the number of seeds is reduced.
However, for the same reason, an axon is more likely to remain unsegmented
when watershed centroids seeds are used.

Figure 4: The workflow of instance segmentation of myelin. (a) Semantic
segmentation of myelin. (b) Instance segmentation of the intra-axonal space of
myelinated axons. (c) Watershed transformation using the segmented intra-
axonal spaces partitions the image space into disjoint subvolumes. (d)
Instance segmentation of myelin. (e) Instance segmentation of myelin and the
intra-axonal spaces. Colors of distinct instances in the figure have been
randomly generated and they are not comparable between different panels.

After the initial segmentation of the intra-axonal space of myelinated axons
using BVG, gACSON refines the segmentation of large volumes $V_{i}$ (with more
than $5000$ voxels) using supervoxels generated by the simple linear iterative
clustering (SLIC) method [27]. In more detail, suppose that we have $Q$
supervoxels $S_{q},\;q=1,\ldots,Q$. Then, we re-define $V_{i}$ as
$V_{i}^{\prime}=\bigcup\limits_{q\in I_{i}}S_{q}$ where

$I_{i}=\bigg{\\{}q^{\prime}\bigg{|}\dfrac{|S_{q^{\prime}}\cap
V_{i}|}{|S_{q^{\prime}}|}\geq 0.8\land\forall j<i:q^{\prime}\notin
I_{j}\bigg{\\}}.$

The second condition states that a supervoxel can belong to only a single
axon. The volumes larger than 5000 voxels are traversed in the order they are
indexed in the segmentation. Refining the segmentation with SLIC technique
eliminates most of the small volumes by relabeling and attaching them to the
larger segments. The segmentation of myelinated axons may include unmyelinated
axons or fragments of other cell processes. To decide whether an axon is
myelinated or unmyelinated, investigates the thick hollow cylinder enclosing
the axon of interest. The enclosing cylinder is formed by those supervoxels
having a common face with the axon. If the enclosing cylinder contains more
myelin than a threshold value, e.g., when more than $70\%$ of the cylinder
volume is myelin, then the axon is considered as myelinated.

Figure 5: Measurement of myelin thickness. (a) A cross section of a segmented
myelin. (b)The Euclidean distance transformation of segmented myelin. (c) The
myelin curve skeleton (red) and the boundaries of the myelin (yellow) are
superimposed on b. The myelin curve skeleton remains in the middle of the
myelin cross section. The value of the distance transform of the myelin cross-
section at any point on the myelin skeleton curve is the half-thickness of the
myelin at that point.

####  2.1.4 Instance segmentation of myelin

gACSON provides instance segmentation of myelin by mapping myelin voxels to
the corresponding intra-axonal space, with a one-to-one correspondence between
myelinated axons and myelin instances, see Fig. 4. To this end, the distance
transform of the complement of the union of segmented intra-axonal spaces,
followed by the watershed transform. Figure 4c shows how the watershed
transform divides the image space into $n$ disjoint subvolumes, where $n$ is
the number of intra-axonal spaces. To achieve the instance segmentation of the
myelin, the labels from the watershed segmentation of the intra-axonal spaces
to the myelin voxels, as shown in Figure 4d, and assigns the same label to the
myelin voxels as to the corresponding intra-axonal space, as shown in Figure
4e. Note that the watershed-based segmentation divides myelin content equally
between two adjacent axons.

Figure 6: The user-interface of the gACSON software.

###  2.2 gACSON morphometry analysis

The gACSON morphometry algorithm extracts and quantifies morphological
features such as equivalent diameter, eccentricity, or g-ratio of segmented
myelinated axons. automatically finds the axonal skeleton based on a distance
transformation-based algorithm as described in [28, 29] and extracts the
cross-sections of the axons using a perpendicular plane to the axonal skeleton
at each point along the axonal skeletons. For an EM volume with an anisotropic
voxel size, first resample the bounding box of an axon to an isotropic voxel
size, and then quantifies the morphology parameter on the resampled volume.
For each 2D cross-section, quantifies the equivalent diameter of the intra-
axonal space, i.e., the diameter of a circle with the same area as the cross-
section, the length of the minor and major axes of a fitted ellipse, and its
eccentricity.

extracts the myelin thickness and g-ratio for each myelinated axon. To measure
the myelin thickness in a cross-section (Figure 5a), first computes the
Euclidean distance transform of the myelin cross-section, as shown in (Figure
5b), where higher intensity values imply longer distances to the myelin
boundaries. calculates the curve that remains in the center of the myelin
cross-section, called the myelin skeleton curve, using the Euclidean distance
transform map (Figure 5c). The value of the distance transform at each point
on the myelin skeleton curve is half the myelin thickness at that point. We
use the median of the myelin thickness values to represent the myelin
thickness for a cross-section of the myelin instance. The g-ratio of the
myelinated axon as $g=\frac{d_{e}}{d_{e}+d_{m}}$, where $d_{e}$ is the
equivalent diameter of the intra-axonal space and $d_{m}$ is the myelin
thickness.

###  2.3 Graphical User Interface

The gACSON user interface is shown in Figure 6. The program menu panel is used
to load image files or import images from the Matlab workspace using file
menu. Images from the file list panel are displayed using the display popup
menu of the image view panel. The user can select the overlay pop-up menu to
visualize an overlay image, such as a segmentation overlaid with the
corresponding intensity image, or change the transparency of the overlay image
by moving the Transparency slider. The toolbar contains buttons to zoom in/out
or display the value of a selected voxel. The toolbar’s Ruler button measures
the Euclidean distance between any two points selected by the user on the
displayed image. The toolbar also contains a multiplanar view for 3D image
files. The software displays 2D slices, in the x-y plane, of 3D-EM data,
allowing the user to browse between 2D slices along the z-axis.

####  2.3.1 ACSON segmentation panel

The ACSON segmentation panel sets the parameter values of the gACSON pipeline
to segment the myelin and intra-axonal space of myelinated axons in 3D- EM
data, as described in Section 2.1.

####  2.3.2 Manual segmentation panel

the user to manually segment images and correct the results of an automatic
segmentation method, The user can manually draw a 2D region to either create a
new label or append the region to an existing label, draw a 2D region to
subtract from an existing label, e.g., to correct segmentation boundaries. .
The user can merge multiple components into a single label to correct the
over-segmentation error, i.e., a single object segmented into multiple
components. It is also possible to display the 3D representation of the
surface of a segment selected by the user, which allows 3D evaluation of the
shape of the segment to detect possible under- or over-segmentation errors.

####  2.3.3 Evaluation panel

The evaluation panel quantifies the accuracy of instance segmentation at the
object-level using precision (positive predictive value), recall
(sensitivity), and F1 score (harmonic mean of precision and recall measures).
To calculate these metrics, the user manually marks all true positives
(correctly detected segments), false positives (incorrectly detected
segments), and false negatives (undetected segments) on a cross-section of a
segmentation overlaid on the EM image. This allows users to assess the
validity of a segmentation as a whole without a need for manually segmenting
complete images which would be burdensome and time consuming. Figure 7
demonstrates an expert the performance of an EM segmentation algorithm.

Figure 7: gACSON evaluation panel to assess a segmentation. The expert
evaluates the segmentation at the object level, quantifying the number of
true-positives (TP), false-positives (FP), and false-negatives (FN) to
calculate the precision, recall, and F1 score. We have applied the evaluation
panel of gACSON to validate our DeepACSON pipeline for the instance
segmentation of axons in corpus callosum and cingulum in [3].

####  2.3.4 Morphometry panel

to set the parameter values of the proposed morphometry algorithm of Section
2.2. requires the 3D instance segmentation of myelinated axons, both the
intra-axonal space and myelin, to calculate the true cross-sectional
morphology of myelinated axons in 3D. The length of the minor and major axes
of a fitted ellipse, equivalent diameter, eccentricity, myelin thickness, and
g-ratio are measured for each axonal cross-section.

Figure 8: 3D representation of the intra-axonal space of myelinated axons in the somatosensory cortex of a TBI animal. The field of view of the segmented volume is $14\times 14\times 15$ $\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{3}$, which corresponds to $1024\times 1024\times 300$ voxels in $x$, $y$, and $z$ directions, respectively. Note the variability of axonal diameter between myelinated axons. Table 2: The size and resolution of the SBEM volumes. The column Axons provides the number of myelinated axons segmented by gACSON that were longer than $4\text{\,}\mathrm{\SIUnitSymbolMicro m}$.  Rat | Hemisphere | Volume ($\mu m^{3}$) | Resolution ($nm^{3}$) | Voxels | Axons
---|---|---|---|---|---
TBI1 | Contralateral | $14\times 14\times 15$ | $13.75\times 13.75\times 50$ | $1024\times 1024\times 300$ | 41
Ipsilateral | $19.1\times 19.3\times 15$ | $18.33\times 18.33\times 50$ | $1044\times 1054\times 300$ | 52
TBI2 | Contralateral | $14\times 14\times 15$ | $13.75\times 13.75\times 50$ | $1024\times 1024\times 300$ | 67
Ipsilateral | $16.9\times 16.9\times 15$ | $16.5\times 16.5\times 50$ | $1024\times 1024\times 300$ | 18
Sham1 | Contralateral | $17.2\times 17.2\times 15$ | $16.81\times 16.81\times 50$ | $1024\times 1024\times 300$ | 133
Ipsilateral | $14\times 14\times 15$ | $13.75\times 13.75\times 50$ | $1024\times 1024\times 300$ | 57

##  3 Results

In this section, we use gACSON to study axonal morphometry in gray matter
using a novel dataset consisting of six 3D SBEM volumes. Tissue samples were
obtained from the somatosensory cortex of one sham-operated and two TBI rats,
both ipsi- and contralateral to the sham operation or injury. Our aim is to
investigate the persistent morphological changes of myelinated axons after
brain injury.

###  3.1 Datasets

Figure 9: Morphology analysis of a myelinated axon. (a) A myelinated axon with cross-sections perpendicular to its skeleton; myelin is blue and the intra-axonal space is red. (b) The diameter of the intra-axonal space : mean = $0.37\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and std = $0.15\text{\,}\mathrm{\SIUnitSymbolMicro m}$. (c) The myelin thickness along the axonal skeleton: mean = $0.15\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and std = $0.06\text{\,}\mathrm{\SIUnitSymbolMicro m}$. (d) The g-ratio along the axonal skeleton: mean = $0.56$ and std = $0.1$. Table 3: Comparison of myelinated axon morphology between sham-operated and TBI rats. The mean and 95$\%$ confidence interval of each group are denoted as $\mu$ and 95$\%$ CI, respectively. Eq diameter and myl thickness represent equivalent diameter and myelin thickness, respectively. Confidence intervals were calculated with the bias-corrected and accelerated percentile bootstrap method using the Matlab function bootci. The tstat column indicates the Welch t-statistic value for the group difference. The asterisk indicates a statistically significant difference between groups.  Hemisphere |  | Morphology
---
parameter
$\mu_{\text{Sham}}$ | 95$\%$ CI Sham | $\mu_{\text{TBI}}$ | 95$\%$ CI TBI | tstat | p-value
Ipsilateral | Eq diameter | 0.48 | [0.44, 0.54] | 0.42 | [0.39, 0.46] | 2.12 | 0.032*
Minor axis | 0.39 | [0.36, 0.43] | 0.35 | [0.33, 0.38] | 1.68 | 0.088
Major axis | 0.65 | [0.58, 0.75] | 0.56 | [0.51, 0.65] | 1.71 | 0.082
Eccentricity | 0.74 | [0.71, 0.77] | 0.74 | [0.72, 0.75] | 0.17 | 0.854
Myl thickness | 0.15 | [0.14, 0.16] | 0.14 | [0.13, 0.15] | 1.18 | 0.250
G-ratio | 0.61 | [0.59, 0.63] | 0.61 | [0.58, 0.65] | 0.00 | 0.999
Contralateral | Eq diameter | 0.59 | [0.55, 0.65] | 0.55 | [0.51, 0.61] | 1.08 | 0.280
Minor axis | 0.48 | [0.44, 0.52] | 0.44 | [0.40, 0.48] | 1.50 | 0.143
Major axis | 0.77 | [0.71, 0.85] | 0.75 | [0.68, 0.83] | 0.49 | 0.631
Eccentricity | 0.75 | [0.74, 0.77] | 0.77 | [0.76, 0.79] | -1.87 | 0.065
Myl thickness | 0.16 | [0.15, 0.17] | 0.16 | [0.15, 0.17] | 0.78 | 0.436
G-ratio | 0.63 | [0.62, 0.65] | 0.63 | [0.61, 0.65] | 0.20 | 0.843
Table 4: Comparison of the morphology of myelinated axons between contralateral and ipsilateral hemispheres in sham-operated and in TBI rats. See Table 3 for notation. Group |  | Morphology
---
parameter
tstat | p-value
Sham | Eq diameter | 3.17 | 0.005*
Minor axis | 3.38 | 0.003*
Major axis | 2.16 | 0.056
Eccentricity | 0.82 | 0.414
Myl thickness | 2.41 | 0.021*
G-ratio | 1.86 | 0.065
TBI | Eq diameter | 4.43 | 0.001*
Minor axis | 3.59 | 0.001*
Major axis | 3.71 | 0.001*
Eccentricity | 3.00 | 0.002*
Myl thickness | 2.35 | 0.018*
G-ratio | 1.08 | 0.294

###  3.2 Segmentation pipeline

We used BM4D to denoise the six pre-processed 3D SBEM data of the primary
somatosensory cortex (section 2.1.1). For the semantic segmentation of myelin,
we applied the bounded volume growing algorithm, where the initial seeds are
generated randomly (section 2.1.2). We segmented the intra-axonal space of
myelinated axons using the bounded volume growing algorithm, where the initial
seeds were generated using the watershed centroids method (section 2.1.3). We
set the $H$ for the H-maxima transform equal to two, the similarity threshold
of the volume growing method to 0.1, and excluded those connected components
with an area smaller than 100 voxels and larger than one million voxels. After
the semantic segmentation of myelin and the instance segmentation of the
intra-axonal spaces of myelinated axons, we performed the instance
segmentation of myelin (section 2.1.4) to assign myelin to its corresponding
intra-axonal space. Figure 8 shows the 3D rendering of the intra-axonal space
of myelinated axons from a TBI dataset. From the six SBEM volumes, the
pipeline segmented 368 axons that were longer than
$4\text{\,}\mathrm{\SIUnitSymbolMicro m}$ as shown in Table 2.

###  3.3 Morphometry pipeline

We analyzed the morphology of the segmented myelinated axons using our
morphometry algorithm described in section 2.2. We extracted the minor and
major diameters of fitted ellipses, equivalent diameter, eccentricity, myelin
thickness, and g-ratio for every cross-section along all myelinated axons in
all 3D-EM data. Figure 9a shows a myelinated axon with three cross-sections
perpendicular to its axonal skeleton. The cross-sectional diameter of the
intra-axonal space of the myelinated axon is shown along its length, where the
mean is equal to $0.37\text{\,}\mathrm{\SIUnitSymbolMicro m}$, and the
standard deviation is equal to $0.15\text{\,}\mathrm{\SIUnitSymbolMicro m}$
(Figure 9b), indicating considerable variation of axonal diameter along the
axon. Figure 9c shows the myelin thickness along the length of the myelinated
axon where the mean value equals $0.15\text{\,}\mathrm{\SIUnitSymbolMicro m}$
and the standard deviation equals $0.06\text{\,}\mathrm{\SIUnitSymbolMicro
m}$. The g-ratio along the length of the myelinated axon (Figure 9d), where
the mean value is equal to 0.56, and the standard deviation is equal to 0.1.

###  3.4 Morphology of myelinated axons in the somatosensory cortex of sham-
operated and TBI rats

We quantified the morphology of myelinated axons in the somatosensory cortex
to investigate the persistent changes in myelinated axon morphology after TBI.
To do this, we first established a threshold for the length of myelinated
axons and selected for analysis those myelinated axons that were long enough
to traverse approximately one-third of the SBEM volume (longer than
$4\text{\,}\mathrm{\SIUnitSymbolMicro m}$). Thresholding myelinated axons
based on their length excluded those axons that were incorrectly split or
components that were incorrectly labelled as myelinated axons. We represented
an axonal parameter by the median of its cross-sectional values.

To compare the morphology of myelinated axons in sham-operated and TBI
animals, we performed the two-sample Studentized permutation test using the
Welch t-statistic as the test statistic [32].The test was performed separately
for both hemispheres. Because we had only one animal in the sham-operated
group, we pooled the axonal measurements from the two animals in the TBI
group. The alpha threshold for statistical significance was set at 0.05 for
all analyses. The analysis was performed using custom code in Matlab R2018b.
Similarly, we assessed differences between axons in the ipsi and contralateral
hemispheres..

The results of the analysis are shown in Figure 1 and Table 3. For example,
the average equivalent diameter of ipsilateral axons was
$0.48\text{\,}\mathrm{\SIUnitSymbolMicro m}$ (95% CI
[$0.43\text{\,}\mathrm{\SIUnitSymbolMicro m}$
$0.53\text{\,}\mathrm{\SIUnitSymbolMicro m}$] in the sham-operated animal and
$0.48\text{\,}\mathrm{\SIUnitSymbolMicro m}$ (95% CI
[$0.38\text{\,}\mathrm{\SIUnitSymbolMicro m}$
$0.45\text{\,}\mathrm{\SIUnitSymbolMicro m}$] in the injured animals. Our
analysis shows that the equivalent diameter of the intra-axonal space of
myelinated axons was significantly smaller in the TBI animals when the groups
were compared ipsilaterally (Table 3, Figure 1) . We did not find the same
effect on equivalent diameter contralaterally (Table 3, Figure 1). We found no
significant difference between minor and major axis length, eccentricity,
myelin thickness, or g-ratio when comparing sham-operated and TBI rats ipsi-
and contralaterally (Table 3, Figure 1). The comparison between ipsi- and
contralateral hemispheres is shown in Table 4. Equivalent diameter, minor
axis, and myelin thickness were greater in the contralateral hemisphere in
both sham-operated and TBI rats. In addition, eccentricity and major axis were
greater in the contralateral hemisphere of animals after TBI.

###  3.5 Computation time and memory consumption

Table 5 shows the computation time of the image segmentation and morphometry
for the SBEM images used in this study. Using a Hp EliteDesk 800 G4 TWRWS
(Windows 10 (64-bit) operating system, an Intel(R) Core(TM) i7-8700K CPU @3.70
GHz processor, and 64 GB RAM), gACSON averaged 5.7 hours for segmentation and
2.4 hours for morphometry for our SBEM datasets. The differences in
computation time across the datasets are due to differences in the number of
seeds evaluated for BVG during the segmentation and differences in the number
and lengths of axons quantified during the morphometry.

Table 5: The computation time of gACSON segmentation (Seg), seed location are defined using the watershed centroids, and the gACSON morphology (Morph) analysis. Rat | Hemisphere | Seg (h) | Morph (h)
---|---|---|---
TBI1 | Contralateral | 5 | 1.5
Ipsilateral | 3 | 1.6
TBI2 | Contralateral | 4 | 3.3
Ipsilateral | 6 | 0.8
Sham1 | Contralateral | 8 | 4.4
Ipsilateral | 8 | 2.7

##  4 Discussion

In this paper, we present a GUI-based gACSON software for automatic
visualization, segmentation and analysis of myelinated axon morphology in 3D-
EM data. We developed gACSON to provide an extensible open-source platform for
both automatic and semi-automatic segmentation of brain tissue ultrastructure
in 3D- EM data. The source code is available under the MIT license at
https://github.com/AndreaBehan/g-ACSON. To demonstrate the applicability of
the , we performed instance segmentation of myelin and the intra-axonal space
of myelinated axons in six SBEM datasets from the somatosensory cortex of one
sham-operated and two TBI rats.

gACSON provides a framework for automatic segmentation of brain tissue
ultrastructures and a means for manual segmentation and correction of
automatic segmentations. The segmentation pipeline is easy to use even for
inexperienced users. Only a few intuitive parameters need to be set, such as
the similarity threshold or the maximum and minimum volume of a segmented
component. provides several methods for semantic segmentation of myelin that
are applicable in different situations. enables instance segmentation of
myelinated axons and provides several methods for generating the initial seeds
for the BVG algorithm [7], which extend our previous algorithms for axon
segmentation [7, 3].

Axon instance segmentation makes no assumptions about the shape of an axon,
unlike axon and myelin segmentation tools such as AxonSeg [12] or MyelTracer
[19]. AxonSeg and MyelTracer work on 2D- EM data with a tissue section
perpendicular to the main axonal orientation, so that the cross sections of
the axons are approximately circular. Thus, these tools assume that axons are
highly parallel and that there is a main axonal orientation against which the
tissue is cut perpendicularly. However, this assumption is not generally true:
axons can orient, cross, or bend in any direction [3]. Thus, while some axons
can be cut perpendicular to their skeleton, others can be cut obliquely and
along their length. Unlike MyelTracer [19], which separates adjacent myelin
sheaths based on user input, gACSON performs this task automatically. The
gACSON segmentation pipeline enables automatic instance segmentation of
myelin, linking the intra-axonal space to its surrounding myelin sheath. This
is a unique feature compared to existing software tools such as [8, 14, 10,
15, 6, 18].

Several functionalities of gACSON, such as browsing 3D data, proofreading
segmentations, or evaluating and correcting a segmentation result, are not
limited to EM data, but can be used for any 2D or 3D image data. However, the
segmentation and morphometry pipelines are specifically designed for the
analysis of myelinated axons in 3D-EM of gray or white matter. extract the
morphology of the segmented myelinated axons in 3D, allowing accurate analysis
of the shape of the myelinated axons. These measurements are important for
studying the conduction velocity of myelinated axons in electrophysiology [34]
or for studying the characterization of axonal diameter with magnetic
resonance imaging [35]. We illustrated the use of gACSON by analyzing an
unpublished dataset of 3D SBEM volumes from rat somatosensory cortex.

Our morphology analyses showed that the equivalent diameter of ipsilateral
axons was significantly smaller in TBI rats than in sham-operated rats.
Obviously, this result should be interpreted with caution because of the small
number of animals available. However, the result is consistent with the
decreased diameter of myelinated axons in the same animals in both the corpus
callosum and cingulum [3].

## Acknowledgments

We would like to thank Maarit Pulkkinen for her assistance with the animal
handling, Mervi Lindman and Antti Salminen for their help preparing the
samples for SBEM.

## Funding

This research has been funded in part by the Academy of Finland grant
$\\#$316258 (JT) and $\\#$323385 (AS), European Social Fund project S21770
(JT), and Biocenter Finland and University of Helsinki (IB, EJ, and SBEM
imaging).

## Competing interest statement

The authors declare that they have no competing interests.

## Availability of data and material

The source code of gACSON is available at
https://github.com/AndreaBehan/g-ACSON under the MIT License. The code is also
archived at Zenodo DOI: 10.5281/zenodo.3693563. Quantitative axon morphology
measures are attached as a supplementary file to this manuscript. The 3D-EM
data and their segmentations will be made available at
https://doi.org/10.23729/bad417ca-553f-4fa6-ae0a-22eddd29a230.

## References

  * [1] J. S. Phelps, D. G. C. Hildebrand, B. J. Graham, A. T. Kuan, L. A. Thomas, T. M. Nguyen, J. Buhmann, A. W. Azevedo, A. Sustar, S. Agrawal, M. Liu, B. L. Shanny, J. Funke, J. C. Tuthill, W. C. A. Lee, Reconstruction of motor control circuits in adult Drosophila using automated transmission electron microscopy, Cell 184 (3) (2021) 759–774.  doi:10.1016/j.cell.2020.12.013.
URL https://doi.org/10.1016/j.cell.2020.12.013

  * [2] P. C. Nahirney, M.-E. Tremblay, Brain ultrastructure: Putting the pieces together, Frontiers in Cell and Developmental Biology 9 (2021) 187.  doi:10.3389/fcell.2021.629503.
URL https://www.frontiersin.org/article/10.3389/fcell.2021.629503

  * [3] A. Abdollahzadeh, I. Belevich, E. Jokitalo, A. Sierra, J. Tohka, DeepACSON automated segmentation of white matter in 3D electron microscopy, Communications Biology 4 (1) (2021) 179.  doi:10.1038/s42003-021-01699-w.
URL
http://dx.doi.org/10.1038/s42003-021-01699-whttp://www.nature.com/articles/s42003-021-01699-w

  * [4] M. Kleinnijenhuis, E. Johnson, J. Mollink, S. Jbabdi, K. L. Miller, A semi-automated approach to dense segmentation of 3D white matter electron microscopy, bioRxiv (2020).  doi:10.1101/2020.03.19.979393.
  * [5] Z. Yuan, X. Ma, J. Yi, Z. Luo, J. Peng, HIVE-Net: Centerline-aware hierarchical view-ensemble convolutional network for mitochondria segmentation in EM images, Computer Methods and Programs in Biomedicine 200 (2021) 105925.  doi:10.1016/j.cmpb.2020.105925.
URL https://doi.org/10.1016/j.cmpb.2020.105925

  * [6] M. Berning, K. Boergens, M. Helmstaedter, SegEM: Efficient image analysis for high-resolution connectomics, Neuron 87 (6) (2015) 1193–1206.  doi:10.1016/j.neuron.2015.09.003.
URL https://doi.org/10.1016/j.neuron.2015.09.003

  * [7] A. Abdollahzadeh, I. Belevich, E. Jokitalo, J. Tohka, A. Sierra, Automated 3D Axonal Morphometry of White Matter, Scientific Reports 9 (1) (2019) 6084.  doi:10.1038/s41598-019-42648-2.
  * [8] I. Belevich, M. Joensuu, D. Kumar, H. Vihinen, E. Jokitalo, Microscopy image browser: A platform for segmentation and analysis of multidimensional datasets, PLOS Biology 14 (1) (2016) e1002340.  doi:10.1371/journal.pbio.1002340.
URL https://doi.org/10.1371/journal.pbio.1002340

  * [9] I. Belevich, E. Jokitalo, DeepMIB: User-friendly and open-source software for training of deep learning network for biological image segmentation, PLOS Computational Biology 17 (3) (2021) e1008374.  doi:10.1371/JOURNAL.PCBI.1008374.
URL
https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008374

  * [10] M. Helmstaedter, K. Briggman, W. Denk, High-accuracy neurite reconstruction for high-throughput neuroanatomy., Nature Neuroscience (2011) 1081–1088doi:10.1038/nn.2868.
  * [11] K. M. Boergens, M. Berning, T. Bocklisch, D. Bräunlein, F. Drawitsch, J. Frohnhofen, T. Herold, P. Otto, N. Rzepka, T. Werkmeister, D. Werner, G. Wiese, H. Wissler, M. Helmstaedter, WebKnossos: Efficient online 3D data annotation for connectomics, Nature Methods 14 (7) (2017) 691–694.  doi:10.1038/nmeth.4331.
URL http://dx.doi.org/10.1038/nmeth.4331

  * [12] A. Zaimi, T. Duval, A. Gasecka, D. Côté, N. Stikov, J. Cohen-Adad, AxonSeg: Open source software for axon and myelin segmentation and morphometric analysis, Frontiers in Neuroinformatics 10 (AUG) (2016) 1–13.  doi:10.3389/fninf.2016.00037.
  * [13] A. Zaimi, M. Wabartha, V. Herman, P.-L. Antonsanti, C. S. Perone, J. Cohen-Adad, AxonDeepSeg: automatic axon and myelin segmentation from microscopy data using convolutional neural networks, Scientific Reports 8 (1) (2018) 3816.  doi:10.1038/s41598-018-22181-4.
URL http://www.nature.com/articles/s41598-018-22181-4

  * [14] A. Cardona, S. Saalfeld, J. Schindelin, I. Arganda-Carreras, S. Preibisch, M. Longair, P. Tomancak, V. Hartenstein, R. Douglas, TrakEM2 software for neural circuit reconstruction, PLoS ONE 7 (6) (2012) e38011.  doi:10.1371/journal.pone.0038011.
URL https://doi.org/10.1371/journal.pone.0038011

  * [15] S. Saalfeld, A. Cardona, V. Hartenstein, P. Tomancak, CATMAID: collaborative annotation toolkit for massive amounts of image data, Bioinformatics 25 (15) (2009) 1984–1986.  doi:10.1093/bioinformatics/btp266.
URL https://doi.org/10.1093/bioinformatics/btp266

  * [16] D. R. Berger, H. S. Seung, J. W. Lichtman, VAST (Volume Annotation and Segmentation Tool): Efficient manual and semi-automatic labeling of large 3D image stacks, Frontiers in Neural Circuits 12 (October) (2018).  doi:10.3389/fncir.2018.00088.
  * [17] A. Jorstad, J. Blanc, G. Knott, NeuroMorph: A Software Toolset for 3D Analysis of Neurite Morphology and Connectivity, Frontiers in Neuroanatomy 12 (July) (2018) 1–12.  doi:10.3389/fnana.2018.00059.
  * [18] C. Sommer, C. Straehle, U. Kothe, F. Hamprecht, ilastik: Interactive learning and segmentation toolkit, in: Eighth IEEE International Symposium on Biomedical Imaging (ISBI 2011).Proceedings, 2011, pp. 230–233.  doi:10.1109/ISBI.2011.5872394.
  * [19] T. Kaiser, H. M. Allen, O. Kwon, B. Barak, J. Wang, Z. He, M. Jiang, G. Feng, Myeltracer: A semi-automated software for myelin g-ratio quantification, eNeuro 8 (4) (2021) 1–9.  doi:10.1523/ENEURO.0558-20.2021.
  * [20] M. Maggioni, V. Katkovnik, K. Egiazarian, A. Foi, Nonlocal transform-domain filter for volumetric data denoising and reconstruction, IEEE Transactions on Image Processing 22 (1) (2013) 119–133.  doi:10.1109/TIP.2012.2210725.
  * [21] K. S. Sim, J. T. L. Thong, J. C. H. Phang, Effect of shot noise and secondary emission noise in scanning electron microscope images, Scanning 26 (1) (2006) 36–40.  doi:10.1002/sca.4950260106.
  * [22] L. Breiman, Random forests, Machine Learning 45 (1) (2001) 5–32.  doi:10.1023/a:1010933404324.
URL https://doi.org/10.1023/a:1010933404324

  * [23] S. H. Lee, H. I. Koo, N. I. Cho, Image segmentation algorithms based on the machine learning of features, Pattern Recognition Letters 31 (14) (2010) 2325–2336.
  * [24] R. Adams, L. Bischof, Seeded region growing, IEEE Transactions on Pattern Analysis and Machine Intelligence 16 (6) (1994) 641–647.  doi:10.1109/34.295913.
  * [25] J. Canny, A computational approach to edge detection, IEEE Transactions on Pattern Analysis and Machine Intelligence PAMI-8 (6) (1986) 679–698.  doi:10.1109/TPAMI.1986.4767851.
  * [26] L. Najman, M. Schmitt, Watershed of a continuous function, Signal Processing 38 (1) (1994) 99–112.  doi:10.1016/0165-1684(94)90059-0.
URL https://doi.org/10.1016/0165-1684(94)90059-0

  * [27] R. Achanta, A. Shaji, K. Smith, A. Lucchi, P. Fua, S. Süsstrunk, Slic superpixels compared to state-of-the-art superpixel methods, IEEE transactions on pattern analysis and machine intelligence 34 (11) (2012) 2274–2282.
  * [28] M. Hassouna, A. Farag, Robust Centerline Extraction Framework Using Level Sets, in: 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), Vol. 1, IEEE, 2005, pp. 458–465.  doi:10.1109/CVPR.2005.306.
URL http://ieeexplore.ieee.org/document/1467303/

  * [29] A. Abdollahzadeh, A. Sierra, J. Tohka, Cylindrical Shape Decomposition for 3D Segmentation of Tubular Objects, IEEE Access 9 (2021) 23979–23995.  doi:10.1109/ACCESS.2021.3056958.
URL
https://ieeexplore.ieee.org/document/9345673/http://dx.doi.org/10.1109/ACCESS.2021.3056958

  * [30] I. Kharatishvili, J. Nissinen, T. McIntosh, A. Pitkänen, A model of posttraumatic epilepsy induced by lateral fluid-percussion brain injury in rats, Neuroscience 140 (2) (2006) 685–697.  doi:10.1016/j.neuroscience.2006.03.012.
URL https://doi.org/10.1016/j.neuroscience.2006.03.012

  * [31] T. Deerinck, E. Bushong, V. Lev-Ram, X. Shu, R. Tsien, M. Ellisman, Enhancing serial block-face scanning electron microscopy to enable high resolution 3-d nanohistology of cells and tissues, Microscopy and Microanalysis 16 (S2) (2010) 1138–1139.  doi:10.1017/s1431927610055170.
URL https://doi.org/10.1017/s1431927610055170

  * [32] A. Janssen, Studentized permutation tests for non-iid hypotheses and the generalized behrens-fisher problem, Statistics & probability letters 36 (1) (1997) 9–21.
  * [33] A. Gelman, J. Hill, M. Yajima, Why we (usually) don’t have to worry about multiple comparisons, Journal of research on educational effectiveness 5 (2) (2012) 189–211.
  * [34] T. Chomiak, B. Hu, What is the optimal value of the g-ratio for myelinated fibers in the rat CNS? A theoretical approach, PLOS ONE 4 (11) (2009) 1–7.  doi:10.1371/journal.pone.0007754.
URL https://doi.org/10.1371/journal.pone.0007754

  * [35] R. A. Salo, I. Belevich, E. Manninen, E. Jokitalo, O. Gröhn, A. Sierra, Quantification of anisotropy and orientation in 3D electron microscopy and diffusion tensor imaging in injured rat brain, NeuroImage 172 (February) (2018) 404–414.  doi:10.1016/j.neuroimage.2018.01.087.

##

Supplementary Figure 1: The morphology analysis of myelinated axons in the
ipsilateral (a-f) and contralateral (g-l) somatosensory cortex of sham-
operated and TBI rats: (a, g) equivalent diameter, (b, h) minor axis length
(c, i) major axis length, (d, j) axonal eccentricity, (e, k) myelin thickness,
and (f, l) g-ratio. On each box, the central mark indicates the median, and
the bottom and top edges of the box indicate the $25^{th}$ and $75^{th}$
percentiles, respectively. The whiskers extend to the most extreme data points
not considered outliers, and the outliers are plotted individually using the
‘o’ symbol. The asterisk in (a) shows that there is a significant reduction in
the equivalent diameter of the intra-axonal space of myelinated axons,
comparing the sham-operated and TBI groups ipsilaterally.

Supplementary Table 1: Description of the parameters of the segmentation pipeline. Parameter | Section |  Description | Default | Fixed
---|---|---|---|---
$\sigma$ | 2.1.2 |  STD for calculating Gaussian smoothed image | 2 | Yes
$\sigma$ | 2.1.2 |  STD for calculating Laplacian of Gaussian image | 0.5 | Yes
$\sigma$ | 2.1.2 |  STD for calculating gradient magnitude image | 2 | Yes
$\sigma$ | 2.1.2 |  STD for calculating difference of Gaussian image | 5 | Yes
H | 2.1.3 |  H-maxima transform suppresses all local maxima whose height is less than H | 2 | Yes
Similarity threshold | 3.2 |  Similarity threshold in BVG to determine whether a voxel can be appended to the growing region. | 0.1 | No
Min volume | 3.2 |  In BVG to exclude volumes below this threshold | 100 voxels | No
Max volume | 3.2 |  In BVG to exclude volumes above this threshold | 1,000,000 voxels | No
Length threshold | 3.4 |  In morphology analysis to exclude myelinated axons that are shorter than this threshold | 4 $\mu m$ | Yes
Myelin min volume |  In semantic segmentation of myelin to exclude volumes under this threshold | 500 voxels | No