# Inferring and Learning Multi-Robot Policies from Observations

Pietro Pierpaoli1, Harish Ravichandar2, Nicholas Waytowich3, Anqi Li4,
Derrik Asher3, and Magnus Egerstedt1 *This work was supported by the Army
Research Lab under Grant DCIST CRA W911NF-17-2-0181.1 School of Electrical and
Computer Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA.
<EMAIL_ADDRESS>School of Interactive Computing, Georgia
Institute of Technology, Atlanta, GA 30332, USA.
<EMAIL_ADDRESS>US Army Research Laboratory, Aberdeen Proving
Grounds, MD, 21005 USA {nicholas.r.waytowich<EMAIL_ADDRESS>Department of Computer Science, University of Washington, Seattle, WA 98195,
USA<EMAIL_ADDRESS>
###### Abstract

We present a technique for learning how to solve a multi-robot mission that
requires interaction with an external environment by observing an expert
system executing the same mission. We define the expert system as a team of
robots equipped with a library of controllers, each designed to solve a
specific task, supervised by an expert policy that appropriately selects
controllers based on the states of robots and environment. The objective is
for an un-trained team of robots (i.e., imitator system) equipped with the
same library of controllers, but agnostic to the expert policy, to execute the
mission, with performances comparable to those of the expert system. From un-
annotated observations of the expert system, a multi-hypothesis filtering
technique is used to estimate individual controllers executed by the expert
policy. Then, the history of estimated controllers and environmental states is
used to train a neural network policy for the imitator system. Considering a
perimeter protection scenario on a team of differential-drive robots, we show
that the learned policy endows the imitator system with performances
comparable to those of the expert system.

##  I Introduction

Motivated by the benefits of robots capable to collaborate (e.g., resilience
to faults [1, 2], spatial distribution in the environment [3]), significant
amount of research has investigated the study of hand-specified controllers
designed to solve multi-robot tasks, while overcoming the complexity of
coordination [4]. Prior work has demonstrated that local interaction rules
between the robots can be exploited for the design of task-oriented behaviors,
e.g. [4, 5]. However, these task-oriented controllers cannot be directly
implemented towards the solution of complex, real-world mission, which usually
require an assortment of them (see, for instance [6, 7], for implementations
of this idea).

Figure 1: An overview of the proposed approach to learning multi-robot
policies from observations.

While it might be possible to enhance the expressiveness of individual
controllers through an hand-specified sequential combination of them [8], a
significant amount of human-level intelligence and expertise is required in
order to design such controller compositions that are suitable for complex and
unknown environments [9]. One way to avoid the burden of requiring significant
expertise and design effort is to learn multi-agent policies in an end-to-end
manner (e.g., [10, 11]). However, unlike hand-specified controllers, end-to-
end learned policies require vast amounts of data, are not interpretable, or
lack low-level control guarantees.

In this work, leveraging the relative merits of both well-understood multi-
robot coordinated controllers and learning, we explore a hybrid approach.
Specifically, inspired by the paradigm of learning from demonstration [12,
13], we propose a solution to the problem of sequencing multi-robot, task-
oriented controllers by leveraging observations of an expert system. We
consider an expert system composed by a team of robots performing a
collaborative mission that requires interaction with the environment. Robots
belonging to the expert systems are capable of executing a finite set of low-
level tasks (e.g., rendezvous, flocking, cyclic pursuit), which are selected
by a coordinating, high-level expert policy based on the state of robots and
environment. We consider the expert policy to be a deterministic automated
strategy capable of performing the mission with a satisfactory level of
performance. Note that our approach is agnostic to whether the expert policy
was hand-specified or learned.

We also consider an imitator system, composed by a team of robots, capable of
executing the same list of task-oriented primitive controllers used by the
expert system, but unaware of how to compose them in order to solve the
mission. The imitator system is capable of collecting noisy observations of
both the expert system and the environment as it performances the mission.
Note that we do not assume access to demonstrations with annotations of
controllers being used. Our objective is a procedure that allows the imitator
system to construct a policy capable of achieving expert-like performances in
the same mission.

We propose a two-part framework to learn a controller selection policy from
un-annotated demonstrations (see Fig. 1). First, based on noisy observations
of the expert system’s state, we present a state-estimation inference engine
that estimates the controller executed by the expert system at any given
instant. Second, based on the history of estimated controllers and observed
environment variables, we train a neural network to approximate the expert
policy. We evaluate the performance of our framework on a perimeter defense
task on a team of differential-drive robots. Our results indicate that our
approach is capable of successfully learning the expert policy from only noisy
and un-annotated demonstrations.

###  I-A Related Work

Robotic systems capable of inferring expert behavior have been studied in
variety of contexts, such as human-robot interaction (HRI) [14, 15],
autonomous navigation [16, 17], and interactive learning [18, 19]. For
instance, within the context of HRI, inferring human behavior is often
referred to as intention inference, wherein the robot attempts to identify the
underlying intention or goal of the human partner [20, 21]. In context of
interactive learning, the robot aims to infer what the human partner is
teaching [18]. In these contexts, it is either assumed that labelled training
data can be acquired or that it is sufficient to classify the unlabelled data
into abstract clusters. In our work, we are interested not only in inferring
the expert’s behavior, but also learning from such inferences.

Learning from demonstrations (LfD) is a popular paradigm in robotics that
provides a plethora of techniques aimed at learning from and imitating expert
behavior [12]. A vast majority of techniques, however, rely on labelled
demonstrations and solve a supervised learning problem. A few exceptions to
this assumption include learning reward functions from unlabelled
demonstrations (e.g., [22, 23]), and learning manipulation skills from
unlabelled videos (e.g., [24]).

The methods discussed thus far, however, are limited to single-robot tasks and
scenarios. In the context of multi-robot systems, prior work has explored
identifying joint-intention and shared plans of a group of agents [25, 26].
However, similar to examples in single-robot systems, the inferred information
is not utilized to learn control policies. On the other hand, algorithms that
learn multi-agent policies (e.g., [27, 28, 29]) rely on labelled expert
demonstrations.

As mentioned earlier, a rich body of works exists on task-oriented controllers
for multi agent teams. Coordinated controllers based on weighted consensus
protocol have been used to achieve, for example, flocking, coverage, formation
control, and cyclic pursuit (see [4] and references therein). Examples of
alternative methods include the Null Space Methods [30], Navigation Functions
[31], and Model Predictive Control [32]. Related to our contribution,
solutions to the problem of composing sequences of controllers include formal
methods [33], path planning [6], Finite State Machines [34], Petri Nets [35],
Behavior Trees [36], and Reinforcement Learning [37]. We take inspiration from
a variety of communities in order to develop a unified framework to learn
multi-robot controller selection policies from unlabelled and noisy
observations of an expert system.

We note that the idea of learning how to create sequences from a pre-defined
list of controllers, rather than learning entirely new controllers, is
motivated by the rich literature on coordinated control for multi-agent system
[38]. In addition, although beyond the scope of the current work, the use of
coordinated controllers makes our framework naturally predisposed for fully
distributed implementations. Finally, by focusing on high-level strategies
that compose well-understood controllers, we can interpret the actions
executed by the robots.

##  II Problem Formulation

In this section, we introduce the systems and algorithms relevant to our
approach and provide a formal definition of our problem.

Expert system: Let us consider a team of $N$ robots interacting with an
external environment in order to complete of a mission. We denote by
$x(k)\in\mathcal{X}\subseteq\mathbb{R}^{nN}$ the state of the robots, where
$n$ is number of dimensions necessary to describe the state of each robot. In
addition, we assume the external environment to be a deterministic system with
unknown dynamics, whose configuration is described by a measurable, discrete-
time variable $e(k)\in\mathcal{E}\subseteq\mathbb{R}^{m}$, for all $k>0$. At
each time step, the motion of the robots is governed by a controller, selected
from a library of controllers $\mathcal{L}=\\{f_{1}(x),\dots,f_{M}(x)\\}$,
where each element $f_{j}(\cdot):\mathcal{X}\mapsto\mathcal{X}$ is a
continuous function of the state of the robots. At any given time, all robots
are driven by the same controller, which excludes the possibility of sub-teams
acting in parallel.

Imitator system: In addition to the expert system just described, we assume
the existence of a second team of $N$ robots that are identical to those in
the expert system, whose dynamics are governed by the same library of
controllers $\mathcal{L}$. We refer to this second system as imitator system
and we denote its state by $x_{I}\in\mathcal{X}$.

Measurement model: Both the expert system and imitator system make
observations of the environment under the following observation model:

$\hat{e}(k)=e(k)+w_{e}(k)$ (1)

where $\hat{e}(k)$ is the state of the environment affected by measurement
noise, and $w_{e}(k)$ is a zero-mean Gaussian white-noise process with known
variance. We denote the sequences of environment measurements by
$E^{k}=\\{\hat{e}(k),\hat{e}(k-1),\cdots,\hat{e}(0)\\}$.

In order to learn from the expert, the imitator also collects noisy
observations of the expert’s states as follows

$z(k)=h(x(k))+w_{x}(k)$ (2)

where $h:\mathcal{X}\rightarrow\mathcal{Z}$ is the imitator’s measurement
function, and $w_{x}(k)$ is a zero-mean Gaussian white-noise process with
known variance. We denote the sequences of observations by
$Z^{k}=\\{z(k),z(k-1),\cdots,z(0)\\}$.

Expert policy: A supervising deterministic policy
$F:\mathcal{X}\times\mathcal{E}\mapsto\mathcal{L}$ selects controllers in
$\mathcal{L}$ as function of the states of robots and environment. Then, the
discrete-time dynamics of the robots is

$x(k+1)=F(x(k),\hat{e}(k))+v(k),$ (3)

where $v(k)$ is a zero-mean Gaussian white-noise processes, with known
variance. The policy $F$ selects controllers so that the mission is executed
with a measurable mission performance $\Pi^{F}$. We refer to the system of
robots with state $x$ as expert system, to $F$ as expert’s policy, and to
$\Pi^{F}$ as the expert’s mission performance under the policy $F$.

Imitator policy: We denote the imitator’s learned policy as
$\Phi:\mathcal{X}\times\mathcal{E}\mapsto\mathcal{L}$, and thus the dynamics
of the imitator are

$x_{I}(k+1)=\Phi(x_{I}(k),\hat{e}(k))+v(k).$ (4)

Given access to noisy observations of the expert system and its interaction
with the environment, the objective of our work is to learn that imitator
policy $\Phi$ such that it imitates the expert policy $F$ and achieves
comparable mission performance.

In the next two sections, we describe the main contribution of this paper,
namely a two-stage procedure for the imitator to learn $\Phi$. First, from
noisy observation of the expert’s state and exact state of the environment,
the imitator creates an online estimate of the controllers being executed by
the expert. Then, based on these estimates, we propose a policy approximation
procedure that allows the imitator to solve the task with performance
comparable to the expert’s performance on the same task.

##  III Inferring Controllers from Observations

In this section, we describe the Interactive Multiple Model (IMM) technique
[39], implemented by the imitator to estimate the sequence of controllers
$f_{j}(\cdot)\in\mathcal{L}$ the expert executes, given a history of
observations $Z^{k}$. The choice of this technique is motivated by its
popularity and computational efficiency and we insist that the same objective
could be achieved with alternative estimation schemes.

###  III-A Interactive Multiple Models Estimator

The IMM estimator requires a finite list of controllers (or modes) whose
transitions are described by a Markov process. We denote by $p_{ij}$ the
process probability of switching from mode $f_{j}$ to $f_{i}$, i.e.

$p_{ij}=\text{Prob}(F_{k}=f_{i}\,|\,F_{k-1}=f_{j}).$ (5)

In addition, we define a bank of filters (e.g. Kalman filters, EKF) each
corresponding to a controller in the library $\mathcal{L}$. By combining
expert’s state observations with modes’ transition probabilities, the IMM
computes 1) the probability of $f_{j}$ being active at time $k$, which we
denote by $\mu_{j}(k)$ and 2) the estimate of the expert’s state
$\hat{x}(k|k)$. The estimation technique follows three main steps, which we
briefly discuss here for completeness. We refer the reader to [40] for a
review of IMM implementations and its possible variations. In the following we
denote with $[u]_{\otimes}$ the outer product of vector $u$, defined as
$u\,u^{T}$ and we use the shorthand notation $x(k^{-})$ to denote $x(k-1)$.

##### Interaction

First, mixing probability are computed by propagating each modes’ probability
through the Markov process as

$\mu_{i|j}(k^{-}|k^{-})=\frac{1}{\bar{c}_{j}}p_{ij}\mu_{i}(k^{-})\quad\forall
i,j\in M$

where $\bar{c}_{j}$ is a normalizing factor. Then, we compute the effect of
the modes probabilities on previous state estimates and covariance matrix:

$\displaystyle\hat{x}_{0j}(k^{-}|k^{-})$
$\displaystyle=\sum_{i=1}^{M}\hat{x}_{i}(k^{-}|k^{-})\,\mu_{i|j}(k^{-}|k^{-})$
$\displaystyle P_{0j}(k^{-}|k^{-})$
$\displaystyle=\sum_{i=1}^{M}\Big{[}P_{i}(k^{-}|k^{-})+$
$\displaystyle[\hat{x}_{i}(k^{-}|k^{-})-\hat{x}_{0j}(k^{-}|k^{-})]_{\otimes}\Big{]}\mu_{i|j}(k^{-}|k^{-})$

##### Filtering

From noisy measurements $z(k)$ of the expert, the posterior estimates
$\hat{x}_{j}(k|k)$, covariance matrix $P_{j}(k|k)$, and likelihood
$\Lambda_{j}(k)$ for each mode are computed by applying EKF iterations for
each of the known modes. Then, individual mode probabilities are computed as
follows:

$\mu_{j}(k)=\frac{1}{c}\Lambda_{j}(k)\sum_{i=1}^{M}p_{ij}\mu_{i}(k^{-})$

where $c$ is a normalizing factor.

##### Combination

Final state estimates and convariance matrix are then obtained by combining
values from each filter weighted by the probability of the corresponding mode

$\displaystyle\hat{x}(k|k)$
$\displaystyle=\sum_{j=1}^{M}\hat{x}_{j}(k|k)\mu_{j}(k)$ $\displaystyle
P(k|k)$
$\displaystyle=\sum_{j=1}^{M}\Big{[}P_{j}(k|k)+[\hat{x}_{j}(k|k)-\hat{x}(k|k)]_{\otimes}\Big{]}\mu_{j}(k)$

Every time a new observation is collected, the imitator performs an iteration
of the procedure just described, from which we obtain the probability
distribution for all the controllers in the library. Then, estimated
controller being executed by the expert system a time step $k$ is obtained
from a maximum a posteriori (MAP) estimate

$f^{*}_{k}=\operatorname*{arg\,max}_{j\in\\{1,\dots,M\\}}\,\mu_{j}(k).$ (6)

###  III-B Multi-Robot Coordinated Controllers Inference

We demonstrate the implementation of the IMM estimator to the identification
of $5$ controllers on a team of $5$ robots. Taking inspiration from the
literature on coordinated behaviors for multi-agent systems [4], we assume
each of the controllers in $\mathcal{L}$ to be composed by two terms. The
first term corresponds to a coordinated effect between the robots, which is
represented by a weighted consensus term [38]. The second term is the agent’s
individual objective, used to represent, for example, one or more leaders in
the team. Each controller can then be described as:

$f_{j}(x):=-L_{w,j}(x)+v_{j}(x)$ (7)

where $L_{w,j}:\mathcal{X}\mapsto\mathcal{X}$ is the weighted Laplacian of the
graph $\mathcal{G}_{j}$ describing the interactions between robots required by
controller $j$ [38], while $v_{j}:\mathcal{X}\mapsto\mathcal{X}$, represents a
leader-like controller. We consider $\mathcal{L}_{b}$ to be composed by the
following $5$ behaviors.

Cyclic pursuit:

$u_{i}=\sum_{j\in\mathcal{N}_{i}}R(\theta)\,(x_{j}-x_{i})+(x_{c}-x_{i})$

where $\theta=2r\,\sin\frac{\pi}{N}$, $r$ is the radius of the cycle formed by
the robots, and $R(\theta)\in SO(2)$. The point $x_{c}\in\mathcal{X}$ is the
center of the cycle and $\mathcal{N}_{i}$ is the set of neighbors to agent $i$
specified by $\mathcal{G}_{j}$.

Leader-follower:

$\displaystyle u_{i}$
$\displaystyle=\sum_{j\in\mathcal{N}_{i}}(\|x_{i}-x_{j}\|^{2}-\delta_{f}^{2})(x_{j}-x_{i})\quad
i=2,\dots,4$ $\displaystyle u_{1}$
$\displaystyle=\sum_{j\in\mathcal{N}_{1}}\left((\|x_{1}-x_{j}\|^{2}-\delta_{f}^{2})(x_{j}-x_{1})\right)+(x_{g}-x_{1})$

where $\delta_{f}$ is the desired separation between the robots,
$x_{g}\in\mathcal{X}$ represents the leader’s goal, and subscript $1$ denotes
the leader.

Circular / Wedge / Star formations:

$\displaystyle u_{i}$
$\displaystyle=\sum_{j\in\mathcal{N}_{i}}(\|x_{i}-x_{j}\|^{2}-\bar{\delta}_{ij}^{2})(x_{j}-x_{i})\quad
i=2,\dots,4$ $\displaystyle u_{1}$
$\displaystyle=\sum_{j\in\mathcal{N}_{i}}(\|x_{i}-x_{j}\|^{2}-\bar{\delta}_{ij}^{2})(x_{j}-x_{i})+(x_{g}-x_{1})$

where $\bar{\delta}_{ij}\in\mathbb{R}_{+}$ is the desired separation between
robots $i$ and $j$ designed to achieve circle, wedge, and star shapes, while
$x_{g}$ is the formation leader’s goal.

To illustrate the inference approach, we consider an expert system that
randomly selects one of the controllers described above at Poisson instants of
time. Results from simulations are reported in Fig. 2. As we can observe from
results in Fig. 2, the IMM correctly estimates the expert’s behaviors,
provided a minimum sojourn time on each behavior.

Figure 2: IMM estimation performance for a library of $5$ behaviors. Solid
line represents the behavior being executed by the expert system, while the
dashed line is the IMM estimate. Controllers order: 1) cyclic pursuit, 2)
leader-follower, 3) circle, 4) wedge, and 5) star formations

##  IV Learning to Compose Behaviors

In this section, we introduce our approach to learn the expert policy $F$. We
propose to learn the imitator’s policy, denoted by $\Phi$, that approximates
the expert policy in mapping the current environmental and estimated robot
states onto the appropriate behavior from the controller library
$\mathcal{L}$. As noted in (4), the policy is attempting to capture the
general strategy of choosing the behaviors from the expert as opposed to
encoding a deterministic sequence of behaviors. In order to train such a
policy, we collect the training data from $N$ episodes. Thus, the training
data is given by
$\mathcal{D}=\\{\mathcal{D}^{1},\mathcal{D}^{2},\cdots,\mathcal{D}^{N}\\}$,
where $\mathcal{D}^{i}=\\{{f^{*}_{k}}^{i},\hat{e}^{i}(k),\hat{x}^{i}(k)\\},\
\forall k$ denotes the data associated with the $i^{\text{th}}$ episode.

As can be seen from the notation, the training process does not assume access
to either the true behavior sequence of the expert systems or the true state
of the expert system. That is, the system is trained using inferred quantities
$f^{*}$ and $\hat{x}$ provided by the IMM filter. Once the imitator policy
$\Phi$ is learned, it can be utilized to compose the individual behaviors of
the imitator system in order to solve the mission.

To illustrate our approach, we utilize a neural network (NN) to represent the
imitator’s policy. The parameters of the network (i.e. weights) are trained
using the standard back propagation algorithm. Note that the choice of neural
networks over other models is motivated by its universal approximation
capability [41] and is not central to our framework. Any function
approximation method capable of sufficiently capturing the expert policy would
be applicable.

##  V Experimental Validation

In this section, we implement the policy approximation technique described in
the previous sections, on the observations of an expert system performing a
perimeter protection mission. During expert’s execution of the mission, the
imitator collects observations of the expert’s state and environment’s state,
in order to estimate what controllers are used by the experts for different
configurations of the environment. Then, using data from the estimation
process, we train the imitator’s policy, which is finally executed by the
imitator on the same scenario.

###  V-A Perimeter Defense Scenario

Perimeter protection can be considered both as a stand-alone scenario or as an
abstraction of sub-tasks common to many multi-robot missions. As such, it is
an ideal testing scenarios for multi-agent control protocols [42]. In this
scenario, the expert team is tasked with defending a region while an
adversarial team, representing the environment, tries to intrude. The
intruders’ objective is to reach the innermost region (dark blue), which we
denote by $\mathcal{B}$. The objective of the defenders is to prevent
intruders from reaching the protected area $\mathcal{B}$. To this end, the
policy governing the defending agents reactively chooses the appropriate
controller in $\mathcal{L}$.

Figure 3: Perimeter protection scenario executed on the Robotarium. Intruders
are distinguishable by the yellow disc representing their capture area. The
three circular regions in different shades of blue, from dark to light,
represent $\mathcal{B}$, $r_{r}$, and $r_{R}$ respectively. Text in figure
denotes the controller from the expert policy (left) and the controller from
the learned policy (right).

####  V-A1 Intruders’ strategy

We assume three intruders, with collective state
$e(k)=[e_{1}(k),e_{2}(k),e_{3}(k)]^{T}$, each described at any given time by
three possible modes, namely loiter, attack, retreat. In addition, each
intruder is equipped with a circular protected area of radius
$\varepsilon_{d}$ around it. The intruder’s strategy is illustrated in Fig. 4.
At initial time, intruders start at random points inside region $\mathcal{S}$,
in loiter mode. At Poisson distributed instants of time, a uniformly
distributed random number of intruders are selected and switched to attack
mode. Intruders in attack mode proceed towards the center of region
$\mathcal{B}$. During attacks, if any of the intruders encounter a defending
robot within distance less than $\varepsilon_{d}$, the intruder switches its
state to retreat and moves back to uniformly selected points in region
$\mathcal{S}$.

Figure 4: Expert’s (top) and intruders’ (bottom) policies.

####  V-A2 Expert Policy

The policy employed by the defending robots acts by selecting controllers in
the library $\mathcal{L}$ (see Fig. 4). Here, we assume $\mathcal{L}$ to be
composed of the same controllers described in Section III-B, namely leader-
follower, cyclic-pursuit, and three different formation assembling protocols
wedge, circular, and star. When none of the intruders are within distance
$r_{R}$, defending robots execute the circle formation centered inside region
$\mathcal{B}$. If one intruder is within distance $r_{r}$, defenders switch to
the leader-follower controller. If two or more intruders are within distance
$r_{r}$, defenders first move closer to $\mathcal{B}$ following the circle
formation controller then switch to cyclic pursuit. Then, whether one or more
intruders are between distances $r_{R}$ and $r_{r}$ robots switch to either
star or wedge formations respectively.

###  V-B Experimental Design

In this section, we describe the details and design choices associated with
our experiments.

####  V-B1 Simulation and Experimentation Platform

We realized both a simulated and a physical implementation of the perimeter
protection scenario on the Robotarium [43], which includes a number of
features of the real robots, e.g., differential drive kinematics and collision
avoidance through barrier control functions. In our experiments, we define an
episode as the time between initiations of two consecutive intruders’ attacks.
In total, the data collected in simulation produced the demonstration set
$\mathcal{D}$ consisting of $1000$ episodes.

####  V-B2 Imitator variations

In order to evaluate the performance of the proposed approach, we learn the
imitator policy under two distinct conditions:

  * $\bullet$

$\Phi^{GT}:\mathcal{X}\times\mathcal{E}\mapsto\mathcal{L}$; Policy learned
from the controller sequence generated by the expert policy, ground-truth
robot state, and noisy environmental state
($\mathcal{D}^{i}_{GT}=\\{{f_{k}}^{i},\hat{e}^{i}(k),x^{i}(k)\\},\ \forall
i,k$)

  * $\bullet$

$\Phi^{IMM}:\mathcal{X}\times\mathcal{E}\mapsto\mathcal{L}$; Policy learned
from the sequence of estimated controllers, estimated robot states, and noisy
environmental states
($\mathcal{D}^{i}_{IMM}=\\{{f^{*}_{k}}^{i},\hat{e}^{i}(k),\hat{x}^{i}(k)\\},\
\forall i,k$)

Note that the imitator variation $\Phi^{GT}$ represents the “best-case”
imitation learning scenario in which ground-truth labeled demonstrations are
available.

####  V-B3 Training details

The two policies described above are represented using a multi-layer neural
network with two 32 neuron, fully connected layers with hyperbolic tangent
activation and an output classification layer with softmax activation to
classify between the $5$ control modes. For each policy, the corresponding
demonstration data were split into $80\%$ training ($800$ episodes) and $20\%$
validation data ($200$ episodes).

####  V-B4 Metrics

We measure the performance of each imitator variation described above in terms
of the following metrics:

Inference accuracy: To evaluate the performance of the inference algorithm, we
report the fraction of instants in which our approach successfully inferred
the correct expert behavior.

Learning accuracy: We report training and validation accuracy of the neural
networks associated with both $\Phi^{GT}$ and $\Phi^{IMM}$. The accuracy for
$\Phi^{GT}$ ($\Phi^{IMM}$) is computed by comparing the predicted behaviors
and the behaviors found in the validation set of $\mathcal{D}^{i}_{GT}$
($\mathcal{D}^{i}_{IMM}$).

Imitation performance: To investigate how well the imitator has learned the
expert system’s policy from noisy data, we compute the imitation accuracy of
the learned policy by comparing their predicted behavior to that of the expert
policy.

Mission performance: To evaluate the performance on the actual mission, we
compute the ratio of the number of attacks that were successfully thwarted by
the defending robots to the total number of initiated attacks.

###  V-C Results and Discussion

####  V-C1 How well can we infer expert behavior?

The performance of the IMM in estimating controllers executed by the expert
during the perimeter defense scenario is reported in Fig. 5. We note that, in
this scenario, accurately estimating circle-formation, wedge-formation, and
leader-follower is particularly challenging since, by definition of the expert
policy, these controllers are executed for a limited period of time, and as a
result, the estimator is unlikely to acquire sufficient evidence before the
change in underlying behavior. Nevertheless, the estimator is able to achieve
an overall accuracy of $75\%$.

####  V-C2 How well do the NNs approximate the available data?

Training results for the two different policies are reported in Fig. 6. As one
would expect, training with ground truth data ($\Phi^{GT}$) results in high
learning accuracy in the validation set. In support of our framework, the IMM-
based policy ($\Phi^{IMM}$) is also shown to exhibit similarly high learning
accuracy. Thus, in both cases, the corresponding neural network is capable of
accurately learning the available data.

####  V-C3 How well can we learn the expert policy from noisy data?

In addition to training and validating of the learned policies, we evaluate
how effectively the imitators learn to approximate the expert’s policies
($F$). In particular, as shown in Fig. 7, over a sequence of $120$ episodes,
the behavior chosen by the policy $\Phi^{IMM}$ matched what would have been
chosen by the expert policy $75\%$ of the time on average. Indeed, the
imitation performance is affected by both the inference accuracy and learning
accuracy. We note that the imitation performance mirrors that of inference.
This is likely due to the high learning performance as depicted in Fig. 6.

Figure 5: Inference accuracy in the perimeter protection scenario. Translucent
bars represent the total number of behavior occurrences while each solid bar
represents the number of times the corresponding behavior was correctly
estimated. Figure 6: Learning accuracy on the validation set for the IMM and
GT imitator variations.  Figure 7: Imitation performance in the perimeter
protection scenario. Translucent bars represent the total number of behavior
occurrences in the ground truth data while each solid bar represents the
number of times $\Phi^{IMM}$ predicted the correct behavior.

####  V-C4 What is the effect of the learned policies on mission performance?

Finally, the mission performance of the learned policies was evaluated in the
perimeter defense scenario on a sequence of $120$ episodes. As shown in Fig.
8, we compare the performance $\Pi$ of the expert policy $F$, the imitator’s
policies $\Phi^{GT}$ and $\Phi^{IMM}$, and a random policy. For the random
policy, behaviors are switched uniformly at random, at exponentially
distributed instants of times111The switching intervals obtained by fitting an
exponential probability density function to the switching intervals observed
from the expert policy. Despite the limited performance in estimating three of
the behaviors as shown in Fig. 5, we observe that the imitator’s mission
performance does not deviate significantly from that of the expert. This is
explained by the limited impact of the misclassified behaviors on the mission
performance.

Figure 8:  Mission performance of the expert policy ($F$), imitator policies
($\Phi^{GT}$ and $\Phi^{IMM}$), and the random policy ($\Phi^{RAND}$).

##  VI Conclusions

In this paper we investigated the problem of learning a multi-agent behaviors
selection policy from un-annotated observations of an expert executing a
mission. Continuous observations of the expert’s positions are first converted
by the imitator into a list of estimated controllers. Then, the history of
estimated controllers, robots’ and environment’s configurations are used to
learn the approximated policy. The work presented in this paper is a step
towards the development of observation-based learning frameworks for multi-
robot systems. In particular, our learning-from-observations framework could
be extended to the solution of other forms of multi-robot selection processes,
such as task-assignments or heterogeneous robots team-composition, which are
notoriously computationally hard to solve with traditional techniques.

## References

  * [1] R. K. Ramachandran, J. A. Preiss, and G. S. Sukhatme, “Resilience by reconfiguration: Exploiting heterogeneity in robot teams,” _arXiv preprint arXiv:1903.04856_ , 2019.
  * [2] P. Pierpaoli, D. Sauter, and M. Egerstedt, “Fault tolerant control for networked mobile robots,” in _2018 IEEE Conference on Control Technology and Applications (CCTA)_. IEEE, 2018, pp. 374–379.
  * [3] P. Culbertson and M. Schwager, “Decentralized adaptive control for collaborative manipulation,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 278–285.
  * [4] J. Cortés and M. Egerstedt, “Coordinated control of multi-robot systems: A survey,” _SICE Journal of Control, Measurement, and System Integration_ , vol. 10, no. 6, pp. 495–503, 2017.
  * [5] M. Schwager, D. Rus, and J.-J. Slotine, “Unifying geometric, probabilistic, and potential field approaches to multi-robot deployment,” _The International Journal of Robotics Research_ , vol. 30, no. 3, pp. 371–383, 2011\.
  * [6] S. Nagavalli, N. Chakraborty, and K. Sycara, “Automated sequencing of swarm behaviors for supervisory control of robotic swarms,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 2674–2681.
  * [7] P. Pierpaoli, A. Li, M. Srinivasan, X. Cai, S. Coogan, and M. Egerstedt, “A sequential composition framework for coordinating multi-robot behaviors,” _arXiv preprint arXiv:1907.07718_ , 2019.
  * [8] P. Twu, P. Martin, and M. Egerstedt, “Graph process specifications for hybrid networked systems,” _IFAC Proceedings Volumes_ , vol. 43, no. 12, pp. 65–70, 2010.
  * [9] S. Nikolaidis and J. Shah, “Human-robot cross-training: computational formulation, modeling and evaluation of a human team training strategy,” in _Proceedings of the 8th ACM/IEEE international conference on Human-robot interaction_. IEEE Press, 2013, pp. 33–40.
  * [10] C. Devin, A. Gupta, T. Darrell, P. Abbeel, and S. Levine, “Learning modular neural network policies for multi-task and multi-robot transfer,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 2169–2176.
  * [11] Y. Xiao, J. Hoffman, T. Xia, and C. Amato, “Multi-robot deep reinforcement learning with macro-actions,” _arXiv preprint arXiv:1909.08776_ , 2019.
  * [12] B. D. Argall, S. Chernova, M. Veloso, and B. Browning, “A survey of robot learning from demonstration,” _Robotics and autonomous systems_ , vol. 57, no. 5, pp. 469–483, 2009.
  * [13] H. Ravichandar, A. S. Polydoros, S. Chernova, and A. Billard, “Recent advances in robot learning from demonstration,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 3, 2020.
  * [14] Z. Wang, K. Mülling, M. P. Deisenroth, H. Ben Amor, D. Vogt, B. Schölkopf, and J. Peters, “Probabilistic movement modeling for intention inference in human–robot interaction,” _The International Journal of Robotics Research_ , vol. 32, no. 7, pp. 841–858, 2013.
  * [15] R. Luo, R. Hayne, and D. Berenson, “Unsupervised early prediction of human reaching for human–robot collaboration in shared workspaces,” _Autonomous Robots_ , vol. 42, no. 3, pp. 631–648, 2018.
  * [16] R. Kelley, A. Tavakkoli, C. King, M. Nicolescu, M. Nicolescu, and G. Bebis, “Understanding human intentions via hidden markov models in autonomous mobile robots,” in _Proceedings of the 3rd ACM/IEEE international conference on Human robot interaction_. ACM, 2008, pp. 367–374.
  * [17] A. Jain, A. Singh, H. S. Koppula, S. Soh, and A. Saxena, “Recurrent neural networks for driver activity anticipation via sensory-fusion architecture,” in _2016 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2016, pp. 3118–3125.
  * [18] M. Cakmak and A. L. Thomaz, “Designing robot learners that ask good questions,” in _Proceedings of the seventh annual ACM/IEEE international conference on Human-Robot Interaction_. ACM, 2012, pp. 17–24.
  * [19] R. Dillmann, “Teaching and learning of robot tasks via observation of human performance,” _Robotics and Autonomous Systems_ , vol. 47, no. 2-3, pp. 109–116, 2004.
  * [20] H. C. Ravichandar and A. P. Dani, “Human intention inference using expectation-maximization algorithm with online model learning,” _IEEE Transactions on Automation Science and Engineering_ , vol. 14, no. 2, pp. 855–868, 2016.
  * [21] H. C. Ravichandar, A. Kumar, and A. Dani, “Gaze and motion information fusion for human intention inference,” _International Journal of Intelligent Robotics and Applications_ , vol. 2, no. 2, pp. 136–148, 2018.
  * [22] M. Babes, V. Marivate, K. Subramanian, and M. L. Littman, “Apprenticeship learning about multiple intentions,” in _Proceedings of the 28th International Conference on Machine Learning (ICML-11)_ , 2011, pp. 897–904.
  * [23] Y. Liu, A. Gupta, P. Abbeel, and S. Levine, “Imitation from observation: Learning to imitate behaviors from raw video via context translation,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 1118–1125.
  * [24] P. Sermanet, C. Lynch, Y. Chebotar, J. Hsu, E. Jang, S. Schaal, S. Levine, and G. Brain, “Time-contrastive networks: Self-supervised learning from video,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 1134–1141.
  * [25] Y. Demiris, “Prediction of intent in robotics and multi-agent systems,” _Cognitive processing_ , vol. 8, no. 3, pp. 151–158, 2007.
  * [26] M. Pei, Y. Jia, and S.-C. Zhu, “Parsing video events with goal inference and intent prediction,” in _2011 International Conference on Computer Vision_. IEEE, 2011, pp. 487–494.
  * [27] S. Natarajan, G. Kunapuli, K. Judah, P. Tadepalli, K. Kersting, and J. Shavlik, “Multi-agent inverse reinforcement learning,” in _2010 Ninth International Conference on Machine Learning and Applications_. IEEE, 2010, pp. 395–400.
  * [28] A. Šošić, W. R. KhudaBukhsh, A. M. Zoubir, and H. Koeppl, “Inverse reinforcement learning in swarm systems,” in _Proceedings of the 16th Conference on Autonomous Agents and MultiAgent Systems_. International Foundation for Autonomous Agents and Multiagent Systems, 2017, pp. 1413–1421.
  * [29] R. P. Bhattacharyya, D. J. Phillips, B. Wulfe, J. Morton, A. Kuefler, and M. J. Kochenderfer, “Multi-agent imitation learning for driving simulation,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 1534–1539.
  * [30] G. Antonelli, F. Arrichiello, and S. Chiaverini, “The null-space-based behavioral control for autonomous robotic systems,” _Intelligent Service Robotics_ , vol. 1, no. 1, pp. 27–39, 2008.
  * [31] H. G. Tanner and A. Kumar, “Towards decentralization of multi-robot navigation functions,” in _Proceedings of the 2005 IEEE International Conference on Robotics and Automation_. IEEE, 2005, pp. 4132–4137.
  * [32] A. Richards and J. How, “Decentralized model predictive control of cooperating UAVs,” in _2004 43rd IEEE Conference on Decision and Control (CDC)(IEEE Cat. No. 04CH37601)_ , vol. 4. IEEE, 2004, pp. 4286–4291.
  * [33] H. Kress-Gazit, M. Lahijanian, and V. Raman, “Synthesis for robots: Guarantees and feedback for robot behavior,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 1, pp. 211–236, 2018.
  * [34] A. Marino, L. Parker, G. Antonelli, and F. Caccavale, “Behavioral control for multi-robot perimeter patrol: A finite state automata approach,” in _2009 IEEE International Conference on Robotics and Automation_. IEEE, 2009, pp. 831–836.
  * [35] E. Klavins and D. E. Koditschek, “A formalism for the composition of concurrent robot behaviors,” in _Proceedings 2000 ICRA. Millennium Conference. IEEE International Conference on Robotics and Automation. Symposia Proceedings (Cat. No. 00CH37065)_ , vol. 4. IEEE, 2000, pp. 3395–3402.
  * [36] M. Colledanchise and P. Ögren, “How behavior trees modularize robustness and safety in hybrid systems,” in _2014 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2014, pp. 1482–1488.
  * [37] T. R. Mehta and M. Egerstedt, “An optimal control approach to mode generation in hybrid systems,” _Nonlinear Analysis: Theory, Methods & Applications_, vol. 65, no. 5, pp. 963–983, 2006.
  * [38] M. Mesbahi and M. Egerstedt, _Graph theoretic methods in multiagent networks_. Princeton University Press, 2010, vol. 33.
  * [39] H. A. Blom and Y. Bar-Shalom, “The interacting multiple model algorithm for systems with markovian switching coefficients,” _IEEE transactions on Automatic Control_ , vol. 33, no. 8, pp. 780–783, 1988.
  * [40] E. Mazor, A. Averbuch, Y. Bar-Shalom, and J. Dayan, “Interacting multiple model methods in target tracking: a survey,” _IEEE Transactions on aerospace and electronic systems_ , vol. 34, no. 1, pp. 103–123, 1998.
  * [41] G. Cybenko, “Approximation by superpositions of a sigmoidal function,” _Mathematics of control, signals and systems_ , vol. 2, no. 4, pp. 303–314, 1989.
  * [42] D. Shishika and V. Kumar, “Local-game decomposition for multiplayer perimeter-defense problem,” in _2018 IEEE Conference on Decision and Control (CDC)_. IEEE, 2018, pp. 2093–2100.
  * [43] S. Wilson, P. Glotfelter, L. Wang, S. Mayya, G. Notomista, M. Mote, and M. Egerstedt, “The robotarium: Globally impactful opportunities, challenges, and lessons learned in remote-access, distributed control of multirobot systems,” _IEEE Control Systems Magazine_ , vol. 40, no. 1, pp. 26–44, 2020\.


# Robust universal approach to identify travelling chimeras and synchronized
clusters in spiking networks

Olesia Dogonasheva  Centre for Cognition and Decision Making, National
Research University Higher School of Economics, Moscow, Russia  Dmitry
Kasatkin  Department of Nonlinear dynamics, Institute of Applied Physics RAS,
Nizhny Novgorod, Russia  Boris Gutkin  Centre for Cognition and Decision
Making, National Research University Higher School of Economics, Moscow,
Russia  Group of Neural Theory and LNC2 INSERM U960, École Normale Supérieure
PSL* University, Paris, France  Denis Zakharov<EMAIL_ADDRESS>Centre for
Cognition and Decision Making, National Research University Higher School of
Economics, Moscow, Russia

###### Abstract

We propose a robust universal approach to identify multiple dynamical states,
including stationary and travelling chimera states based on an adaptive
coherence measure. Our approach allows automatic disambiguation of
synchronized clusters, travelling waves, chimera states, and asynchronous
regimes. In addition, our method can determine the number of clusters in the
case of cluster synchronization. We further couple our approach with a new
speed calculation method for travelling chimeras. We validate our approach by
an example of a ring network of type II Morris-Lecar neurons with asymmetrical
nonlocal inhibitory connections where we identify a rich repertoire of
coherent and wave states. We propose that the method is robust for the
networks of phase oscillators and extends to a general class of relaxation
oscillator networks.

Spiking neuronal networks, synchronization, corder parameter, luster
synchronization, chimera state, travelling chimera state, travelling wave

††preprint: APS/123-QED

##  I Introduction

Chimera states in homogeneous networks, characterized by heterogeneous spatio-
temporal patterns containing both synchronous and asynchronous activity
regions, have been a concerted focus of modern nonlinear science and
synchronization theory [1]. Since their recent dyscovery [2] they have been
identified in multiple systems, such as phase oscillator networks (e.g. [3]),
the Stuart-Landau systems [4, 5], and found experimentally in multiple systems
[6, 7, 8, 9, 10, 11]. Studies of chimera states in the neuronal networks,
despite their interest for computational functional importance, are rapidly
developing yet are less numerous ([12]).

In fact, coexistence of synchronous and asynchronous neuronal activity
(chimera states) have been found experimentally accross brain regions and
experimental paradigms [13, 14, 15, 16, 17]. The chimera/multichimera states
in the neuronal networks have been mostly studied for the nonlocal ring
networks with electrical and chemical synapses as well as different neuronal
models (e.g. [18, 19, 20, 21, 22, 23]; reviewed in [12]).

One of the more interesting cases of chimera states are travelling chimeras
where the incoherent domain moves across the network. Such states have been
recently identified across a wide variety of dynamical systems: Kuramoto
networks [24, 25], mechanical systems [26], electromagnetic oscillators [27],
hierarchical topology systems [28], and neuronal networks [29, 30].

Despite the increased interest in chimera states, ability to robustly and
automatically identify the structure of such complex spatio-temporal dynamics
correctly, remains a key challenge. Arguably, without such methodology,
results remain exemplary. Over the decades several approaches have been used.
The simplest – is the visual method that uses rasterplots, snapshots of
instantaneous distributions of the state variable (e.g. neuronal membrane
potentials) and frequency distribution diagrams (e.g. [2, 20]). This approach
gives the network state for an individual set of network parameters and
initial conditions. For a more automatic and generalizable parametric search
for the chimera states, coherence measures have been proposed, like the
Kuramoto order parameter [2], strength of incoherence (together with a
discontinuity measure) [31], and the $\chi^{2}$-parameter [32]. All these
methods have significant drawbacks for relaxation systems and, in particular,
for neuronal networks. The Kuramoto order parameter calculations require well
characterized phases of the network elements; these are difficult to define
correctly given the fast-slow relaxation dynamics of the spiking and bursting
neurons and phase perturbations caused by neuronal interactions. The strength
of incoherence is exquisitely sensitive to two method parameters (the number
of bins and the coherence threshold) that are difficult to select correctly,
especially when the state space includes of a wide variety of synchronous
regimes. The $\chi^{2}$-parameter is more straightforward to compute, yet it
suffer from disadvantages inherent in all of these methods. However, since it
also has some drawbacks, which we will discuss below, one cannot automatically
determine state space regions distinguishing basic dynamical regimes in
neuronal networks.

We introduce a new robust universal approach that allows to identify semi-
automatically chimera/multichimera states, both stationary and travelling,
based on an extension of the $\chi^{2}$-parameter. Our approach is free from
the above-mentioned drawbacks and can be used to identify correctly and
comprehensively dynamical regimes in neuronal networks: global
synchronization, cluster synchronization, travelling waves, chimera states, as
well as asynchronous behavior.

We validate the approach, by giving as a paradigmatic example, calculations of
dynamical regimes of a ring network of type II-excitable Morris-Lecar neurons
[33] with the asymmetrical nonlocal inhibitory connectivity:

$\begin{cases}C\dot{V}_{i}&=I_{app}-g_{Ca}m^{inf}_{i}(V_{i}-E_{Ca})-\\\
&-g_{K}w_{i}(V_{i}-E_{K})-g_{L}(V_{i}-E_{L})+I^{syn}_{i},\\\
\dot{w}_{i}&=\phi(w^{inf}_{i}-w_{i})\cosh\frac{V_{i}-V_{3}}{2V_{4}}\\\
w^{inf}_{i}&=0.5\Bigl{(}1+\tanh\frac{V_{i}-V_{3}}{V_{4}}\Bigr{)}\\\
m^{inf}_{i}&=0.5\Bigl{(}1+\tanh\frac{V_{i}-V_{1}}{V_{2}}\Bigr{)}\end{cases}$
(1)

where $i\in{1..N}$, $N$ is number of neurons. $V_{i}$ is a membrane potential
of $i$-th neuron. $w_{i}$ and $m_{i}$ are the fraction of open K+ and Ca2+
channels, respectively. $E_{K}$, $E_{Ca}$, $E_{L}$ are the reversible
potentials for potassium, calcium and leak channels, respectively. $g_{K}$,
$g_{Ca}$, $g_{L}$ are corresponding conductances. $I_{app}$ is an applied
current. $C$ is a membrane capacity. The neurons are of II type of
excitability.

$I^{syn}_{i}$ is a synaptic current supplied to $i$-th neuron. The inhibitory
GABAA ($\gamma$-aminobutyric acid-A) synapses are given by first order
kinetics [34]:

$\begin{cases}I^{syn}_{i}&=\frac{g_{syn}}{N}\sum^{i+R}_{j=0}{x_{ij}}(V_{R}-V_{i})\\\
\dot{x}_{ij}&=\alpha\frac{1}{1+\exp(-\frac{V_{j}-V_{syn}}{K_{p}})}(1-x_{ij})-\beta
x_{ij}\end{cases}$ (2)

Here $V_{R}$ is a reversal potential, $V_{syn}$ is a threshold, $K_{p}$ is the
synaptic activation, $g_{syn}$ is a synaptic strength, $r=\frac{R}{N}$ is a
connectivity parameter, $R$ is a number of connections, $x_{i}$ is the
proportion of open receptors on the postsynaptic membrane of $i$-th neuron.
The topology of the network is a ring. All connections are unidirectional. To
study network dynamical regimes, we fix all parameters (Table 1) except the
control parameters: the external current $I_{app}$, synaptic strength
$g_{syn}$ and connectivity parameter $r$.

$g_{K}=8$ mS/cm2 |  $E_{K}=-80$ mV |  $V_{1}=-1.2$ mV | $\alpha=1.1$
---|---|---|---
$g_{Ca}=4.4$ mS/cm2 |  $E_{Ca}=120$ mV |  $V_{2}=18$ mV | $\beta=0.19$
$g_{L}=2$ mS/cm2 |  $E_{L}=-60$ mV |  $V_{3}=2$ mV | $N=500$
$\phi=1/25$ |  $C=20$ $\mu$F/cm2 |  $V_{4}=30$ mV | $K_{p}=5$
$V_{syn}=2$ mV |  $V_{R}=-60$ mV |  |
Table 1: The fixed neuronal network parameters

Numerical integration was performed (1-2) using the fixed step Euler method
($100$ $\mu$s). Simulation time was $30$ s.

##  II Identification of travelling chimera states

###  II.1 Adaptive coherence measure

To estimate coherence of a network, one can use $\chi^{2}$-parameter [32]:

$\chi^{2}=\frac{\sigma^{2}_{V}}{\frac{1}{N}\sum^{N}_{i=1}{\sigma^{2}_{V_{i}}}},$
(3)

where $\sigma^{2}_{V}$ is variance of average membrane potential of the
network $V(t)=\frac{1}{N}\sum_{i=1}^{N}V_{i}(t)$:

$\sigma^{2}_{V}=\langle V^{2}(t)\rangle_{t}-\langle V(t)\rangle_{t}^{2},$ (4)

and $\sigma^{2}_{V_{i}}$ is variance of membrane potential of the i-th neuron:

$\sigma^{2}_{V_{i}}=\langle V_{i}^{2}(t)\rangle_{t}-\langle
V_{i}(t)\rangle_{t}^{2}.$ (5)

The parameter $\chi^{2}$ takes the values in [0,1]: $\chi^{2}=0$ corresponds
to asynchronous state and $\chi^{2}=1$ for full synchrony. For the values
$0<\chi^{2}<1$ the chimera states may exist. Compared with the Kuramoto order
parameter and the strength of incoherence, the advantage of this parameter is
that there is no need to calculate the neuronal phases and fine tune the
method’s intrinsic parameters. However, it cannot distinguish travelling
waves, cluster synchronization and chimera states. To solve this problem, we
suggest an extended approach, that assumes the identification coherent states
with multiple synchronous clusters firing with the time lags. We introduce an
Adaptive Coherence Measure (ACM) that based on the solution of the
optimization problem:

$R^{2}=\max_{\mathbf{\Delta t}=(\Delta t_{1},\Delta t_{2},..,\Delta
t_{L})}\chi^{2}(\\{V_{i}(t-\Delta t_{i})\\}_{i=1}^{N}),$ (6)

where $\mathbf{\Delta t}=(\Delta t_{1},\Delta t_{2},..,\Delta t_{L})$ is a
vector of the unique time lags.

An analogous approach was used for the similarity function in [35].

Direct solution of the equation (6) involves multidimensional optimization and
may require significant computational capacity. To avoid it, let us align the
first spikes of all neurons in the network and in this way find the vector
$\mathbf{\Delta t}$. This removes the phase shift between the neuronal
clusters. For bursters, it is sufficient to align the first spikes of the
bursts. Thus, if $R^{2}=1$ and $1<L\ll N$, there exists cluster
synchronization and quantity of synchronous clusters is equal to $L)$. For
example, in the Fig. 1 there are two synchronous antiphase clusters,
$R^{2}=0.9999\approx 1$ and $L=2$ that verify our approach. For comparison,
$\chi^{2}=0.2254$ and does not correctly identify the cluster synchronization
regime. Note that in the trivial case, when the network is globally
synchronized, both $R^{2}$ and $\chi^{2}$ takes the equal value of 1 and $L=1$
which corresponds to a single synchronous cluster.

Figure 1: Rasterplot (A), frequency distribution (B), and snapshot at $t=1000$
ms (C) for the case of two synchronous antiphase clusters ($I_{app}=95$
$\mu$A/cm2, $g_{syn}=1$ mS/cm2, $r=0.9$)

In the case of travelling wave dynamics (Fig. 2), $\chi^{2}=0.0006$ (like for
an asynchronous state) but $R^{2}=1$ and $L=N$. When the voltage traces are
aligned, neurons are fully synchronous, and all of them have a phase shift
relative to each other.

Figure 2: Rasterplot (A), frequency distribution (B), and snapshot at $t=1000$
ms (C) for a travelling wave ($I_{app}=95$ $\mu$A/cm2, $g_{syn}=5$ mS/cm2,
$r=0.2$)

For a static multichimera state (Fig. 3) both approaches give the values of
coherence parameters in $(0,1)$: $\chi^{2}=0.1591$ and $R^{2}=0.8642$. Marking
the neurons by their time lags, we can determine that there are two large
synchronous groups of neurons ($L_{lsg}=2$) and a large population of
asynchronous neurons in the network. To make sure that these clusters are
synchronous one can calculate the ACM only for the neurons that belong to the
clusters. One can distinguish the chimera and multichimera states by searching
for alternating coherent and incoherent domains marked by the time lags.

Figure 3: Rasterplot (A), frequency distribution (B), and snapshot at $t=1000$
ms (C) for the case of a multichimera state ($I_{app}=95$ $\mu$A/cm2,
$g_{syn}=3$ mS/cm2, $r=0.84$)

There also exist a travelling chimera state (Fig. 4) containing two
synchronous clusters for which we can determine the chimera speed
$v_{chi}=0.00878$ neurons per time unit (the method of chimera speed
calculation is considered in the next subsection). Changing to the appropriate
travelling coordinates, we can “stop” the travelling chimera and plot the
frequency distribution diagram (Fig. 4B). Note that our approach gives an
accurate estimation of this chimera state: $R^{2}=0.7779$, $L_{lsg}=2$ in
comparison with $\chi^{2}=0.1269$ that improperly corresponds to a low level
of coherence.

Figure 4: Rasterplot (A), frequency distribution (B), and snapshot at $t=1000$
ms (C) for the case of a travelling multichimera state ($I_{app}=95$
$\mu$A/cm2, $g_{syn}=3$ mS/cm2, $r=0.78$)

Another interesting regime is a travelling chimera state containing domains of
synchronous firing and domains of subthreshold oscillations (Fig. 5) which has
the speed $v_{chi}=0.07121$ neurons per time unit. As it was for the case of
the travelling chimera in Fig. 4 our coherence parameter ACM correctly
characterizes this state $R^{2}=0.5945$ in comparison with $\chi^{2}=0.0908$
that incorrectly shows almost asynchronous behavior.

Figure 5: Rasterplot (A), frequency distribution (B), and snapshot at $t=1000$ ms (C) for the case of the travelling chimera state with domains of subthreshold oscillations ($I_{app}=95$ $\mu$A/cm2, $g_{syn}=5.5$ mS/cm2, $r=0.73$) Regime | ACM | dimension | number
---|---|---|---
|  | of $\mathbf{\Delta t}$ | of clusters
Asynchronous | $R^{2}=0$ | $-$ | $-$
state |  |  |
Global | $R^{2}=1$ | $L=1$ | $L$
synchronization |  |  |
Cluster | $R^{2}=1$ | $1<L\ll N$ | $L$
synchronization |  |  |
Travelling wave | $R^{2}=1$ | $L=N$ | $-$
Chimera state | $0<R^{2}<1$ | - | $L_{lsg}$
Table 2: Table for parameter values for identification of the different
dynamical regimes

###  II.2 Frequency distributions for travelling chimeras and calculation of
their speed

One of the visual criteria to identify a chimera state is a frequency
distribution diagram. We define the average spiking frequency (average
interspike interval) of the $i$-th neuron as

$f[V(i,t)]=\frac{1}{T_{total}}\left|\left\\{t:V(i,t)=V_{thr}\text{ }\&\text{
}\dot{V}(i,t)>0\right\\}\right|,$ (7)

where $V(i,t)$ is the membrane potential of the $i$-th neuron, $\dot{V}(i,t)$
is its time derivative, $V_{thr}$ is the value for Poincare section of the
neuronal voltage, and $T_{total}$ is the duration of the time series.

For the travelling chimera state that moves across the ring, it is impossible
to plot such diagram directly. Let us define the chimera speed $v_{chi}$ as
the number of network elements, the chimera moves per unit of time. Thus,
using the travelling wave coordinates ($V(i,t)\rightarrow V(i+v_{chi}t,t)$) we
can “stop” the chimera (like in Fig. 6B,C). In this case, all neurons can be
split into the coherent and incoherent domains. In other words, there are some
neurons that fire periodically with the same frequencies and the others that
show irregular activity (Fig. 6C). Actually, if the number of synchronous
neurons is maximized than we can distinctly split the coherent and incoherent
chimera domains (Fig. 6E) and find properly the value of chimera speed
$v_{chi}$. If we also have some clusters of subthreshold oscillations (belong
to a chimera like in Fig. 5) it would be better to take them into account and
also “align” them horizontally in the raster plot. According to this idea, let
us consider heuristic functional $H[v]$:

$H[v]=\frac{N_{coh}[v]}{||N_{coh}[v]||}+\frac{N_{inh}[v]}{||N_{inh}[v]||},$
(8)

which is the weighted sum of the heuristics $N_{coh}[v]$ and $N_{inh}[v]$
that, for the fixed value $v$, show the maximum numbers of neurons oscillating
with the same spiking frequency and of neurons demonstrating subthreshold
oscillations, respectively. Now, one can calculate the chimera speed by
solving the optimization problem. In the one-dimensional case it has a form:

$v_{chi}=\operatorname*{arg\,max}_{v}Hf[V\left((i+tv)\text{ mod }N,t\right)].$
(9)

Note that the method works successfully for $N_{coh}[v]$ only, but taking into
account inhibited neurons (if they exist) makes it more precise. The example
of the typical form of $H[v]$ is shown in Fig. 6A. It has a global maximum
that can be reliably found by the automated algorithms.

Figure 6: Heuristic functional $H[v]$ for the case of travelling chimera state
(A). Initial (B) and traveling chimera coordinates (C) rasterplots and
corresponding frequency diagrams (D) and (E). Parameters: $I_{app}=95$
$\mu$A/cm2, $g_{syn}=4.5$ mS/cm2, $r=0.82$

Alternatively, frequency variance can be used as an heuristic functional:

$H[v]=D[v]=\sum_{j=1}^{N}\left(f_{j}[v]-\langle
f_{i}[v]\rangle_{i}\right)^{2}$ (10)

where $f_{i}[v]=f[V\left((i+tv)\text{ mod }N,t\right)]$ and $\langle
f_{i}[v]\rangle_{i}$ is the spatially average. The only drawback of this
method is that in some cases a several maxima exist. Therefore, extra effort
is needed to select the correct one.

##  III Conclusion

In this paper, we introduced a robust and universal approach for studying
synchronization processes in neuronal networks based on the ACM-parameter. In
comparison with other methods, it has the following advantages: absence of
internal method parameters, simplicity of calculations and the ability to
automatically distinguish between stationary and travelling chimera states,
global and cluster synchronization, and travelling waves (Table 2). In
addition, our approach can find the number of synchronous clusters and chimera
speed if it travels.

Of course, there are situations where previously proposed methods are
appropriate. For example, if the elements of a network have well defined
phases, the Kuramoto order parameter [2] is highly recommended and may also
apply to neural networks under specific circumstances (see, for example, [20,
36, 37]). We also point out that although our approach was created for spiking
neuronal networks, it also gives excellent results for multiple network
classes, and, in particular, for the networks of phase oscillators
(simulations not shown). Therefore, we would suggest that our methodology is
more universal and is able to differentiate the cluster and chimera states, as
well as determine the number and size of synchronous clusters, and can be very
useful in networks where it is necessary to separate these states in the
parameter space like [38].

We should also like to comment that the strength of incoherence [31] together
with the discontinuity measure, while being able to differentiate chimera and
multichimera states, cannot do so for coherent states consisting of multiple
synchronous clusters and chimeras. In this case, to avoid confusion, the
authors suggested studying the distribution of membrane potentials differences
in the network to eliminate the removable discontinuity. Most importantly,
this method is exquisitely sensitive to the choice of method parameters (bin
number and the coherence threshold). For example, for travelling waves the
strength of incoherence may be equal to either 1 or 0 (with a jump between
these values) depending on the coherence threshold. Almost same situation is
observed for the case of travelling chimeras and it cannot distinguish between
cluster synchronization and asynchronous state.

In addition, we suggest the new method of chimera speed calculation that
allows to stop the chimera and properly separate coherent and incoherent
chimera domains. Also, there exists a method [28] based on Fourier spectra
that has been suggested for calculation of the chimera speed for the lattice
limit cycle (LLC) model and shown to be useful for neuronal networks (see, for
example, [29]). Our method has some significant advantages. First, it is
simpler for implementation and calculations can be easily automated. Second,
one does not need to select the “right” maximum in the Fourier spectra, but to
find the global maximum of $H[v]$. Also, our method is less demanding on the
simulation time and, in fact, the neuronal time series should contain only one
non-coherent cluster of a chimera. Increasing of a simulation time only leads
to a slight growth in accuracy.

To sum up, we propose a universal and robust methodology to characterize
multiple complex spatio-temporal patterns in a wide range of coupled networks,
notably those with relaxation dynamics of the constituent elements.

##  IV Acknowledgements

The research leading to these results has received funding from the Basic
Research Program at the National Research University Higher School of
Economics. This research was supported in part through computational resources
of HPC facilities at NRU HSE. BG acknowledge support from the ANR Project
ERMUNDY (Grant No ANR-18-CE37-0014)

## References

  * Parastesh  _et al._ [2020] F. Parastesh, S. Jafari, H. Azarnoush, Z. Shahriari, Z. Wang, S. Boccaletti, and M. Perc, Chimeras, Physics Reports (2020).
  * Abrams and Strogatz [2004] D. M. Abrams and S. H. Strogatz, Chimera states for coupled oscillators, Physical review letters 93, 174102 (2004).
  * Panaggio and Abrams [2015] M. J. Panaggio and D. M. Abrams, Chimera states: coexistence of coherence and incoherence in networks of coupled oscillators, Nonlinearity 28, R67 (2015).
  * Zakharova  _et al._ [2014] A. Zakharova, M. Kapeller, and E. Schöll, Chimera death: Symmetry breaking in dynamical networks, Physical review letters 112, 154101 (2014).
  * Laing [2010] C. R. Laing, Chimeras in networks of planar oscillators, Physical Review E 81, 066221 (2010).
  * Hagerstrom  _et al._ [2012] A. M. Hagerstrom, T. E. Murphy, R. Roy, P. Hövel, I. Omelchenko, and E. Schöll, Experimental observation of chimeras in coupled-map lattices, Nature Physics 8, 658 (2012).
  * Tinsley  _et al._ [2012] M. R. Tinsley, S. Nkomo, and K. Showalter, Chimera and phase-cluster states in populations of coupled chemical oscillators, Nature Physics 8, 662 (2012).
  * Martens  _et al._ [2013] E. A. Martens, S. Thutupalli, A. Fourrière, and O. Hallatschek, Chimera states in mechanical oscillator networks, Proceedings of the National Academy of Sciences 110, 10563 (2013).
  * Viktorov  _et al._ [2014] E. A. Viktorov, T. Habruseva, S. P. Hegarty, G. Huyet, and B. Kelleher, Coherence and incoherence in an optical comb, Physical review letters 112, 224101 (2014).
  * Larger  _et al._ [2015] L. Larger, B. Penkovsky, and Y. Maistrenko, Laser chimeras as a paradigm for multistable patterns in complex systems, Nature communications 6, 1 (2015).
  * Uy  _et al._ [2019] C.-H. Uy, L. Weicker, D. Rontani, and M. Sciamanna, Optical chimera in light polarization, APL Photonics 4, 056104 (2019).
  * Wang and Liu [2020] Z. Wang and Z. Liu, A brief review of chimera state in empirical brain networks, Frontiers in Physiology 11 (2020).
  * Mukhametov  _et al._ [1977] L. Mukhametov, A. Y. Supin, and I. Polyakova, Interhemispheric asymmetry of the electroencephalographic sleep patterns in dolphins., Brain research (1977).
  * Rattenborg  _et al._ [2000] N. C. Rattenborg, C. Amlaner, and S. Lima, Behavioral, neurophysiological and evolutionary perspectives on unihemispheric sleep, Neuroscience & Biobehavioral Reviews 24, 817 (2000).
  * Compte  _et al._ [2000] A. Compte, N. Brunel, P. S. Goldman-Rakic, and X.-J. Wang, Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model, Cerebral cortex 10, 910 (2000).
  * Lainscsek  _et al._ [2019] C. Lainscsek, N. Rungratsameetaweemana, S. S. Cash, and T. J. Sejnowski, Cortical chimera states predict epileptic seizures, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 121106 (2019).
  * Andrzejak  _et al._ [2016] R. G. Andrzejak, C. Rummel, F. Mormann, and K. Schindler, All together now: Analogies between chimera state collapses and epileptic seizures, Scientific reports 6, 1 (2016).
  * Sakaguchi [2006] H. Sakaguchi, Instability of synchronized motion in nonlocally coupled neural oscillators, Physical Review E 73, 031907 (2006).
  * Calim  _et al._ [2018] A. Calim, P. Hövel, M. Ozer, and M. Uzuntarla, Chimera states in networks of type-i morris-lecar neurons, Physical Review E 98, 062217 (2018).
  * Omelchenko  _et al._ [2013] I. Omelchenko, E. Omel’chenko, P. Hövel, and E. Schöll, When nonlocal coupling between oscillators becomes stronger: patched synchrony or multichimera states, Physical review letters 110, 224101 (2013).
  * Shepelev  _et al._ [2017] I. A. Shepelev, T. E. Vadivasova, A. Bukh, G. Strelkova, and V. Anishchenko, New type of chimera structures in a ring of bistable fitzhugh–nagumo oscillators with nonlocal interaction, Physics Letters A 381, 1398 (2017).
  * Wang  _et al._ [2020] Z. Wang, S. Baruni, F. Parastesh, S. Jafari, D. Ghosh, M. Perc, and I. Hussain, Chimeras in an adaptive neuronal network with burst-timing-dependent plasticity, Neurocomputing (2020).
  * Bera  _et al._ [2016a] B. K. Bera, D. Ghosh, and M. Lakshmanan, Chimera states in bursting neurons, Physical Review E 93, 012205 (2016a).
  * Xie  _et al._ [2014] J. Xie, E. Knobloch, and H.-C. Kao, Multicluster and traveling chimera states in nonlocal phase-coupled oscillators, Physical Review E 90, 022919 (2014).
  * Omel’chenko [2019] O. Omel’chenko, Travelling chimera states in systems of phase oscillators with asymmetric nonlocal coupling, Nonlinearity 33, 611 (2019).
  * Dudkowski  _et al._ [2019] D. Dudkowski, K. Czołczyński, and T. Kapitaniak, Traveling chimera states for coupled pendula, Nonlinear Dynamics 95, 1859 (2019).
  * Parastesh  _et al._ [2019] F. Parastesh, S. Jafari, and H. Azarnoush, Traveling patterns in a network of memristor-based oscillators with extreme multistability, The European Physical Journal Special Topics 228, 2123 (2019).
  * Hizanidis  _et al._ [2015] J. Hizanidis, E. Panagakou, I. Omelchenko, E. Schöll, P. Hövel, and A. Provata, Chimera states in population dynamics: networks with fragmented and hierarchical connectivities, Physical Review E 92, 012915 (2015).
  * Bera  _et al._ [2016b] B. K. Bera, D. Ghosh, and T. Banerjee, Imperfect traveling chimera states induced by local synaptic gradient coupling, Physical Review E 94, 012215 (2016b).
  * Mishra  _et al._ [2017] A. Mishra, S. Saha, D. Ghosh, G. V. Osipov, and S. K. Dana, Traveling chimera pattern in a neuronal network under local gap junctional and nonlocal chemical synaptic interactions, Opera Medica et Physiologica (2017).
  * Gopal  _et al._ [2014] R. Gopal, V. Chandrasekar, A. Venkatesan, and M. Lakshmanan, Observation and characterization of chimera states in coupled dynamical systems with nonlocal coupling, Physical Review E 89, 052914 (2014).
  * Golomb  _et al._ [2001] D. Golomb, D. Hansel, and G. Mato, Mechanisms of synchrony of neural activity in large networks, in  _Handbook of biological physics_ , Vol. 4 (Elsevier, 2001) pp. 887–968.
  * Morris and Lecar [1981] C. Morris and H. Lecar, Voltage oscillations in the barnacle giant muscle fiber, Biophysical journal 35, 193 (1981).
  * Destexhe  _et al._ [1998] A. Destexhe, Z. F. Mainen, and T. J. Sejnowski, Kinetic models of synaptic transmission, Methods in neuronal modeling 2, 1 (1998).
  * Rosenblum  _et al._ [1997] M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, From phase to lag synchronization in coupled chaotic oscillators, Physical Review Letters 78, 4193 (1997).
  * Omelchenko  _et al._ [2015] I. Omelchenko, A. Provata, J. Hizanidis, E. Schöll, and P. Hövel, Robustness of chimera states for coupled fitzhugh-nagumo oscillators, Physical Review E 91, 022917 (2015).
  * Chouzouris  _et al._ [2018] T. Chouzouris, I. Omelchenko, A. Zakharova, J. Hlinka, P. Jiruska, and E. Schöll, Chimera states in brain networks: Empirical neural vs. modular fractal connectivity, Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 045112 (2018).
  * Kasatkin  _et al._ [2017] D. Kasatkin, S. Yanchuk, E. Schöll, and V. Nekorkin, Self-organized emergence of multilayer structure and chimera states in dynamical networks with adaptive couplings, Physical Review E 96, 062211 (2017).