# Inclination-type resonance in two-planet systems

Sotiris Sotiriadis,1 Anne-Sophie Libert,1
1naXys, Department of Mathematics, University of Namur, 8 Rempart de la
Vierge, 5000 Namur, Belgium E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We investigate the inclination-growth mechanisms for two-planet systems during
the late protoplanetary disc phase. In previous works, much attention has been
directed to the inclination-type resonance, and it has been shown that it asks
for high eccentricities to be acquired during the migration of the giant
planets. By adopting eccentricity and inclination damping formulae based on
hydrodynamical simulations (instead of the $K$-prescription), we have carried
out 20000 numerical simulations, where we vary the initial planetary
eccentricities, the migration rate, and the dispersal time of the gas disc.
Our results confirm that highly mutually inclined systems are unlikely to be
produced by an inclination-type resonance of two migrating giant planets.
However, in $\sim 1~{}\%$ of the simulations, inclination-type resonance is
observed, and a dynamical study of the evolutions reveals that the
inclination-type resonance mechanism operates in three cases: i) when the
inner planet reaches the inner cavity of the disc, ii) at moderate to high
eccentricities for faster migration rates, and iii) at low to moderate
eccentricities after a phase of orbital destabilization and re-arrangement.

###### keywords:

planet and satellites: formation – planet-disc interactions – planets and
satellites: dynamical evolution and stability

††pubyear: 2019††pagerange: Inclination-type resonance in two-planet
systems–Inclination-type resonance in two-planet systems

##  1 Introduction

The diversity in the orbits of extrasolar planets can be generated by planet-
planet interactions during the migration in the protoplanetary disc (e.g.,
Moorhead & Adams (2005); Moeckel et al. (2008); Marzari et al. (2010);
Matsumura et al. (2010); Libert & Tsiganis (2011b); Moeckel & Armitage
(2012)). When the migrating planets are captured in mean-motion resonance
(MMR), a growth of their eccentricities is produced, and the system can also
possibly enters an inclination-type resonance (whose resonant angle is a
combination of the mean longitudes and the longitudes of the ascending nodes),
which induces a rapid increase of the inclinations. The inclination-type
resonance was discussed by Thommes & Lissauer (2003) for the 2/1 MMR in two-
planet systems and by Libert & Tsiganis (2009) for captures into other MMRs
(e.g., 3/1, 4/1 and 5/1). Both studies concluded that an inclination-type
resonance is observed in two-planet systems as long as the inner planet is not
very massive, one of the planets acquires high eccentricity ($e>0.4$) and the
eccentricity damping is not very efficient. In the two studies, inclination
growth were observed only for $\tau_{e}/\tau_{II}>0.2$, where $\tau_{e}$ is
the timescale for eccentricity damping and $\tau_{II}$ is the timescale for
Type-II migration.

More recently, Teyssandier & Terquem (2014) investigated, both analytically
and numerically, the conditions for the onset of an inclination-type resonance
for two migrating giant planets trapped in 2:1 MMR. They confirmed that, for a
less massive inner planet, the system cannot enter an inclination-type
resonance if the ratio $\tau_{e}/\tau_{II}$ is smaller than 0.2, even if only
the eccentricity of the outer planet is damped. Therefore, they concluded that
if eccentricity damping rate is much more faster than the migration rate in
protoplanetary disc, i.e., $\tau_{e}/\tau_{II}\sim\\!10^{-2}$ as currently
believed, inclination-type resonances are unlikely to occur during the
convergent migration of two gas giants.

All the previous works made use of a $K$-prescription for the eccentricity
(and inclination) damping, consisting in a first-order approximation
($\dot{e}/e=-K|\,\dot{a}/a|$). The $K$-factor to be considered is currently
not well determined but usually estimated in the wide range 1–100. On the
other hand, formulae for the eccentricity and inclination damping were derived
from hydrodynamical simulations in Bitsch et al. (2013). These explicit
formulae, suitable for $n$-body codes, depend on the planetary eccentricity,
inclination, and mass, as well as on the disc mass.

In Sotiriadis et al. (2017), extensive $n$-body simulations of three giant
planets on quasi-circular and quasi-coplanar orbits in the late stage of the
disc were performed, using the damping formulae for eccentricity and
inclination of Bitsch et al. (2013). It was shown that while the semi-major
axis and eccentricity distributions matched the observed ones for the detected
giant planets, a significant proportion of the systems were found in a three-
dimensional (3D) configuration with a mutual inclination higher than
$10^{\circ}$, a large number of them being strongly affected by a MMR or even
captured in a three-body inclination-type resonance during the migration phase
(see e.g. Libert & Tsiganis (2011a) for an analytical study of the three-body
inclination-type resonance).

Let us note that Libert et al. (2018) have recently shown that the
establishment of an inclination-type resonance can also occur at small to
moderate eccentricities in two-planet systems that have previously experienced
an ejection or a merging and are near a resonant commensurability at the end
of the gas phase. The inclination excitation is produced by a (temporary)
capture in an inclination-type resonance (associated with a vertical critical
orbit along a planar elliptic family of periodic orbits), and the proximity to
a spatial periodic orbit can maintain the system in a 3D configuration for a
long period of time.

The goal of the present work is to study the conditions for the onset of an
inclination-type resonance in two-planet systems during the disc phase,
pursuing the efforts of Thommes & Lissauer (2003), Libert & Tsiganis (2009)
and Teyssandier & Terquem (2014), but making use of the eccentricity and
inclination damping formulae of Bitsch et al. (2013) instead of the simplistic
$K$-prescription. As stated above, a variety of dynamical evolutions were
reported in $n$-body simulations of three planets when adopting these damping
formulae. In doing so, the present study extents the previous works of Thommes
& Lissauer (2003) and Libert & Tsiganis (2009) by also considering the
influence of inclination damping on the establishment of the inclination-type
resonance. Moreover, unlike Teyssandier & Terquem (2014), the study presented
here will not be restricted to the 2:1 MMR only.

The set-up of our numerical explorations is described in Section 2. In Section
3, we present several typical evolutions of two-planet systems during the disc
phase, and carefully describe the inclination-growth mechanisms that operate
to form the 3D system architectures observed in the simulations. An analysis
of the occurrence of the different mechanisms is presented in Section 4.
Finally, a comparison with previous studies is provided in Section 5, where
our results are summarized.

Table 1: Parameters of the seven sets of simulations. For each set, the second column describes the main features of the set, the third one the number of simulations, the fourth one the lifetime of the disc and the last column the integration time of the simulations. See text for further details. Set | Parameters | $N_{Sys}$ |  $T_{disc}$ $[yr]$ |  $t_{stop}$ $[yr]$
---|---|---|---|---
1 | Standard set | 11000 | $\sim\\!1\times 10^{6}$ | $1.5\times 10^{6}$
2 | $e_{0}\in[10^{-3},0.5]$ | 1500 | $\sim\\!1\times 10^{6}$ | $1.5\times 10^{6}$
3 | $2~{}T_{disc}$ | 1500 | $\sim\\!2\times 10^{6}$ | $3.5\times 10^{6}$
4 | $T_{disc}/2$ | 1500 | $\sim\\!0.5\times 10^{6}$ | $1.5\times 10^{6}$
5 |  $2~{}T_{disc}$ and $10~{}\tau_{II}$ | 1500 | $\sim\\!2\times 10^{6}$ | $3.5\times 10^{6}$
6 |  $T_{disc}/2$ and $\tau_{II}/5$ | 1500 | $\sim\\!0.5\times 10^{6}$ | $1.5\times 10^{6}$
7 |  $T_{disc}/2$ and $\tau_{II}/10$ | 1500 | $\sim\\!0.5\times 10^{6}$ | $1.5\times 10^{6}$

##  2 Set-up of the n-body simulations

In this study we follow the orbital evolution of two-planet systems and adopt
the same approach as in Sotiriadis et al. (2017), combining the speed and
efficiency of a $n$-body integrator with a symplectic scheme (SyMBA code), and
the improved modelization of the gas effect promoted by the hydrodynamical
simulations of Bitsch et al. (2013). Our simulations include two fully formed
giant planets that evolve around a Solar-type star in the late stage of the
protoplanetary disc phase. We consider that there is enough gas in the disc to
affect the orbital eccentricity and inclination of both planets, while the
inward migration is applied to the outer planet only.

We have realized 20000 numerical simulations, exploring a variety of initial
conditions for the planetary orbital parameters and planetary mass ratios, as
well as disc masses. Seven different sets of simulations have been performed,
and are presented in Table 1. In each ensemble, six initial disc mass values,
namely 1, 2, 4, 8, 16 and 32 $M_{\rm Jup}$, were considered with equal
proportion. The numerical parameters for the viscosity and the structure of
the gas disc are set to the classical values $\alpha=0.005$ for the Shakura-
Sunyaev parameter and $h=0.05$ for the disc aspect ratio. The initial surface
density profile is $\Sigma\propto r^{-0.5}$ and the disc’s inner and outer
edges are set to $R_{\rm in}=0.05$ AU and $R_{\rm out}=30$ AU. A smooth
transition in the gas-free inner cavity is applied with the hyperbolic tangent
function $\tanh{\frac{(r-R_{\rm in})}{\Delta r}}$, where $\Delta r=0.001$ AU
(Matsumoto et al., 2012).

The first set of simulations (Set 1) is referred to as the standard set. For
this set, the dispersal time for the gas is fixed to $\sim\\!10^{6}$ yr and we
neglect the planet-disc interactions when $dM_{\rm Disc}/dt<10^{-9}$ $M_{\rm
star}/yr$. The rate for Type-II migration is given by the equation:

$\tau_{II}=\frac{2}{3}\alpha^{-1}h^{-2}\Omega_{pl}^{-1}\times\max\left(1,\frac{M_{p}}{(4\pi/3)\Sigma(r_{\rm
pl})r_{\rm pl}^{2}}\right),$ (1)

where $\Omega_{pl}$ is the orbital frequency of the planet and
$(4\pi/3)\Sigma(r_{\rm pl})r_{\rm pl}^{2}$ the local mass of the disc, thus
the robustness of the mechanism is defined by the local mass of the disc in
the vicinity of the planet. For the purpose of studying as many different MMR
captures as possible, the initial semi-major axes for the inner and the outer
planet are chosen randomly from the intervals $a_{1}\in[3,5]$ AU and
$a_{2}\in[8,16]$ AU, respectively. We assume that the random initial
eccentricities and inclinations follow a uniform distribution in the ranges
$[0.001,0.01]$ and $[0.01^{\circ},0.1^{\circ}]$, and the same stands for the
initial phase angles in the range $[0.001^{\circ},360^{\circ}]$. Planetary
masses are chosen also randomly from a log-uniform distribution in the
interval $[0.65,5]$ $M_{\rm Jup}$.

In the second set of simulations (Set 2), we examine whether planets initially
on non-circular orbits result in different 3D configurations after the gas
phase, compared to the standard set. Thus, we keep the same modelization as in
the first set, except that the initial eccentricities of the bodies are chosen
randomly in the range $[0.001,0.05]$. The effect of the lifetime of the
protoplanetary disc is studied in Set 3 and Set 4, where the dispersal time of
the disc is fixed to $\sim\\!2\times 10^{6}$ yr and $\sim\\!0.5\times 10^{6}$
yr, respectively. In the last three sets, we investigate the impact of
different Type-II rates on the final orbital architectures of the planetary
systems. Set 5 is identical to Set 3, apart from the migration rate which is
now 10 times slower. In Set 6 and Set 7, we assume that the lifetime of the
disc is reduced by a factor of 2, and the migration timescales are 5 and 10
times smaller, respectively. For all the different sets, we use an integration
timestep of $dt=0.001$ yr, and the systems are followed until $1.5\times
10^{6}$ yr, except for Set 3 and Set 5 where the simulations are carried out
up to $3.5\times 10^{6}$ yr.

Figure 1: The most common outcome. The giant planets emerge from the gas phase
in a resonant coplanar configuration. Eccentricity damping is faster than the
migration rate for the outer giant planet and $K={\tau_{e}}/{\tau_{II}}$ ratio
remains roughly constant throughout the evolution of the system (bottom left
panel). The planetary masses are $m_{1}=1.17$ $M_{\rm Jup}$ and $m_{2}=2.4$
$M_{\rm Jup}$. Simulation from Set 1. Figure 2: Formation of a 3D system by an
inclination-type resonance occuring just after that the planet closer to the
star reaches the inner cavity of the disc. The inner edge of the disc at 0.5
AU is indicated with the dashed black line. The planetary masses are
$m_{1}=2.66$ $M_{\rm Jup}$ and $m_{2}=2.23$ $M_{\rm Jup}$, and the initial
mass of the disc is 32 $M_{\rm Jup}$ (Set 3). Figure 3: Formation of a 3D
system by inclination-type resonance with ten times faster migration rate (Set
7). The planetary masses are $m_{1}=2.23$ $M_{\rm Jup}$ and $m_{2}=0.68$
$M_{\rm Jup}$. Figure 4: Formation of a 3D system by an inclination-type
resonance at low to moderate eccentricities. The gas is dispersed quite
rapidly (Set 4) and when inclination excitation occurs, the damping effects
are not efficient enough to maintain the planets in coplanar orbits. The
planetary masses are $m_{1}=1.25$ $M_{\rm Jup}$ and $m_{2}=2.20$ $M_{\rm
Jup}$.

##  3 Dynamical evolutions

In this section we present four typical evolutions observed in the
simulations, and focus on the possible dynamical mechanisms producing
inclination increase in migrating two-planet systems.

Fig. 1 illustrates the most common outcome of our simulations, where the two
giant planets emerge from the gas phase in a stable and coplanar resonant
configuration. The system belongs to the standard set (Set 1), with an initial
disc mass of 8 $M_{\rm Jup}$. The planetary masses are $m_{1}=1.17$ $M_{\rm
Jup}$ and $m_{2}=2.4$ $M_{\rm Jup}$, thus the mass ratio is $q\sim\\!0.5$. As
the outer planet migrates, the system enters the 2:1 MMR at approximately
$2\times 10^{5}$ yr and stays deep into the resonance until the end of the
integration. The two resonant angles,
$\theta_{1}=\lambda_{1}-2\lambda_{2}+\varpi_{1}$ and
$\theta_{2}=\lambda_{1}-2\lambda_{2}+\varpi_{2}$, start to librate just after
the resonant capture, as shown by the top panels of Fig. 1 (right column).
Whereas the eccentricity of the inner less massive body increases and reaches
the moderate value $\sim\\!0.12$, the outer planet remains on a quasi-circular
orbit. Although the system is captured in a low-order MMR, the eccentricity
values reached during the evolution will remain low due to the efficiency of
the damping exerted by the gas disc. The bottom panel of the left column shows
the evolution of ${\tau_{e}}/{\tau_{II}}$, the ratio of the eccentricity
damping timescale and the Type-II timescale, for the outer giant planet. We
see that the damping on eccentricity induced by the disc is strong since
${\tau_{e}}$ is approximately four times less than the migration timescale
throughout the disc’s lifetime (i.e., $K\sim\\!4$). In this case, no
inclination-type resonance is observed (see the evolution of the inclination-
type resonant angles in the right column of Fig. 1).

In the next examples, we illustrate the formation of 3D systems through
inclination-type resonance. Fig. 2 shows a simulation from Set 3, where an
inclination-type resonance gets underway due to the presence of a planet
inside of the inner cavity of the disc. While the massive disc (32 $M_{\rm
Jup}$) drives the planets ($q\simeq 1.15$) rapidly towards the central star,
the system enters very early the 2:1 MMR, and the two angles
$\theta_{1}=\lambda_{1}-2\lambda_{2}+\varpi_{1}$ and
$\theta_{2}=\lambda_{1}-2\lambda_{2}+\varpi_{2}$ start to librate about
$0^{\circ}$ (top right panels of Fig. 2). However, the efficient eccentricity
damping ($K\sim 5$, see the bottom left panel of Fig. 2) does not allow the
bodies to acquire very eccentric orbits. Once the giant planet closer to the
star reaches the disc’s inner cavity, the interaction with the gas stops and
its eccentricity starts to increase more. After $\sim\\!5\times 10^{5}$ yr,
the outer planet also reaches the inner edge of the disc where its orbit is
unaffected by the disc. The eccentricities of both bodies are rapidly high
enough to trigger a capture into an inclination-type resonance. The two
inclination-type resonant angles,
$\theta_{i_{1}^{2}}=2\lambda_{1}-4\lambda_{2}+2\Omega_{1}$ and
$\theta_{i_{2}^{2}}=2\lambda_{1}-4\lambda_{2}+2\Omega_{2}$, librate about
$180^{\circ}$, and at $\sim 2\times 10^{6}$ yr the planets reach their maximum
inclinations, both above $10^{\circ}$. The bodies continue to evolve in a
stable and non-coplanar resonant configuration until the end of the
simulation.

Non-coplanar systems are also formed in Set 6 and Set 7, where we artificially
decrease the Type-II migration timescale by a factor of 5 and 10,
respectively. An example with $\tau_{II}/10$ is presented in Fig. 3 (Set 7).
The outer planet, initially at $a_{2}=15.3$ AU, migrates rapidly inwards and,
just before $5\times 10^{4}$ yr, the 2:1 MMR is established with
$\theta_{1}=\lambda_{1}-2\lambda_{2}+\varpi_{1}$ and
$\theta_{2}=\lambda_{1}-2\lambda_{2}+\varpi_{2}$ librating around $0^{\circ}$
and $180^{\circ}$, respectively (top right panels of Fig. 3). When the
eccentricity of the less massive outer giant planet ($m_{1}=2.23$ and
$m_{2}=0.68$ $M_{\rm Jup}$) exceeds $\sim 0.4$, the system enters an
inclination-type resonance, as
$\theta_{i_{1}^{2}}=2\lambda_{1}-4\lambda_{2}+2\Omega_{1}$ and
$\theta_{i_{2}^{2}}=2\lambda_{1}-4\lambda_{2}+2\Omega_{2}$ librate about
$0^{\circ}$ and $180^{\circ}$, respectively. Let us note that the ratio
${\tau_{e}}/{\tau_{II}}$ is kept above $\sim 0.5$ throughout the lifetime of
the disc (i.e., $K<2$, see bottom left panel of Fig. 3). The final orbits of
the planets are eventually mutually inclined with $I_{mut}>20^{\circ}$.

Finally, in Fig. 4, we show an example of the establishment of an inclination-
type resonance at low to moderate eccentricities. The system belongs to Set 4,
wherein the dispersal time of the disc is reduced by a factor of 2. The
initial mass of the disc is 32 $M_{\rm Jup}$ and the planetary masses are
$m_{1}=1.25$ $M_{\rm Jup}$ and $m_{2}=2.20$ $M_{\rm Jup}$. The planets are
driven rapidly to the 2:1 MMR, but the system is destabilized at $\sim 2\times
10^{5}$ yr and orbital re-arrangement occurs. Rapidly the system enters the
5:1 MMR (see the libration of the resonant angle
$\theta_{3}=\lambda_{1}-5\lambda_{2}+2\varpi_{1}+2\varpi_{2}$ in the second
right panel), and subsequent inclination growth occurs at $\sim 4\times
10^{5}$ yr due to an inclination-type resonance (at low to moderate
eccentricities), although no clear libration of the inclination-type resonant
angles can be seen in Fig. 4 (right column). To describe more deeply the
mechanism, we display, in the bottom panel of Fig. 4, the moving average every
$3\times 10^{3}$ yr of the angle $2\omega_{1}$ (black curve) and the moving
average of the inclination of the inner planet $i_{1}$ (red curve, in
logarithmic scale), during the period of inclination growth. A correlation
between the secular libration of $2\omega_{1}$ and the inclination increase is
clearly visible. Since $\theta_{i_{3}^{2}}=\theta_{3}-2\omega_{1}$, an
inclination-type resonance can be seen, in the planetary case, as a Lidov-
Kozai resonance (Lidov, 1962; Kozai, 1962) embedded in a MMR. Thus, the
libration of the inclination-type resonant angles is hidden, in Fig. 4, by the
short period oscillations. This mechanism was previously described in Libert
et al. (2018), to which we refer for more details. The planets eventually come
out from the disc phase with $I_{\rm mut}\sim 30^{\circ}$ and eccentricities
close to 0.25.

In the next section, we discuss the occurence of the dynamical mechanisms
hereabove described.

Figure 5: Eccentricities of the 153 two-planet systems with
$I_{mut}>10^{\circ}$ at the end of the simulation. In the top panel, the color
code refers to the different sets, while in the bottom panel it indicates the
mutual inclination value at the end of the simulation. Three groups of systems
are clearly visible: GR-I at high eccentricity of the inner planet, GR-II at
high eccentricity of the outer planet and GR-III at low to moderate
eccentricities for both planets. See text for a detailed description.

##  4 Orbital configurations after the dispersal of the disc

In this section, we focus on the orbital characteristics of the systems
observed in the simulations with final $I_{mut}$ higher than $10^{\circ}$, and
investigate the different formation mechanisms leading to the inclination
excitation.

The vast majority of the systems formed in the simulations end up in a
coplanar configuration as a result of the strong eccentricity and inclination
damping of the gas disc. Nonetheless, as shown in the previous section,
instabilities and MMR captures can lead the migrating planets away from the
disc midplane. Thereafter, the inclined bodies either return to coplanar
orbits due to the strong damping, or maintain their 3D configuration until the
end of the protoplanetary disc phase.

Table 2: Frequency of the different inclination excitation mechanisms producing the highly mutually inclined systems ($I_{mut}>10^{\circ}$), for each set of simulations. | Set 1 | Set 2 | Set 3 | Set 4 | Set 5 | Set 6 | Set 7
---|---|---|---|---|---|---|---
Number | 11000 | 1500 | 1500 | 1500 | 1500 | 1500 | 1500
GR-I | 30 | 3 | 9 | 0 | 0 | 38 | 11
GR-II | 0 | 0 | 0 | 0 | 0 | 0 | 28
GR-III | 12 | 6 | 2 | 2 | 0 | 6 | 0
Other/Instab | 3 | 0 | 0 | 0 | 0 | 1 | 2

More specifically, the huge majority of the 19413 systems that consist of the
two planets at the end of simulation111In the simulations, 587 planetary
systems suffer from merging or ejection, especially the systems issued from
simulations with reduced migration timescales (Set 6 and Set 7)., are found in
a coplanar configuration. However, for $\sim 5\%$ of the simulations (1050
systems), the giant planets are found with mutually inclined orbits
($I_{mut}\\!>$ $1^{\circ}$), while significantly mutually inclined
architectures with $I_{mut}>10^{\circ}$ represent $\sim 1\%$ of the
simulations (153 systems). Hereafter, we focus on the analysis of the latter
category of systems and try to perceive the different formation scenarios.

The final eccentricities of the 153 highly mutually inclined systems are
examined in Fig. 5. The color code gives additional information on the
different sets in the top panel and the final mutual inclinations in the
bottom panel. We clearly observe three distinct regions of systems in the
$(e_{1},e_{2})$ plane. The first region, denoted GR-I, covers the area for
$e_{inner}>~{}0.25$ and $e_{outer}<~{}0.35$. It mainly consists of systems
that have been captured in an inclination-type resonance (at moderate to high
eccentricities), when the inner planet is interior to the inner edge of the
disc (Set 6), similarly to the evolution shown in Fig. 2. The mean mass ratio
$<q>$ for this region is $\sim$ 1.

The second region, GR-II, is found at $e_{inner}<~{}0.2$ and
$e_{outer}>~{}0.4$. It only consists of systems from the simulations with ten
times faster Type-II migration rate, for which the inclination-type resonance
occurs at moderate to high eccentricities during the evolution in the disc. In
this region, $<q>\sim$ 3, i.e., the inner planet is always more massive than
the outer one, which explains the higher value of the eccentricity of the
outer planet. The evolutionary path of a system that belongs to this group has
been displayed in Fig. 3.

The GR-III region consists of 3D systems that have acquired their inclinations
through an inclination-type resonance at low to moderate eccentricities (see
the example of Fig. 4), after a phase of destabilization of the system and
subsequent orbital re-arrangement.

Figure 6: Ratio of the eccentricity damping timescale over the migration
timescale ${\tau_{e}}/{\tau_{II}}$ vs. the mass ratio, for the inclination-
type resonant systems of GR-II. The ratio is evaluated for the outer planet,
just before the end of the interaction with the disc.

Let us note that a large majority of the 3D planetary systems are found in the
2:1 MMR. It is especially the case for the systems from Set 6 and Set 7, since
a fast migration makes difficult a capture in a high-order MMR (see e.g.
Libert & Tsiganis (2009)). For these systems, higher eccentricity and mutual
inclination values (see the second panel of Fig. 5) are reached during the
evolution in the gas disc. However, as observed in our simulations, an
inclination-type resonance can also possibly follow a capture in 3:1, 4:1, 5:1
or 6:1 MMR (see for instance the evolution in Fig. 4).

Table 2 shows the frequency of the above-mentioned formation mechanisms, for
each set of simulations. We observe that around $59\%$ of the highly mutually
inclined systems acquire their 3D configuration by an inclination-type
resonance with planets inside of the inner cavity of the disc (GR-I). They
principally originate from the standard set as well as Set 6 and Set 7 (faster
migration rates). Around a fifth of the systems belong to GR-II characterized
by inclination-type resonance captures associated with a fast migration rate,
and they are all issued from Set 7. Libration of the inclination-type resonant
angles at low to moderate eccentricities (GR-III) is observed in $18\%$ of the
systems. The GR-III systems emanate from nearly all the simulation sets, but
to a large extent from the standard set. The remaining $4\%$ of simulations
generate inclination increase due to orbital instability (no inclination-type
resonance observed).

According to Table 2, about $0.5\%$ of the planetary systems formed in the
standard set (Set 1) are highly mutually inclined, and $\sim 25\%$ of them do
experience an inclination-type resonance at low to moderate eccentricities.
While modifications in the initial planetary eccentricies (Set 2) and the
lifetime of the disc (Set 3, Set 4 and Set 5) do not significantly increase
the number of systems formed in highly mutually inclined configuration, faster
migration rates clearly boost inclination excitation. Approximately $3\%$ of
the planetary systems are found highly non-coplanar in Set 6 and Set 7.

##  5 Discussion and conclusions

In this work, we investigated the inclination-growth mechanisms for systems
with two giant planets during the late protoplanetary disc phase. We performed
20000 simulations considering Type-II migration and the eccentricity and
inclination damping modelling of Bitsch et al. (2013). Besides adopting
different initial orbital parameters of the planets and disc masses, seven
sets of simulations were carried out with different disc lifetimes and
migration rates. Let us recall that the migration of the planets in Type-II
regime due to angular momentum exchange with the disc is on a similar
timescale as the viscous accretion time. However, the actual Type-II migration
rate is still under debate and may depart from the viscous time (Hasegawa &
Ida, 2013; Duffell & Chiang, 2015; Dürmann & Kley, 2015; Kanagawa et al.,
2018; Scardoni et al., 2020). In particular, the prescription for the
viscosity depends on the $\alpha$ parameter whose appropriate value is
unknown. Considering different migration rates includes this uncertainty.

The vast majority of the two-planet systems formed in our simulations end up
in coplanar orbits, as inclination-type resonance capture is rarely observed.
This is in agreement with previous works. In particular, Thommes & Lissauer
(2003) found that inclination-type resonance cannot be established for a
planetary system in the 2:1 MMR if $K>5$ and the mass ratio $q<2$, since, due
to the strong damping, the planetary eccentricities do not reach sufficiently
high values for the inclination-type resonance to be established. Lee &
Thommes (2009) showed that an inclination-type resonance capture is also
possible for $q>2$ if the Type-II migration rate is slower than the one
adopted by Thommes & Lissauer (2003). This result was extended for other MMRs
by Libert & Tsiganis (2009), showing that inclination-type resonance is only
observed when the inner planet is not very massive, one of the planets
acquires high eccentricity ($e>0.4$) and the eccentricity damping is not too
efficient ($K<5$). Regarding the 2:1 MMR, Teyssandier & Terquem (2014)
confirmed that, for a less massive inner planet, inclination-type resonance
only occurs if ${\tau_{e}}/{\tau_{II}}>0.2$ (i.e., $K<5$). All these works
conclude that, given the current estimation of the $K$-factor to 100,
inclination-type resonance capture is unlikely to occur in protoplanetary
discs.

Using a different approach for the eccentricity and inclination damping based
on hydrodynamical simulations, we have found that the evolution presented in
Fig. 1 showing a planetary system remaining coplanar during the protoplanetary
disc phase is also the typical outcome of our simulations. For comparison with
previous studies using the $K$-factor, our simulations are characterized by
the mean ratio $<{\tau_{e}}/{\tau_{II}}>\sim\\!0.12$. The strong eccentricity
damping explains that the onset of inclination-type resonance at moderate to
high eccentricities is quite uncommon in our study.

However, around $1\%$ of the systems formed in the simulations are non-
coplanar, with a mutual inclination higher than $10^{\circ}$. We deeply
studied their dynamical evolution and identified three groups of systems,
depending on the values of their final eccentricities. GR-I involves almost
half of the systems, and consists of the simulations where an inclination-type
resonance has been triggered interior to the disc’s inner edge. The reason for
inclination excitation is obvious, since any interaction between the planets
and the gas is halted in that region and no more damping is applied to the
planets (Fig. 2 is a typical case from this category).

GR-II consists of the systems where the inclination-type resonance occurs
inside the protoplanetary disc at moderate to high eccentricities ($\sim 20\%$
of the 3D systems, all from Set 7 with ten times faster migration rate, see
example in Fig. 3). To show how the ratio of the eccentricity damping
timescale over the migration timescale affects the onset of the inclination-
type resonance for these systems, we present, in Fig. 6, the ratio
${\tau_{e}}/{\tau_{II}}$, for the outer planet, as a function of the mass
ratio of the planets. We observe that the establishment of an inclination-type
resonance for these systems requires ${\tau_{e}}/{\tau_{II}}>0.5$ (the mean
value is 2.14 in this subset), and a more massive inner giant ($q>2$). Whereas
our simulations adopt a more complex prescription for eccentricity and
inclination damping, our results are in agreement with the work of Teyssandier
& Terquem (2014). For this reason, no system from Set 1 (standard set) was
found in this category.

However, our work highlights that inclination increase can be driven in two-
planet systems by a (temporary) capture in an inclination-type resonance at
low to moderate eccentricities. We reported that, for a third group of nearly
$20\%$ of the 3D systems, GR-III, an inclination-type resonance at low to
moderate eccentricities can excite the inclinations, especially after a phase
of orbital instability and re-arrangement, as observed in Fig. 4. In this
case, the ratio ${\tau_{e}}/{\tau_{II}}$ can be lower than 0.5. This dynamical
mechanism can also operate in three-planet systems after a phase of ejection
or merging, as previously observed in Libert et al. (2018).

Although only observed in $\sim 1\%$ of the simulations, inclination-type
resonance is a possible mechanism for the formation of highly non-coplanar
systems. This percentage is strongly dependent on the parameters of the
damping formulae for eccentricity and inclination (e.g., in our case, the disc
mass considered) and on the Type-II migration timescale which is still in
debate.

## Acknowledgements

This work was supported by the Fonds de la Recherche Scientifique-FNRS under
Grant No. T.0029.13 (ExtraOrDynHa research project). Computational resources
have been provided by the PTCI (Consortium des Équipements de Calcul Intensif
CECI), funded by the FNRS-FRFC, the Walloon Region, and the University of
Namur (Conventions No. 2.5020.11, GEQ U.G006.15, 1610468 et RW/GEQ2016).

## References

  * Bitsch et al. (2013) Bitsch B., Crida A., Libert A.-S., Lega E., 2013, Astronomy & Astrophysics, 555, A124
  * Duffell & Chiang (2015) Duffell P. C., Chiang E., 2015, The Astrophysical Journal, 812, id. 94
  * Dürmann & Kley (2015) Dürmann C., Kley W., 2015, Astronomy & Astrophysics, 574, A52
  * Hasegawa & Ida (2013) Hasegawa Y., Ida S., 2013, The Astrophysical Journal, 774, 146
  * Kanagawa et al. (2018) Kanagawa K. D., Tanaka H., Szuszkiewicz E., 2018, ApJ, 861, 140
  * Kozai (1962) Kozai Y., 1962, Astronomical Journal, 67, 591
  * Lee & Thommes (2009) Lee M., Thommes E., 2009, The Astrophysical Journal, 702, 1662
  * Libert & Tsiganis (2009) Libert A.-S., Tsiganis K., 2009, Monthly Notices of the Royal Astronomical Society, 400, 1373
  * Libert & Tsiganis (2011a) Libert A.-S., Tsiganis K., 2011a, Celestial Mechanics and Dynamical Astronomy, 111, 201
  * Libert & Tsiganis (2011b) Libert A.-S., Tsiganis K., 2011b, Monthly Notices of the Royal Astronomical Society, 412, 2353
  * Libert et al. (2018) Libert A.-S., Sotiriadis S., Antoniadou K. I., 2018, Celestial Mechanics and Dynamical Astronomy, 130, 19
  * Lidov (1962) Lidov M., 1962, Planetary and Space Science, 9, 719
  * Marzari et al. (2010) Marzari F., Baruteau C., Scholl H., 2010, Astronomy & Astrophysics, 514, L4
  * Matsumoto et al. (2012) Matsumoto Y., Nagasawa M., Ida S., 2012, Icarus, 221, 624
  * Matsumura et al. (2010) Matsumura S., Thommes E., Chatterjee S., Rasio F., 2010, The Astrophysical Journal, 714, 194
  * Moeckel & Armitage (2012) Moeckel N., Armitage P., 2012, Monthly Notices of the Royal Astronomical Society, 419, 366
  * Moeckel et al. (2008) Moeckel N., Raymond S. N., Armitage P. J., 2008, The Astrophysical Journal, 688, 1361
  * Moorhead & Adams (2005) Moorhead A., Adams F., 2005, Icarus, 178, 517
  * Scardoni et al. (2020) Scardoni C. E., Rosotti G. P., Lodato G., Clarke C. J., 2020, MNRAS, 492, 1318
  * Sotiriadis et al. (2017) Sotiriadis S., Libert A.-S., Bitsch B., Crida A., 2017, Astronomy & Astrophysics, 598, A70
  * Teyssandier & Terquem (2014) Teyssandier J., Terquem C., 2014, Monthly Notices of the Royal Astronomical Society, 443, 568
  * Thommes & Lissauer (2003) Thommes E. W., Lissauer J. J., 2003, The Astrophysical Journal, 597, 566


#
Imposing Connectome-Derived Topology on an Echo State Network

Jacob Morra1 Department of Computer Science
The Brain and Mind Institute
The University of Western Ontario
London, ON Canada
<EMAIL_ADDRESS>Mark Daley2 Department of Computer Science
The Brain and Mind Institute
The University of Western Ontario
London, ON Canada
<EMAIL_ADDRESS>
###### Abstract

Can connectome-derived constraints inform computation? In this paper we
investigate the contribution of a fruit fly connectome’s topology on the
performance of an Echo State Network (ESN) — a subset of Reservoir Computing
which is state of the art in chaotic time series prediction. Specifically, we
replace the reservoir layer of a classical ESN — normally a fixed, random
graph represented as a 2-d matrix — with a particular (female) fruit fly
connectome-derived connectivity matrix. We refer to this experimental class of
models (with connectome-derived reservoirs) as “Fruit Fly ESNs” (FFESNs). We
train and validate the FFESN on a chaotic time series prediction task; here we
consider four sets of trials with different training input sizes (small,
large) and train-validate splits (two variants). We compare the validation
performance (Mean-Squared Error) of all of the best FFESN models to a class of
control model ESNs (simply referred to as “ESNs”). Overall, for all four sets
of trials we find that the FFESN either significantly outperforms (and has
lower variance than) the ESN; or simply has lower variance than the ESN.

##  I Introduction

Reservoir Computers (RC) were first introduced as easily-trainable Recurrent
Neural Networks; however, they have since grown to occupy a niche in chaotic
time-series prediction. For this task type, RCs yield best-in-class
performance whilst requiring less training time and memory [1]. In the space
of Machine Learning models, RCs can be thought of as a class of predictors
which resemble biological brains [2]; notably, this resemblance requires some
squinting.

Much of the recent work on Brain-Inspired Artificial Intelligence has been
“top-down”, applying some conceptual understanding of biological intelligence
to a machine learning approach [3, 4, 5]. Here we instead wish to consider a
“mesoscale” approach from early principles; specifically in that we will
derive a particular aspect of our RC architecture from a connectome (a “map”
of brain connections). We call our approach one from “early principles” in
that we opt to select the simplest biological brain which retains a set of
desirable capabilities (for a predictor) — in particular, efficiency and
discriminatory ability. Intuitively, this leads us to consider the olfactory
system neural network of the fruit fly.

In considering olfaction alone, the fruit fly (Drosophila Melanogaster) is
both effective and efficient. For the former, it has to be effective at odor
discrimination for survival and mating: among other behaviours, it can
selectively target odors from fruit metabolism [6]; it can also classify
pheromone odors of potential mates to determine their age [7]. For the latter,
it also needs to be efficient at odor discrimination (again, for its own
survival): as an example, it identifies new odors in one-shot, “tagging” new
scents with an accuracy similar to those it has observed hundreds of times
([8, 9]).

How can we insert a connectome into a Reservoir Computer? The key bridging
factor here is the RC’s reservoir: a random, fixed graph of nodes which
transforms inputs into a non-linear higher dimensional space. Since the
reservoir can be represented as a 2-d matrix, here we can generate and insert
a connectivity matrix (see Methodology) from the connectome in its place. For
clarity, since we are only concerned with imposing connectome topology, we
will only consider the non-spiking class of RCs: “Echo State Networks” (ESNs).
The basic ESN architecture is illustrated below (Figure 1). To train an ESN,
discrete time series input is fed into the network and allowed to “echo”
around the reservoir; ultimately, a simple linear regression (with an L2
regularization penalty) minimizes the error between the target and predicted
outputs. For a more formal description of the ESN, please see “Methodology”.

Figure 1: The Echo State Network architecture (adapted from [10], Fig. 1).

In our previous work, we created three connectivity matrices: one for each of
the major ROIs (regions of interest) of the olfactory system — the antenna
lobe, the mushroom body calyx, and the lateral horn. We then implanted these
matrices into a Multi-layer Perceptron (MLP) with three hidden layers (one for
each ROI); here the ROI neuron counts were used to specify the number of nodes
in each corresponding hidden layer, and the connectivity matrices were used as
the weight vectors between layers [11]. This approach was limiting for two
reasons: First, the discrete structure of the hidden layers prevented within-
layer connections (thus many weights were not included); Second, because
hidden layer weights would change during training, the connectome-derived
weights were not preserved. The ESN addresses both of these issues: First, the
entire connectivity matrix (i.e. as opposed to three separate connectivity
matrices) is preserved within the reservoir; Second, the reservoir weights are
fixed (only the readout layer of the network is trainable).

Our contribution in this work is a set of performance comparisons between a
control model ESN and an experimental model ESN — one which has had its
reservoir topology constrained by that of the fruit fly connectome; we call
this latter model the “Fruit Fly ESN” (FFESN). Ultimately, we predict that the
FFESN will achieve superior performance to the control model on one or more
variants of a particular time-series prediction task (see Methodology).

In the next section (“Methodology”) we will further describe the process of
generating a connectivity matrix from the fruit fly connectome; we will also
formally describe the ESN, the time series prediction task (and variants)
considered, our particular model classes, and our process for selecting
optimal hyperparameters. In the “Results and Discussion” section, we will
identify the best selected models (from hyperparameter optimization) and
present and discuss the performance of these models. Finally, in the
“Conclusion” section we will summarize the results, present the project
limitations, and highlight directions for future work.

##  II Methodology

###  II-A Generating the Connectivity Matrix

Here we describe the process of creating our connectivity matrix, which is
used to impose a topological constraint on our ESN reservoir by forcing a
particular, connectome-derived structure, as opposed to a randomly-generated
one. For the connectome we consider Janelia’s hemibrain [12]. The hemibrain is
the largest and most complete fruit fly connectome reconstruction to date. It
is an Electron-Microscopy-derived 3-d volume (from a particular female
subject) which comprises over 25,000 neurons, each grouped by a specific cell
type and in their respective brain regions of interest. To retrieve the
hemibrain data we first access the API through a generated key; next, we query
the hemibrain database to select all neurons and connections within the
olfactory system (the antenna lobe, mushroom body calyx, and lateral horn); we
then store all neuron body IDs, all connections (i.e. from each neuron body ID
to all other neuron body IDs it is connected to), and all “weights” between
neurons. Here, weight is the number of synapse-synapse connections between two
neurons (as originally described in [12]). From this point, we select only
those neurons which reside in the right lateral horn: this is a simplification
which we have opted for in order to reduce training times of our control and
experimental models.

With all connectivity information (neuron body IDs plus the weights between
all other neuron body IDs) for the right lateral horn in hand, we generate the
connectivity matrix (see Figure 2). This $n$x$n$ matrix — where $n$ is the
number of right lateral horn neurons — identifies all connections from a
particular body ID (indicated by the row index, as in Figure 2) to a
particular body ID (indicated by the column index). As previously mentioned,
the cell values in this matrix correspond to the number of synapse-synapse
connections between neurons (the weights). Under an assumption that each
neuron is “connected to itself”, the diagonal of the weight matrix is filled
with ones. In Figure 2 we illustrate a small subsection of the (normalized)
connectivity matrix.

###  II-B The Echo State Network

Here let us formally describe the network dynamics of the ESN. Given a
discrete-time input signal $\mathbf{u}(n)\in\mathbb{R}^{N_{u}}$ and target
signal $\mathbf{y}(n)\in\mathbb{R}^{N_{y}}$, the goal is to learn
$\hat{\mathbf{y}}(n)\in\mathbb{R}^{N_{y}}$ such that an error measure
$E\bigg{(}\hat{\mathbf{y}}(n),\mathbf{y}(n)\bigg{)}$ is minimized [2]. The
recurrent reservoir dynamics (as previously specified in [2]) can be described
by the following:

$\tilde{\mathbf{x}}(n)=\mathit{tanh}\bigg{(}\mathbf{W}^{in}[1;\mathbf{u}(n)]+\mathbf{Wx}(n-1)\bigg{)}$
(1) $\mathbf{x}(n)=(1-\alpha)\mathbf{x}(n-1)+\alpha\tilde{\mathbf{x}}(n)$ (2)

Here $[\cdot;\cdot]$ is a column vector; $\mathbf{W}^{in}$ is a weight vector
between the input and reservoir layers
($\mathbf{W}^{in}\in\mathbb{R}^{N_{x}\times(1+N_{u})}$);

Figure 2: A subset of the connectivity matrix used for this work which
illustrates the number of synapse-synapse connections (normalized) between
neurons in the fruit fly connectome (the hemibrain). Blank entries are zero-
valued.

$\mathbf{W}$ is the reservoir weight vector
($\mathbf{W}\in\mathbb{R}^{N_{x}\times(N_{x})}$); $x(n)$ contains the
reservoir activations ($x(n)\in\mathbb{R}^{N_{x}}$); $\tilde{\mathbf{x}}(n)$
contains the updated activations; $\alpha$ is the leaking rate,
($\alpha\in(0,1]$).

The ESN is trained only at the readout (output) layer. The readout layer is
defined as:

$\mathbf{y}(n)=\mathbf{W}^{out}[1;\mathbf{u}(n);\mathbf{x}(n)]$ (3)

We can rewrite this in matrix notation:

$\mathbf{Y}=\mathbf{W}^{out}\mathbf{X}$ (4)

A standard approach is to solve this system with L2-penalized regression
(Ridge Regression):

$\mathbf{W}^{out}=\mathbf{Y}\mathbf{X}^{T}\bigg{(}\mathbf{X}\mathbf{X}^{T}+\lambda\mathbf{I}\bigg{)}^{-1}$
(5)

Here $\lambda$ is the regularization coefficient [2].

###  II-C Time Series Prediction

For a time series prediction task we consider the Mackey-Glass dataset (Figure
3) [13]. This chaotic, 1-d time series (describing physiological control
systems) is a conventional benchmarking dataset for comparison of ESNs [14].
We generate this data from the following ODE [15]:

$\frac{dx}{dt}=\beta x(t)+\frac{\gamma x(t-\tau)}{1+x(t-\tau)^{10}}$ (6)

Here, $x(t)$ is the value of the Mackey-Glass time series at time $t$, and
$[\gamma,\beta,\tau]$ are some real-valued constants. The objective of the
time series prediction task for a series of size $T$, given $x(t)$ for $t$
$\in$ $[0,1,...,m-1]$, is the correct prediction of $x(t)$ for $t$ $\in$
$[m+r,m+r+1,...,T]$. Here $r$ controls how far into the future our model
predicts; $m\in\mathbb{Z^{+}}$ and $r\in\mathbb{Z^{+}}$.

Consider a subset of the Mackey-Glass dataset: specifically $x(t)$ for
$t\in[0,T]$. We define the training input for this dataset (on a time series
prediction task) as the subset $[0,A],A<<T$; here we consider $A=300$ and
$A=900$. For these two training input sizes we will also select from two
train-validate split variants: we will refer to these as “Variant A” and
“Variant B”. The Variant A (as in [16]) train-validate split follows a train-
validate-train-validate ordering whereby the training input time series is
followed by the validation input time series, and then the training output
time series is followed by the validation output time series. Here the inputs
and outputs for each set must be the same sizes; moreover, we select the
majority of the data for training (input and output). The Variant B train-
validate split follows a train-train-validate-validate ordering, where
training output follows training input and validation output follows
validation input. Figure 4 illustrates these split variants. A distinction
between the two variants lies in their perceived difficulty; specifically, for
a particular training input size $A$, we note that Variant A requires
prediction of values which are far into the future as compared to Variant B —
we therefore would expect poorer performance from both models for the former
variant.

###  II-D Two Classes of Models

We consider two classes of models: the Fruit Fly Echo State Network (FFESN)
and the Echo State Network (ESN). Both are identical with the exception of the
topology of their reservoirs. The ESN (as previously defined in subsection B)
serves as a control model: it has a fixed, random reservoir (a 2-d matrix)
with its values drawn from a uniform distribution over $[0,1)$. As in [16], a
mask of 0’s is applied to impose a sparsity of 20%; the locations of 0
elements in the mask are determined by comparing the designated reservoir
density ($0.2$) with a 2-d array filled with values sampled from a uniform
distribution over $[0,1)$. To ensure that our ESN reservoir is fixed through
all iterations, we enforce a random seed; for added measure, we also store,
re-assign, and display this array after each iteration (i.e. for all
predictions). The FFESN (our experimental model) has a fixed reservoir derived
from a connectivity matrix, as described in subsection A. To implement all
models we use the easyesn package [16].

###  II-E Hyperparameter Selection

Given the FFESN and ESN classes of models, for each of the training input
sizes and train-validate split variants (described in subsection C of this
section) we wish to identify the best-performing models (of each class) for
each task. We therefore optimize the hyperparameters $[\lambda,\alpha]$ for
each task-specific model class using a grid search with $\lambda\in[0,1]$,
$\alpha\in(0,1]$. The leaking rate, $\alpha$, determines the contribution of
the previous and current reservoir activations on the current reservoir
activation (see Equation 2). The regularization coefficient, $\lambda$,
determines the magnitude of the Ridge Regression penalty (see Equation 5).
Notably, we have opted out of selecting for other hyperparameters such as the
spectral radius as these would alter our experimental and control reservoir
topologies.

Figure 3: The “Mackey-Glass” chaotic time series.

For each training input size (900, 300) and train-validate split variant (A,
B), we train and validate a population of models with varying hyperparameters.
From this population of models we then select the best performers and choose
their specific combination of hyperparameters for each class of model (FFESN,
ESN). For a given split variant and training input size, we then repeatedly
train and validate our chosen models 50 times in order to generate 50 time
series predictions. Finally, for a performance metric we consider the Mean-
Squared Error (MSE) between each best-model prediction time series and the
ground truth time series; we report a 95% Confidence Interval (CI) around each
MSE (see Table I). The Mean-Squared Error is defined as the following:

$\text{MSE}=\frac{1}{n}\sum_{i=1}^{n}(y(t_{i})-\hat{y}(t_{i}))^{2},$ (7)

where $n$ is the number of time steps; $y(t_{i})$ is the labelled output at
time $t_{i}$, and $\hat{y}(t_{i})$ is the predicted output at time $t_{i}$.

##  III Results and Discussion

###  III-A Time Series Prediction (Variant A)

For Variant A we consider training input sizes $A=900$ and 300. For $A=900$ we
select the first 2000 discrete time series points from the Mackey-Glass
dataset. We assign training input and output sets as $x(t)$ for $t$
$\in[0:900]$ and $[1000:1900]$, respectively; validation input and output sets
are $x(t)$ for $t$ $\in[900:1000]$ and $[1900:2000]$. We retain this train-
validate sequence for $A=300$ and again use 200 time series points in the
validation set (the rest is used for training). After performing
hyperparameter optimization, for $A=900$ we select
$[\alpha,\lambda]=[0.1,0.0024]$ for the FFESN and
$[\alpha,\lambda]=[0.2,4.53$E-$05]$ for the ESN; for $A=300$ we select
$[\alpha,\lambda]=[0.1,0.0024]$ for the FFESN and
$[\alpha,\lambda]=[0.2,0.0024]$ for the ESN.

Figures 5-8 illustrate our best-model predictions and prediction errors on the
validation set for 50 train-validate iterations; here we consider training
input sizes of 900 and 300 (Variant A). In Figure 5 (top subfigure) we observe
(for $A=900$) that the FFESN better captures the two peaks in the original
data. In the bottom subfigure of Figure 5 we observe that the ESN CI is much
wider than the FFESN CI. Intuitively, these results are congruent with the
error-based plots in Figures 6 and 7: In Figure 6, it is apparent that for
$A=900$ the FFESN MSE is much lower than the ESN MSE; the distribution of
errors is also more narrow. In Figure 7, for $A=300$ again we see a tighter
distribution of errors with the FFESN, but also report similar MSEs. Figure 8
elucidates our findings from the directly preceding figures: in summary, for
50 train-validate iterations the FFESN significantly outperforms the ESN for
$A=900$; importantly, the FFESN also yields a fraction of the ESN’s variance
(in MSE) for both training input sizes ($A=900$ and $A=300$).

Figure 4: The two variants of train-validate split considered in this work (in
this case, for a training input size of $A=300$).

Figure 5: Time series predictions (Variant A) for the last 100 steps of a
portion of the Mackey-Glass dataset (for $A=900$ and 300, respectively from
top to bottom). The green and orange dotted lines around the prediction curves
(blue and red lines) denote their 95% Confidence Interval margins.

Figure 6: A histogram illustrating the Mean-Squared Error for 50 iterations of
best-model (FFESN and ESN) predictions on the validation set (Variant A). Here
we consider a training input size of 900 time steps. TABLE I: Summary of
Comparisons between Experimental (FFESN) and Control (ESN) Model Classes

|  $\mathbf{A=900}$ (training input size) |  $\mathbf{A=300}$ (training input size)
---|---|---
| Variant A | Variant B | Variant A | Variant B
| FFESN | ESN | FFESN | ESN | FFESN | ESN | FFESN | ESN
Best $\lambda$ | 0.0024 | 4.53E-05 | 0.00127 | 0.0357 | 0.0024 | 0.0024 | 0.189 | 1
Best $\alpha$ | 0.1 | 0.2 | 0.2 | 0.6 | 0.1 | 0.2 | 0.1 | 0.2
MSE 95% CI | [0.042, 0.043] | [0.046, 0.051] | [0.016, 0.017] | [0.013,0.025] | [0.0048,0.0057] | [0.0022, 0.0067] | [0.0041, 0.0042] | [0.0049, 0.0062]
MSE Variance | 3.76E-08 | 1.46E-06 | 4.89E-08 | 9.03E-06 | 5.21E-08 | 1.24E-06 | 1.36E-09 | 1.03E-07

222

Bold values indicate the (significantly) best-performing or lowest-variance
model for particular train-validate split variants.

###  III-B Time Series Prediction (Variant B)

For the Variant B train-validate split again we consider large ($A=900$) and
small ($A=300$) training input sizes. As in Variant A, for $A=900$ we select
the first 2000 Mackey-Glass data points. In contrast to Variant A, here the
first 900 discrete time steps are for training input and the next 900 steps
are for training output; following this, the next 100 steps and remaining 100
steps are reserved for the validation input and output, respectively. This
train-train-validate-validate split is maintained for $A=300$, with the only
difference being that there are 800 time series points in total (instead of
2000). After performing hyperparameter optimization, for a training input size
of 900 we select $[\alpha,\lambda]=[0.2,0.00127]$ for the FFESN and
$[\alpha,\lambda]=[0.6,0.0357]$ for the ESN; for 300 training inputs we select
$[\alpha,\lambda]=[0.1,0.189]$ for the FFESN and $[\alpha,\lambda]=[0.2,1]$
for the ESN.

Figure 9 provides an illustration of the prediction performance of the best
FFESN and ESN models for 50 train-validate iterations on the validation set
for Variant B. In the figure we observe that for 900 training inputs, whilst
there is no significant difference in prediction performance, the variance of
the FFESN is a fraction of that of the ESN. For a training input size of 300,
we observe that the FFESN outperforms the ESN significantly and also has a
fraction of the MSE variance. Overall, this complements the results for the
Variant A train-validate split, suggesting that the FFESN either outperforms
and has lower variance than the ESN, or simply has lower variance than the
ESN.

Table I summarizes the results from Figures 5-9. It also includes the optimal
hyperparameters for each model class on a particular train-validate split
(i.e. A or B) with a certain training input size (900 or 300). Importantly,
from this table it is evident that for all of the train-validate split
variants and training input sizes we observe a hundred-fold decrease in MSE
variance for the FFESN class of models (as compared to the ESN class of
models). We also observe that the FFESN significantly outperforms the ESN
(outside of two standard deviations of the mean) in the $A=900$, Variant A
trials and $A=300$, Variant B trials.

From our results, the observed reduction in variance from the FFESN is
surprising. It is particularly of importance as

Figure 7: A histogram showing the MSE for 50 iterations of best-model (FFESN
and ESN) predictions on the validation set (Variant A). For this plot we have
used a training input size of $A=300$ time steps.

Figure 8: For Variant A, these boxplots show the Mean-Squared Error for 50
trials of predictions on the validation set for the FFESN and ESN (training
input sizes are $A=900$ and $A=300$).

Figure 9: For Variant B, these boxplots show the MSE for 50 trials of
validation predictions for the FFESN and ESN with training input sizes of 900
and 300.

high variance is a known challenge for ESNs ([17, 18]). To consider a
potential explanation for this outcome, we intuit that the fruit fly brain
retains a topological structure which (likely as a result of selective
evolutionary pressures) favours consistent output responses (i.e. low
variance); this idea is not unlike current theories regarding brain modularity
as a result of evolutionary pressure [19]; however, at this stage it is purely
speculative and remains an open question.

##  IV Conclusion

In summary, we have provided experimental evidence supporting a significant
(100-fold) reduction in performance variance in conjunction with similar or
improved performance metrics after imposing brain-derived topological
restrictions on an Echo State Network. Specifically, after hyperparameter
optimization of the experimental (FFESN) and control (ESN) classes of models
on two train-validate split variants — each with two training input sizes
considered — we have observed that the FFESN either significantly outperforms
the ESN with greatly reduced variance; or yields comparable performance (not
significantly different) to the ESN while maintaining a significant reduction
in variance. With our results we have validated the initial proposition that
the FFESN would outperform the ESN in a time series prediction task — that is,
for a particular dataset, training input size, and split variant.

To extend our previous discussion above, it is apparent that some particular
aspect of the topological structure of the fruit fly olfactory brain is
capable of contributing to a reduction in ESN variance (and a performance
increase); however, we do not yet know which specific feature is responsible.
As a next step, we therefore propose an investigation into three structural
components of the fruit fly neural network and their contributions to ESN
performance (and variance): small-worldness (i.e. the clustering coefficient)
[20], sparsity, and network weights. Here we would create three experimental
models (one for each aspect). As another direction for future work, a study on
the impact of connectome topology on ESN training efficiency would be worth
considering; primarily because the fruit fly is an efficient learner ([8, 9]);
but also in light of recent work in the Reservoir Computing space on “Next
Generation Reservoir Computers”, where [1] reports similar performance to
conventional Reservoir Computers but with reduced computational requirements
(shorter training time, lower data requirements).

Speaking to one of the inherent limitations of this work, it should be noted
here that whilst all of the right lateral horn neurons are present in the
hemibrain, future iterations of this work — wherein we would hope to use the
complete olfactory neural system — would be missing a portion of neurons which
make up the network topology; for example, only 83% of right antenna lobe
neurons are mapped within the hemibrain [12]. To mitigate this we will strive
to only consider ROIs with some tolerated level of trace completion.

Overall, with this work we have presented a particular set of results which
can serve as motivation for looking towards biological brains to yield
performance improvements in a given class of machine learning models;
moreover, from our results we wish to facilitate further lines of comparison
between biologically-motivated, topologically-constrained machine learning
models and their conventional counterparts.

## References

  * [1] D. J. Gauthier, E. Bollt, A. Griffith, and W. A. S. Barbosa, “Next generation reservoir computing,” Nat. Commun., vol. 12, no. 1, p. 5564, Sep. 2021.
  * [2] M. Lukosevicius, “A Practical Guide to Applying Echo State Networks”, [Online]. Available: https://www.ai.rug.nl/minds/uploads/PracticalESN.pdf
  * [3] B. M. Lake, T. D. Ullman, J. B. Tenenbaum, and S. J. Gershman, “Building machines that learn and think like people,” Behav. Brain Sci., vol. 40, p. e253, Jan. 2017.
  * [4] T. P. Lillicrap, A. Santoro, L. Marris, C. J. Akerman, and G. Hinton, “Backpropagation and the brain,” Nat. Rev. Neurosci., vol. 21, no. 6, pp. 335–346, Jun. 2020.
  * [5] G. Bellec et al., “A solution to the learning dilemma for recurrent networks of spiking neurons,” Nat. Commun., vol. 11, no. 1, p. 3625, Jul. 2020.
  * [6] G.A. Wright, “Olfaction: Smells Like Fly Food,” Current Biology, vol. 25, no. 4, PR144-R146, Feb. 2015, doi: 10.1016/j.cub.2014.12.052.
  * [7] G. Smit, “How fruit flies find your food (and mates!),” Noldus. https://www.noldus.com/blog/how-fruit-flies-find-your-food (accessed May 20, 2020).
  * [8] “Creating enoses using fruit flies.” Kurzweil ai. https://www.kurzweilai.net/creating-e-noses-using-fruit-flies (accessed May 20, 2020).
  * [9] S. Dasgupta, C.F. Stevens, and S. Navlakha, “A neural algorithm for a fundamental computing problem,” Science, vol. 358, no. 6364, pp. 793-796, Nov. 2017, doi:10.1126/science.aam9868.
  * [10] H. Jaeger, “The “echo state” approach to analysing and training recurrent neural networks-with an erratum note”. Bonn, Germany: German National Research Center for Information Technology GMD Technical Report, vol. 148.
  * [11] J. Morra and M. Daley, “A Fully-Connected Neural Network Derived from an Electron Microscopy Map of Olfactory Neurons in Drosophila Melanogaster for Odor Classification,” in 2020 IEEE International Conference on Systems, Man, and Cybernetics (SMC), Oct. 2020, pp. 4504–4509.
  * [12] C. Shan Xu et al., “A Connectome of the Adult Drosophila Central Brain”, bioRxiv, Jan. 2020. Accessed: Feb. 01, 2020, doi: 10.1101/2020.01.21.911859. [Online]. Available: https://www.biorxiv.org/content/early/2020/01/21/2020.01.21.911859
  * [13] M. C. Mackey and L. Glass, “Oscillation and chaos in physiological control systems,” Science, vol. 197, no. 4300, pp. 287–289, Jul. 1977.
  * [14] J. R. Maat, N. Gianniotis, and P. Protopapas, “Efficient Optimization of Echo State Networks for Time Series Datasets,” arXiv [cs.NE], Mar. 12, 2019. [Online]. Available: http://arxiv.org/abs/1903.05071
  * [15] C. H. López-Caraballo, I. Salfate, J. A. Lazzús, P. Rojas, M. Rivera, and L. Palma-Chilla, “Mackey-Glass noisy chaotic time series prediction by a swarm-optimized neural network,” J. Phys. Conf. Ser., vol. 720, no. 1, p. 012002, May 2016.
  * [16] R. Zimmermann and L. Thiede. easyesn (2020). Available at: https://github.com/kalekiu/easyesn.
  * [17] Q. Wu, E. Fokoué, and D. Kudithipudi, “On the Statistical Challenges of Echo State Networks and Some Potential Remedies,” Feb. 20, 2018. Accessed: Jan. 23, 2022. [Online]. Available: http://dx.doi.org/
  * [18] A. Millea, “Explorations in Echo State Networks,” 2016, Accessed: Jan. 23, 2022. [Online]. Available: https://www.ai.rug.nl/m̃wiering/Thesis_Adrian_Millea.pdf
  * [19] A. Kelkar and J. D. Medaglia, “Evidence of Brain Modularity,” in Encyclopedia of Evolutionary Psychological Science, T. K. Shackelford and V. A. Weekes-Shackelford, Eds. Cham: Springer International Publishing, 2018, pp. 1–10.
  * [20] M. D. Humphries and K. Gurney, “Network “small-world-ness”: a quantitative method for determining canonical network equivalence,” PLoS One, vol. 3, no. 4, p. e0002051, Apr. 2008.