††thanks: Corresponding author. E-mail<EMAIL_ADDRESS>
# Designing two-dimensional limit-cycle oscillators with prescribed
trajectories and phase-response characteristics

Norihisa Namura  Department of Systems and Control Engineering, Tokyo
Institute of Technology, Tokyo 152-8552, Japan  Tsubasa Ishii  Department of
Systems and Control Engineering, Tokyo Institute of Technology, Tokyo
152-8552, Japan  Hiroya Nakao  Department of Systems and Control Engineering,
Tokyo Institute of Technology, Tokyo 152-8552, Japan

###### Abstract

We propose a method for designing two-dimensional limit-cycle oscillators with
prescribed periodic trajectories and phase response properties based on the
phase reduction theory, which gives a concise description of weakly-perturbed
limit-cycle oscillators and is widely used in the analysis of synchronization
dynamics. We develop an algorithm for designing the vector field with a stable
limit cycle, which possesses a given shape and also a given phase sensitivity
function. The vector field of the limit-cycle oscillator is approximated by
polynomials whose coefficients are estimated by convex optimization. Linear
stability of the limit cycle is ensured by introducing an upper bound to the
Floquet exponent. The validity of the proposed method is verified numerically
by designing several types of two-dimensional existing and artificial
oscillators. As applications, we first design a limit-cycle oscillator with an
artificial star-shaped periodic trajectory and demonstrate global entrainment.
We then design a limit-cycle oscillator with an artificial high-harmonic phase
sensitivity function and demonstrate multistable entrainment caused by a high-
frequency periodic input.

##  I INTRODUCTION

Synchronization of rhythmic systems have found wide applications in various
fields of engineering in recent years. Some examples include human-robot
interactions mortl2014rhythm ; Jouaiti2018hebbian , frequency tuning or
stabilization in electrical oscillators millimeter1 ; millimeter2 ; lockrange1
, power networks Dorfler2012synchronization , suppression of pulsus alternans
in the heart wilson2017spatiotemporal ; Monga2019optimal , and adjustment of
circadian rhythms in shift works Stone2019application . Various methods for
controlling synchronization dynamics of rhythmic systems have also been
developed Bai2019asymptotic ; Oviedo2018synchronization .

Rhythmic systems are commonly modeled as limit-cycle oscillators, namely,
nonlinear dynamical systems with stable limit-cycle trajectories Strogatz .
Some examples include brain waves stankovski2015coupling ;
stankovski2017neural , animal gaits gait1 ; gait2 ; kobayashi2016 ;
funato2016evaluation , heartbeats and breathing heart , and passive walking
hobbelen2007limit ; garcia1998simplest . Synchronization of limit-cycle
oscillators typically occur when they are periodically perturbed or mutually
coupled. For example, when a periodic input is given to a limit-cycle
oscillator, entrainment or phase locking Kuramoto1984 ;
pikovsky2001synchronization , where the oscillator synchronizes with the
periodic input, can be observed. When two or more oscillators are coupled
together, mutual synchronization can be observed, in which their rhythms are
aligned with each other. In the case of living organisms, the circadian
rhythms, internal body clocks that exist in many organisms from bacteria to
mammals, are entrained to the environmental light cycles. Fireflies flicker
synchronously in unison due to the mutual interaction, namely, by sensing each
other’s luminescence.

The phase reduction theory Kuramoto1984 ; Hoppensteadt1997 ; Winfree2001 ;
Ermentrout2010 ; Nakao2016 ; monga2019phase ; Kuramoto2019 ; Ermentrout2019
is useful for analyzing synchronization dynamics of limit-cycle oscillators
subjected to weak inputs. It represents the multidimensional state of a limit-
cycle oscillator by using a single phase value defined along its limit cycle
and approximately describes the oscillator dynamics by a one-dimensional phase
equation. The phase equation, determined by the natural frequency and phase-
response characteristics of the oscillator, has been extensively used for
analyzing various types of synchronization dynamics of limit-cycle oscillators
Kuramoto1984 ; Winfree2001 ; Nakao2016 . In particular, various methods for
optimizing or designing the periodic waveforms of the external inputs to
achieve desired synchronization dynamics have been proposed on the basis of
the phase reduction theory Harada2010optimal ; Zlotnik2013optimal ;
Tanaka2014optimal ; Tanaka2015optimal ; Pikovsky2015maximizing ;
Zlotnik2016phase ; Qiao2017entrainment ; Kato2021optimization ; Takata2021fast
.

Studies on the design of oscillators with desirable periodic trajectories have
also been conducted righetti2009adaptive ; isjpeert2013dynamical ;
ajalooeian2013general ; Khoramshahi2017adaptive ; hongu2022nonlinear . For
instance, rhythmic dynamics have been designed by constructing central pattern
generators in the field of robotics Ijspeert2001connectionist ;
Righetti2006programmable ; ajalooeian2012design ; Jouaiti2018hebbian ;
pasandi2022integrated . Oscillatory trajectories have also been designed by
using neural networks Ruiz1998existence ; Townley2000existence ;
Zegars2003trajectory ; kuroe2005method ; Jouffroy2008design . Several studies
have proposed methods to design dynamical systems with stable limit cycles of
given shapes by constructing the vector fields Okada2002Polynomial ;
okada2003hierachical ; pasandi2019data ; pasandi2020programmable . In
applications, it is necessary to control the synchronization dynamics of the
oscillators by appropriately choosing their phase-response characteristics and
the applied external inputs or mutual interactions. It is also important to
ensure the stability of the periodic trajectory. However, the phase-response
characteristics of the oscillators, which are essentially important in
determining their synchronization properties, have not been considered in most
of the conventional methods.

In this study, we propose a method for designing two-dimensional limit-cycle
oscillators with given phase-response characteristics in addition to given
stable periodic trajectories from the viewpoint of synchronization control. To
ensure the stability of the limit cycle, we use the Floquet theory
hartman1964ordinary ; guckenheimer1983nonlinear  and introduce an inequality
constraint on the Floquet exponent characterizing the linear stability of the
limit-cycle trajectory. We approximate the vector field of the oscillator by
polynomials and optimize their coefficients to satisfy the given conditions.
Our optimization problem is convex and can be numerically solved without
difficulty to yield the target vector field. As applications, we design an
oscillator with an artificial star-shaped periodic trajectory and an
oscillator with an artificial phase sensitivity function that leads to
multistable entrainment with multiple stable phase differences.

This paper is organized as follows. We first explain the phase reduction
theory and synchronization in Sec. II, then describe the method for designing
the vector field of two-dimensional oscillators in Sec. III. In Sec. IV, we
verify the proposed method by numerical simulations for several types of two-
dimensional oscillators. We conclude this study in Sec. V.

##  II PHASE REDUCTION AND SYNCHRONIZATION

In this section, we briefly explain the phase reduction theory and
synchronization Kuramoto1984 ; Hoppensteadt1997 ; Winfree2001 ; Ermentrout2010
; Nakao2016 ; monga2019phase ; Kuramoto2019 ; Ermentrout2019 .

###  II.1 Phase reduction theory

####  II.1.1 Asymptotic phase

We consider a limit-cycle oscillator described by

$\displaystyle\frac{d}{dt}\bm{X}(t)=\bm{F}(\bm{X}(t)),$ (1)

where $\bm{X}(t)\in\mathbb{R}^{N}$ is the system state at time $t$. We assume
that the system has an exponentially stable limit-cycle trajectory
$\bm{X}_{0}(t)$ with a natural period $T$ and frequency $\omega=2\pi/T$, which
is a $T$-periodic function of $t$ satisfying $\bm{X}_{0}(t+T)=\bm{X}_{0}(t)$.

First, we assign a phase $\theta\in[0,2\pi)$ for each state on the limit
cycle, where $0$ and $2\pi$ are considered identical, by choosing a state
$\bm{X}_{0}(0)$ on the limit cycle as the phase origin, $\theta=0$, and
defining the phase of the state $\bm{X}_{0}(t)$ at $t>0$ as $\theta=\omega
t\mod{2\pi}$. We denote the state on the limit cycle with phase $\theta$ by
$\bm{\mathcal{X}}_{0}(\theta)=\bm{X}_{0}(\theta/\omega)$.

Next, we extend the definition of the phase to the basin of the limit cycle.
The phase of a system state $\bm{X}_{A}$ in the basin is defined as $\theta$
if
$\lim_{t\to\infty}\left\|\Psi_{t,t_{0}}\bm{X}_{A}-\Psi_{t,t_{0}}\bm{\mathcal{X}}_{0}(\theta)\right\|=0$
holds, where $\Psi_{t,t_{0}}$ is the flow of Eq. (1) and $\|\cdot\|$ is the
$L_{2}$ norm. That is, if the state $\Psi_{t,t_{0}}\bm{X}_{A}$ started from
$\bm{X}_{A}$ at $t_{0}$ converges to the state
$\Psi_{t,t_{0}}\bm{\mathcal{X}}_{0}(\theta)$ started from
$\bm{\mathcal{X}}_{0}(\theta)$ at $t_{0}$ as $t\to\infty$, we consider that
the state $\bm{X}_{A}$ has the same phase $\theta$ as
$\bm{\mathcal{X}}_{0}(\theta)$. This defines a phase function $\Theta$ for all
states $\bm{X}$ in the basin, which assigns a phase value
$\theta(t)=\Theta(\bm{X}(t))$ to the state $\bm{X}(t)$ and satisfies

$\displaystyle\frac{d}{dt}\theta(t)=\frac{d}{dt}\Theta(\bm{X}(t))=\omega.$ (2)

The phase defined in this way is called the asymptotic phase, and the level
sets of the phase function are called “isochrons” Kuramoto1984 ;
Hoppensteadt1997 ; Winfree2001 ; Ermentrout2010 ; Nakao2016 .

####  II.1.2 Phase sensitivity function

The phase sensitivity function (PSF, also known as the infinitesimal phase
resetting curve, iPRC) is a key quantity in analyzing synchronization of
limit-cycle oscillators by weak inputs. It is defined by the gradient of the
phase function at $\bm{\mathcal{X}}_{0}(\theta)$ on the limit cycle as

$\displaystyle\bm{Z}(\theta)=\nabla\Theta(\bm{X})|_{\bm{X}=\bm{\mathcal{X}}_{0}(\theta)}.$
(3)

The PSF describes linear phase-response characteristics of the oscillator
caused by weak external inputs Kuramoto1984 ; Winfree2001 ; Brown2004 ;
Ermentrout2010 .

The PSF cannot be obtained analytically in general, so it should be calculated
numerically or measured experimentally. If the mathematical model of the
oscillator is known, the PSF $\bm{Z}(\theta)$ can be obtained by numerically
calculating the $2\pi$-periodic solutions to the following adjoint equation
Brown2004 ; Ermentrout2010 ; Kuramoto2019 :

$\displaystyle\omega\frac{d}{d\theta}\bm{Z}(\theta)$
$\displaystyle=-\tilde{\bm{J}}^{\top}(\theta)\bm{Z}(\theta),$ (4)

where $\tilde{\bm{J}}^{\top}(\theta)$ denotes the transposed Jacobian matrix
of the vector field $\bm{F}$ at $\bm{X}=\bm{\mathcal{X}}_{0}(\theta)$. The PSF
$\bm{Z}(\theta)$ should satisfy
$\bm{Z}(\theta)\cdot\bm{F}\left(\bm{\mathcal{X}}_{0}(\theta)\right)=\omega$ as
a normalization condition.

In this study, we consider the PSF as a function with a time argument,

$\displaystyle\bm{Z}(\theta)=\bm{Z}(\omega t)\coloneqq\tilde{\bm{Z}}(t).$ (5)

For this $\tilde{\bm{Z}}(t)$, the adjoint equation (6) and normalization
condition (7) are given by

$\displaystyle\frac{d}{dt}\tilde{\bm{Z}}(t)=-\bm{J}^{\top}(t)\tilde{\bm{Z}}(t),$
(6)
$\displaystyle\tilde{\bm{Z}}(t)\cdot\bm{F}\left(\bm{\mathcal{X}}_{0}(\theta(t))\right)=\omega,$
(7)

where $\bm{J}^{\top}(t)=\tilde{\bm{J}}^{\top}(\theta(t))$ is the transposed
Jacobian matrix of $\bm{F}$ evaluated at $\bm{\mathcal{X}}_{0}(\theta(t))$ on
the limit cycle.

####  II.1.3 Phase reduction

Weakly-perturbed limit-cycle oscillators can generally be reduced to a one-
dimensional phase equation by using the PSF $\bm{Z}(\theta)$. Let us consider
the dynamics of a weakly-perturbed limit-cycle oscillator described by

$\displaystyle\frac{d}{dt}\bm{X}(t)=\bm{F}(\bm{X}(t))+\varepsilon\tilde{\bm{\eta}}(\bm{X}(t),t),$
(8)

where $\varepsilon\tilde{\bm{\eta}}(\bm{X},t)$ denotes the weak input and
$0<\varepsilon\ll 1$ is a small parameter representing the intensity of the
input. Since $\varepsilon$ is small, the deviation of the oscillator state
$\bm{X}(t)$ from the state $\bm{\mathcal{X}}_{0}(\theta(t))=\bm{X}_{0}(t)$ on
the limit cycle with the same phase $\theta(t)$ is $O(\varepsilon)$, i.e.,
$\bm{X}(t)=\bm{\mathcal{X}}_{0}(\theta(t))+O(\varepsilon)$. The time evolution
of the phase $\theta(t)=\Theta(\bm{X}(t))$ of the limit-cycle oscillator is
then described by

$\displaystyle\begin{aligned}
\frac{d}{dt}\theta(t)&=\frac{d}{dt}\Theta(\bm{X}(t))\\\
&=\omega+\varepsilon\nabla\Theta(\bm{X}(t))\cdot\tilde{\bm{\eta}}(\bm{X}(t),t)\\\
&=\omega+\varepsilon\nabla\Theta(\bm{\mathcal{X}}_{0}(\theta(t)))\cdot\tilde{\bm{\eta}}(\bm{\mathcal{X}}_{0}(\theta(t)),t)+O(\varepsilon^{2})\\\
&=\omega+\varepsilon\bm{Z}(\theta(t))\cdot\bm{\eta}(\theta(t),t)+O(\varepsilon^{2}).\end{aligned}$
(9)

Therefore, the phase $\theta(t)$ obeys the following approximate phase
equation up to $O(\varepsilon)$:

$\displaystyle\frac{d}{dt}\theta(t)=\omega+\varepsilon\bm{Z}(\theta(t))\cdot\bm{\eta}(\theta(t),t).$
(10)

The simplicity of this phase equation has facilitated detailed and extensive
analysis of synchronization caused by weak inputs Kuramoto1984 ; Winfree2001 ;
Brown2004 ; Ermentrout2010 .

###  II.2 Synchronization

By using phase reduction, we can analyze the synchronization of limit-cycle
oscillators with a weak periodic external input. We consider a limit-cycle
oscillator subjected to weak periodic external input described by

$\displaystyle\frac{d}{dt}\bm{X}(t)=\bm{F}(\bm{X}(t))+\varepsilon\tilde{\bm{q}}(\bm{X},t),$
(11)

where $\tilde{\bm{q}}(\bm{X},t)$ denotes a weak periodic external input with
the period $\tau$ and the frequency $\Omega=2\pi/\tau$ satisfying
$\tilde{\bm{q}}(\bm{X},t+\tau)=\tilde{\bm{q}}(\bm{X},t)$. When the intensity
$\varepsilon$ is sufficiently small, the time evolution of the phase
$\theta(t)$ is given by

$\displaystyle\frac{d}{dt}\theta(t)=\omega+\varepsilon\bm{Z}(\theta(t))\cdot\bm{q}(\theta(t),t)$
(12)

by phase reduction, where the external input is approximately evaluated on the
limit cycle, namely,
$\bm{q}(\theta(t),t)=\bm{q}(\bm{\mathcal{X}}(\theta(t)),t)$.

Assuming that the natural frequency $\omega$ of the oscillator and the
frequency $\Omega$ of the periodic external input are sufficiently close, we
denote the frequency difference by $\varepsilon\Delta=\omega-\Omega$, where
“$\varepsilon\Delta$” indicates that the frequency difference is
$O(\varepsilon)$. The oscillator phase relative to the external input,
$\phi(t)=\theta(t)-\Omega t$, obeys

$\displaystyle\frac{d}{dt}{\phi(t)}=\varepsilon\left[\Delta+\bm{Z}(\phi{(t)}+\Omega
t)\cdot\bm{q}(\phi{(t)}+\Omega t,t)\right],$ (13)

where the range of $\phi$ is extended outside $[0,2\pi)$ and $\bm{Z}$ and
$\bm{q}$ are regarded as periodic functions of period $2\pi$. Since the right-
hand side is $O(\varepsilon)$ and $\phi$ is a slow variable, this equation can
further be simplified by the averaging approximation Kuramoto1984 ;
Hoppensteadt1997 ; Nakao2016 , namely, by averaging the right-hand side over
one period of the external input while keeping $\phi(t)$ fixed. This yields a
single-variable autonomous system,

$\displaystyle\frac{d}{dt}{\phi(t)}=\varepsilon\left[\Delta+\Gamma(\phi{(t)})\right],$
(14)

where we defined a $2\pi$-periodic phase coupling function

$\displaystyle\begin{aligned}
\Gamma(\phi)&=\frac{1}{\tau}\int_{0}^{\tau}\bm{Z}(\phi+\Omega
t)\cdot\bm{q}(\phi+\Omega t,t)dt\\\
&=\frac{1}{2\pi}\int_{0}^{2\pi}\bm{Z}(\phi+\psi)\cdot\bm{q}\left(\phi+\psi,\psi/\Omega\right)d\psi,\end{aligned}$
(15)

where $\psi=\Omega t$ is the phase of the external input.

We can analyze the entrainment of the oscillator to the weak periodic input by
using this simplified equation. Specifically, the stable fixed point of Eq.
(14) corresponds to a phase-locking point, where the oscillator phase is
entrained to the periodic input.

##  III DESIGN OF OSCILLATOR DYNAMICS

In this section, we propose a method for designing a stable two-dimensional
limit-cycle oscillator with a prescribed periodic trajectory and PSF.

###  III.1 Conditions for the periodic trajectory and PSF

Our aim is to design a vector field with a prescribed stable limit cycle
$\bm{p}(t)$ and PSF $\tilde{\bm{Z}}(t)$ of period $T$ and frequency
$\omega=2\pi/T$. It is noted that the periodic trajectory and the PSFs are not
completely independent. Specifically, the PSFs should satisfy the
normalization condition $\dot{\bm{p}}(t)\cdot\tilde{\bm{Z}}(t)=\omega$.

We approximate the vector field of the oscillator by using polynomials as
namura2022estimating

$\displaystyle\bm{F}(\bm{X})=\begin{bmatrix}F_{1}(\bm{X})\\\
F_{2}(\bm{X})\end{bmatrix}\simeq\begin{bmatrix}\bm{U}^{\top}(\bm{X})\bm{\zeta}_{1}\\\
\bm{U}^{\top}(\bm{X})\bm{\zeta}_{2}\end{bmatrix},$ (16)

where $\bm{X}=[x_{1}\ x_{2}]^{\top}$,

$\displaystyle\bm{U}(\bm{X})=\begin{bmatrix}1&\overline{x_{1}}&\overline{x_{2}}&\overline{x_{1}^{2}}&\overline{x_{1}x_{2}}&\overline{x_{2}^{2}}&\cdots&\overline{x_{2}^{n}}\end{bmatrix}^{\top}\in\mathbb{R}^{P}$
(17)

is a vector of monomials, $n$ is the maximum degree of the polynomial, and
$\bm{\zeta}_{1},\bm{\zeta}_{2}\in\mathbb{R}^{P}$ are the coefficient vectors.
Here, the overline denotes standardization
($\overline{z}\coloneqq(z-\mu)/\sigma$), where $\mu$ and $\sigma$ are the mean
and the standard deviation of $z$, respectively.

We assume that the system has a non-intersecting differentiable periodic
trajectory $\bm{p}(t)=[p_{1}(t)\ p_{2}(t)]^{\top}$ and impose the following
condition on the vector field:

$\displaystyle\bm{F}(\bm{p}(t))=\frac{d}{dt}\bm{p}(t).$ (18)

This equation can be expressed in the polynomial approximation as follows:

$\displaystyle\begin{aligned}
\bm{U}^{\top}(\bm{p}(t))\bm{\zeta}_{1}\simeq\frac{d}{dt}p_{1}(t),\\\
\bm{U}^{\top}(\bm{p}(t))\bm{\zeta}_{2}\simeq\frac{d}{dt}p_{2}(t).\end{aligned}$
(19)

Introducing the coefficient $\bm{\xi}=\left[\bm{\zeta}_{1}^{\top}\
\bm{\zeta}_{2}^{\top}\right]^{\top}\in\mathbb{R}^{2P}$ and discretizing the
time $t$ into $L$ points, $t_{k}=(k-1)\Delta t,\;(k=1,\dots,L)$ with $\Delta
t=T/L$, we express the left-hand side as

$\displaystyle\begin{bmatrix}\bm{U}^{\top}(\bm{p}(t_{k}))&\bm{0}\\\
\bm{0}&\bm{U}^{\top}(\bm{p}(t_{k}))\end{bmatrix}\bm{\xi}\coloneqq\bm{A}_{\bm{p},k}\bm{\xi}$
(20)

for each $t_{k}$, where $\bm{A}_{\bm{p},k}\in\mathbb{R}^{2\times 2P}$. We also
denote the right-hand side by

$\displaystyle\bm{b}_{\bm{p},k}\coloneqq\begin{bmatrix}\dot{p}_{1}(t_{k})\\\
\dot{p}_{2}(t_{k})\end{bmatrix},$ (21)

where $\bm{b}_{\bm{p},k}\in\mathbb{R}^{2}$. The difference between both sides
of Eq. (18) under polynomial approximation is then expressed as

$\displaystyle\bm{A}_{\bm{p},k}\bm{\xi}-\bm{b}_{\bm{p},k}.$ (22)

Next, we introduce the conditions for the PSF into the vector field. Assuming
that the oscillator has a PSF $\tilde{\bm{Z}}(t)=[\tilde{Z}_{1}(t)\
\tilde{Z}_{2}(t)]^{\top}$, we require the adjoint equation (6), represented as
a function of $t$, is satisfied. This equation is expressed in the polynomial
approximation as

$\displaystyle\begin{aligned}
\tilde{Z}_{1}(t)\bm{U}_{1}^{\top}(\bm{p}(t))\bm{\zeta}_{1}+\tilde{Z}_{2}(t)\bm{U}_{1}^{\top}(\bm{p}(t))\bm{\zeta}_{2}\simeq-\frac{d}{dt}\tilde{Z}_{1}(t),\\\
\tilde{Z}_{1}(t)\bm{U}_{2}^{\top}(\bm{p}(t))\bm{\zeta}_{1}+\tilde{Z}_{2}(t)\bm{U}_{2}^{\top}(\bm{p}(t))\bm{\zeta}_{2}\simeq-\frac{d}{dt}\tilde{Z}_{2}(t),\end{aligned}$
(23)

where $\bm{U}_{i}(\bm{X})=\left(\nabla\bm{U}(\bm{X})\right)_{i}$.

Introducing the coefficient $\bm{\xi}$ and $\\{t_{k}\\}_{k=1}^{L}$, we express
the left-hand side as

$\displaystyle\begin{bmatrix}\tilde{Z}_{1}(t_{k})\bm{U}_{1}^{\top}(\bm{p}(t_{k}))&\tilde{Z}_{2}(t_{k})\bm{U}_{1}^{\top}(\bm{p}(t_{k}))\\\
\tilde{Z}_{1}(t_{k})\bm{U}_{2}^{\top}(\bm{p}(t_{k}))&\tilde{Z}_{2}(t_{k})\bm{U}_{2}^{\top}(\bm{p}(t_{k}))\end{bmatrix}\bm{\xi}\coloneqq\bm{A}_{\bm{Z},k}\bm{\xi}$
(24)

for each $t_{k}$, where $\bm{A}_{\bm{Z},k}\in\mathbb{R}^{2\times 2P}$, and
also denote the right-hand side as

$\displaystyle\bm{b}_{\bm{Z},k}\coloneqq\begin{bmatrix}-\dot{\tilde{Z}}_{1}(t_{k})\\\
-\dot{\tilde{Z}}_{2}(t_{k})\end{bmatrix},$ (25)

where $\bm{b}_{\bm{Z},k}\in\mathbb{R}^{2}$. The difference between both sides
of Eq. (6) under polynomial approximation is then expressed as

$\displaystyle\bm{A}_{\bm{Z},k}\bm{\xi}-\bm{b}_{\bm{Z},k}.$ (26)

We seek the polynomial vector field that satisfies Eq. (18) and (6) as much as
possible by minimizing the overall approximation errors for the periodic
trajectory and PSF over one period $T$. By introducing matrices

$\displaystyle\bm{A}$
$\displaystyle=\begin{bmatrix}\bm{A}_{\bm{p},1}^{\top}&\cdots&\bm{A}_{\bm{p},L}^{\top}&\bm{A}_{\bm{Z},1}^{\top}&\cdots&\bm{A}_{\bm{Z},L}^{\top}\end{bmatrix}^{\top}\in\mathbb{R}^{4L\times
2P}$ (27)

and

$\displaystyle\bm{b}$
$\displaystyle=\begin{bmatrix}\bm{b}_{\bm{p},1}^{\top}&\cdots&\bm{b}_{\bm{p},L}^{\top}&\bm{b}_{\bm{Z},1}^{\top}&\cdots&\bm{b}_{\bm{Z},L}^{\top}\end{bmatrix}^{\top}\in\mathbb{R}^{4L},$
(28)

the sum of the squared errors in Eqs. (22) and (26) can be expressed as

$\displaystyle E=\frac{1}{2}\|\bm{A}\bm{\xi}-\bm{b}\|^{2}.$ (29)

This error is the objective function to be minimized in our method.

In addition to the above conditions, we introduce regularization of the vector
field by adding the squared norm $\|\bm{\xi}\|^{2}$ of the coefficient
$\bm{\xi}$ with an appropriate weight $\gamma$ to the objective function in
order to obtain a better model. This regularization prevents the coefficients
of the vector field from becoming excessively large, which is expected to
reduce the complexity of the vector field and ensures the uniqueness of the
solution of the optimization problem.

###  III.2 Condition for the linear stability

In addition to the conditions for the trajectory and PSF, we also impose the
condition for the stability of the periodic trajectory on the vector field.
From the Floquet theory hartman1964ordinary ; guckenheimer1983nonlinear , the
two Floquet exponents $\lambda_{1}$ and $\lambda_{2}$ characterizing the
linear stability of the periodic trajectory satisfy

$\displaystyle\exp(\lambda_{1}T)\exp(\lambda_{2}T)=\exp\left(\int_{0}^{T}\mathrm{tr}(\bm{J}(t))dt\right),$
(30)

where $\bm{J}(t)$ is the Jacobian matrix at $\bm{X}_{0}(t)$ on the periodic
trajectory and $\mathrm{tr}$ denotes the trace of a matrix. Taking the
logarithm, we obtain

$\displaystyle\lambda_{1}+\lambda_{2}=\frac{1}{T}\int_{0}^{T}\mathrm{tr}(\bm{J}(t))dt$
(31)

since $\lambda_{1}$ and $\lambda_{2}$ are real for two-dimensional limit-cycle
oscillators.

From the condition (18) for the periodic trajectory, the first Floquet
exponent $\lambda_{1}=0$, thus we only need to evaluate $\lambda=\lambda_{2}$,
which should be negative for the stability. Discretizing the time $t$ into
$\\{t_{k}\\}_{k=1}^{L}$, we can approximate the above integral as

$\displaystyle\lambda\simeq\frac{1}{L}\sum_{k=1}^{L}\left(\bm{J}_{11}(t_{k})+\bm{J}_{22}(t_{k})\right),$
(32)

which is expressed in the polynomial approximation as

$\displaystyle\begin{aligned}
\lambda=\lambda_{2}&\simeq\frac{1}{L}\sum_{k=1}^{L}\left(\bm{U}_{1}^{\top}(\bm{p}(t_{k}))\bm{\zeta}_{1}+\bm{U}_{2}^{\top}(\bm{p}(t_{k}))\bm{\zeta}_{2}\right)\\\
&=\frac{1}{L}\sum_{k=1}^{L}\begin{bmatrix}\bm{U}_{1}^{\top}(\bm{p}(t_{k}))&\bm{U}_{2}^{\top}(\bm{p}(t_{k}))\end{bmatrix}\bm{\xi}\coloneqq\bm{C}\bm{\xi}.\end{aligned}$
(33)

We impose the condition

$\displaystyle\lambda=\bm{C}\bm{\xi}\leq\lambda_{\mathrm{tol}}$ (34)

for some maximum tolerance value $\lambda_{\mathrm{tol}}<0$ so that the
periodic trajectory is sufficiently stable.

###  III.3 Optimization problem

Summarizing the conditions for the periodic trajectory $\bm{p}(t)$, PSF
$\tilde{\bm{Z}}(t)$, and maximal tolerance $\lambda_{\mathrm{tol}}$ of the
linear stability explained in the previous subsections, our optimization
problem for the coefficient vector $\bm{\xi}$ is formulated as

$\displaystyle\begin{gathered}\bm{\xi}^{*}=\mathop{\rm
argmin}\limits_{\bm{\xi}}\frac{1}{2}\|\bm{A}\bm{\xi}-\bm{b}\|^{2}+\gamma\|\bm{\xi}\|^{2}\\\
\mathrm{s.t.}\quad\bm{C}\bm{\xi}\leq\lambda_{\mathrm{tol}}.\end{gathered}$
(37)

Since this is a quadratic programming problem for the parameter $\bm{\xi}$
with a linear inequality constraint, which is convex, it can be easily solved
uniquely and globally.

In the numerical implementation, the weight parameter $\gamma$ is chosen so
that the designed oscillator has a stable limit cycle with the prescribed
properties and no unnecessary attractors arises around the limit cycle. Also,
to match the scales of the errors for the periodic trajectory and PSF in Eqs.
(22) and (26), the conditions for the PSF are multiplied by an appropriate
normalizing constant.

##  IV RESULTS

###  IV.1 Reconstruction of existing oscillators

In this subsection, we test the validity of the proposed method by
reconstructing the periodic trajectories and PSFs of two well-known
oscillators.

####  IV.1.1 van der Pol oscillator

As the first example, we design an oscillator with a periodic trajectory and
PSF of the van der Pol (vdP) oscillator van1927vii ; van1927frequency ;
van1926lxxxviii . The vdP oscillator is described by

$\displaystyle\frac{d}{dt}\begin{bmatrix}x_{1}\\\
x_{2}\end{bmatrix}=\begin{bmatrix}x_{2}\\\
\nu\left(1-x_{1}^{2}\right)x_{2}-x_{1}\end{bmatrix},$ (38)

where we assume $\nu=3$. This vdP oscillator has a period $T=8.860$, natural
frequency $\omega=0.7092$, and the second Floquet exponent
$\lambda_{2}=-3.9396$. The PSF can be calculated by numerically solving the
adjoint equation (4).

We designed the vector field using the time interval $\Delta t=0.005$, the
degree of the polynomial as $n=10$, the weight parameter $\gamma=1.0$, and the
maximum tolerance value of the second Floquet exponent
$\lambda_{\mathrm{tol}}=-0.5$. By optimization, we obtained an oscillator with
a period $T=8.860$, natural frequency $\omega=0.7092$, and the second Floquet
exponent $\lambda_{2}=-3.7260$. The period and natural frequency are almost
identical to those of the original vdP oscillator, while the second Floquet
exponent is much smaller than the assumed maximal tolerance value and close to
the value of the original vdP oscillator.

The periodic trajectory on the $(x_{1},x_{2})$ plane and the velocity
$\dot{p}_{1}(t)$ and $\dot{p}_{2}(t)$ of the designed oscillator are compared
with those of the original vdP oscillator in Figs. 1 (a), (b), and (c); both
of them are almost identical to those of the original oscillator. The vector
field of the designed oscillator is compared with that of the original vdP
oscillator in Fig. 2. Both vector fields are similar in the whole region,
though the scales of the vector fields are slightly different. Figures 1 (d)
and (e) show the PSFs of the designed oscillator and the original vdP
oscillator, which are almost identical to each other. Thus, our method can
reproduce the functional forms of the periodic orbit and PSF of the vdP
oscillator successfully.

Figure 1:  Periodic trajectories and PSFs of the designed oscillator and the
original vdP oscillator. (a) Periodic trajectories. (b,c) Velocities on the
periodic trajectories. (b) $x_{1}$ component, (c) $x_{2}$ component. (d,e)
PSFs. (d) $x_{1}$ component, (e) $x_{2}$ component. In each graph, the red
line shows the designed functional form and the black dotted line shows the
original one, respectively.  Figure 2:  Vector fields of the designed
oscillator and the original vdP oscillator. (a) designed, (b) original.

####  IV.1.2 FitzHugh-Nagumo oscillator

As the second example, we design an oscillator with a periodic trajectory and
PSF of the FitzHugh-Nagumo (FHN) oscillator Nagumo1962Active  described by

$\displaystyle\frac{d}{dt}\begin{bmatrix}x_{1}\\\
x_{2}\end{bmatrix}=\begin{bmatrix}x_{1}(x_{1}-a)(1-x_{1})-x_{2}\\\
c(x_{1}-bx_{2})\end{bmatrix},$ (39)

where we consider $a=-0.1$, $b=0.5$, and $c=0.01$. The FHN oscillator is a
fast-slow system with a large timescale difference between the two variables.
This FHN oscillator has a period $T=126.5$, natural frequency $\omega=0.0497$,
and the second Floquet exponent $\lambda_{2}=-0.4586$. The PSF is calculated
by solving the adjoint equation numerically.

We designed the vector field by setting the time interval as $\Delta t=0.05$,
the degree of the polynomial as $n=10$, the weight parameter $\gamma=1.0\times
10^{-3}$, and the maximum tolerance value of the second Floquet exponent as
$\lambda_{\mathrm{tol}}=-0.5$. The designed oscillator had a period $T=126.5$,
natural frequency $\omega=0.0497$, and the second Floquet exponent
$\lambda_{2}=-0.5000$, where the period and the natural frequency are almost
identical to the original FHN oscillator. We note that the maximum tolerance
value of the second Floquet exponent $\lambda_{\mathrm{tol}}=-0.5$ is below
the original value $\lambda_{2}=-0.4586$ of the FHN oscillator. This shows
that the stability of the periodic trajectory can be varied from that of the
original oscillator to a certain extent.

As shown in Figs. 3 (a), (b), and (c), the original and designed oscillators
have almost identical limit cycles. The vector field of the designed
oscillator is compared to that of the original FHN oscillator in Fig. 4. In
the vicinity of the periodic trajectory, the vector fields are almost similar
and the nullclines of the fast component look alike. They are considerably
different in the regions far away from the periodic trajectory, though it does
not affect the resulting limit cycles largely. The PSFs of the original and
designed oscillators are also almost identical as shown in Figs. 3 (d) and
(e).

Figure 3:  Periodic trajectories and PSFs of the designed oscillator and the
original FHN oscillator. (a) Periodic trajectories. (b,c) Velocities on the
periodic trajectories. (b) $x_{1}$ component, (c) $x_{2}$ component. (d,e)
PSFs. (d) $x_{1}$ component, (e) $x_{2}$ component. In each graph, the red
line shows the designed functional form and the black dotted line shows the
original one, respectively.  Figure 4:  Vector fields of the designed
oscillator and the original FHN oscillator. (a) designed, (b) original.

###  IV.2 Design of an oscillator with an artificial periodic trajectory and
PSF

####  IV.2.1 Design of an oscillator with a star shape

Next, we design a limit-cycle oscillator with an artificial periodic
trajectory and PSF. We assume that the oscillator has the following periodic
trajectory (40) and PSF (41):

$\displaystyle\bm{p}(t)$
$\displaystyle=\begin{bmatrix}\sqrt{2}\cos(t)+\frac{1}{4}\sin(4t)\\\
\sqrt{2}\sin(t)+\frac{1}{4}\cos(4t)\end{bmatrix},$ (40)
$\displaystyle\tilde{\bm{Z}}(t)$
$\displaystyle=\begin{bmatrix}-\sqrt{2}\sin(t)-\cos(4t)\\\
\sqrt{2}\cos(t)+\sin(4t)\end{bmatrix},$ (41)

where the given PSF satisfies the normalization condition
$\dot{\bm{p}}(t)\cdot\tilde{\bm{Z}}(t)=\omega=1$. The periodic trajectory has
a star shape with a $5$-fold symmetry, which is not often seen in real-world
oscillators. The oscillation period is $T=2\pi$.

We designed the vector field that possesses the given artificial periodic
trajectory and PSF by using the number of data points $L=1000$, the degree of
the polynomial $n=10$, the weight parameter $\gamma=1.0\times 10^{-3}$, and
the maximum tolerance value of the second Floquet exponent
$\lambda_{\mathrm{tol}}=-1$. We obtained an oscillator with a period
$T=6.2832$, natural frequency $\omega=1.0000$, which are almost identical to
the assumed values, and the second Floquet exponent $\lambda_{2}=-1.0001$,
which is almost equal to $\lambda_{\mathrm{tol}}$.

The limit cycle of the designed oscillator and the given periodic trajectory
are compared in Figs. 5 (a), (b), and (c). Both the periodic trajectory and
the velocities on the periodic trajectory are almost identical to the assumed
ones. The PSF of the designed oscillator is shown in Figs. 5 (d) and (e),
which is also almost identical to the assumed one. The vector field of the
designed oscillator is shown in Fig. 6.

Figure 5:  Periodic trajectory and PSFs of the designed oscillator and the
prescribed functional forms. (a) Periodic trajectories. (b,c) Velocities on
the periodic trajectories. (b) $x_{1}$ component, (c) $x_{2}$ component. (d,e)
PSFs. (d) $x_{1}$ component, (e) $x_{2}$ component. In each graph, the red
line shows the designed functional form and the black dotted line shows the
original one, respectively.  Figure 6: Vector field and limit cycle of the
designed oscillator with an artificial star-shaped periodic trajectory.

####  IV.2.2 Global entrainment

Here, we demonstrate global entrainment of the star-shaped oscillator designed
above. We apply an external periodic input with a frequency $\omega$, which is
equal to the natural frequency. Assuming that the external periodic input has
the form of $\bm{q}(t)=\left[-\sin(\omega t),\cos(\omega t)\right]^{\top}$, we
can write the reduced phase equation of the oscillator as

$\displaystyle\frac{d}{dt}{\theta(t)}=\omega+\varepsilon\bm{Z}(\theta{(t)})\cdot\bm{q}(t).$
(42)

From Eq. (14), the time evolution of the relative phase
$\phi{(t)}=\theta(t)-\omega t$ obeys

$\displaystyle\begin{aligned}
\frac{d}{dt}{\phi(t)}&=\varepsilon\Gamma(\phi{(t)})\\\
&=\sqrt{2}\varepsilon\cos(\phi{(t)}).\end{aligned}$ (43)

The phase coupling function $\Gamma(\phi)$ is shown in Fig. 7. There are two
fixed points (one stable and one unstable) satisfying $\dot{\phi}=0$ within
$\phi\in[0,2\pi)$; the only stable fixed point is at $\phi{{}^{*}}=\pi/2$,
indicating that the relative phase converges to this point irrespective of the
initial value (except for the unstable fixed point $\phi=3\pi/2$), namely, the
oscillator exhibits global $1:1$ entrainment.

We performed numerical simulations to confirm that global entrainment occurs.
We assumed $\varepsilon=0.01$ and evolved $100$ independent oscillators from
$100$ random initial points on the star-shape periodic trajectory. The
oscillator states at time $t=0,\;20\pi,\;40\pi$, and $100\pi$ are shown in
Fig. 8. As time passes, the phases of the oscillators gradually form a
cluster, and at $t=100\pi$, the oscillators almost converge to the unique
fixed point.

Figure 7:  Phase coupling function $\Gamma(\phi)$ for global entrainment. The
black solid point shows the stable fixed point and the black circle point
shows the unstable fixed point. If $\dot{\phi}=0$ and
$\Gamma^{\prime}(\phi)<0$ at $\phi=\phi^{*}$, where $\Gamma^{\prime}(\phi)$
denotes the derivative of $\Gamma$, $\phi^{*}$ is a stable fixed point.
Figure 8:  Global entrainment of $N=100$ designed artificial oscillators. (a)
$t=0$, (b) $t=20\pi$, (c) $t=40\pi$, and (d) $t=100\pi$. In (d), we can see
the convergence to the unique stable fixed point.

###  IV.3 Design of an oscillator exhibiting multistable entrainment

####  IV.3.1 Design of an oscillator with a high-harmonic PSF

Finally, we design an oscillator with an artificial PSF that leads to
nontrivial multistable entrainment. We assume that the oscillator has the
following periodic trajectory and PSF:

$\displaystyle\bm{p}(t)$ $\displaystyle=\begin{bmatrix}\cos(t)\\\
\sin(t)\end{bmatrix},$ (44) $\displaystyle\tilde{\bm{Z}}(t)$
$\displaystyle=\begin{bmatrix}-\sin(5t)\\\
2\cos(t)-2\cos(3t)+\cos(5t)\end{bmatrix},$ (45)

where the given PSF satisfies the normalization condition
$\dot{\bm{p}}(t)\cdot\tilde{\bm{Z}}(t)=\omega=1$. This periodic trajectory is
simply a unit circle and the system state rotates with a constant frequency
$\omega=1$, that is, the period is $T=2\pi$.

We designed the vector field that possesses the given artificial periodic
trajectory and PSF by setting the number of data points as $L=1000$, the
degree of the polynomial as $n=7$, the weight parameter $\gamma=1.0\times
10^{-2}$, and the maximum tolerance value of the second Floquet exponent as
$\lambda_{\mathrm{tol}}=-1$. We obtained an oscillator with a period
$T=6.2832$, natural frequency $\omega=1.0000$, which are almost identical to
the assumed values, and the second Floquet exponent $\lambda_{2}=-0.9998$,
which is almost equal to the maximum tolerance value $\lambda_{\mathrm{tol}}$.

The limit cycle of the designed oscillator is compared with the given periodic
trajectory in Figs. 9 (a), (b), and (c). Both the periodic trajectory and
velocity are almost identical to the assumed ones. The PSF of the designed
oscillator is shown in Figs. 9 (d) and (e), which is also almost identical to
the prescribed one. The vector field of the designed oscillator is shown in
Fig. 10.

Figure 9:  Periodic trajectory and PSFs of the designed oscillator and the
prescribed functional forms. (a) Periodic trajectories. (b,c) Velocities on
the periodic trajectories. (b) $x_{1}$ component, (c) $x_{2}$ component. (d,e)
PSFs. (d) $x_{1}$ component, (e) $x_{2}$ component. In each graph, the red
line shows the designed functional form and the black dotted line shows the
original one, respectively.  Figure 10: Vector field and limit cycle of the
designed oscillator with an artificial PSF.

####  IV.3.2 Multistable entrainment

Here, we demonstrate multistable entrainment of the oscillator designed above.
We apply an external periodic input with a frequency $\Omega=5\omega$, which
is five times the natural frequency $\omega$. Assuming that the external
periodic input has the form of $\bm{q}(t)=\left[\sin(\Omega
t),0\right]^{\top}$, we can write the reduced phase equation of the oscillator
as

$\displaystyle\frac{d}{dt}{\theta(t)}=\omega+\varepsilon\bm{Z}(\theta{(t)})\cdot\bm{q}(t).$
(46)

Since the ratio of the frequency of the external periodic input to that of the
oscillator is $5$, $5:1$ entrainment is expected to occur. Generalizing the
analysis in Sec. II and defining the phase of the external input to be $\Omega
t/5$, we introduce the phase difference

$\displaystyle\phi{(t)}=\theta{(t)}-\frac{1}{5}\Omega t.$ (47)

The time evolution of $\phi{(t)}$ is then

$\displaystyle\begin{aligned}
\frac{d}{dt}{\phi(t)}&=\frac{d}{dt}{\theta(t)}-\frac{1}{5}\Omega\\\
&=\omega-\frac{1}{5}\Omega+\varepsilon\bm{Z}\left(\phi{(t)}+\frac{1}{5}\Omega
t\right)\cdot\bm{q}(t)\\\ &=\varepsilon\bm{Z}\left(\phi{(t)}+\frac{1}{5}\Omega
t\right)\cdot\bm{q}(t),\end{aligned}$ (48)

and by applying the averaging approximation Kuramoto1984 ; Hoppensteadt1997 ;
Nakao2016 , $\phi{(t)}$ obeys

$\displaystyle\begin{aligned}
\frac{d}{dt}{\phi(t)}&=\varepsilon\Gamma(\phi{(t)})\\\
&=\frac{\varepsilon}{T}\int_{0}^{T}\bm{Z}\left(\phi{(t)}+\frac{1}{5}\Omega
t{{}^{\prime}}\right)\cdot\bm{q}(t{{}^{\prime}})dt{{}^{\prime}}\\\
&=-\frac{\varepsilon}{2}\cos(5\phi{(t)}).\end{aligned}$ (49)

The phase coupling function $\Gamma(\phi)$ is shown in Fig. 11. There are ten
fixed points (five stable and five unstable) satisfying $\dot{\phi}=0$ within
$\phi\in[0,2\pi)$; the five stable fixed points are
$\phi{{}^{*}}=3\pi/10,\;7\pi/10,\;11\pi/10,\;3\pi/2$, and $19\pi/10$.

We performed numerical simulations to confirm that entrainment with multiple
phase differences occurs. We assumed $\varepsilon=0.05$ and evolved $100$
independent oscillators from $100$ random initial points on the unit circle.
The oscillator states at time $t=0,\;2\pi,\;4\pi$, and $16\pi$ are shown in
Fig. 12. As time passes, the phases of the oscillators gradually form
clusters, and at $t=16\pi$, the oscillators almost converge to either of the
five points.

We note that synchronization with such many fixed points is rarely observed in
natural limit-cycle oscillators. Thus, by using the present method, we can
design artificial oscillators with desirable synchronization properties.

Figure 11:  Phase coupling function $\Gamma(\phi)$ for multistable
entrainment. The black solid points show the stable fixed points and the black
circle points show the unstable fixed points.  Figure 12:  Multistable
entrainment of $N=100$ designed artificial oscillators. (a) $t=0$, (b)
$t=2\pi$, (c) $t=4\pi$, and (d) $t=16\pi$. In (d), we can see the convergence
to the five stable fixed points.

##  V CONCLUDING REMARKS

We proposed a method for designing two-dimensional limit-cycle oscillators
that possess stable prescribed periodic trajectories and PSFs. Using this
method, we could design the vector fields that exhibit the periodic
trajectories and PSFs of several types of existing or artificial oscillators.
Furthermore, we were able to design an artificial oscillator with a high-
harmonic PSF that rarely exists in the real world, and were able to
demonstrate multistable entrainment caused by a high-frequency periodic input.
The proposed method is simple and generally applicable to design the vector
fields of various two-dimensional oscillators. Generalizations of the proposed
method to higher-dimensional oscillators will be reported in our future study.

###### Acknowledgements.

H. N. thanks financial support from JSPS KAKENHI (Nos. JP22K11919, JP22H00516,
JPJSBP120202201) and JST CREST (No. JP-MJCR1913).

## References

  * (1) Alexander MÃ¶rtl, Tamara Lorenz, and Sandra Hirche.  Rhythm patterns interaction - synchronization behavior for human-robot joint action.  PLOS ONE, 9(4):1–17, 04 2014.
  * (2) Melanie Jouaiti, Lancelot Caron, and Patrick HÃ©naff.  Hebbian plasticity in CPG controllers facilitates self-synchronization for human-robot handshaking.  Frontiers in Neurorobotics, 12, 2018.
  * (3) A.S. Daryoush.  Optical synchronization of millimeter-wave oscillators for distributed architecture.  IEEE Transactions on Microwave Theory and Techniques, 38(5):467–476, 1990.
  * (4) Kenichi Kawasaki, Yoshiyuki Akiyama, Kenji Komori, Masahiro Uno, Hidenori Takeuchi, Tomoari Itagaki, Yasufumi Hino, Yoshinobu Kawasaki, Katsuhisa Ito, and Ali Hajimiri.  A millimeter-wave intra-connect solution.  IEEE Journal of Solid-State Circuits, 45(12):2655–2666, 2010.
  * (5) Tomoharu Nagashima, Xiuqin Wei, Hisa-Aki Tanaka, and Hiroo Sekiya.  Locking range derivations for injection-locked class-e oscillator applying phase reduction theory.  IEEE Transactions on Circuits and Systems I: Regular Papers, 61(10):2904–2911, 2014.
  * (6) Florian Dörfler and Francesco Bullo.  Synchronization and transient stability in power networks and nonuniform Kuramoto oscillators.  SIAM Journal on Control and Optimization, 50(3):1616–1642, 2012\.
  * (7) Dan Wilson and Jeff Moehlis.  Spatiotemporal control to eliminate cardiac alternans using isostable reduction.  Physica D: Nonlinear Phenomena, 342:32–44, 2017.
  * (8) Bharat Monga and Jeff Moehlis.  Optimal phase control of biological oscillators using augmented phase reduction.  Biological Cybernetics, 113(1):161–178, Apr 2019.
  * (9) Julia E. Stone, Xavier L. Aubert, Henning Maass, Andrew J. K. Phillips, Michelle Magee, Mark E. Howard, Steven W. Lockley, Shantha M. W. Rajaratnam, and Tracey L. Sletten.  Application of a limit-cycle oscillator model for prediction of circadian phase in rotating night shift workers.  Scientific Reports, 9(1):11032, Jul 2019.
  * (10) He Bai and John T. Wen.  Asymptotic synchronization of phase oscillators with a single input.  IEEE Transactions on Automatic Control, 64(4):1611–1618, 2019.
  * (11) Erick Oviedo, Nimrod VÃ¡zquez, and Ricardo Femat.  Synchronization technique of grid-connected power converters based on a limit cycle oscillator.  IEEE Transactions on Industrial Electronics, 65(1):709–717, 2018\.
  * (12) Steven H. Strogatz.  Nonlinear Dynamics and Chaos.  CRC Press, Boca Raton, 2015.
  * (13) T. Stankovski, V. Ticcinelli, P. V. McClintock, and A. Stefanovska.  Coupling functions in networks of oscillators.  New Journal of Physics, 17(3):035002, Mar. 2015.
  * (14) Tomislav Stankovski, Valentina Ticcinelli, Peter VE McClintock, and Aneta Stefanovska.  Neural cross-frequency coupling functions.  Frontiers in Systems Neuroscience, 11:33, 2017.
  * (15) Lotta Borgius, Hiroshi Nishimaru, Vanessa Caldeira, Yuka Kunugise, Peter Löw, Ramon Reig, Shigeyoshi Itohara, Takuji Iwasato, and Ole Kiehn.  Spinal glutamatergic neurons defined by EphA4 signaling are essential components of normal locomotor circuits.  Journal of Neuroscience, 34(11):3841–3853, Mar. 2014.
  * (16) J. J. Collins and I. N. Stewart.  Coupled nonlinear oscillators and the symmetries of animal gaits.  Journal of Nonlinear Science, 3(1):349–392, Dec. 1993.
  * (17) Ryota Kobayashi, Hiroshi Nishimaru, and Hisao Nishijo.  Estimation of excitatory and inhibitory synaptic conductance variations in motoneurons during locomotor-like rhythmic activity.  Neuroscience, 335:72–81, 2016.
  * (18) Tetsuro Funato, Yuki Yamamoto, Shinya Aoi, Takashi Imai, Toshio Aoyagi, Nozomi Tomita, and Kazuo Tsuchiya.  Evaluation of the phase-dependent rhythm control of human walking using phase response curves.  PLoS Computational Biology, 12(5):e1004950, 2016.
  * (19) Björn Kralemann, Matthias Frühwirth, Arkady Pikovsky, Michael Rosenblum, Thomas Kenner, Jochen Schaefer, and Maximilian Moser.  In vivo cardiac phase response curve elucidates human respiratory heart rate variability.  Nature Communications, 4(1):2418, Sep. 2013.
  * (20) Matthias Hackel.  Humanoid Robots.  IntechOpen, Rijeka, 2007.
  * (21) Mariano Garcia, Anindya Chatterjee, Andy Ruina, and Michael Coleman.  The simplest walking model: Stability, complexity, and scaling.  Journal of Biomechanical Engineering, 120(2):281–288, 04 1998.
  * (22) Yoshiki Kuramoto.  Chemical Oscillations, Waves, and Turbulence.  Springer, Berlin, 1984.
  * (23) Arkady Pikovsky, Michael Rosenblum, and Jürgen Kurths.  Synchronization: A universal concept in nonlinear science.  Cambridge University Press, 2001.
  * (24) Frank C. Hoppensteadt and Eugene M. Izhikevich.  Weakly Connected Neural Networks.  Springer, New York, 1997.
  * (25) Arthur T. Winfree.  The Geometry of Biological Time.  Springer, New York, 2001.
  * (26) G. Bard Ermentrout and David H. Terman.  Mathematical Foundations of Neuroscience.  Springer, New York, 2010.
  * (27) Hiroya Nakao.  Phase reduction approach to synchronisation of nonlinear oscillators.  Contemporary Physics, 57(2):188–214, Oct. 2016.
  * (28) Bharat Monga, Dan Wilson, Tim Matchen, and Jeff Moehlis.  Phase reduction and phase-based optimal control for biological systems: a tutorial.  Biological cybernetics, 113(1):11–46, 2019.
  * (29) Yoshiki Kuramoto and Hiroya Nakao.  On the concept of dynamical reduction: the case of coupled oscillators.  Philosophical Transactions of the Royal Society A, 377(2160):20190041, Oct. 2019.
  * (30) Bard Ermentrout, Youngmin Park, and Dan Wilson.  Recent advances in coupled oscillator theory.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 377(2160):20190092, Oct. 2019.
  * (31) Takahiro Harada, Hisa-Aki Tanaka, Michael J. Hankins, and István Z. Kiss.  Optimal waveform for the entrainment of a weakly forced oscillator.  Phys. Rev. Lett., 105:088301, Aug 2010.
  * (32) Anatoly Zlotnik, Yifei Chen, István Z. Kiss, Hisa-Aki Tanaka, and Jr-Shin Li.  Optimal waveform for fast entrainment of weakly forced nonlinear oscillators.  Phys. Rev. Lett., 111:024102, Jul 2013.
  * (33) Hisa-Aki Tanaka.  Optimal entrainment with smooth, pulse, and square signals in weakly forced nonlinear oscillators.  Physica D: Nonlinear Phenomena, 288:1–22, 2014.
  * (34) Hisa-Aki Tanaka, Isao Nishikawa, JÃ¼rgen Kurths, Yifei Chen, and István Z. Kiss.  Optimal synchronization of oscillatory chemical reactions with complex pulse, square, and smooth waveforms signals maximizes tsallis entropy.  EPL (Europhysics Letters), 111(5):50007, sep 2015.
  * (35) Arkady Pikovsky.  Maximizing coherence of oscillations by external locking.  Phys. Rev. Lett., 115:070602, Aug 2015.
  * (36) Anatoly Zlotnik, Raphael Nagao, István Z. Kiss, and Jr-Shin Li.  Phase-selective entrainment of nonlinear oscillator ensembles.  Nature Communications, 7:10788, Mar 2016.
  * (37) Wei Qiao, John T. Wen, and Agung Julius.  Entrainment control of phase dynamics.  IEEE Transactions on Automatic Control, 62(1):445–450, 2017.
  * (38) Yuzuru Kato, Anatoly Zlotnik, Jr-Shin Li, and Hiroya Nakao.  Optimization of periodic input waveforms for global entrainment of weakly forced limit-cycle oscillators.  Nonlinear Dynamics, 105(3):2247–2263, Aug 2021.
  * (39) Shohei Takata, Yuzuru Kato, and Hiroya Nakao.  Fast optimal entrainment of limit-cycle oscillators by strong periodic inputs via phase-amplitude reduction and Floquet theory.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(9):093124, Aug. 2021.
  * (40) Ludovic Righetti, Jonas Buchli, and Auke Jan Ijspeert.  Adaptive frequency oscillators and applications.  The Open Cybernetics & Systemics Journal, 3(1):64–69, 2009.
  * (41) Auke Jan Ijspeert, Jun Nakanishi, Heiko Hoffmann, Peter Pastor, and Stefan Schaal.  Dynamical movement primitives: Learning attractor models for motor behaviors.  Neural Computation, 25(2):328–373, 02 2013.
  * (42) Mostafa Ajallooeian, Jesse van den Kieboom, Albert Mukovskiy, Martin A. Giese, and Auke J. Ijspeert.  A general family of morphed nonlinear phase oscillators with arbitrary limit cycle shape.  Physica D: Nonlinear Phenomena, 263:41–56, 2013.
  * (43) Mahdi Khoramshahi, Rezvan Nasiri, Mohammad Shushtari, Auke Jan Ijspeert, and Majid Nili Ahmadabadi.  Adaptive natural oscillator to exploit natural dynamics for energy efficiency.  Robotics and Autonomous Systems, 97:51–60, 2017.
  * (44) Junichi Hongu and Daisuke Iba.  Nonlinear oscillator design based on phase reduction method for closed-loop system.  International Journal of Adaptive Control and Signal Processing, 36(8):1958–1997, 2022.
  * (45) Auke Jan Ijspeert.  A connectionist central pattern generator for the aquatic and terrestrial gaits of a simulated salamander.  Biological Cybernetics, 84(5):331–348, Apr 2001.
  * (46) L. Righetti and Auke Jan Ijspeert.  Programmable central pattern generators: an application to biped locomotion control.  In Proceedings 2006 IEEE International Conference on Robotics and Automation, 2006. ICRA 2006., pages 1585–1590, 2006.
  * (47) Mostafa Ajallooeian, Majid Nili Ahmadabadi, Babak Nadjar Araabi, and Hadi Moradi.  Design, implementation and analysis of an alternation-based central pattern generator for multidimensional trajectory generation.  Robotics and Autonomous Systems, 60(2):182–198, 2012.
  * (48) Venus Pasandi, Hamid Sadeghian, Mahdi Keshmiri, and Daniele Pucci.  An integrated programmable CPG with bounded output.  IEEE Transactions on Automatic Control, pages 4658–4673, 2022.
  * (49) A. Ruiz, D.H. Owens, and S. Townley.  Existence, learning, and replication of periodic motions in recurrent neural networks.  IEEE Transactions on Neural Networks, 9(4):651–661, 1998.
  * (50) S. Townley, A. Ilchmann, M.G. Weiss, W. Mcclements, A.C. Ruiz, D.H. Owens, and D. Pratzel-Wolters.  Existence and learning of oscillations in recurrent neural networks.  IEEE Transactions on Neural Networks, 11(1):205–214, 2000.
  * (51) P. Zegers and M.K. Sundareshan.  Trajectory generation and modulation using dynamic neural networks.  IEEE Transactions on Neural Networks, 14(3):520–533, 2003.
  * (52) Y. Kuroe and K. Miura.  A method of oscillatory trajectory generation using recurrent hybrid neural networks.  In Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005., volume 2, pages 706–711 vol. 2, 2005.
  * (53) Guillaume Jouffroy.  Design of oscillatory recurrent neural network controllers with gradient based algorithms.  In 16th European Symposium on Artificial Neural Networks, pages 7–12, 2008.
  * (54) M. Okada, K. Tatani, and Y. Nakamura.  Polynomial design of the nonlinear dynamics for the brain-like information processing of whole body motion.  In Proceedings 2002 IEEE International Conference on Robotics and Automation (Cat. No.02CH37292), volume 2, pages 1410–1415, 2002.
  * (55) Masafumi Okada, Daisuke Nakamura, and Yoshihiko Nakamura.  Hierachical design of dynamics based information processing system for humanoid motion generation.  In Proceedings of the 2nd International Symposium on Adaptive Motion of Animals and Machines, Kyoto, March, pages 4–8, 2003.
  * (56) Venus Pasandi, Aiko Dinale, Mehdi Keshmiri, and Daniele Pucci.  A data driven vector field oscillator with arbitrary limit cycle shape.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 8007–8012, 2019.
  * (57) V. Pasandi, A. Dinale, M. Keshmiri, and D. Pucci.  A programmable central pattern generator with bounded output.  Robotics and Autonomous Systems, 125:103423, 2020.
  * (58) Philip Hartman.  Ordinary Differential Equations.  Willy, New York, 1964.
  * (59) John Guckenheimer and Philip Holmes.  Nonlinear oscillations, dynamical systems, and bifurcations of vector fields.  Springer, New York, 1983.
  * (60) Eric Brown, Jeff Moehlis, and Philip Holmes.  On the phase reduction and response dynamics of neural oscillator populations.  Neural Computation, 16(4):673–715, Apr. 2004.
  * (61) Norihisa Namura, Shohei Takata, Katsunori Yamaguchi, Ryota Kobayashi, and Hiroya Nakao.  Estimating asymptotic phase and amplitude functions of limit-cycle oscillators from time series data.  Phys. Rev. E, 106:014204, Jul 2022.
  * (62) Balth van der Pol.  Forced oscillations in a circuit with non-linear resistance.(Reception with reactive triode).  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 3(13):65–80, 1927.
  * (63) Balth van der Pol and Jan van der Mark.  Frequency demultiplication.  Nature, 120(3019):363–364, 1927.
  * (64) Balth van der Pol.  On relaxation-oscillations.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11):978–992, 1926.
  * (65) J. Nagumo, S. Arimoto, and S. Yoshizawa.  An active pulse transmission line simulating nerve axon.  Proceedings of the IRE, 50(10):2061–2070, 1962.


# Sophisticated Students in Boston Mechanism and Gale-Shapley Algorithm for
School Choice Problem

Abhishek Paudel  George Mason UniversityFairfaxVirginiaUSA

###### Abstract.

We present our experimental results of simulating the school choice problem
which deals with the assignment of students to schools based on each group’s
complete preference list for the other group using two algorithms: Boston
mechanism and student-proposing Gale-Shapley algorithm. We compare the effects
of sophisticated students altering their preference lists with regards to
these two algorithms. Our simulation results show that sophisticated students
can benefit more in Boston mechanism compared to Gale-Shapley algorithm based
on multiple evaluation metrics.

##  1\. Introduction

Finding a set of matches between students and schools is an interesting
problem commonly known as the school choice problem. Various algorithms have
been used to find best matching between students and schools in real world
scenario. In this paper, we study two variants of algorithms that are used to
solve this problem: Boston mechanism also known as immediate acceptance
algorithm and student-proposing Gale-Shapley algorithm also known as deferred
acceptance algorithm. We introduce the school choice problem as follows.

###  1.1. School Choice Problem

The school choice problem(Abdulkadiroğlu and Sönmez, 2003) consists of the
following.

  * •

A set of $n$ students $A=\\{a_{1},a_{2},...a_{n}\\}$.

  * •

A set of $m$ schools $B=\\{b_{1},b_{2},...b_{m}\\}$.

  * •

Each school $b_{j}\in B$ has a capacity of $q_{j}$ which is the maximum number
of students it can admit.

  * •

Each student $a_{i}\in A$ has a strict preference list $P_{i}$ that ranks all
schools in $B$. The set of all students’ preference lists is
$P=\\{P_{1},P_{2},...P_{n}\\}$

  * •

Each school $b_{j}\in B$ has a strict preference list $Q_{j}$ that ranks all
students in $A$. The set of all schools’ preference lists is
$Q=\\{Q_{1},Q_{2},...Q_{m}\\}$

The task is to find a matching $M:A\rightarrow B$ such that each student
$a_{i}$ is matched to one and only one school, and each school $b_{j}$ is
matched to no more than $q_{j}$ students. By taking into account the
preference lists of schools and students, different algorithms can find
different sets of matches with different characteristics.

###  1.2. Boston Mechanism

Boston mechanism is one of the widely used matching algorithms used in the
United States(Miralles, 2009). The Boston mechanism tries to match as many
students to their most preferred school as possible. In each round $k$, we
take all unassigned students and take the schools with $k$th priority in these
students’ priority lists. Then, these schools are assigned to the students
based on the school’s priority over the students. If the school hasn’t reached
its maximum capacity, the student is admitted, else the student remains
unassigned in this round and will be attempted to match again in the next
round. This is repeated until no student remains unassigned. Note that the
maximum number of rounds we need to iterate over until no students remain
unassigned is equal to the number of schools assuming that the cumulative
capacity of all schools is greater than or equal to the total number of
students.

###  1.3. Gale-Shapley Algorithm

Gale-Shapley algorithm(Gale and Shapley, 1962) is known for its stable
matching properties which is desired in many of such matching problems. In the
student-proposing variant of the Gale-Shapley algorithm, each student picks
their most preferred school that hasn’t already rejected them previously. If
the school hasn’t already reached its capacity, the student is provisionally
admitted to the school. However, if the school has reached its capacity, the
school looks at its priority list over the students. If this student has
higher priority than any of the provisionally admitted students, then the
admitted student with lowest priority is unassigned from the school, the new
student with higher priority is admitted in its place. Once there remains no
unassigned students, the obtained provisional assignments become the final
assignments. Note that we assume that the cumulative capacity of all schools
is greater than or equal to the total number of students.

###  1.4. Sophisticated and Sincere Students

The two algorithms described above result in different sets of matches between
students and schools. And depending upon the preference lists of students and
schools, the matches generated by each of these algorithms can also be
different. This suggests that by altering the preference lists over schools, a
student can get admitted to a different school than they would otherwise be
admitted to if they hadn’t altered preferences(Pathak and Sönmez, 2008). We
call such students sophisticated students who will alter their preferences so
as to get admitted to a different school which might otherwise have not been
possible if they didn’t alter their preference. Sincere students are those who
will report their true preferences.

##  2\. Methodology

We aim to study the effects of sophisticated students altering their
preference so as to benefit in the process by experimentally simulating the
school choice problem. We generate a set of students, a set of schools with
certain capacities, and their preferences. Then we use Boston mechanism and
student-proposing Gale-Shapley algorithm to find the matches. Next, we
randomly sample certain fraction of the students as sophisticated students who
will alter their preference lists using various strategies. We then obtain new
sets of matches corresponding to these new preferences using Boston mechanism
and student-proposing Gale-Shapley algorithm. Finally, we evaluate the
properties of resulting matches obtained after preference alteration and
compare them between both algorithms.

###  2.1. Preference List Generation

We generate preference lists for students and schools by taking into
consideration various factors. The preference lists of students for schools
are generated based on the following four features: distance to schools,
whether they have a sibling studying at a school, school tier and random
factor. The preference lists of schools are generated based only upon the
distance to students. For schools, we discretize the generated preference
values based on distance into four values 1, 2, 3 and 4 where lower values
denotes higher preference of a student in a school. For generating the
distances, we first randomly generate students and schools in a 2D plane where
x coordinates and y coordinates are uniformly sampled between 0 and 1. Then,
we calculate the pairwise distance between them. Next, we sample whether a
student has a sibling or not based on Bernoulli distribution, and assign this
sibling to random school. The school tiers are sampled with varying
probabilities from tier-1, tier-2, tier-3 or tier-4 where tier-1 represents a
top tier school and tier-4 represents a bottom tier school. Random factor is
uniformly sampled between 0 and 1.

After generating the four features as described above, we combine all these
features to generate the preference values of students for each school. The
preference value is generated by normalizing the above features appropriately
and taking weighted sum of these features. A lower preference value for a
school means that the student prefers that school more compared to another
school with a higher preference value.

###  2.2. Preference List Alteration

We sample a varying fraction of sophisticated and sincere students from the
student pool. We choose sophisticated students by randomly choosing students
from the student pool without replacement until desired number of
sophisticated students are sampled. Sophisticated students will then have
their preference lists altered using various strategies. We employ three
strategies that sophisticated students use to alter their preferences.

  1. (1)

Strategy A: Find a school that most students have on rank 1. And if you also
have that school on rank 1, swap schools on rank 1 and rank 2.

  2. (2)

Strategy B: Find a school that most students have on rank 1. And if you also
have that school on rank 1, move that school to the last in your preference
list.

  3. (3)

Strategy C: Find a school that most students have on the bottom half of their
preference lists (a proxy for less popular school). And if you have that
school on top half of your preference list, move that school to rank 1.

###  2.3. Evaluation Metrics

We use the following evaluation metrics as appropriate based on the preference
alteration strategies used.

  1. (1)

EM-Top3: The percentage of sophisticated students that got admitted to a
school that is on top 3 in their true preference list after they altered their
preference list. Before altering preference, they would not be admitted to
their top 3 schools.

  2. (2)

EM-Higher: The percentage of sophisticated students that got admitted to a
school that is ranked higher in their true preference list after they altered
their preference list.

  3. (3)

EM-Selected: The percentage of sophisticated students that got admitted to a
specific selected school after they altered their preference list. This metric
is specifically useful for Strategy C with the less popular school as selected
school.

##  3\. Experiments

We run our experiments for 100 simulations with the number of students fixed
to 2000 and number of schools fixed to 20 for all our simulation instances
with different instances of preference lists of students and schools and
aggregate the results from all simulations. For each instance of simulation,
school capacities are generated randomly for all schools such that cumulative
capacities of all schools are equal to the total number of students. We then
randomly generate preference lists for students and schools as described in
the previous section. The preference values of schools are discretized into 1,
2, 3 and 4 corresponding to the bins 0, 0.3, 0.5, 0.7 and 1 respectively so
that each student has a priority of 1, 2, 3 or 4 where 1 denotes highest
priority and 4 denotes lowest priority. Each feature for preference list
generation of students are weighted as follows: distance to schools (0.5),
sibling (0.2), school tier (0.2) and random factor (0.1). The Bernoulli
distribution for sampling siblings uses the success probability of 0.5 meaning
that approximately half of the students have a sibling in one of the schools.
For sampling the school tiers, we sample tier-1, tier-2, tier-3 and tier-4
schools with probabilities 0.1, 0.2, 0.3 and 0.4 respectively.

###  3.1. Strategy A

The results for preference alteration using Strategy A are shown in Figure 1.
We observe similar trends for both evaluation metrics EM-Higher and EM-Top3
showing that more percentage of sophisticated students are able to benefit in
Boston mechanism compared to Gale-Shapley algorithm. This remains true for any
number of sophisticated students in the student pool. We also observe that as
the number of sophisticated students increases, less percentage of them are
able to benefit using this strategy. However, after the number of
sophisticated students increases to more than 1400, more fraction of them are
benefited. This could be because when a large number of students prefer one
school over other, there are going to be some sophisticated students who will
have preferred more a school that other students preferred less, and hence are
likely to get admitted to that school.

Figure 1. Percentage of sophisticated students who got admitted to a higher
ranked school (left) and their top-3 school (right) after altering their
preferences using Strategy A (Results are averaged over 100 trials)

###  3.2. Strategy B

Results for preference alteration using Strategy B are shown in Figure 2. We
observe that more percentage of sophisticated are benefited when there the
number of sophisticated students is small. We also see that Boston mechanism
is more beneficial to sophisticated students compared to Gale-Shapley
algorithm. As the number of sophisticated students starts to increase in the
student pool, less and less fraction of them are able to benefit until except
all students become sophisticated at which point there is slight increase in
benefited fraction of sophisticated students likely for the same reason as
with Strategy A explained earlier.

Figure 2. Percentage of sophisticated students who got admitted to a higher
ranked school (left) and their top-3 school (right) after altering their
preferences using Strategy B (Results are averaged over 100 trials)

###  3.3. Strategy C

This strategy is somewhat of a proxy for finding a less popular school that
one likes and putting that school on rank 1. We evaluate this strategy based
on what fraction of sophisticated students are able to get admitted to this
less popular school that they like as represented by the evaluation metric EM-
Selected with the selected school being the less popular school that they
like. The results are shown in 3. With this strategy, we see a similar trend
where as the number of sophisticated students increase, lesser fraction of
them are benefited with this strategy. However, this strategy benefits the
most compared to other two strategies as seen in 3 which starts at around 18%
to 20% sophisticated students being benefited when there are around 200
sophisticated students. Interestingly, the gap between Boston mechanism and
Gale-Shapley algorithm closes immediately when number of sophisticated
students starts to increase with both algorithms being equally likely to be
beneficial for sophisticated students using Strategy C.

Figure 3. Percentage of sophisticated students who got admitted to a selected
less popular school (EM-Selected) after altering their preferences using
Strategy C (Results are averaged over 100 trials)

##  4\. Conclusion and Future Direction

Overall, we see that Boston mechanism is more beneficial to sophisticated
students than Gale-Shapley algorithm in almost all of our experiments. This
suggests Gale-Shapley algorithm is more robust when it comes to sophisticated
students trying to game the matching mechanism. We also find that in general,
lesser the number of sophisticated students in the student pool, more likely
they will benefit by altering their preferences. Additionally , this is also
more true with Boston mechanism than with Gale-Shapley algorithm. The future
direction of this research could be headed to study the effects of preference
alteration by sophisticated students in sincere students. It would also be
interesting study various other types of matching algorithms that are derived
from Boston mechanism and Gale-Shapley algorithm for various matching settings
in real world scenario.

## References

  * (1)
  * Abdulkadiroğlu and Sönmez (2003) Atila Abdulkadiroğlu and Tayfun Sönmez. 2003.  School Choice: A Mechanism Design Approach.  _American Economic Review_ 93, 3 (June 2003), 729–747.  https://doi.org/10.1257/000282803322157061
  * Gale and Shapley (1962) D. Gale and L. S. Shapley. 1962.  College Admissions and the Stability of Marriage.  _The American Mathematical Monthly_ 69, 1 (Jan. 1962), 9.  https://doi.org/10.2307/2312726
  * Miralles (2009) Antonio Miralles. 2009\.  School choice: The case for the Boston mechanism. In _International conference on auctions, market mechanisms and their applications_. Springer, 58–60.
  * Pathak and Sönmez (2008) Parag A. Pathak and Tayfun Sönmez. 2008.  Leveling the Playing Field: Sincere and Sophisticated Players in the Boston Mechanism.  _American Economic Review_ 98, 4 (September 2008), 1636–52.  https://doi.org/10.1257/aer.98.4.1636