sufficiently small (depending on $c,C_{0}$ and $\varepsilon_{3}$), then

$|\dot{f}_{y,r}(x_{1})-\dot{f}_{y^{\prime},r^{\prime}}(x_{1})|\geq\frac{\varepsilon_{3}}{10}|y-y^{\prime}|,$
(A.13)

and thus by (A.10), (A.5) holds with
$\varepsilon=\frac{\varepsilon_{3}}{30C}$.

To establish (A.13), we will make use of the following inequality

$\Big{|}\frac{A}{B}-\frac{A^{\prime}}{B^{\prime}}\Big{|}\geq|A-A^{\prime}|\frac{|B|}{|BB^{\prime}|}-|A|\frac{|B-B^{\prime}|}{|BB^{\prime}|}$
(A.14)

In practice, when we apply Inequality (A.14), $B$ and $B^{\prime}$ will always
have $\frac{1}{2}\leq|B|\leq 2$ and similarly for $B^{\prime}$;
$|A-A^{\prime}|\gtrsim|y-y^{\prime}|$; $|B-B^{\prime}|\lesssim|y-y^{\prime}|$
(with implicit constants that might depend on $c$ and $C_{0}$), and
$|A|\leq\varepsilon_{1}$. If $A,A^{\prime},B,B^{\prime}$ satisfy these
conditions, then we will have
$|A/B-A^{\prime}/B^{\prime}|\geq\frac{1}{8}|A-A^{\prime}|$, provided
$\varepsilon_{1}$ is chosen sufficiently small.

We have

$\begin{split}|\dot{f}_{y,r}(x_{1})-\dot{f}_{y^{\prime},r^{\prime}}(x_{1})|&=\Big{|}\frac{\Phi_{1}(x,y)}{\Phi_{2}(x,y)}-\frac{\Phi_{1}(x^{\prime},y^{\prime})}{\Phi_{2}(x^{\prime},y^{\prime})}\Big{|}\\\
&\geq\Big{|}\frac{\Phi_{1}(x^{\prime},y)}{\Phi_{2}(x^{\prime},y)}-\frac{\Phi_{1}(x^{\prime},y^{\prime})}{\Phi_{2}(x^{\prime},y^{\prime})}\Big{|}-\Big{|}\frac{\Phi_{1}(x,y)}{\Phi_{2}(x,y)}-\frac{\Phi_{1}(x^{\prime},y)}{\Phi_{2}(x^{\prime},y)}\Big{|}.\end{split}$
(A.15)

By (A.1), (A.6), and (A.11), we have

$\Big{|}\frac{\Phi_{1}(x,y)}{\Phi_{2}(x,y)}-\frac{\Phi_{1}(x^{\prime},y)}{\Phi_{2}(x^{\prime},y)}\Big{|}\leq
4C|x-x^{\prime}|\leq 4C_{0}\varepsilon_{2}|y-y^{\prime}|.$

Applying Inequality (A.14) and using assumption (A.12), we conclude that

$\Big{|}\frac{\Phi_{1}(x^{\prime},y)}{\Phi_{2}(x^{\prime},y)}-\frac{\Phi_{1}(x^{\prime},y^{\prime})}{\Phi_{2}(x^{\prime},y^{\prime})}\Big{|}\geq\frac{1}{8}|\Phi_{1}(x^{\prime},y)-\Phi_{1}(x^{\prime},y^{\prime})|\geq\frac{\varepsilon_{3}}{8}|y-y^{\prime}|.$

Thus if $\varepsilon_{2}$ is chosen sufficiently small depending on $C_{0}$
and $\varepsilon_{3}$, then (A.13) follows from (A.15).

Step 3
Suppose now that (A.12) is false, i.e.

$|\Phi_{1}(x^{\prime},y)-\Phi_{1}(x,y^{\prime})|<\varepsilon_{3}|y-y^{\prime}|.$
(A.16)

Then by (A.7), we have

$\Big{|}\frac{\Phi_{11}(x^{\prime},y)}{|\nabla\Phi(x,y)|}-\frac{\Phi_{11}(x^{\prime},y^{\prime})}{|\nabla\Phi(x^{\prime},y^{\prime})|}\Big{|}\geq\frac{c}{4}|y-y^{\prime}|.$
(A.17)

We will show that if $\varepsilon_{1},\varepsilon_{2},\varepsilon_{3}$ were
chosen sufficiently small, then

$|\ddot{f}_{y,r}(x_{1})-\ddot{f}_{y^{\prime},r^{\prime}}(x_{1})|\geq\frac{c}{10}|y-y^{\prime}|,$
(A.18)

and thus by (A.10), (A.5) holds with $\varepsilon=\frac{c}{30C}.$ Thus (A.18)
implies that set $M$ defined in (A.4) is a cinematic family of functions. It
remains to establish (A.18). The proof of this inequality uses the same ideas
as the proof of (A.13), though the steps are slightly more complicated.

Define

$F(a,b)=\frac{\Phi_{11}(a,b)\Phi_{2}(a,b)^{2}-2\Phi_{12}(a,b)\Phi_{1}(a,b)\Phi_{2}(a,b)+\Phi_{22}(a,b)\Phi_{1}(a,b)^{2}}{\Phi_{2}(a,b)^{3}}$

Then

$|\ddot{f}_{y,r}(x_{1})-\ddot{f}_{y^{\prime},r^{\prime}}(x_{1})|=|F(x,y)-F(x^{\prime},y^{\prime})|\geq|F(x^{\prime},y)-F(x,y^{\prime})|-|F(x,y)-F(x^{\prime},y)|.$

By (A.1), (A.6), and (A.11), we have

$|F(x,y)-F(x^{\prime},y)|\lesssim C|x-x^{\prime}|\lesssim
C\varepsilon_{2}|y-y^{\prime}|.$ (A.19)

Thus our goal is to establish a lower bound for
$|F(x^{\prime},y)-F(x,y^{\prime})|$. We will estimate

$\begin{split}|F(x^{\prime},y)-F(x,y^{\prime})|\geq\Big{|}&\frac{\Phi_{11}(x^{\prime},y)}{\Phi_{2}(x^{\prime},y)}-\frac{\Phi_{11}(x^{\prime},y^{\prime})}{\Phi_{2}(x^{\prime},y^{\prime})}\Big{|}\\\
&-2\Big{|}\frac{\Phi_{12}(x^{\prime},y)\Phi_{1}(x^{\prime},y)}{\Phi_{2}(x^{\prime},y)^{2}}-\frac{\Phi_{12}(x^{\prime},y^{\prime})\Phi_{1}(x^{\prime},y^{\prime})}{\Phi_{2}(x^{\prime},y^{\prime})^{2}}\Big{|}\\\
&-\Big{|}\frac{\Phi_{22}(x^{\prime},y)\Phi_{1}(x^{\prime},y)^{2}}{\Phi_{2}(x^{\prime},y)^{3}}-\frac{\Phi_{22}(x^{\prime},y^{\prime})\Phi_{1}(x^{\prime},y^{\prime})^{2}}{\Phi_{2}(x^{\prime},y^{\prime})^{3}}\Big{|}.\end{split}$
(A.20)

For the first term, we will use the inequality
$AB-A^{\prime}B^{\prime}=(A-A^{\prime})B^{\prime}+(B-B^{\prime})A$ to write

$\begin{split}&\frac{\Phi_{11}(x^{\prime},y)}{\Phi_{2}(x^{\prime},y)}-\frac{\Phi_{11}(x^{\prime},y^{\prime})}{\Phi_{2}(x^{\prime},y^{\prime})}=\frac{\Phi_{11}(x^{\prime},y)}{|\nabla\Phi(x^{\prime},y)|}\Big{(}\frac{\Phi_{1}^{2}(x^{\prime},y)}{\Phi_{2}^{2}(x^{\prime},y)}+1\Big{)}^{1/2}-\frac{\Phi_{11}(x^{\prime},y^{\prime})}{|\nabla\Phi(x^{\prime},y^{\prime})|}\Big{(}\frac{\Phi_{1}^{2}(x^{\prime},y^{\prime})}{\Phi_{2}^{2}(x^{\prime},y^{\prime})}+1\Big{)}^{1/2}\\\
&=\Big{(}\frac{\Phi_{11}(x^{\prime},y)}{|\nabla\Phi(x^{\prime},y)|}-\frac{\Phi_{11}(x^{\prime},y^{\prime})}{|\nabla\Phi(x^{\prime},y^{\prime})|}\Big{)}\Big{(}\frac{\Phi_{1}^{2}(x^{\prime},y^{\prime})}{\Phi_{2}^{2}(x^{\prime},y^{\prime})}+1\Big{)}^{1/2}+\Big{[}\Big{(}\frac{\Phi_{1}^{2}(x^{\prime},y)}{\Phi_{2}^{2}(x^{\prime},y)}+1\Big{)}^{1/2}-\Big{(}\frac{\Phi_{1}^{2}(x^{\prime},y^{\prime})}{\Phi_{2}^{2}(x^{\prime},y^{\prime})}+1\Big{)}^{1/2}\Big{]}\frac{\Phi_{11}(x^{\prime},y)}{|\nabla\Phi(x^{\prime},y)|}\end{split}$

Using (A.8) and (A.6), and then (A.17), We have

$\begin{split}\Big{|}\Big{(}\frac{\Phi_{11}(x^{\prime},y)}{|\nabla\Phi(x^{\prime},y)|}-\frac{\Phi_{11}(x^{\prime},y^{\prime})}{|\nabla\Phi(x^{\prime},y^{\prime})|}\Big{)}\Big{(}\frac{\Phi_{1}^{2}(x^{\prime},y^{\prime})}{\Phi_{2}^{2}(x^{\prime},y^{\prime})}+1\Big{)}^{1/2}\Big{|}&\geq\frac{1}{2}\Big{|}\frac{\Phi_{11}(x^{\prime},y)}{|\nabla\Phi(x^{\prime},y)|}-\frac{\Phi_{11}(x^{\prime},y^{\prime})}{|\nabla\Phi(x^{\prime},y^{\prime})|}\Big{|}\\\
&\geq\frac{c}{4}|y-y^{\prime}|,\end{split}$

and by (A.8),

$\begin{split}&\Big{[}\Big{(}\frac{\Phi_{1}^{2}(x^{\prime},y)}{\Phi_{2}^{2}(x^{\prime},y)}+1\Big{)}^{1/2}-\Big{(}\frac{\Phi_{1}^{2}(x^{\prime},y^{\prime})}{\Phi_{2}^{2}(x^{\prime},y^{\prime})}+1\Big{)}^{1/2}\Big{]}\frac{\Phi_{11}(x^{\prime},y)}{|\nabla\Phi(x^{\prime},y)|}\\\
&\leq\Big{[}\Big{(}\frac{\Phi_{1}^{2}(x^{\prime},y)}{\Phi_{2}^{2}(x^{\prime},y)}+1\Big{)}^{1/2}-\Big{(}\frac{\Phi_{1}^{2}(x^{\prime},y^{\prime})}{\Phi_{2}^{2}(x^{\prime},y^{\prime})}+1\Big{)}^{1/2}\Big{]}(2\varepsilon_{1})\\\
&\leq(2C)|y-y^{\prime}|(2\varepsilon_{1}).\end{split}$

Thus if $\varepsilon_{1}$ is selected sufficiently small, we have

$\Big{|}\frac{\Phi_{11}(x^{\prime},y)}{\Phi_{2}(x^{\prime},y)}-\frac{\Phi_{11}(x^{\prime},y^{\prime})}{\Phi_{2}(x^{\prime},y^{\prime})}\Big{|}\geq\frac{c}{5}|y-y^{\prime}|.$

Bounding the second and third term on the RHS of (A.20) is straightforward,
using the estimates (A.6), (A.8), and (A.16). We conclude that

$|F(x^{\prime},y)-F(x,y^{\prime})|\geq\frac{c}{6}|y-y^{\prime}|.$

Combining this with (A.19), we obtain (A.18).

## References

  * [1] P. Agarwal and M. Sharir.  Davenport-Schinzel sequences and their geometric applications.  Cambridge University Press, Cambridge, 1995.
  * [2] J. Bourgain.  The discretized sum-product and projection theorems.  J. Anal. Math., 112:193–236, 2010.
  * [3] A. Chang, G. Dosidis, and J. Kim.  Nikodym sets and maximal functions associated with spheres.  In preparation, 2022+.
  * [4] K.L. Clarkson, H Edelsbrunner, L.J. Guibas, M Sharir, and E. Welzl.  Combinatorial complexity bounds for arrangements of curves and spheres.  Discrete Comput. Geom., 5:99–160, 1990.
  * [5] K. Falconer, J. Fraser, and X. Jin.  Sixty years of fractal projections.  In Christoph Bandt, Kenneth Falconer, and Martina Zähle, editors, Fractal Geometry and Stochastics V, pages 3–25. Springer International Publishing, 2015.
  * [6] K. Fässler and T. Orponen.  On restricted families of projections in $\mathbb{R}^{3}$.  Proc. Lond. Math. Soc. (3), 109(2):353–381, 2014.
  * [7] S. Gan, S. Guo, L. Guth, T. L. J. Harris, D. Maldague, and H. Wang.  On restricted projections to planes in $\mathbb{R}^{3}$.  arXiv:2207.13844, 2022.
  * [8] L. Guth and N.H. Katz.  On the Erdős distinct distances problem in the plane.  Ann. of Math., 181(1):155–190, 2015.
  * [9] T. L. J. Harris.  Length of sets under restricted families of projections onto lines.  arXiv:2208.06896, 2022.
  * [10] W. He.  Orthogonal projections of discretized sets.  J. Fractal Geom., 7(3):271–317, 2020.
  * [11] A. Käenmäki, T. Orponen, and L. Venieri.  A marstrand-type restricted projection theorem in $\mathbb{R}^{3}$.  arXiv: Classical Analysis and ODEs, 2021.  arXiv:1708.04859.
  * [12] N. H. Katz and T. Tao.  Some connections between Falconer’s distance set conjecture and sets of Furstenburg type.  New York J. Math., 7:149–187, 2001.
  * [13] N.H. Katz, S. Wu, and J. Zahl.  The $SL_{2}$ Kakeya conjecture.  In preparation, 2022+.
  * [14] R. Kaufman.  On Hausdorff dimension of projections.  Mathematika, 15:153–155, 1968.
  * [15] L. Kolasa and T. Wolff.  On some variants of the Kakeya problem.  Pacific J. Math., 190(1):111–154, 1999.
  * [16] J. Liu.  Dimension estimates on circular $(s,t)$-Furstenberg sets.  arXiv:2204.01770, 2022.
  * [17] A. Marcus and G. Tardos.  Intersection reverse sequences and geometric applications.  J. Combin. Theory Ser. A, 113(4):675–691, 2006.
  * [18] J.M. Marstrand.  Some fundamental geometrical properties of plane sets of fractional dimensions.  Proc. Lond. Math. Soc., 4(3):257–302, 1954.
  * [19] P. Mattila.  Hausdorff dimension, orthogonal projections and intersections with planes.  Ann. Acad. Sci. Fenn. A Math., 1:227–244, 1975.
  * [20] P. Mattila.  Geometry of Sets and Measure in Euclidean Spaces.  Cambridge Studies in Advanced Mathematics. Cambridge University Press, Cambridge, 1995.
  * [21] W. Schlag.  On continuum incidence problems related to harmonic analysis.  J. Funct. Anal., 201(2):480–521, 2003.
  * [22] C. Sogge.  Propagation of singularities and maximal functions in the plane.  Inv. Math., 104:349–376, 1991.
  * [23] T. Tao.  The two-ends reduction for the Kakeya maximal conjecture, 2009.  https://terrytao.wordpress.com/2009/05/15/the-two-ends-reduction-for-the-kakeya-maximal-conjecture/.
  * [24] H. Wang and J. Zahl.  Sticky Kakeya sets in $\mathbb{R}^{3}$.  In preparation, 2022+.
  * [25] T. Wolff.  A Kakeya-type problem for circles.  Amer. J. Math., 119(5):985–1026, 1997.
  * [26] T. Wolff.  Recent work connected with the Kakeya problem.  In Prospects in mathematics (Princeton, NJ, 1996), pages 129–162. Amer. Math. Soc., Providence, RI, 1999.
  * [27] T. Wolff.  Local smoothing type estimates on $L^{p}$ for large $p$.  Geom. Funct. Anal., 10(5):1237–1288, 2000.
  * [28] J. Zahl.  $L^{3}$ estimates for an algebraic variable coefficient Wolff circular maximal function.  Rev. Mat. Iberoam., 28(4):1281–1295, 2012.
  * [29] J. Zahl.  On the Wolff circular maximal function.  Illinois J. Math., 56(4):1281–1295, 2012.


# Loss-driven topological transitions in lasing

Grazia Salerno<EMAIL_ADDRESS>Rebecca Heilmann  Kristian Arjas
Kerttu Aronen  Jani-Petri Martikainen  Päivi Törmä<EMAIL_ADDRESS>Department of Applied Physics, Aalto University School of Science, P.O. Box
15100, Aalto, FI-00076, Finland

###### Abstract

We experimentally observe lasing in a hexamer plasmonic lattice and find that
when tuning the scale of the unit cell, the polarization properties of the
emission change. By a theoretical analysis we identify the lasing modes as
quasi bound states in continuum (quasi-BICs) of topological charges of zero,
one or two. A T-matrix simulation of the structure reveals that the mode
quality(Q)-factors depend on the scale of the unit cell, with highest-Q modes
favored by lasing. The system thus shows a loss-driven transition between
lasing in modes of trivial and high-order topological charge.

Bound states in the continuum (BICs) are peculiar eigenstates of quantum or
classical wave systems: their energies lie within a continuum of other states
but they are nevertheless completely uncoupled from those states. For optical
structures such as gratings, cavities and photonic crystals this means that
their BIC modes are decoupled from free-space radiation and appear as dark
states: light cannot in- or out-couple from/to the far-field at the energy and
momentum of the BIC Hsu  _et al._ (2016); Hwang  _et al._ (2022). From the
fundamental perspective, photonic BICs can be viewed as vortex centers of
light polarization. When one measures the polarization of light radiated by
the structure from different angles, corresponding to different momentum, the
polarization winds around the BIC momentum. Due to this winding, it is
impossible to define what would be the polarization exactly at the BIC
momentum, which therefore must appear dark, similarly as the core of a vortex
is empty. This description has naturally connected the BICs to the realm of
more general topological phenomena, including topological robustness and
topological transitions Ozawa  _et al._ (2019); Rider  _et al._ (2019); Liu
_et al._ (2021); Rider  _et al._ (2022). The polarization vortex is associated
with the existence of a protected and quantized topological charge which tells
how many times the polarization winds in a path circulating the BIC momentum
Doeleman  _et al._ (2018); Zhen  _et al._ (2014). This charge is conserved and
cannot be removed with small perturbations. In order to observe BICs, they can
be weakly coupled to the radiative continuum via an intentionally designed
leakage mechanism or due to material loss Azzam  _et al._ (2021); Kodigala
_et al._ (2017); Ha  _et al._ (2018); the presence of an edge can also act as
a leakage channel Heilmann  _et al._ (2022); Asamoah  _et al._ (2022). Such
leaky or lossy BICs are called quasi-BICs. Quasi-BICs can have extremely large
Q-factors Iwahashi  _et al._ (2011); Zito  _et al._ (2019); Wang  _et al._
(2020); Huang  _et al._ (2020); Dyakov  _et al._ (2021); Zhang  _et al._
(2021); Hwang  _et al._ (2021); Wu  _et al._ (2022); Kang  _et al._ (2022),
which is promising for numerous applications in lasing, filtering, or sensing
Gentry and Popović (2014); Kodigala  _et al._ (2017); Ha  _et al._ (2018);
Guan  _et al._ (2020a, b); Silva  _et al._ (2020).

BICs of topological charge $|q|=1$ have been observed in many systems, from
photonic crystals to nanoplasmonic structures Doeleman  _et al._ (2018); Zhen
_et al._ (2014); De Angelis  _et al._ (2019); Jin  _et al._ (2019); Heilmann
_et al._ (2022); Mohamed  _et al._ (2022); Asamoah  _et al._ (2022); Kodigala
_et al._ (2017). Switching between non-degenerate BICs with different
topological charge of $|q|\leq 1$ was recently achieved Mohamed  _et al._
(2022). The interest in BICs with large topological charges ($|q|>1$) and
their generation, evolution, or annihilation is only just starting. BICs of
high-order charge are predicted to exist in systems with a six-fold rotational
symmetry Zhen  _et al._ (2014); Bai  _et al._ (2021) and it has been
theoretically proposed that a BIC of charge $|q|=2$ can undergo a topological
transition and break into two BICs of charge $|q|=1$ Yoda and Notomi (2020).
The first experimental studies demonstrating vortex beams generated by BIC of
higher topological charge in photonic crystals have recently emerged Zhang
_et al._ (2018); Wang  _et al._ (2020), while observation of lasing in BICs of
higher topological charges still lacks in plasmonic lattices.

In this Letter, we observe lasing in a plasmonic lattice from a BIC of
topological charge $q=-2$ and demonstrate several topological transitions
driven by losses as a structural parameter of the unit cell is changed. We use
plasmonic nanoparticle arrays due to their strong dipole moments in the
nanoparticles Moerland  _et al._ (2017); Wang  _et al._ (2018a, b); Väkeväinen
_et al._ (2014); Hakala  _et al._ (2018); Kravets  _et al._ (2018). The lasing
mode is identified by using theory of coupled dipoles together with group
theory and T-matrix calculations Nečada and Törmä (2021); Nec (2022). We show
that the energy of the modes as well as their loss properties change as the
structural parameter is varied. The Q-factor of the modes changes, and the
mode with the highest Q-factor lases. By comparing the theory with the
experimental measurements, we characterize the polarization vorticity of the
lasing mode and confirm topological transitions between a bright mode and
quasi-BIC modes of topological charges $q=-2$, $q=-1$ and $q=+1$.

Figure 1: (a) The hexamer is composed of six particles (black disks) arranged
in a regular hexagon of side length $b$. Purple and orange arrows represent
the longitudinal $\mathbf{e}_{L}^{n}$ and transverse $\mathbf{e}_{T}^{n}$ unit
vectors used in the simple hexamer model, having indicated only nearest-
neighbor. The hexamer constitutes the unit cell of a triangular lattice of
periodicity $a$ (red points). Due to a finite size of the particles, whose
diameter is $d$, the edge-to-edge distance between nanoparticles is $t=b-d$.
(b) Modes of an isolated hexamer, where the arrows correspond to the electric
dipole moment. The modes are labeled according to the corresponding IRs of the
symmetry group $C_{6}$. (c-f) Amplitude and phase of the four IRs, filtered
with a polarizer oriented according to the arrow in each panel, while the
phase in the last column is unfiltered. The modes and their respective
topological charges are (c) $A^{\prime}$: $q=+1$; (d) $E_{1}^{\prime}$: $q=0$;
(e) $E_{2}^{\prime}$: $q=-1$, and (f) $B^{\prime}$: $q=-2$. (d-e) show
balanced superpositions of the degenerate modes with zero relative phase.

The topological charge of a BIC is defined as Zhen  _et al._ (2014)

$q=\frac{1}{2\pi}\oint_{\mathcal{C}}d\mathbf{k}\cdot\nabla_{\mathbf{k}}\Phi(\mathbf{k}),$
(1)

where
$\Phi(\mathbf{k})=\arg[\mathbf{p}(\mathbf{k})\cdot\hat{x}+i\mathbf{p}(\mathbf{k})\cdot\hat{y}]$
is the phase of the polarization vector projected onto the $xy$ plane in
momentum space
$\mathbf{p}(\mathbf{k})=(\hat{x}\cdot\langle\mathbf{u}_{\mathbf{k}}(\mathbf{r},z)\rangle)\hat{x}+(\hat{y}\cdot\langle\mathbf{u}_{\mathbf{k}}(\mathbf{r},z)\rangle)\hat{y}$,
and $\langle\mathbf{u}_{\mathbf{k}}\rangle$ indicates the spatial average of
the periodic part of the electric field
$\mathbf{E}_{\mathbf{k}}(\mathbf{r},z)=e^{i\mathbf{k}\cdot\mathbf{r}}\mathbf{u}_{\mathbf{k}}(\mathbf{r},z)$
at constant vertical distance $z$ from the lattice plane much larger than the
unit cell size. The charge $q$ counts the number of times the polarization
vector $\mathbf{p}(\mathbf{k})$ goes around a specific $\mathbf{k}$ point,
e.g. the $\Gamma$-point. The sign of the topological charge gives the
orientation of the winding.

We start by theoretically analyzing the modes of an isolated hexamer, Fig.
1(a). Each particle has an in-plane dipole moment
$\mathbf{p}_{i}=(p_{i}^{x},p_{i}^{y})$, coupled to all other dipoles in a
polarization-dependent way, as in Ref. Heilmann  _et al._ (2022). We set
$\Omega_{L}$ ($\Omega_{T}$) to be the coupling between dipoles oriented
longitudinally (transversely) to the link
$\mathbf{e}_{L}^{j}=(\mathbf{r}_{i}-\mathbf{r}_{i+j})/|\mathbf{r}_{i}-\mathbf{r}_{i+j}|$
(such that $\mathbf{e}_{T}^{j}\cdot\mathbf{e}_{L}^{j}=0$) connecting particles
$i$ and $i+j$ (typically for dipoles $\Omega_{T}\gg\Omega_{L}$). The bare
dipole oscillation frequency is $\omega_{0}\gg\Omega_{L,T}$. The modes are
found from

$\ddot{\mathbf{p}}_{i}\\!=\\!\omega_{0}\mathbf{p}_{i}+\\!\sum_{j\neq
i}\left[\frac{\Omega_{L}}{R_{j}^{3}}(\mathbf{p}_{i+j}\\!\cdot\\!\mathbf{e}_{L}^{j})\mathbf{e}_{L}^{j}\\!+\\!\frac{\Omega_{T}}{R_{j}^{3}}(\mathbf{p}_{i+j}\\!\cdot\\!\mathbf{e}_{T}^{j})\mathbf{e}_{T}^{j}\right],$

where $R_{j}=|\mathbf{r}_{i}-\mathbf{r}_{i+j}|$. In the basis
$(p_{i}^{x},p_{i}^{y})$, we find twelve modes centered around $\omega_{0}$.
The spatial structure of the lower energy set of eigenmodes is shown in Fig.
1(b) where arrows correspond to the electric dipole moment orientation. The
modes belong to the symmetry class of the $C_{6}$ group; the singlets to the
$A$ and $B$ irreducible representations (IRs), the doublets to the
$E_{1}=(E_{1}^{1},E_{1}^{2})$ and the $E_{2}=(E_{2}^{1},E_{2}^{2})$ IRs. The
ordering of these modes does not depend on the interparticle distance, in the
isolated hexamer case as losses are not included. From Finite-Element-Method
simulations on arrays with periodic boundary conditions, we find the same
modes, and also another set of six modes Sup (2022) that have a larger out-of-
plane character. Those are not excited in our setup and therefore we will
focus on the $A^{\prime}$, $E_{1}^{\prime}$, $E_{2}^{\prime}$ and $B^{\prime}$
modes.

From the electric dipole moment orientation in each mode, we calculate the
real-space diffraction pattern generated by the nanoparticle array, together
with the far field in momentum space. BICs have distinct momentum-space
features, as visible in Fig. 1(c-f), where the far field is obtained by
Fourier transforming the real space field with a polarization filter oriented
according to the arrow in each panel; these calculations are done as in Refs.
Guo  _et al._ (2019); Heilmann  _et al._ (2022). The last column in Fig.
1(c-f) shows the phase of each mode’s polarization calculated from the Stokes
parameters ($\Phi=\frac{1}{2}\arctan(S_{1}/S_{2})$) without any polarization
filter. Fig. 1(d) shows that the $E_{1}^{\prime}$ doublet modes are bright
modes, since their emission in momentum space is featureless along all
polarization directions while the phase is uniform. These modes have no
topological charge ($q=0$), and correspond to dipoles oscillating in phase
along the same direction. All other modes are quasi-BICs, and the dipoles
oscillate in such a way that the emission from $\mathbf{k}=0$ is dark. With a
polarization filter, we see a non-zero emission in the vicinity of the
$\Gamma$-point, on two or four opposite lobes of a donut-shaped amplitude
pattern Iwahashi  _et al._ (2011); Wu  _et al._ (2020); Hwang  _et al._
(2021); Mohamed  _et al._ (2022). As the polarization filter is rotated from
vertical polarization in a clockwise direction, the lobes rotate as well,
following the non-trivial phase polarization (see last column). In particular,
the polarized amplitude of mode $A^{\prime}$ in Fig. 1(c) has two lobes that
rotate in the clockwise direction, which correspond to $q=1$. Accordingly, the
phase in the last column of Fig. 1(c) shows a $2\pi$ rotation, corresponding
to a unit winding. Similarly, the doublet modes $E_{2}^{\prime}$ have $q=-1$
in Fig. 1(e), since the two lobes rotate counter-clockwise; this is also
visible in the different phase winding. Finally, the charge $q=-2$ of mode
$B^{\prime}$ in Fig. 1(f) is evidenced by a four-lobe rotation in a counter-
clockwise direction. Consequently, the phase in the last column has an overall
$4\pi$ winding. The topological charge $q_{\text{IR}}$ of a mode corresponding
to an IR can be derived Sup (2022) from the n-fold rotational symmetry
$C_{n}$:

$q_{\text{IR}}=1-\frac{n}{2\pi}\arg(\epsilon_{\text{IR}}),$ (2)

where $\epsilon_{\text{IR}}$ can be found in character tables. Values given by
Eq. (2) for $C_{6}$ equal those given by our analysis of coupled dipoles in
Fig. 1.

Figure 2: (a) Energy dispersion of an array of cylindrical gold nanoparticles
of diameter $d=80$ nm, height $h=50$ nm, arranged in a triangular array of
periodicity $a=680$ nm, with $t=50$ nm, as a function of the $y$-component of
the in-plane momentum. (b) Peak intensity as a function of the pump fluence.
(c) Spectra of the lasing emission at a pump fluence of 1.78 mJ/cm2. (d)
Energies of the $\Gamma$-point modes in the complex plane calculated with the
T-matrix approach, vertically shifted for clarity (without the shift in the
Supplemental Material Sup (2022)). The color scale indicates the value of $t$
[nm]. (e) SEM images of arrays with $t=$ 50, 120, 170, 200 nm from left to
right. (f) Q-factor of the modes as a function of $t$ [nm]. The light-colored
region indicates the highest Q-factor mode, which is expected to be the lasing
mode. The corresponding experimental observations are marked with the matching
dark-colored regions. (g-l) Angle-resolved images of the lasing emission with
different polarization filters in the detection path, indicated by the double-
headed arrow in the lower left corner of each image. The array is pumped with
LCP polarization and pump fluences of 1 mJ/cm2 for t = 50 nm (g), 1.58 mJ/cm2
for t = 120 (h) and t= 170 nm (i), and 0.63 mJ/cm2 for t= 200 nm (l).

We now consider an extended array composed of hexamers arranged in a finite-
size triangular lattice. The edge-to-edge distance between nanoparticles of
the same hexamer is $t$. The geometry of the system is depicted in Fig. 1(a).
We fabricated arrays of cylindrical gold nanoparticles with an edge-to-edge
distance that is varied between $45<t<200$ nm. The system hosts dispersive
plasmonic-photonic modes Moerland  _et al._ (2017); Wang  _et al._ (2018a, b);
Kravets  _et al._ (2018). The dispersion of such an array for $t=50$ nm is
shown in Fig. 2(a). The arrays were combined with a solution of fluorescent
dye molecules (IR-140) with a concentration of 10 mM. Under optical pumping
with a left-circularly polarized (LCP) 100 fs laser pulse (center wavelength
$800$ nm (1.55 eV)), we observe a non-linear increase in the emission at
$k_{y}=0$ with increasing pump fluence as shown in Fig. 2(b): arrays with
different $t$ show different behavior. Furthermore, the peak intensity differ
as well: the lasing peak of the $t=120$ nm array is about one order of
magnitude larger than the one for the $t=200$ nm array. More details of the
experiments are given in the Supplemental Material Sup (2022), where we also
show the real space emission from the arrays in the lasing regime.

The different thresholds, peak intensities, and lasing energies in Fig. 2(b-c)
hint towards having a qualitatively different mode lasing for different array
geometries. To understand this, we recorded angle- and polarization-resolved
measurements of the lasing emission from arrays with a wide range of different
hexamer sizes. In Fig. 2(g-l) we show angle-resolved images of arrays with
$t=$ 50, 120, 170, and 200 nm with different polarization filters. For $t=50$
nm, we see that the momentum-space emission in Fig. 2(g) mainly concentrates
in four lobes which rotate in a counter-clockwise direction, compatibly with
the $B^{\prime}$ mode of $q=-2$. For $t=120$ nm in Fig. 2(h), there is only a
single lobe with no visible change as the polarization filter is rotated; this
mode is the doublet mode $E_{1}^{\prime}$, which has no topological charge. As
we further increase the size of the unit cell, for $t=170$ nm, in Fig. 2(i) we
see two lobes rotating in a counter-clockwise direction, compatibly with the
doublet $E_{2}^{\prime}$ of $q=-1$. For even higher $t=200$ nm in Fig. 2(l),
the momentum-space emission comes from two lobes that rotate in a clockwise
direction, as for mode $A^{\prime}$ with $q=+1$. These observations directly
correspond to the amplitudes and phases theoretically calculated in Fig. 1(c).

We have employed T-matrix calculations Nečada and Törmä (2021); Nec (2022)
(see also Sup (2022)) to explain why there are four regimes in which the
lasing mode corresponds to a different IR. We obtain the energy of the modes
at the $\Gamma$-point by varying $t$ in the experimentally relevant parameter
range, and plot their real and imaginary parts in the complex plane, in Fig.
2(d). We see that the mode energies mostly lie on a linear trajectory. The
imaginary part of the energy accounts for the losses, which in our system are
radiative and/or Ohmic losses. The modes $A^{\prime}$, $B^{\prime}$, and
$E_{2}^{\prime}$ are dark modes that do not radiate to the continuum: for this
reason, they only experience Ohmic losses, and the complex energy curve has
the same slope for all three. The bright mode $E_{1}^{\prime}$ is subjected to
both Ohmic and radiative losses, thus its energy slope for $t<150$ nm is
steeper than for the three dark modes. However, for $t>150$ nm the slope of
$E_{1}^{\prime}$ gets smaller, indicating that the mode experiences less
radiative losses. In general, all modes’ losses behave differently for small
($t<100$ nm) and large ($t>210$ nm) hexamers; for example, at small $t$, mode
$B^{\prime}$ seems to have a very small imaginary part, and the same is true
for mode $A^{\prime}$ for large $t$ instead. As the modes have different
radiative properties when the hexamer size $t$ is varied, a loss-induced
transition between the modes could occur in lasing (see Sup (2022) for results
also on a quadrumer array). The change of losses for the various modes with
different $t$ is supported by FEM simulations where we see that the modes’
electric field intensity distribution in the unit cell changes significantly
with changing $t$, despite the symmetry of the modes being the same Sup
(2022). The field hot-spots away from the nanoparticles can contribute to how
lossy a given mode and hence its quality factor Q, and also its coupling to
the gain medium.

From the complex energies of the modes in Fig. 2(d), we calculated the
Q-factor of each IR as a function of the hexamer size, see Fig. 2(f). For
small hexamers (light green region $t<100$ nm) the singlet $B^{\prime}$ with
$q=-2$ has the largest Q-factor; the experiment shows lasing from a quasi-BIC
of charge $q=-2$ in the dark green region. For intermediate hexamers (light
violet region for $100<t<145$ nm) the largest Q-factor mode is the
$E_{1}^{\prime}$ doublet with $q=0$, which in the experiment lases for the
dark violet region. For large hexamers (light pink region $t>150$ nm) the
other doublet $E_{2}^{\prime}$ with $q=-1$ has the largest Q-factor in
T-matrix simulations and the experiment evidences lasing with this topological
charge in the dark pink region. Lastly, the experiment shows lasing in the
$A^{\prime}$ IR with $q=+1$ for $t\approx 200$ nm, although this mode has a
lower Q-factor than $E_{2}^{\prime}$. We attribute this inconsistency to the
different spatial overlap of the mode with the gain medium mentioned above,
causing mode competition between $A^{\prime}$ and $E_{2}^{\prime}$ Sup (2022).

In an equilibrium topological transition, the topological character of the
ground state changes. In our non-equilibrium lasing system, modes of different
topological charge coexist with nearly identical energies, but the mutual
ordering of their Q-factors changes with a structural parameter. Therefore
lower losses (effectively higher gain) can drive a transition of lasing
between modes of different topological charge. We have demonstrated lasing in
a topologically trivial bright mode and quasi-BIC states of topological
charges $q=\pm 1$ and $q=2$. Since topological charges take only discrete
values, there must be a transition, not a smooth crossover between the
regimes. Our results inspire future studies of the transition boundaries: is
there multimode lasing, no lasing, bistability, fractional topological
charges? Furthermore, the edge states between areas of different structural
parameter may have properties distinct from the edge states of equilibrium
systems. Other possible outlooks include the study of splitting the doubly-
degenerate BICs at high-symmetry points of the Brillouin zone and the fate of
their topological charge when degeneracy is removed. In general, our results
introduce a versatile new platform for studies of non-Hermitian, out-of-
equilibrium topological physics Yu  _et al._ (2016); Gong  _et al._ (2018);
Bergholtz  _et al._ (2021).

###### Acknowledgements.

We thank Javier Cuerda for useful discussions on the finite element method
simulations. This work was supported by the Academy of Finland under project
number 318937 (PROFI), the Academy of Finland Flagship Programme, Photonics
Research and Innovation (PREIN), project numbers 320167 and by Centre for
Quantum Engineering (CQE) at Aalto University. Part of the research was
performed at the OtaNano Nanofab cleanroom (Micronova Nanofabrication Centre),
supported by Aalto University. We acknowledge the computational resources
provided by the Aalto Science-IT project. GS has received funding from the
European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement No 101025211 (TEBLA).

## References

  * Hsu  _et al._ (2016) C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, and M. Soljačić, Nature Review Materials 1, 1 (2016).
  * Hwang  _et al._ (2022) M.-S. Hwang, K.-Y. Jeong, J.-P. So, K.-H. Kim, and H.-G. Park, Communications Physics 5, 1 (2022).
  * Ozawa  _et al._ (2019) T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, Review of Modern Physics 91, 015006 (2019).
  * Rider  _et al._ (2019) M. S. Rider, S. J. Palmer, S. R. Pocock, X. Xiao, P. Arroyo Huidobro, and V. Giannini, Journal of Applied Physics 125, 120901 (2019).
  * Liu  _et al._ (2021) W. Liu, W. Liu, L. Shi, and Y. Kivshar, Nanophotonics 10, 1469 (2021).
  * Rider  _et al._ (2022) M. S. Rider, A. Buendía, D. R. Abujetas, P. A. Huidobro, J. A. Sánchez-Gil, and V. Giannini, ACS Photonics 9, 1483–1499 (2022).
  * Doeleman  _et al._ (2018) H. M. Doeleman, F. Monticone, W. den Hollander, A. Alù, and A. F. Koenderink, Nature Photonics 12, 397 (2018).
  * Zhen  _et al._ (2014) B. Zhen, C. W. Hsu, L. Lu, A. D. Stone, and M. Soljačić, Physical Review Letters 113, 257401 (2014).
  * Azzam  _et al._ (2021) S. I. Azzam, K. Chaudhuri, A. Lagutchev, Z. Jacob, Y. L. Kim, V. M. Shalaev, A. Boltasseva, and A. V. Kildishev, Laser & Photonics Reviews 15, 2000411 (2021).
  * Kodigala  _et al._ (2017) A. Kodigala, T. Lepetit, Q. Gu, B. Bahari, Y. Fainman, and B. Kanté, Nature 541, 196 (2017).
  * Ha  _et al._ (2018) S. T. Ha, Y. H. Fu, N. K. Emani, Z. Pan, R. M. Bakker, R. Paniagua-Domínguez, and A. I. Kuznetsov, Nature Nanotechnology 13, 1042 (2018).
  * Heilmann  _et al._ (2022) R. Heilmann, G. Salerno, J. Cuerda, T. K. Hakala, and P. Törmä, ACS Photonics 9, 224–232 (2022).
  * Asamoah  _et al._ (2022) B. O. Asamoah, M. Nečada, W. Liu, J. Heikkinen, S. Mohamed, A. Halder, H. T. Rekola, M. Koivurova, A. I. Väkeväinen, P. Törmä, J. Turunen, T. Setälä, A. T. Friberg, L. Shi, and T. K. Hakala, arXiv:2206.05011 10.48550/arXiv.2206.05011 (2022).
  * Iwahashi  _et al._ (2011) S. Iwahashi, Y. Kurosaka, K. Sakai, K. Kitamura, N. Takayama, and S. Noda, Optics Express 19, 11963 (2011).
  * Zito  _et al._ (2019) G. Zito, S. Romano, S. Cabrini, G. Calafiore, A. C. De Luca, E. Penzo, and V. Mocella, Optica 6, 1305 (2019).
  * Wang  _et al._ (2020) B. Wang, W. Liu, M. Zhao, J. Wang, Y. Zhang, A. Chen, F. Guan, X. Liu, L. Shi, and J. Zi, Nature Photonics 14, 623 (2020).
  * Huang  _et al._ (2020) C. Huang, C. Zhang, S. Xiao, Y. Wang, Y. Fan, Y. Liu, N. Zhang, G. Qu, H. Ji, J. Han, _et al._ , Science 367, 1018 (2020).
  * Dyakov  _et al._ (2021) S. A. Dyakov, M. V. Stepikhova, A. A. Bogdanov, A. V. Novikov, D. V. Yurasov, M. V. Shaleev, Z. F. Krasilnik, S. G. Tikhodeev, and N. A. Gippius, Laser & Photonics Reviews 15, 2000242 (2021).
  * Zhang  _et al._ (2021) Z. Zhang, Z. Lan, Y. Xie, M. Chen, W. Sha, and Y. Xu, Physical Review Applied 16, 064036 (2021).
  * Hwang  _et al._ (2021) M.-S. Hwang, H.-C. Lee, K.-H. Kim, K.-Y. Jeong, S.-H. Kwon, K. Koshelev, Y. Kivshar, and H.-G. Park, Nature Communications 12, 1 (2021).
  * Wu  _et al._ (2022) Y. Wu, L. Kang, and D. H. Werner, New Journal of Physics 24, 033002 (2022).
  * Kang  _et al._ (2022) L. Kang, Y. Wu, X. Ma, S. Lan, and D. H. Werner, Advanced Optical Materials 10, 2101497 (2022).
  * Gentry and Popović (2014) C. M. Gentry and M. A. Popović, Optics Letters 39, 4136 (2014).
  * Guan  _et al._ (2020a) J. Guan, L. K. Sagar, R. Li, D. Wang, G. Bappi, N. E. Watkins, M. R. Bourgeois, L. Levina, F. Fan, S. Hoogland, _et al._ , Nano Letters 20, 1468 (2020a).
  * Guan  _et al._ (2020b) J. Guan, L. K. Sagar, R. Li, D. Wang, G. Bappi, W. Wang, N. Watkins, M. R. Bourgeois, L. Levina, F. Fan, _et al._ , ACS Nano 14, 3426 (2020b).
  * Silva  _et al._ (2020) S. V. Silva, T. A. Morgado, and M. G. Silveirinha, Physical Review B 101, 041106(R) (2020).
  * De Angelis  _et al._ (2019) L. De Angelis, T. Bauer, F. Alpeggiani, and L. Kuipers, Optica 6, 1237 (2019).
  * Jin  _et al._ (2019) J. Jin, X. Yin, L. Ni, M. Soljačić, B. Zhen, and C. Peng, Nature 574, 501 (2019).
  * Mohamed  _et al._ (2022) S. Mohamed, J. Wang, H. Rekola, J. Heikkinen, B. Asamoah, L. Shi, and T. K. Hakala, Laser & Photonics Reviews , 2100574 (2022).
  * Bai  _et al._ (2021) T. Bai, Q. Li, Y. Wang, Y. Chen, Z.-D. Hu, and J. Wang, Optics Express 29, 25270 (2021).
  * Yoda and Notomi (2020) T. Yoda and M. Notomi, Physical Review Letters 125, 053902 (2020).
  * Zhang  _et al._ (2018) Y. Zhang, A. Chen, W. Liu, C. W. Hsu, B. Wang, F. Guan, X. Liu, L. Shi, L. Lu, and J. Zi, Physical Review Letters 120, 186103 (2018).
  * Moerland  _et al._ (2017) R. J. Moerland, T. K. Hakala, J.-P. Martikainen, H. T. Rekola, A. I. Väkeväinen, and P. Törmä,  _Quantum Plasmonics_ , edited by S. I. Bozhevolnyi, L. Martin-Moreno, and F. J. Garcia-Vidal (Springer, 2017) Chap. Strong Coupling Between Organic Molecules and Plasmonic Nanostructures.
  * Wang  _et al._ (2018a) W. Wang, M. Ramezani, A. I. Väkeväinen, P. Törmä, J. G. Rivas, and T. W. Odom, Materials Today 21, 303 (2018a).
  * Wang  _et al._ (2018b) D. Wang, W. Wang, M. P. Knudson, G. C. Schatz, and T. W. Odom, Chemical Reviews 118, 2865 (2018b).
  * Väkeväinen  _et al._ (2014) A. I. Väkeväinen, R. J. Moerland, H. T. Rekola, A.-P. Eskelinen, J.-P. Martikainen, D.-H. Kim, and P. Törmä, Nano Letters 14, 1721 (2014).
  * Hakala  _et al._ (2018) T. K. Hakala, A. J. Moilanen, A. I. Väkeväinen, R. Guo, J.-P. Martikainen, K. S. Daskalakis, H. T. Rekola, A. Julku, and P. Törmä, Nature Physics 14, 739 (2018).
  * Kravets  _et al._ (2018) V. G. Kravets, A. V. Kabashin, W. L. Barnes, and A. N. Grigorenko, Chemical Reviews 118, 5912 (2018).
  * Nečada and Törmä (2021) M. Nečada and P. Törmä, Communications in Computational Physics 30, 357 (2021).
  * Nec (2022) Marek Nečada, QPMS Photonic Multiple Scattering suite, (2020), source code available at https://repo.or.cz/qpms.git. (2022).
  * Sup (2022) See Supplemental Material at URLwillbeinsertedbypublisher. (2022).
  * Guo  _et al._ (2019) R. Guo, M. Nečada, T. K. Hakala, A. I. Väkeväinen, and P. Törmä, Physical Review Letters 122, 013901 (2019).
  * Wu  _et al._ (2020) M. Wu, S. T. Ha, S. Shendre, E. G. Durmusoglu, W.-K. Koh, D. R. Abujetas, J. A. Sánchez-Gil, R. Paniagua-Dominguez, H. V. Demir, and A. I. Kuznetsov, Nano Letters 20, 6005 (2020).
  * Yu  _et al._ (2016) K. Yu, Z. Guo, H. Jiang, and H. Chen, Journal of Applied Physics 119, 203102 (2016).
  * Gong  _et al._ (2018) Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Physical Review X 8, 031079 (2018).
  * Bergholtz  _et al._ (2021) E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Reviews of Modern Physics 93, 015005 (2021).

See Supplemental_Material.pdf


# Spike-based causal inference for weight alignment

###### Abstract

In artificial neural networks trained with gradient descent, the weights used
for processing stimuli are also used during backward passes to calculate
gradients. For the real brain to approximate gradients, gradient information
would have to be propagated separately, such that one set of synaptic weights
is used for processing and another set is used for backward passes. This
produces the so-called “weight transport problem” for biological models of
learning, where the backward weights used to calculate gradients need to
mirror the forward weights used to process stimuli. This weight transport
problem has been considered so hard that popular proposals for biological
learning assume that the backward weights are simply random, as in the
feedback alignment algorithm. However, such random weights do not appear to
work well for large networks. Here we show how the discontinuity introduced in
a spiking system can lead to a solution to this problem. The resulting
algorithm is a special case of an estimator used for causal inference in
econometrics, regression discontinuity design. We show empirically that this
algorithm rapidly makes the backward weights approximate the forward weights.
As the backward weights become correct, this improves learning performance
over feedback alignment on tasks such as Fashion-MNIST, SVHN, CIFAR-10 and
VOC. Our results demonstrate that a simple learning rule in a spiking network
can allow neurons to produce the right backward connections and thus solve the
weight transport problem.

Jordan Guerguiev1,2 Konrad P. Kording3 and Blake A. Richards4,5,6,7,*

1 Department of Biological Sciences, University of Toronto Scarborough,
Toronto, ON, Canada

2 Department of Cell and Systems Biology, University of Toronto, Toronto, ON,
Canada

3 Department of Bioengineering, University of Pennsylvania, PA, United States

4 Mila, Montreal, QC, Canada

5 Department of Neurology & Neurosurgery, McGill University, Montreal, QC,
Canada

6 School of Computer Science, McGill University, Montreal, QC, Canada

7 Canadian Institute for Advanced Research, Toronto, ON, Canada

* Corresponding author, email<EMAIL_ADDRESS>

##  1 Introduction

Any learning system that makes small changes to its parameters will only
improve if the changes are correlated to the gradient of the loss function.
Given that people and animals can also show clear behavioral improvements on
specific tasks (Shadmehr et al., 2010), however the brain determines its
synaptic updates, on average, the changes in must also correlate with the
gradients of some loss function related to the task (Raman et al., 2019). As
such, the brain may have some way of calculating at least an estimator of
gradients.

To-date, the bulk of models for how the brain may estimate gradients are
framed in terms of setting up a system where there are both bottom-up,
feedforward and top-down, feedback connections. The feedback connections are
used for propagating activity that can be used to estimate a gradient
(Williams, 1992; Lillicrap et al., 2016; Akrout et al., 2019; Roelfsema &
Ooyen, 2005; Lee et al., 2015; Scellier & Bengio, 2017; Sacramento et al.,
2018). In all such models, the gradient estimator is less biased the more the
feedback connections mirror the feedforward weights. For example, in the
REINFORCE algorithm (Williams, 1992), and related algorithms like AGREL
(Roelfsema & Ooyen, 2005), learning is optimal when the feedforward and
feedback connections are perfectly symmetric, such that for any two neurons
$i$ and $j$ the synaptic weight from $i$ to $j$ equals the weight from $j$ to
$i$, e.g. $W_{ji}=W_{ij}$ (Figure 1). Some algorithms simply assume weight
symmetry, such as Equilibrium Propagation (Scellier & Bengio, 2017). The
requirement for synaptic weight symmetry is sometimes referred to as the
“weight transport problem”, since it seems to mandate that the values of the
feedforward synaptic weights are somehow transported into the feedback
weights, which is not biologically realistic (Crick, 1989-01-12; Grossberg,
1987). Solving the weight transport problem is crucial to biologically
realistic gradient estimation algorithms (Lillicrap et al., 2016), and is thus
an important topic of study.

Several solutions to the weight transport problem have been proposed for
biological models, including hard-wired sign symmetry (Moskovitz et al.,
2018), random fixed feedback weights (Lillicrap et al., 2016), and learning to
make the feedback weights symmetric (Lee et al., 2015; Sacramento et al.,
2018; Akrout et al., 2019; Kolen & Pollack, 1994). Learning to make the
weights symmetric is promising because it is both more biologically feasible
than hard-wired sign symmetry (Moskovitz et al., 2018) and it leads to less
bias in the gradient estimator (and thereby, better training results) than
using fixed random feedback weights (Bartunov et al., 2018; Akrout et al.,
2019). However, of the current proposals for learning weight symmetry some do
not actually work well in practice (Bartunov et al., 2018) and others still
rely on some biologically unrealistic assumptions, including scalar value
activation functions (as opposed to all-or-none spikes) and separate error
feedback pathways with one-to-one matching between processing neurons for the
forward pass and error propagation neurons for the backward pass Akrout et al.
(2019); Sacramento et al. (2018).

Interestingly, learning weight symmetry is implicitly a causal inference
problem—the feedback weights need to represent the causal influence of the
upstream neuron on its downstream partners. As such, we may look to the causal
infererence literature to develop better, more biologically realistic
algorithms for learning weight symmetry. In econometrics, which focuses on
quasi-experiments, researchers have developed various means of estimating
causality without the need to actually randomize and control the variables in
question Angrist & Pischke (2008); Marinescu et al. (2018). Among such quasi-
experimental methods, regression discontinuity design (RDD) is particularly
promising. It uses the discontinuity introduced by a threshold to estimate
causal effects. For example, RDD can be used to estimate the causal impact of
getting into a particular school (which is a discontinuous, all-or-none
variable) on later earning power. RDD is also potentially promising for
estimating causal impact in biological neural networks, because real neurons
communicate with discontinuous, all-or-none spikes. Indeed, it has been shown
that the RDD approach can produce unbiased estimators of causal effects in a
system of spiking neurons Lansdell & Kording (2019). Given that learning
weight symmetry is fundamentally a causal estimation problem, we hypothesized
that RDD could be used to solve the weight transport problem in biologically
realistic, spiking neural networks.

Here, we present a learning rule for feedback synaptic weights that is a
special case of the RDD algorithm previously developed for spiking neural
networks (Lansdell & Kording, 2019). Our algorithm takes advantage of a
neuron’s spiking discontinuity to infer the causal effect of its spiking on
the activity of downstream neurons. Since this causal effect is proportional
to the feedforward synaptic weight between the two neurons, by estimating it,
feedback synapses can align their weights to be symmetric with the reciprocal
feedforward weights, thereby overcoming the weight transport problem. We
demonstrate that this leads to the reduction of a cost function which measures
the weight symmetry (or the lack thereof), that it can lead to better weight
symmetry in spiking neural networks than other algorithms for weight alignment
(Akrout et al., 2019) and it leads to better learning in deep neural networks
in comparison to the use of fixed feedback weights (Lillicrap et al., 2016).
Altogether, these results demonstrate a novel algorithm for solving the weight
transport problem that takes advantage of discontinuous spiking, and which
could be used in future models of biologically plausible gradient estimation.

Figure 1: Illustration of weight symmetry in a neural network with feedforward
and feedback connections. Processing of inputs to outputs is mirrored by
backward flow of gradient information. Gradient estimation is best when
feedback synapses have symmetric weights to feedforward synapses (illustrated
with colored circles).

##  2 Related work

Previous work has shown that even when feedback weights in a neural network
are initialized randomly and remain fixed throughout training, the feedforward
weights learn to partially align themselves to the feedback weights, an
algorithm known as feedback alignment (Lillicrap et al., 2016). While feedback
alignment is successful at matching the learning performance of true gradient
descent in relatively shallow networks, it does not scale well to deeper
networks and performs poorly on difficult computer vision tasks (Bartunov et
al., 2018).

The gap in learning performance between feedback alignment and gradient
descent can be overcome if feedback weights are continually updated to match
the sign of the reciprocal feedforward weights (Moskovitz et al., 2018).
Furthermore, learning the feedback weights in order to make them more
symmetric to the feedforward weights has been shown to improve learning over
feedback alignment (Akrout et al., 2019).

To understand the underlying dynamics of learning weight symmetry, Kunin et
al. (2019) define the symmetric alignment cost function,
$\mathcal{R}_{\text{SA}}$, as one possible cost function that, when minimized,
leads to weight symmetry:

$\displaystyle\mathcal{R}_{\text{SA}}$
$\displaystyle\coloneqq\|W-Y^{T}\|^{2}_{F}$ (1)
$\displaystyle=\|W\|^{2}_{F}+\|Y\|^{2}_{F}-2\text{tr}(WY)$

where $W$ are feedforward weights and $Y$ are feedback weights. The first two
terms are simply weight regularization terms that can be minimized using
techniques like weight decay. But, the third term is the critical one for
ensuring weight alignment.

In this paper we present a biologically plausible method of minimizing the
third term. This method is based on the work of Lansdell & Kording (2019), who
demonstrated that neurons can estimate their causal effect on a global reward
signal using the discontinuity introduced by spiking. This is accomplished
using RDD, wherein a piecewise linear model is fit around a discontinuity, and
the differences in the regression intercepts indicates the causal impact of
the discontinuous variable. In Lansdell & Kording (2019), neurons learn a
piece-wise linear model of a reward signal as a function of their input drive,
and estimate the causal effect of spiking by looking at the discontinuity at
the spike threshold. Here, we modify this technique to perform causal
inference on the effect of spiking on downstream neurons, rather than a reward
signal. We leverage this to develop a learning rule for feedback weights that
induces weight symmetry and improves training.

##  3 Our contributions

The primary contributions of this paper are as follows:

  * •

We demonstrate that spiking neurons can accurately estimate the causal effect
of their spiking on downstream neurons by using a piece-wise linear model of
the feedback as a function of the input drive to the neuron.

  * •

We present a learning rule for feedback weights that uses the causal effect
estimator to encourage weight symmetry. We show that when feedback weights
update using this algorithm it minimizes the symmetric alignment cost
function, $\mathcal{R}_{\text{SA}}$.

  * •

We demonstrate that this learning weight symmetry rule improves training and
test accuracy over feedback alignment, approaching gradient-descent-level
performance on Fashion-MNIST, SVHN, CIFAR-10 and VOC in deeper networks.

##  4 Methods

Figure 2: A. Layers of the convolutional network trained on CIFAR-10 and the
corresponding network of LIF neurons that undergoes RDD feedback training.
Fully-connected feedforward weights (blue) and feedback weights (red) are
shared between the two networks. Every training epoch consists of an RDD
feedback training phase where feedback weights in the LIF net are updated
(bold red arrows) and transferred to the convolutional net, and a regular
training phase where feedforward weights in the convolutional net are updated
(bold blue arrows) and transferred back to the LIF net. B. RDD feedback
training protocol. Every 30 s, a different layer in the LIF network receives
driving input and updates its feedback weights (red) using the RDD algorithm.
C. Top: Sample voltage ($v_{i}$, solid line) and input drive ($u_{i}$, dashed
line) traces. Whenever $v_{i}$ approaches the spiking threshold, an RDD window
lasting $T$ ms is triggered. $u_{i}^{\text{max}}$ is the maximum input drive
during this window of time. Bottom: Feedback received at a synapse, $q_{k}$.
$q_{k}^{\text{pre}}$ is the feedback signal at the start of an RDD window,
while $q_{k}^{\text{avg}}$ is the average of the feedback signal during the
time window. D. Samples of $\Delta q_{k}^{\text{avg}}$ vs.
$u_{i}^{\text{max}}$ are used to update a piece-wise linear function of
$u_{i}^{\text{max}}$, and the causal effect $\beta_{ik}$ is defined as the
difference of the left and right limits of the function at the spiking
threshold.

###  4.1 General approach

In this work, we utilize a spiking neural network model for aligning
feedforward and feedback weights. However, due to the intense computational
demands of spiking neural networks, we only use spikes for the RDD algorithm.
We then use the feedback weights learned by the RDD algorithm for training a
non-spiking convolutional neural network. We do this because the goal of our
work here is to develop an algorithm for aligning feedback weights in spiking
networks, not for training feedforward weights in spiking networks on other
tasks. Hence, in the interest of computational expediency, we only used
spiking neurons when learning to align the weights. Additional details on this
procedure are given below.

###  4.2 RDD feedback training phase

At the start of every training epoch of a convolutional neural network, we use
an RDD feedback weight training phase, during which all fully-connected sets
of feedback weights in the network are updated. To perform these updates, we
simulate a separate network of leaky integrate-and-fire (LIF) neurons. LIF
neurons incorporate key elements of real neurons such as voltages, spiking
thresholds and refractory periods. Each epoch, we begin by training the
feedback weights in the LIF network. These weights are then transferred to the
convolutional network, which is used for training the feedforward weights. The
new feedforward weights are then transferred to the LIF net, and another
feedback training phase with the LIF net starts the next epoch (Figure 2A).
During the feedback training phase, the LIF network undergoes a training phase
lasting 90 s of simulated time (30 s per set of feedback weights) (Figure 2B).
We find that the spiking network used for RDD feedback training and the
convolutional neural network are very closely matched in the activity of the
units (Figure S1), which gives us confidence that this approach of using a
separate non-spiking network for training the feedforward weights is
legitimate.

During the feedback training phase, a small subset of neurons in the first
layer receive driving input that causes them to spike, while other neurons in
this layer receive no input (see Appendix A.2). The subset of neurons that
receive driving input is randomly selected every 100 ms of simulated time.
This continues for 30 s in simulated time, after which the same process occurs
for the subsequent hidden layers in the network. This protocol enforces
sparse, de-correlated firing patterns that improve the causal inference
procedure of RDD.

###  4.3 LIF dynamics

During the RDD feedback training phase, each unit in the network is simulated
as a leaky integrate-and-fire neuron. Spiking inputs from the previous layer
arrive at feedforward synapses, where they are convolved with a temporal
exponential kernel to simulate post-synaptic spike responses
$\mathbf{p}=[p_{1},p_{2},...,p_{m}]$ (see Appendix A.1). The neurons can also
receive driving input $\tilde{p}_{i}$, instead of synaptic inputs. The total
feedforward input to neuron $i$ is thus defined as:

$I_{i}:=\begin{cases}\omega\tilde{p}_{i}&\text{if $\tilde{p}_{i}>0$}\\\
\sum_{j=1}^{m}W_{ij}p_{j}&\text{otherwise}\end{cases}$ (2)

where $W_{ij}$ is the feedforward weight to neuron $i$ from neuron $j$ in the
previous layer, and $\omega$ is a hyperparameter. The voltage of the neuron,
$v_{i}$, evolves as:

$\frac{dv_{i}}{dt}=-g_{L}v_{i}+g_{D}(I_{i}-v_{i})$ (3)

where $g_{L}$ and $g_{D}$ are leak and dendritic conductance constants,
respectively. The input drive to the neuron, $u_{i}$, is similarly modeled:

$\frac{du_{i}}{dt}=-g_{L}u_{i}+g_{D}(I_{i}-u_{i})$ (4)

If the voltage $v_{i}$ passes a spiking threshold $\theta$, the neuron spikes
and the voltage is reset to a value $v_{\text{reset}}=-1$ (Figure 2C). Note
that the input drive does not reset. This helps us to perform regressions both
above and below the spike threshold.

In addition to feedforward inputs, spiking inputs from the downstream layer
arrive at feedback synapses, where they create post-synaptic spike responses
$\mathbf{q}=[q_{1},q_{2},...,q_{n}]$. These responses are used in the causal
effect estimation (see below).

###  4.4 RDD algorithm

Whenever the voltage approaches the threshold $\theta$ (ie.
$|v_{i}-\theta|<\alpha$ where $\alpha$ is a constant), an RDD window is
initiated, lasting $T=30$ ms in simulated time (Figure 2C). At the end of this
time window, at each feedback synapse, the maximum input drive during the RDD
window, $u_{i}^{\text{max}}$, and the average change in feedback from
downstream neuron $k$ during the RDD window, $\Delta q_{k}^{\text{avg}}$, are
recorded. $\Delta q_{k}^{\text{avg}}$ is defined as the difference between the
average feedback received during the RDD window, $q_{k}^{\text{avg}}$, and the
feedback at the start of the RDD window, $q_{k}^{\text{pre}}$:

$\Delta q_{k}^{\text{avg}}\coloneqq q_{k}^{\text{avg}}-q_{k}^{\text{pre}}$ (5)

Importantly, $u_{i}^{\text{max}}$ provides a measure of how strongly neuron
$i$ was driven by its inputs (and whether or not it passed the spiking
threshold $\theta$), while $\Delta q_{k}^{\text{avg}}$ is a measure of how the
input received as feedback from neuron $k$ changed after neuron $i$ was driven
close to its spiking threshold. These two values are then used to fit a piece-
wise linear model of $\Delta q_{k}^{\text{avg}}$ as a function of
$u_{i}^{\text{max}}$ (Figure 2D). This piece-wise linear model is defined as:

$f_{ik}(x)\coloneqq\begin{cases}c^{1}_{ik}x+c^{2}_{ik}&\text{if $x<\theta$
}\\\ c^{3}_{ik}x+c^{4}_{ik}&\text{if $x\geq\theta$ }\end{cases}$ (6)

The parameters $c^{1}_{ik}$, $c^{2}_{ik}$, $c^{3}_{ik}$ and $c^{4}_{ik}$ are
updated to perform linear regression using gradient descent:

$\displaystyle L$
$\displaystyle=\frac{1}{2}\|f_{ik}(u_{i}^{\text{max}})-\Delta
q_{k}^{\text{avg}}\|^{2}$ (7) $\displaystyle\Delta c^{l}_{ik}$
$\displaystyle\propto-\frac{\partial L}{\partial c^{l}_{ik}}\quad\text{for
$l\in\\{1,2,3,4\\}$}$ (8)

An estimate of the causal effect of neuron $i$ spiking on the activity of
neuron $k$, $\beta_{ik}$, is then defined as the difference in the two sides
of the piece-wise linear function at the spiking threshold:

$\beta_{ik}\coloneqq\lim_{x\rightarrow\theta^{+}}f_{ik}(x)-\lim_{x\rightarrow\theta^{-}}f_{ik}(x)$
(9)

Finally, the weight at the feedback synapse, $Y_{ik}$, is updated to be a
scaled version of $\beta_{ik}$:

$Y_{ik}=\beta_{ik}\frac{\gamma}{\sigma^{2}_{\beta}}$ (10)

where $\gamma$ is a hyperparameter and ${\sigma^{2}_{\beta}}$ is the standard
deviation of $\beta$ values for all feedback synapses in the layer. This
ensures that the scale of the full set of feedback weights between two layers
in the network remains stable during training.

##  5 Results

Figure 3: A. Evolution of sign alignment (the percent of feedforward and
feedback weights that have the same sign) for each fully-connected layer in
the network when trained on CIFAR-10 using RDD feedback training (blue), using
the algorithm proposed by Akrout et al. (2019) (purple), and using feedback
alignment (red). B. Feedforward vs. feedback weights for each fully-connected
layer at the end of training, with RDD feedback training (blue), the Akrout
algorithm (purple), and feedback alignment (red).

###  5.1 Alignment of feedback and feedforward weights

To measure how well the causal effect estimate at each feedback synapse,
$\beta_{ik}$, and thus the feedback weight $Y_{ik}$, reflects the reciprocal
feedforward weight $W_{ki}$, we can measure the percentage of feedback weights
that have the same sign as the reciprocal feedforward weights (Figure 3A).
When training on CIFAR-10 with no RDD feedback training phase (ie. feedback
weights remain fixed throughout training), the feedback alignment effect
somewhat increases the sign alignment during training, but it is ineffective
at aligning the signs of weights in earlier layers in the network. Compared to
feedback alignment, the addition of an RDD feedback training phase greatly
increases the sign aligmnent between feedback and feedforward weights for all
layers in the network, especially at earlier layers. In addition, the RDD
algorithm increases sign alignment throughout the hierarchy more than the
current state-of-the-art algorithm for weight alignment introduced recently by
Akrout et al. Akrout et al. (2019) (Figure 3A). Furthermore, RDD feedback
training changes feedback weights to not only match the sign but also the
magnitude of the reciprocal feedforward weights (Figure 3B), which makes it
better for weight alignment than hard-wired sign symmetry (Moskovitz et al.,
2018).

###  5.2 Descending the symmetric alignment cost function

The symmetric alignment cost function (Kunin et al., 2019) (Equation 1) can be
broken down as:

$\displaystyle\mathcal{R}_{\text{SA}}$
$\displaystyle=\mathcal{R}_{\text{decay}}+\mathcal{R}_{\text{self}}$ (11)

where we define $\mathcal{R}_{\text{decay}}$ and $\mathcal{R}_{\text{self}}$
as:

$\displaystyle\mathcal{R}_{\text{decay}}$
$\displaystyle\coloneqq\|W\|^{2}_{F}+\|Y\|^{2}_{F}$ (12)
$\displaystyle\mathcal{R}_{\text{self}}$
$\displaystyle\coloneqq-2\text{tr}(WY)$ (13)

$\mathcal{R}_{\text{decay}}$ is simply a weight regularization term that can
be minimized using techniques like weight decay. $\mathcal{R}_{\text{self}}$,
in contrast, measures how well aligned in direction the two matrices are. Our
learning rule for feedback weights minimizes the $\mathcal{R}_{\text{self}}$
term for weights throughout the network (Figure 4). By comparison, feedback
alignment decreases $\mathcal{R}_{\text{self}}$ to a smaller extent, and its
ability to do so diminishes at earlier layers in the network. This helps to
explain why our algorithm induces weight alignment, and can improve training
performance (see below).

Figure 4: Evolution of $\mathcal{R}_{\text{self}}$ for each fully-connected
layer in the network when trained on CIFAR-10 using RDD feedback training
(solid lines) and using feedback alignment (dashed lines). RDD feedback
training dramatically decreases this loss compared to feedback alignment,
especially in earlier layers.

###  5.3 Performance on Fashion-MNIST, SVHN, CIFAR-10 and VOC

We trained the same network architecture (see Appendix A.3) on the Fashion-
MNIST, SVHN, CIFAR-10 and VOC datasets using standard autograd techniques
(backprop), feedback alignment and our RDD feedback training phase. RDD
feedback training substantially improved the network’s performance over
feedback alignment, and led to backprop-level accuracy on the train and test
sets (Figure 5).

Figure 5: Comparison of Fashion-MNIST, SVHN, CIFAR-10 and VOC train error (top
row) and test error (bottom row). RDD feedback training substantially improves
test error performance over feedback alignment in both learning tasks.

##  6 Discussion

In order to understand how the brain learns complex tasks that require
coordinated plasticity across many layers of synaptic connections, it is
important to consider the weight transport problem. Here, we presented an
algorithm for updating feedback weights in a network of spiking neurons that
takes advantage of the spiking discontinuity to estimate the causal effect
between two neurons (Figure 2). We showed that this algorithm enforces weight
alignment (Figure 3), and identified a loss function,
$\mathcal{R}_{\text{self}}$, that is minimized by our algorithm (Figure 4).
Finally, we demonstrated that our algorithm allows deep neural networks to
achieve better learning performance than feedback alignment on Fashion-MNIST
and CIFAR-10 (Figure 5). These results demonstrate the potential power of RDD
as a means for solving the weight transport problem in biologically plausible
deep learning models.

One aspect of our algorithm that is still biologically implausible is that it
does not adhere to Dale’s principle, which states that a neuron performs the
same action on all of its target cells (Strata & Harvey, ). This means that a
neuron’s outgoing connections cannot include both positive and negative
weights. However, even under this constraint, a neuron can have an excitatory
effect on one downstream target and an inhibitory effect on another, by
activating intermediary inhibitory interneurons. Because our algorithm
provides a causal estimate of one neuron’s impact on another, theoretically,
it could capture such polysynaptic effects. Therefore, this algorithm is in
theory compatible with Dale’s principle. Future work should test the effects
of this algorithm when implemented in a network of neurons that are explicitly
excitatory or inhibitory.

## References

  * Akrout et al. (2019) Mohamed Akrout, Collin Wilson, Peter C Humphreys, Timothy Lillicrap, and Douglas Tweed.  Using weight mirrors to improve feedback alignment.  _arXiv preprint arXiv:1904.05391_ , 2019.
  * Angrist & Pischke (2008) Joshua D Angrist and Jörn-Steffen Pischke.  _Mostly harmless econometrics: An empiricist’s companion_.  Princeton university press, 2008.
  * Bartunov et al. (2018) Sergey Bartunov, Adam Santoro, Blake Richards, Luke Marris, Geoffrey E Hinton, and Timothy Lillicrap.  Assessing the scalability of biologically-motivated deep learning algorithms and architectures.  In _Advances in Neural Information Processing Systems_ , pp. 9368–9378, 2018.
  * Crick (1989-01-12) Francis Crick.  The recent excitement about neural networks.  337(6203):129–132, 1989-01-12.  doi: 10.1038/337129a0.  URL http://dx.doi.org/10.1038/337129a0.
  * Grossberg (1987) Stephen Grossberg.  Competitive learning: From interactive activation to adaptive resonance.  11(1):23–63, 1987.  ISSN 1551-6709.
  * Kolen & Pollack (1994) John F Kolen and Jordan B Pollack.  Backpropagation without weight transport.  In _Proceedings of 1994 IEEE International Conference on Neural Networks (ICNN’94)_ , volume 3, pp. 1375–1380. IEEE, 1994.
  * Kunin et al. (2019) Daniel Kunin, Jonathan M Bloom, Aleksandrina Goeva, and Cotton Seed.  Loss landscapes of regularized linear autoencoders.  _arXiv preprint arXiv:1901.08168_ , 2019.
  * Lansdell & Kording (2019) Benjamin James Lansdell and Konrad Paul Kording.  Spiking allows neurons to estimate their causal effect.  _bioRxiv_ , pp. 253351, 2019.
  * Lee et al. (2015) Dong-Hyun Lee, Saizheng Zhang, Asja Fischer, and Yoshua Bengio.  Difference target propagation.  In _Joint European conference on machine learning and knowledge discovery in databases_ , pp. 498–515. Springer, 2015.
  * Lillicrap et al. (2016) Timothy P Lillicrap, Daniel Cownden, Douglas B Tweed, and Colin J Akerman.  Random synaptic feedback weights support error backpropagation for deep learning.  _Nature communications_ , 7:13276, 2016.
  * Marinescu et al. (2018) Ioana E Marinescu, Patrick N Lawlor, and Konrad P Kording.  Quasi-experimental causality in neuroscience and behavioural research.  _Nature human behaviour_ , pp. 1, 2018.
  * Moskovitz et al. (2018) Theodore H Moskovitz, Ashok Litwin-Kumar, and LF Abbott.  Feedback alignment in deep convolutional networks.  _arXiv preprint arXiv:1812.06488_ , 2018.
  * Raman et al. (2019) Dhruva Venkita Raman, Adriana Perez Rotondo, and Timothy O’Leary.  Fundamental bounds on learning performance in neural circuits.  116(21):10537–10546, 2019.  ISSN 0027-8424.
  * Roelfsema & Ooyen (2005) Pieter R Roelfsema and Arjen van Ooyen.  Attention-gated reinforcement learning of internal representations for classification.  _Neural computation_ , 17(10):2176–2214, 2005\.
  * Sacramento et al. (2018) João Sacramento, Rui Ponte Costa, Yoshua Bengio, and Walter Senn.  Dendritic cortical microcircuits approximate the backpropagation algorithm.  pp. 8735–8746, 2018.
  * Scellier & Bengio (2017) Benjamin Scellier and Yoshua Bengio.  Equilibrium propagation: Bridging the gap between energy-based models and backpropagation.  11:24, 2017.  ISSN 1662-5188.  doi: 10.3389/fncom.2017.00024.  URL http://www.ncbi.nlm.nih.gov/pmc/articles/PMC5415673/.
  * Shadmehr et al. (2010) Reza Shadmehr, Maurice A Smith, and John W Krakauer.  Error correction, sensory prediction, and adaptation in motor control.  33:89–108, 2010.  ISSN 0147-006X.
  * (18) Piergiorgio Strata and Robin Harvey.  Dale’s principle.  50(5):349–350.  ISSN 0361-9230.  doi: 10.1016/S0361-9230(99)00100-8.  URL http://www.sciencedirect.com/science/article/pii/S0361923099001008.
  * Williams (1992) Ronald J Williams.  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  _Machine learning_ , 8(3-4):229–256, 1992.

##  Appendix A Appendix

###  A.1 LIF neuron simulation details

Post-synaptic spike responses at feedforward synapses, $\mathbf{p}$, were
calculated from pre-synaptic binary spikes using an exponential kernel
function $\kappa$:

$p_{j}(t)=\sum_{k}\kappa(t-\tilde{t}_{jk})$ (14)

where $\tilde{t}_{jk}$ is the $k^{\text{th}}$ spike time of input neuron $j$
and $\kappa$ is given by:

$\kappa(t)=(e^{-t/\tau_{L}}-e^{-t/\tau_{S}})\Theta(t)/(\tau_{L}-\tau_{s})$
(15)

where $\tau_{s}=0.003$ s and $\tau_{L}=0.01$ s represent short and long time
constants, and $\Theta$ is the Heaviside step function. Post-synaptic spike
responses at feedback synapses, $\mathbf{q}$, were computed in the same way.

###  A.2 RDD feedback training implementation

####  A.2.1 Weight scaling

Weights were shared between the convolutional network and the network of LIF
neurons, but feedforward weights in the LIF network were scaled versions of
the convolutional network weights:

$W_{ij}^{\text{LIF}}=\psi mW_{ij}^{\text{Conv}}/\sigma_{W^{\text{Conv}}}^{2}$
(16)

where $W^{\text{Conv}}$ is a feedforward weight matrix in the convolutional
network, $W^{\text{LIF}}$ is the corresponding weight matrix in the LIF
network, $m$ is the number of units in the upstream layer (ie. the number of
columns in $W^{\text{Conv}}$), $\sigma_{W^{\text{Conv}}}^{2}$ is the standard
deviation of $W^{\text{Conv}}$ and $\psi$ is a hyperparameter. This rescaling
ensures that spike rates in the LIF network stay within an optimal range for
the RDD algorithm to converge quickly, even if the scale of the feedforward
weights in the convolutional network changes during training. This avoids
situations where the scale of feedforward weights is so small that little or
no spiking occurs in the LIF neurons.

####  A.2.2 feedback training paradigm

The RDD feedback training paradigm is implemented as follows. We start by
providing driving input to the first layer in the network of LIF neurons. To
create this driving input, we choose a subset of $20\%$ of the neurons in that
layer, and create a unique input spike train for each of these neurons using a
Poisson process with a rate of $200$ Hz. All other neurons in the layer
receive no driving input. Every $100$ ms, a new set of neurons to receive
driving input is randomly chosen. After $30$ s, this layer stops receiving
driving input, and the process repeats for the next layer in the network.

###  A.3 Network and training details

The network architectures used to train on Fashion-MNIST and CIFAR-10 are
described in Table 1.

Layer | Fashion-MNIST | SVHN & CIFAR-10 | VOC
---|---|---|---
Input | $28\times 28\times 1$ | $32\times 32\times 3$ | $32\times 32\times 3$
1 | Conv2D $5\times 5$, 64 ReLU | Conv2D $5\times 5$, 64 ReLU | Conv2D $5\times 5$, 64 ReLU
2 | MaxPool $2\times 2$, stride 2 | MaxPool $2\times 2$, stride 2 | MaxPool $2\times 2$, stride 2
3 | Conv2D $5\times 5$, 64 ReLU | Conv2D $5\times 5$, 64 ReLU | Conv2D $5\times 5$, 64 ReLU
4 | MaxPool $2\times 2$, stride 2 | MaxPool $2\times 2$, stride 2 | MaxPool $2\times 2$, stride 2
5 | FC 384 ReLU | FC 384 ReLU | FC 384 ReLU
6 | FC 192 ReLU | FC 192 ReLU | FC 192 ReLU
7 | FC 10 ReLU | FC 10 ReLU | FC 21 ReLU
Table 1: Network architectures used to train on Fashion-MNIST, SVHN, CIFAR-10
and VOC.

Inputs were randomly cropped and flipped during training, and batch
normalization was used at each layer. Networks were trained using a minibatch
size of 32.

###  A.4 Akrout et al. (2019) algorithm implementation

In experiments that compared sign alignment using our RDD algorithm with the
Akrout et al. (2019) algorithm, we kept the same RDD feedback training
paradigm (ie. layers were sequentially driven, and a small subset of neurons
in each layer was active at once). However, rather than updating feedback
weights using RDD, we recorded the mean firing rates of the active neurons in
the upstream layer, $\mathbf{r}^{l}$, and the mean firing rates in the
downstream layer, $\mathbf{r}^{l+1}$. We then used the following feedback
weight update rule:

$\Delta\mathbf{Y}=\eta\mathbf{r}^{l}\mathbf{r}^{{(l+1)}^{T}}-\lambda_{\text{WD}}\mathbf{Y}$
(17)

where $Y$ are the feedback weights between layers $l+1$ and $l$, and $\eta$
and $\lambda_{\text{WD}}$ are learning rate and weight decay hyperparameters,
respectively.

Figure S1: Comparison of average spike rates in the fully-connected layers of
the LIF network vs. activities of the same layers in the convolutional
network, when both sets of layers were fed the same input. Spike rates in the
LIF network are largely correlated with activities of units in the
convolutional network.


Conservation laws for a mathematical model of HIV transmission

Winter Sinkala , Andrew Otieno111Corresponding author. Tel.: +27 47 502 2271;
Fax: +27 47 502 2725.
E-mail addresses<EMAIL_ADDRESS><EMAIL_ADDRESS>
Department of Mathematical Sciences and Computing, Walter Sisulu University,
Private Bag X1, Mthatha 5117, South Africa.

###### Abstract

A theorem due to Nail H. Ibragimov (2007) provides a connection between
symmetries and conservation laws for arbitrary differential equations. The
theorem is valid for any system of differential equations provided that the
number of equations is equal to the number of dependent variables. In this
paper we use the theorem to determine conservation laws for a nonlinear system
of ordinary differential equations that represents a mathematical model for
HIV transmission.

KEY WORDS: Lie symmetry analysis, Adjoint equation, Noether symmetry,
Nonlinear equations, Lagrangian, Conservation laws.

##  1 Introduction

Ibragimov [9] provides a general theorem on conservation laws by means of
which conservation laws can be constructed for an arbitrary system of
differential equations admitting Lie symmetries, provided the number of
equations is equal to the number of dependent variables. Ibragimov’s theorem
is an extension of the classical theorem of Noether [3] in that the latter
theorem does not require existence of a Lagrangian. In fact, unlike Noether’s
theorem, Ibragimov’s theorem allows one to associate a conservation law to
every Lie point symmetry admitted by a given arbitrary system of differential
equations. Applications that one comes across involving the use of admitted
Lie point symmetries to construct conservation laws via either the classical
Noether’s theorem or Ibragimov’s theorem are often limited to scalar partial
differential equations [13, 2, 8, 7, 12, 14, 6, 1, 4]. The authors were unable
to find any applications involving systems of first-order ordinary
differential equations, for instance. In this connection the application
reported in this paper adds to the repertoire of nontrivial applications of
Ibragimov’s theorem.

We consider a nonlinear system of three first-order ordinary differential
equations that arise from a model formulated by Anderson [11] which describes
the transmission of HIV/AIDS in male homosexual/bisexual cohorts. A particular
case of the model translates into a coupled nonlinear system of first-order
ordinary differential equations [15]:

$\begin{array}[]{l}\displaystyle F_{1}(t,u,u_{(1)})=u^{1}_{t}+\frac{\beta
cu^{1}u^{2}}{u^{1}+u^{2}+u^{3}}+\mu u^{1}=0\\\ \displaystyle
F_{2}(t,u,u_{(1)})=u^{2}_{t}-\frac{\beta
cu^{1}u^{2}}{u^{1}+u^{2}+u^{3}}+(\nu+\mu)u^{2}=0\\\ \displaystyle
F_{3}(t,u,u_{(1)})=u^{3}_{t}-\nu u^{2}+(\mu+\beta c)u^{3}=0,\end{array}$ (1.1)

where $u_{(1)}$ denotes the first-order partial derivatives $u^{1}_{t}$,
$u^{2}_{t}$ and $u^{3}_{t}$. This system is typical of models that are
formulated to mimic dynamics in epidemiology. Such systems are often highly
nonlinear and difficult to analyze. In their seminal work, Torrisi and Nucci
[15] perform Lie symmetry analysis on (1.1). They determine a solvable Lie
algebra admitted by the system and exploit this to find a solution (albeit in
quadrature form) of the system. In this article we extend the Lie symmetry
analysis of (1.1) by generating conservation laws for the system via
Ibragimov’s new conservation theorem [9].

Conservation laws of physical systems are fundamental to our understanding of
the system being studied. Apart from having a direct physical interpretation,
conservation laws may be essential in studying the integrability of the
system. For example, in the numerical integration of partial differential
equations conservation laws help to control numerical errors in that they
describe the properties of the system that do not change. On the whole
conservation laws play an important role in the analysis of basic properties
of the solutions. Therefore, the construction of conservation laws is one of
the most important applications of symmetries to physical systems [10, 5].

The rest of this paper is organised as follows. We present elements of Lie
symmetry analysis of differential equations in Section 2. An overview of the
theorems of Noether and Ibragimov for constructing conservation laws via
admitted symmetries is provided in Section 3. Conservation Laws of Anderson’s
HIV model are constructed in Section 4. In Section 4 we discuss the results
and give concluding remarks.

##  2 Preliminaries

Let us consider an $r$th-order ($r\geq 1$) system of differential equations
with $m$ dependent variables $u=(u^{1},u^{2},\ldots,u^{m})$ and $n$
independent variable $x=(x^{1},x^{2},\ldots,x^{n})$, $u=u(x)$,

$F_{\alpha}(x,u,\ldots,u_{(r)})=0,\quad\alpha=1,\ldots,m,$ (2.1)

where $u_{(k)}$ denotes the collection $\left\\{u_{k}^{\alpha}\right\\}$ of
$k$th-order derivatives, $k\geq 1$. Suppose that (2.1) admits a one parameter
Lie group of point transformations

$\begin{array}[]{c}\widetilde{x}^{i}=x^{i}+\mbox{$\varepsilon$}\,\xi^{i}(x,u)+O(\mbox{$\varepsilon$}^{2}),\\\
\widetilde{u}^{\alpha}=u^{\alpha}+\mbox{$\varepsilon$}\,\eta^{\alpha}(x,u)+O(\mbox{$\varepsilon$}^{2}),\\\
\end{array}$ (2.2)

where $\varepsilon$ is a real parameter; $\xi^{i}$ and $\eta^{\alpha}$ are
given smooth functions. Invariance of (2.1) under (2.2) is conveniently
expressed in terms of the infinitesimal generator of (2.2) is

$X=\xi^{i}(x,u)\dfrac{\partial}{\partial{x^{i}}}+\eta^{\alpha}(x,u)\dfrac{\partial}{\partial{u^{\alpha}}},$
(2.3)

where the usual convention of summation over repeated indices is adopted [5].
In fact this convention is adopted in all subsequent expressions. We say that
(2.2) is a symmetry of (2.1) if and only if for all $\alpha=1,\ldots,m$

$X^{(r)}[F_{\alpha}(x,u,\ldots,u_{(r)})]\big{|}_{(\ref{3.1})}=0,$ (2.4)

where $X^{(r)}$ is the $r$th extension of (2.3) defined by

$\displaystyle X^{(r)}$ $\displaystyle=$
$\displaystyle\xi^{i}(x,u)\mbox{$\partial$}_{x^{i}}+\eta^{\alpha}(x,u)\mbox{$\partial$}_{u^{\alpha}}+\zeta_{i}^{\alpha}(x,u,u_{(1)})\mbox{$\partial$}_{u_{i}^{\alpha}}+\zeta_{i_{1}i_{2}}^{\alpha}(x,u,u_{(1)},u_{(2)})\mbox{$\partial$}_{u_{i_{1}i_{2}}^{\alpha}}$
(2.5) $\displaystyle+\,\cdots+\zeta_{i_{1}i_{2}\ldots
i_{r}}^{\alpha}(x,u,u_{(1)},\ldots,u_{(r)})\mbox{$\partial$}_{u_{i_{1}i_{2}\ldots
i_{r}}^{\alpha}},$

with the explicit formulas for the extended infinitesimal coefficients given
recursively by

$\displaystyle\zeta_{i}^{\alpha}$ $\displaystyle=$ $\displaystyle
D_{i}\left(\eta^{\alpha}\right)-u_{j}^{\alpha}D_{i}\left(\xi^{j}\right),\quad
i=1,2,\ldots,n,$ (2.6) $\displaystyle\zeta_{i_{1}i_{2}\ldots i_{s}}^{\alpha}$
$\displaystyle=$ $\displaystyle D_{i_{s}}\left(\zeta_{i_{1}i_{2}\ldots
i_{s-1}}^{\alpha}\right)-u_{ji_{1}i_{2}\ldots
i_{s-1}}^{\alpha}D_{i_{s}}\left(\xi^{j}\right),\quad s>1,$ (2.7)

where $D_{i}$ is the total derivative operator with respect to $x^{i}$ defined
by

$D_{i}=\frac{\mbox{$\partial$}}{\mbox{$\partial$}x^{i}}+u_{i}^{\alpha}\frac{\mbox{$\partial$}}{\mbox{$\partial$}u^{\alpha}}+u_{ij}^{\alpha}\frac{\mbox{$\partial$}}{\mbox{$\partial$}u_{j}^{\alpha}}+\cdots+u_{ii_{1}i_{2}\ldots
i_{n}}^{\alpha}\frac{\mbox{$\partial$}}{\mbox{$\partial$}u_{i_{1}i_{2}\ldots
i_{n}}^{\alpha}}+\cdots.$ (2.8)

Note that in terms of the total derivative operator the derivatives of
$u^{\alpha}$ with respect to $x_{i}$ are

$u_{i}^{\alpha}=D_{i}\left(u^{\alpha}\right),\quad
u_{ij}^{\alpha}=D_{i}\left(u_{j}^{\alpha}\right)-D_{i}D_{j}\left(u^{\alpha}\right),\quad\ldots$

A conserved vector of (2.1) is an $n$-tuple

$C=\left(C^{1}(x,u,\ldots,u_{(r-1)}),\ldots,C^{n}(x,u,\ldots,u_{(r-1)})\right)$
(2.9)

such that

$D_{i}\left(C^{i}\right)=0$ (2.10)

on the solution space of (2.1). The expression (2.10) is a conservation law of
(2.1).

##  3 The connection between conservation laws and admitted symmetries

A fundamental relationship between symmetries and conservation laws is
provided by Noether’s theorem [3], which states that for Euler-Lagrange
systems of differential equations, to each _Noether symmetry_ associated with
the Lagrangian there corresponds a conservation law which can be determined
explicitly by a formula [5]. Noether’s theorem therefore reduces the search
for conservation laws to a search for Noether symmetries. However, the
dependence upon knowledge of a suitable Lagrangian to exploit Noether’s
theorem diminishes the applicability of the theorem significantly. Ibragimov’s
theorem [9] extends the application of Noether’s theorem by providing for the
association of a conservation law to every symmetry of a system of
differential equations, albeit with the proviso that the number of equations
in the system equals the number of dependent variables and that the given
system be considered together with the associated adjoint system. The rest of
this section introduces the essential elements of Noether’s theorem and
Ibragimov’s theorem.

###  3.1 Noether’s theorem [3]

Consider a system of differential equations identical with Euler-Lagrange
equations

$\frac{\delta L}{\delta u^{\alpha}}\equiv\frac{\partial L}{\partial
u^{\alpha}}-D_{i}\left(\frac{\partial L}{\partial
u_{i}^{\alpha}}\right)=0,\quad\alpha=1,\ldots,m,$ (3.1)

arising from the variational integral

$\int L(x,u,\ldots,u_{(r)})\mbox{\rm d}x$ (3.2)

taken over an arbitrary $n$-dimensional domain in the space of the independent
variables $x=(x^{1},\ldots,x^{n})$. The Lagrangian $L(x,u,\ldots,u_{(r)})$
involves $x$ and the dependent variables $u=(u^{1},\ldots,u^{n})$, $u=u(x)$,
together with their derivatives $u_{(i)}$.

Let the system of differential equations (2.1) admit a continuous group $G$
with a generator (2.3). Noether’s theorem states that if the variational
integral (3.2) is invariant under the group $G$, then the vector field
$C=\left(C^{1},\ldots,C^{n}\right)$ defined by

$\displaystyle C^{i}$ $\displaystyle=$ $\displaystyle
L\xi^{i}+\left(\eta^{\alpha}-\xi^{i}u_{j}^{\alpha}\right)\left[\frac{\partial
L}{\partial u_{i}^{\alpha}}-D_{j}\left(\frac{\partial L}{\partial
u_{ij}^{\alpha}}\right)+D_{j}D_{k}\left(\frac{\partial L}{\partial
u_{ijk}^{\alpha}}\right)-\cdots\right]$ (3.3)
$\displaystyle+\,\,D_{j}\left(\eta^{\alpha}-\xi^{i}u_{j}^{\alpha}\right)\left[\frac{\partial
L}{\partial u_{ij}^{\alpha}}-D_{k}\left(\frac{\partial L}{\partial
u_{ijk}^{\alpha}}\right)+\cdots\right]$
$\displaystyle+\,\,D_{j}D_{k}\left(\eta^{\alpha}-\xi^{i}u_{j}^{\alpha}\right)\left[\frac{\partial
L}{\partial u_{ijk}^{\alpha}}-\cdots\right]+\cdots$

provides a conservation law for the Euler-Lagrange equations (3.1). Noether’s
theorem states that if the Invariance of the variational integral (3.2) under
the group $G$ is established via the infinitesimal test for invariance,

$X\left(L\right)+LD_{i}\left(\xi^{i}\right)=0,$ (3.4)

where the appropriate prolongation of $X$ is understood.

###  3.2 Extension of Noether’s theorem: Conservation Laws via Ibragimov’s
theorem [9]

Consider a system of $r$th-order differential equations defined in (2.1). We
introduce the differential functions

$F_{\alpha}^{*}(x,u,v,\ldots,u_{(r)},v_{(r)})=\frac{\delta\left(v^{\beta}F_{\beta}\right)}{\delta
u^{\alpha}},\quad\alpha=1,\ldots,m$ (3.5)

where $v=(v^{1},\ldots,v^{m})$ are new dependent variables, $v=v(x)$, and

$\frac{\delta}{\delta u^{\alpha}}=\frac{\partial}{\partial
u^{\alpha}}+\sum(-1)^{r}D_{i_{1}},\ldots,D_{i_{r}}\frac{\partial}{\partial
u_{i_{1},\ldots,i_{r}}^{\alpha}},\quad\alpha=1,\ldots,m.$ (3.6)

The system of adjoint equations to (2.1) is defined by

$F_{\alpha}^{*}(x,u,v,\ldots,u_{(r)},v_{(r)})=0,\quad\alpha=1,\ldots,m.$ (3.7)

We now have that the simultaneous system consisting of the $r$th-order
differential equations (2.1) considered together with its adjoint equation
(3.7) has a Lagrangian defined by

$L=v^{\beta}F_{\beta}(x,u,\ldots,u_{(r)}).$ (3.8)

Furthermore, the adjoint system (3.7) inherits the symmetries of the system
(2.1) in the sense that if the system (2.1) admits a point transformation
group with a generator

$X=\xi^{i}(x,u)\frac{\partial}{\partial
x^{i}}+\eta^{\alpha}(x,u)\frac{\partial}{\partial u^{\alpha}},$ (3.9)

then the adjoint system (3.7) admits the operator (3.9) extended to the
variables $v^{\alpha}$ by the formula

$Y=\xi^{i}\frac{\partial}{\partial
x^{i}}+\eta^{\alpha}\frac{\partial}{\partial
u^{\alpha}}+\eta_{*}^{\alpha}\frac{\partial}{\partial v^{\alpha}},$ (3.10)

with coefficients $\eta_{*}^{\alpha}=\eta_{*}^{\alpha}(x,u,v,\ldots)$ chosen
in such a way that $Y$ satisfies the infinitesimal test for invariance of the
variational integral associated with (3.8), i.e,

$Y\left(L\right)+LD_{i}\left(\xi^{i}\right)=0,$ (3.11)

where the generator $Y$ is prolonged appropriately to the $r$th derivatives
$u_{(r)}$ and $v_{(r)}$. It turns out that

$\eta_{*}^{\alpha}=-\left[\lambda_{\beta}^{\alpha}v^{\beta}+v^{\alpha}D_{i}\left(\xi^{i}\right)\right],\quad\alpha=1,\ldots,m.$
(3.12)

with $\lambda_{\beta}^{\alpha}$ defined by the invariance condition (2.4) in
the form

$X\left(F_{\alpha}\right)=\lambda_{\beta}^{\alpha}F_{\beta},\quad\alpha=1,\ldots,m,$
(3.13)

where the prolongation of $X$ to all derivatives involved in Lagrangian the
system (2.1) is understood.

Noether’s theorem is now employed to furnish the conserved vector
$C^{i}=\left(C^{1},\ldots,C^{n}\right)$, with $C_{i}$ defined by

$\displaystyle C^{i}$ $\displaystyle=$ $\displaystyle
L\xi^{i}+\left(\eta^{\alpha}-\xi^{i}u_{j}^{\alpha}\right)\left[\frac{\partial
L}{\partial u_{i}^{\alpha}}-D_{j}\left(\frac{\partial L}{\partial
u_{ij}^{\alpha}}\right)+D_{j}D_{k}\left(\frac{\partial L}{\partial
u_{ijk}^{\alpha}}\right)-\cdots\right]$ (3.14)
$\displaystyle+\,\,D_{j}\left(\eta^{\alpha}-\xi^{i}u_{j}^{\alpha}\right)\left[\frac{\partial
L}{\partial u_{ij}^{\alpha}}-D_{k}\left(\frac{\partial L}{\partial
u_{ijk}^{\alpha}}\right)+\cdots\right]$
$\displaystyle+\,\,D_{j}D_{k}\left(\eta^{\alpha}-\xi^{i}u_{j}^{\alpha}\right)\left[\frac{\partial
L}{\partial u_{ijk}^{\alpha}}-\cdots\right]+\cdots$

where $u^{\alpha}=\left(u^{1},\ldots,u^{m},v^{1},\ldots,v^{m}\right).$

##  4 Conservation Laws of Anderson’s HIV model

For the problem at hand, (1.1), we have a system of three dependent variables
$u=(u^{1},u^{2},u^{3})$ and one independent variable $t$, where $u=u(t)$. This
system is considered together with the corresponding adjoint system of
equations which is constructed as outlined in Section 3. Let
$v=(v^{1},v^{2},v^{3})$ be the new dependent variables, $v=v(t)$. According to
(3.7) the adjoint system is given by

$F_{\alpha}^{*}(t,u,v,u_{(1)},v_{(1)})=\frac{\delta\left(v^{\beta}F_{\beta}\right)}{\delta
u^{\alpha}}=0,\quad\alpha=1,2,3,$ (4.1)

where $\delta/\delta u^{\alpha}$ is the Euler-Lagrange operator defined by
(3.6), which for the problem being considred reduces to

$\frac{\delta}{\delta u^{\alpha}}=\frac{\partial}{\partial
u^{\alpha}}-D_{t}\frac{\partial}{\partial u_{t}^{\alpha}},\quad\alpha=1,2,3,$
(4.2)

where $D_{t}$ is the total derivative operator with respect to $t$ defined by

$D_{t}=\frac{\mbox{$\partial$}}{\mbox{$\partial$}t}+u^{1}_{t}\frac{\mbox{$\partial$}}{\mbox{$\partial$}u^{1}}+u^{2}_{t}\frac{\mbox{$\partial$}}{\mbox{$\partial$}u^{2}}+u^{3}_{t}\frac{\mbox{$\partial$}}{\mbox{$\partial$}u^{3}}+\cdots.$

Thus (4.1) translates into the following adjoint system of equations:

$\displaystyle\frac{dv^{1}}{dt}$ $\displaystyle=$ $\displaystyle
v^{1}\left(\frac{u^{2}\,\beta\,\delta}{u^{1}+u^{2}+u^{3}}-\frac{u^{1}\,u^{2}\,\beta\,\delta}{{\left(u^{1}+u^{2}+u^{3}\right)}^{2}}+\mu\right)-\frac{v^{2}\,u^{2}\left(u^{2}+u^{3}\right)\,\beta\,\delta}{{\left(u^{1}+u^{2}+u^{3}\right)}^{2}},$
$\displaystyle\frac{dv^{2}}{dt}$ $\displaystyle=$ $\displaystyle
v^{2}\left(\mu+\nu-\frac{u^{1}\left(u^{1}+u^{3}\right)\,\beta\,\delta}{{\left(u^{1}+u^{2}+u^{3}\right)}^{2}}\right)+\frac{v^{1}\,u^{1}\left(u^{1}+u^{3}\right)\,\beta\,\delta}{{\left(u^{1}+u^{2}+u^{3}\right)}^{2}}-v^{3}\,\nu,$
(4.3) $\displaystyle\frac{dv^{3}}{dt}$ $\displaystyle=$ $\displaystyle
v^{3}\,\alpha-\frac{v^{1}\,u^{1}\,u^{2}\,\beta\,\delta}{{\left(u^{1}+u^{2}+u^{3}\right)}^{2}}+\frac{v^{2}\,u^{1}\,u^{2}\,\beta\,\delta}{{\left(u^{1}+u^{2}+u^{3}\right)}^{2}}.$

According to the infinitesimal condition for invariance (2.4) the system (1.1)
admits a symmetry group with the infinitesimal generator

$X=\xi(t,u^{1},u^{2},u^{3})\partial_{t}+\eta^{1}(t,u^{1},u^{2},u^{3})\partial_{u^{1}}+\eta^{2}(t,u^{1},u^{2},u^{3})\partial_{u^{2}}+\eta^{3}(t,u^{1},u^{2},u^{3})\partial_{u^{3}}$

if and only if

$\left.X^{(1)}F_{1}\right|_{(\ref{gptra4abc})}=0,\quad\left.X^{(1)}F_{2}\right|_{(\ref{gptra4abc})}=0,\quad\left.X^{(1)}F_{3}\right|_{(\ref{gptra4abc})}=0,$
(4.4)

where

$X^{(1)}\\!=\\!\xi\partial_{t}+\eta^{1}\partial_{u^{1}}+\eta^{2}\partial_{u^{2}}+\eta^{3}\partial_{u^{3}}+\zeta_{t}^{1}\partial_{u_{t}^{1}}+\zeta_{t}^{2}\partial_{u_{t}^{2}}+\zeta_{t}^{3}\partial_{u_{t}^{3}}$
(4.5)

with

$\displaystyle\zeta_{t}^{1}$ $\displaystyle=$ $\displaystyle
D_{t}\left(\eta^{1}\right)-u^{1}_{t}D_{t}\,\xi$ (4.6)
$\displaystyle\zeta_{t}^{2}$ $\displaystyle=$ $\displaystyle
D_{t}\left(\eta^{2}\right)-u^{2}_{t}D_{t}\,\xi$ (4.7)
$\displaystyle\zeta_{t}^{3}$ $\displaystyle=$ $\displaystyle
D_{t}\left(\eta^{3}\right)-u^{3}_{t}D_{t}\,\xi.$ (4.8)

After making an ansatz on the form of the operator $X$ (assuming that the
functions $\xi$, $\eta^{1}$ and $\eta^{2}$ are polynomials of second degree on
$u^{1}$, $u^{2}$ and $u^{3}$) the solution of the equations (4.4), after
lengthy analysis, leads to a three-dimensional Lie symmetry algebra
$L_{3}=\left<X_{1},X_{2},X_{3}\right>$ admitted by (1.1) with the following
basis operators [15]:

$\displaystyle X_{1}$ $\displaystyle=$ $\displaystyle\partial_{t},$
$\displaystyle X_{2}$ $\displaystyle=$ $\displaystyle
u^{1}\,\partial_{u^{1}}+u^{2}\,\partial_{u^{2}}+u^{3}\,\partial_{u^{3}},$
(4.9) $\displaystyle X_{3}$ $\displaystyle=$ $\displaystyle
e^{-(\mu+\nu)\,t}\partial_{u^{2}}+\frac{u^{1}+u^{3}}{u^{2}}\,\,e^{-(\mu+\nu)\,t}\,\partial_{u^{3}}.$

We shall find conservation laws for the simultaneous system (1.1) and (4)
using each of the symmetries in (4), extended suitably to the adjoint
variables $v^{1}$, $v^{2}$ and $v^{3}$. According to (3.8) the system (1.1)
and (4) considered together has the Lagrangian

$\displaystyle L$ $\displaystyle=$ $\displaystyle
v^{\beta}F_{\beta}=v^{1}\left[u_{t}^{1}-\frac{\beta
cu^{1}u^{2}}{u^{1}+u^{2}+u^{3}}+\mu
u^{1}\right]-v^{2}\left[u_{t}^{2}-\frac{\beta
cu^{1}u^{2}}{u^{1}+u^{2}+u^{3}}+(\nu+\mu)u^{2}\right]$ (4.10)
$\displaystyle\quad\quad\quad+\,v^{3}\left[u_{t}^{3}-\nu u^{2}+(\mu+\beta
c)u^{3}\right].$

Each of the symmetries in (4) has the form

$X=\xi\,\partial_{t}+\eta^{1}\partial_{u^{1}}+\eta^{2}\partial_{u^{2}}+\eta^{3}\partial_{u^{3}}$
(4.11)

and needs to be extended appropriately to the operator

$Y=X+\eta_{*}^{1}\,\partial_{v^{1}}+\eta_{*}^{2}\,\partial_{v^{2}}+\eta_{*}^{3}\,\partial_{v^{3}},$
(4.12)

to cater for the adjoint differential variables. It turns out that (4.12) is
admitted by the adjoint system (4) provided the infinitesimal coefficients
$\eta_{*}^{\alpha}$ in the operator are prescribed as follows [9]:

$\eta_{*}^{\alpha}=-\left[\lambda_{\beta}^{\alpha}v^{\beta}+v^{\alpha}D_{t}\left(\xi\right)\right],\quad\alpha=1,2,3,$
(4.13)

where $\lambda_{\beta}^{\alpha}$ are defined by the invariance condition
(3.13). Taking $X_{1}$ from (4) and extending it once to $X_{1}^{(1)}=X_{1}$ ,
condition (3.13) leads to a set of equations,

$\displaystyle X_{1}^{(1)}F_{1}$ $\displaystyle=$
$\displaystyle\lambda_{1}^{1}\,F_{1}+\lambda_{2}^{1},F_{2}+\lambda_{3}^{1}\,F_{3}$
(4.14) $\displaystyle X_{1}^{(1)}F_{2}$ $\displaystyle=$
$\displaystyle\lambda_{1}^{2}\,F_{1}+\lambda_{2}^{2},F_{2}+\lambda_{3}^{2}\,F_{3}$
(4.15) $\displaystyle X_{1}^{(1)}F_{3}$ $\displaystyle=$
$\displaystyle\lambda_{1}^{3}\,F_{1}+\lambda_{2}^{3},F_{2}+\lambda_{3}^{3}\,F_{3},$
(4.16)

which must be solved for $\lambda_{\beta}^{\alpha}$. Clearly
$X_{1}^{(1)}F_{\alpha}=0$ for all $\alpha=1,2,3$. It follows, therefore, that
$\lambda_{\beta}^{\alpha}=0$ for all $\alpha$ and $\beta$. Furthermore,
$\xi=1$ in the generator $X_{1}$, which leads us to the conclusion that in
this case $\eta_{*}^{\alpha}=0$ for all $\alpha=1,2,3$. Therefore

$\displaystyle Y_{1}$ $\displaystyle=$ $\displaystyle X_{1}=\partial_{t}.$
(4.17)

For $X_{2}$ and $X_{3}$, the coefficient $\xi$ equals zero, which means that
the second term in (3.12) vanishes. Proceeding as we did for $X_{1}$ we
determine $\lambda_{\beta}^{\alpha}$ in each of the two cases and obtain
$Y_{2}$ and $Y_{3}$, via the desired extensions of $X_{2}$ and $X_{3}$
respectively:

$\displaystyle Y_{2}$ $\displaystyle=$ $\displaystyle
u^{1}\,\partial_{u^{1}}+u^{2}\,\partial_{u^{2}}+u^{3}\,\partial_{u^{3}}-v^{1}\,\partial_{v^{1}}-v^{2}\,\partial_{v^{2}}-v^{3}\,\partial_{v^{3}},$
(4.18) $\displaystyle Y_{3}$ $\displaystyle=$ $\displaystyle
e^{-(\mu+\nu)\,t}\left[\partial_{u^{2}}+\frac{u^{1}+u^{3}}{u^{2}}\,\partial_{u^{3}}-c\left(\frac{1}{u^{2}}\,\partial_{v^{1}}-\frac{u^{1}+u^{3}}{(u^{2})^{2}}\,\partial_{v^{2}}+\frac{1}{u^{2}}\,\partial_{v^{3}}\right)\right].$
(4.19)

By applying Noether’s theorem to each of the generators $Y_{i}$ with the
associated Lagragian (4.10) we wish to find the corresponding conservation
laws. Let us rename the dependent variables $v^{1}$, $v^{2}$ and $v^{3}$ of
the adjoint system as $u^{4}$, $u^{5}$ and $u^{6}$, respectively, so that each
of the generators (4.17), (4.18) and (4.19) is in the form

$Y=\xi(t,u)\frac{\partial}{\partial
t}+\eta^{\alpha}(t,u)\frac{\partial}{\partial u^{\alpha}},$ (4.20)

where

$u=\left(u^{1},u^{2},u^{3},v^{1},v^{2},v^{3}\right)=\left(u^{1},u^{2},u^{3},u^{4},u^{5},u^{6}\right).$

According to (2.10) and (3.14), the conservation law corresponding to (4.20)
is

$D_{t}\left\\{\xi L+\left(\eta^{\alpha}-\xi
u^{\alpha}\right)\left[\frac{\partial L}{\partial
u^{\alpha}}-D_{t}\frac{\partial L}{\partial u^{\alpha}}\right]\right\\}=0.$
(4.21)

We therefore obtain the following conservation laws:

$\displaystyle
D_{t}\left\\{\frac{u^{1}\,u^{2}\,\beta\,\delta\left(v^{1}-v^{2}\right)}{u^{1}+u^{2}+u^{3}}+v^{1}\,u^{1}\,\mu+v^{2}\,u^{2}\left(\mu+\nu\right)+v^{3}\left(u^{3}\,\alpha-u^{2}\,\nu\right)\right\\}$
$\displaystyle=$ $\displaystyle 0$ $\displaystyle
D_{t}\left\\{v^{1}\,u^{1}+v^{2}\,u^{2}+v^{3}\,u^{3}\right\\}$ $\displaystyle=$
$\displaystyle 0$ (4.22) $\displaystyle
D_{t}\left\\{e^{-(\mu+\nu)\,t}\left[v^{2}+\frac{v^{3}\left(u^{1}+u^{3}\right)}{(u^{1})^{2}}\right]\right\\}$
$\displaystyle=$ $\displaystyle 0,$

corresponding to each of the extended symmetries $Y_{1}$, $Y_{2}$ and $Y_{3}$,
respectively, where $v=\left(v^{1},v^{2},v^{3}\right)$ is any solution of the
adjoint system of equations (4).

##  5 Concluding remarks

In this paper we have considered a nontrivial nonlinear system of first-order
ordinary differential equations arising from a mathematical model formulated
by Anderson [11] to describe the transmission of HIV/AIDS in male
homosexual/bisexual cohorts. We have applied Ibragimov’s theorem and
constructed conservation laws of the system considered together with the
associated adjoint system. We have, however, not attempted to attach physical
meaning to the adjoint system (4) and/or the constructed conservation laws
(4). We defer this to possible future work on the model (1.1). The application
reported in this paper is an instructive application of Ibragimov’s new
theorem and may be used to construct conservation laws in other settings
involving systems of first-order ordinary differential equations.

##  6 Acknowledgment

Support from the Directorate of Research Development of Walter Sisulu
University is gratefully acknowledged.

## References

  * [1] A.H. Kara, F.M. Mahomed, The relationship between symmetries and conservation laws, Int J Theor Phys 39(2000)23–40.
  * [2] B. Muatjetjeja, C. M. Khalique, Symmetry analysis and conservation laws for a coupled (2 + 1)-dimensional hyperbolic system, Commun Nonlinear Sci Numer Simulat 22(2014)1252–1262.
  * [3] E. Noether, Invariante Variationsprobleme, Nachr. König. Gesell. Wissen. Göttingen, Math.-Phys. Kl., (1918)235–57 English translation in Transport Theory and Statistical Physiscs 1(1971)186-207.
  * [4] E. Yaşar, T. Özer, Invariant solutions and conservation laws to nonconservative FP equation. Comput Math Appl 59(2010)3203–3210.
  * [5] G.W. Bluman, S. Kumei, Symmetries and differential equations, Springer, New York, 1989.
  * [6] G.W.Bluman, Temuerchaolu, S.C. Anco, New conservation laws obtained directly from symmetry action on known conservation laws, J Math Anal Appl 322(2005)233 250
  * [7] I.L. Freire, J.S.C. Sampaio, On the nonlinear self-adjointness and local conservation laws for a class of evolution equations unifying many models, Commun Nonlinear Sci Numer Simulat 19(2014)350–360.
  * [8] L. Wei, J. Zhang, Self-adjointness and conservation laws for Kadomtsev-Petviashvili-Burgers equation, Nonlinear Anal. Real World Appl 23(2015)123–8.
  * [9] N.H. Ibragimov, A new conservation theorem, J Math Anal Appl 333(2007)311–28.
  * [10] P.J. Olver, Application of Lie group to differential equations, Springer, New York, 1993.
  * [11] R.M. Anderson, The Role of Mathematical Models in the study of HIV Transmission and the Epidemiology of AIDS, J Acquir Immune Defic Syndr 1(1988)241–256.
  * [12] R.M. Edelstein, K.S. Govinder, Conservation laws for the Black Scholes equation, Nonlinear Anal. Real World Appl 10(2009)3372–3380
  * [13] R. Naz, F.M. Mahomed,D.P. Mason, Comparison of different approaches to conservation laws for some partial differential equations in fluid mechanics, Appl Math Comput 205(2008)212–230.
  * [14] S.C.Anco, N.M. Ivanova, Conservation laws and symmetries of semilinear radial wave equations, J Math Anal Appl 322(2007)863–876.
  * [15] V. Torrisi, M.C. Nucci, Application of Lie group analysis to a mathematical model which describe HIV transmission, J Amer Math Soc 285(2001)11–20.