##  8\. Proof of the Main Theorem

For clarity, we restate the main theorem of this paper here.

######  Theorem 8.1 (Existence of Entropy Solutions).

Let $(\rho_{0},u_{0})$ be measurable and bounded initial data function
satisfying

$|u_{0}(x)|\leq M_{0}<\frac{1}{\sqrt{\varepsilon}},\quad
0\leq\rho_{0}(x)\leq\rho_{M_{0}}<\rho_{\max}^{\varepsilon}\qquad\,\,\text{ for
a.e. $x\in{\mathbb{R}}$,}$

for some $M_{0}$ and $\rho_{M_{0}}$, and let the pressure function $p(\rho)$
satisfy $p^{\prime}(\rho)>0$, (1.6), and (1.9). Then there exists
$\varepsilon_{0}>0$ such that the following holds: There exists a sequence
$(\rho^{\delta},u^{\delta})$ of the viscosity solutions to the approximate
equations (7.1) such that, if $\varepsilon\leq\varepsilon_{0}$, then sequence
$(\rho^{\delta},u^{\delta})$ converges to an entropy solution $(\rho,u)$ of
(1.1) in the sense of Definition 3.2 such that

$|u(t,x)|\leq M<\frac{1}{\sqrt{\varepsilon}},\quad
0\leq\rho(t,x)\leq\rho_{M}<\rho_{\max}^{\varepsilon}\qquad\,\,\text{ for a.e.
$(t,x)\in{\mathbb{R}}^{2}_{+}$},$

for some $M$ and $\rho_{M}$ depending only on $M_{0}$ and $\rho_{M_{0}}$. The
convergence is a.e. and strong in $L^{p}_{{\text{\rm
loc}}}({\mathbb{R}}^{2}_{+})$ for all $p\in[1,\infty)$.

To prove the theorem, we begin by constructing approximate initial data
satisfying the assumptions of Theorem 7.1 and Lemma 7.3. For each $\delta>0$,
we cut off the initial data outside interval $(-\delta^{-1},\delta^{-1})$,
mollify with a standard Friedrichs mollifier, and add a small positive
constant $\bar{\rho}^{\delta}>0$ to the approximate density. This gives us the
initial data $(\rho_{0}^{\delta},u_{0}^{\delta})$ satisfying the desired
assumptions such that $(\rho_{0}^{\delta},u_{0}^{\delta})\to(\rho_{0},u_{0})$
a.e. and in $L^{p}_{{\text{\rm loc}}}({\mathbb{R}})$. Then Theorem 7.1 gives
the existence of the approximate solutions, and the results of §7 then imply
that these approximate solutions are uniformly bounded and satisfy the
condition:

$\eta(U^{\delta})_{t}+q(U^{\delta})_{x}\qquad\text{ are compact in
$H^{-1}_{{\text{\rm loc}}}({\mathbb{R}}^{2}_{+})$}$

for all weak entropy pairs $(\eta,q)$.

The uniform bound on the approximate solutions gives a constant $\rho_{M}>0$
independent of $\delta$ such that $\rho^{\delta}\leq\rho_{M}$ for all
$\delta$. Taking $\varepsilon_{0}>0$ as in Proposition 5.6 then implies that
the sequence of approximate solutions $(\rho^{\delta},u^{\delta})$ satisfies
the assumptions of the compactness framework of §6. Thus Theorem 6.1 gives the
strong convergence of the approximate solutions
$(\rho^{\delta},u^{\delta})\to(\rho,u)$ in $L^{r}_{{\text{\rm
loc}}}({\mathbb{R}}^{2}_{+})$ for all $r\in[1,\infty)$. We can see that the
obtained limit is an entropy solution of (1.1).

Finally, we remark that the proof of the existence part of Theorem 1.4 is
similar.

##  9\. Newtonian Limit

The final section of this paper is devoted to the proof of our second main
result, Theorem 1.3, concerning the Newtonian limit. We first introduce some
notations to make the explicit dependence of the quantities on $\varepsilon$
in this section. For $\varepsilon\in(0,1)$, we rewrite (1.1) as

$\partial_{t}U^{\varepsilon}+\partial_{x}F^{\varepsilon}(U^{\varepsilon})=0,$
(9.1)

where
$U^{\varepsilon}(\rho^{\varepsilon},u^{\varepsilon})=\Big{(}\frac{n(\rho^{\varepsilon})}{\sqrt{1-\varepsilon(u^{\varepsilon})^{2}}},\frac{(\rho^{\varepsilon}+\varepsilon
p(\rho^{\varepsilon}))u^{\varepsilon}}{1-\varepsilon(u^{\varepsilon})^{2}}\Big{)}^{\top}$
and $F^{\varepsilon}$ is the associated flux. We recall the definition of
$v^{\varepsilon}=v^{\varepsilon}(u^{\varepsilon})$ from §2:

$v^{\varepsilon}(u^{\varepsilon})=\frac{1}{2\sqrt{\varepsilon}}\log\big{(}\frac{1+\sqrt{\varepsilon}u^{\varepsilon}}{1-\sqrt{\varepsilon}u^{\varepsilon}}\big{)}.$

We write the Euler equations (1.3) as

$\partial_{t}U+\partial_{x}F(U)=0,$ (9.2)

where $U=(\rho,m)^{\top}$, $F(U)$ is the associated flux, and
$F=(m,\frac{m^{2}}{\rho}+p(\rho))^{\top}$.

###  9.1. Proof of Theorem 1.3

By assumption, $(\rho^{\varepsilon},v^{\varepsilon})$ is a uniformly bounded
sequence in $L^{\infty}$. Then we may extract a subsequence such that

$(\rho^{\varepsilon},v^{\varepsilon})\stackrel{{\scriptstyle*}}{{\rightharpoonup}}(\rho,v)\qquad\mbox{weakly-
star in $L^{\infty}$}.$

To this sequence, we associate a Young measure $\nu_{t,x}$. Our aim now is to
show that we may apply a reduction argument analogous to that of Theorem 6.2.
However, we observe that Theorem 6.2 holds only for fixed $\varepsilon>0$, not
for a sequence. We therefore recall the following theorem from [4, 5].

######  Theorem 9.1 ([4], Theorem 4.2).

Let $\nu(\rho,v)$ be a probability measure with bounded support in
$\\{\rho\geq 0,\,v\in{\mathbb{R}}\\}$ such that

$\langle\nu,\chi_{*}(s_{1})\sigma_{*}(s_{2})-\chi_{*}(s_{2})\sigma_{*}(s_{1})\rangle=\langle\nu,\chi_{*}(s_{1})\rangle\langle\nu,\sigma_{*}(s_{2})\rangle-\langle\nu,\chi_{*}(s_{2})\rangle\langle\nu,\sigma_{*}(s_{1})\rangle$

for any $s_{1},s_{2}\in{\mathbb{R}}$, where $\chi_{*}$ and $\sigma_{*}$ are
the entropy and entropy-flux kernels of the classical Euler equations (1.3).
Then the support of $\nu$ is either a single point or a subset of the vacuum
line, $\\{\rho=0\\}$.

Defining $m^{\varepsilon}:=\rho^{\varepsilon}v^{\varepsilon}$, we show that

$\partial_{t}\eta(\rho^{\varepsilon},m^{\varepsilon})+\partial_{x}q(\rho^{\varepsilon},m^{\varepsilon})\qquad\text{
is compact in $H^{-1}_{{\text{\rm loc}}}$}$

for any weak entropy pair of $(\eta,q)$ of system (9.2), and hence argue by
using the div-curl lemma as in the proof of Theorem 6.1 to deduce the
commutation relation:

$\langle\nu,\chi_{*}(s_{1})\sigma_{*}(s_{2})-\chi_{*}(s_{2})\sigma_{*}(s_{1})\rangle=\langle\nu,\chi_{*}(s_{1})\rangle\langle\nu,\sigma_{*}(s_{2})\rangle-\langle\nu,\chi_{*}(s_{2})\rangle\langle\nu,\sigma_{*}(s_{1})\rangle.$

We may then apply Theorem 9.1 to deduce the strong convergence of the sequence
$(\rho^{\varepsilon},v^{\varepsilon})$.

It remains to prove the $H^{-1}$ compactness of the entropy dissipation
measures. In a slight abuse of notation, we simply write the argument of
$(\eta,q)$ as $U^{\varepsilon}$. Then

$\partial_{t}\eta(U^{\varepsilon})+\partial_{x}q(U^{\varepsilon})=\big{(}\eta(U^{\varepsilon})-\eta^{\varepsilon}(U^{\varepsilon})\big{)}_{t}+\big{(}q(U^{\varepsilon})-q^{\varepsilon}(U^{\varepsilon})\big{)}_{x}+\eta^{\varepsilon}(U^{\varepsilon})_{t}+q^{\varepsilon}(U^{\varepsilon})_{x},$
(9.3)

where $(\eta,q)$ and $(\eta^{\varepsilon},q^{\varepsilon})$ are weak entropy
pairs generated by the same test function $\psi$ by convolution with the weak
entropy and entropy-flux kernels associated to systems (9.2) and (9.1),
respectively.

Since, for each $\varepsilon$, $U^{\varepsilon}$ is an entropy solution of
(9.1), we know that, for any convex weak entropy pairs
$(\eta^{\varepsilon},q^{\varepsilon})$,

$\eta^{\varepsilon}(U^{\varepsilon})_{t}+q^{\varepsilon}(U^{\varepsilon})_{x}\leq
0.$

Moreover,
$\eta^{\varepsilon}(U^{\varepsilon})_{t}+q^{\varepsilon}(U^{\varepsilon})_{x}$
is uniformly bounded in $H^{-1}_{{\text{\rm loc}}}$ since $U^{\varepsilon}$ is
uniformly bounded. Then Murat’s lemma [22] indicates that the injection of the
positive cone in $H^{1}$ into $W^{-1,q},1\leq q<2$, is compact, which implies
that the convex entropy dissipation measure sequence
$\eta^{\varepsilon}(U^{\varepsilon})_{t}+q^{\varepsilon}(U^{\varepsilon})_{x}$
is compact in $W^{-1,q}_{{\text{\rm loc}}}$ for any $q\in[1,2)$. Likewise, for
any concave entropy pair $(\eta,q)$, we have

$\eta^{\varepsilon}(U^{\varepsilon})_{t}+q^{\varepsilon}(U^{\varepsilon})_{x}\qquad\text{
is compact in $W^{-1,q}_{{\text{\rm loc}}}$ for any $q<2$.}$

On the other hand, this expression is clearly bounded in
$W^{-1,\infty}_{{\text{\rm loc}}}$ by the uniform bounds on $U^{\varepsilon}$.
Then the interpolation compactness theorem of [3, 11] yields that

$\eta^{\varepsilon}(U^{\varepsilon})_{t}+q^{\varepsilon}(U^{\varepsilon})_{x}\qquad\text{
is compact in $W^{-1,2}_{{\text{\rm loc}}}$}$

for any convex or concave entropy.

Considering now the first term on the right-hand side of (9.3),

$\big{(}\eta(U^{\varepsilon})-\eta^{\varepsilon}(U^{\varepsilon})\big{)}_{t}+\big{(}q(U^{\varepsilon})-q^{\varepsilon}(U^{\varepsilon})\big{)}_{x},$

we apply the bounds of Theorems 4.8 and 5.3 to estimate the difference between
the relativistic and classical entropies generated by the same test function
$\psi(s)$ by

$\begin{split}|\eta(\rho,v)-\eta^{\varepsilon}(\rho,v)|\leq&\,\int_{|v-s|\leq
k(\rho)}|\psi(s)||\chi_{*}(\rho,v-s)-\chi^{\varepsilon}(\rho,v-s)|\,{\rm
d}s\\\ \leq&\,C\varepsilon\int_{|v-s|\leq k(\rho)}|\psi(s)|\,{\rm
d}s,\end{split}$ (9.4)

where $\chi^{\varepsilon}$ is the relativistic entropy kernel of Theorem 3.4.
Similarly, $|q(\rho,v)-q^{\varepsilon}(\rho,v)|\leq C\varepsilon$.

Thus, for any compact set $K\subset\mathbb{R}^{2}_{+}$, we have

$(\eta(U^{\varepsilon})-\eta^{\varepsilon}(U^{\varepsilon}))_{t}+(q(U^{\varepsilon})-q^{\varepsilon}(U^{\varepsilon}))_{x}\to
0\qquad\text{ in $H^{-1}(K)$ as $\varepsilon\rightarrow 0$.}$ (9.5)

Therefore, the weak entropy dissipation measures (9.3) are compact in
$H^{-1}_{{\text{\rm loc}}}$. We have therefore shown the compactness of the
weak entropy dissipation measures (9.3) such that the entropy functions are
either convex or concave. Applying the div-curl lemma in the standard way, as
in the proof of Theorem 6.1, we deduce the commutation relation

$\langle\nu,\eta_{1}q_{2}-\eta_{2}q_{1}\rangle=\langle\nu,\eta_{1}\rangle\langle\nu,q_{2}\rangle-\langle\nu,\eta_{2}\rangle\langle\nu,q_{1}\rangle$

for all such weak entropy pairs.

We recall now that a weak entropy for the classical Euler equations (9.2) is
convex (respectively concave) if the generating test function is also convex
(respectively concave). As the linear span of convex functions and concave
functions is dense in our space of test functions, we argue by density to
conclude that the commutation relation also holds for the kernels themselves:

$\langle\nu,\chi_{*}(s_{1})\sigma_{*}(s_{2})-\chi_{*}(s_{2})\sigma_{*}(s_{1})\rangle=\langle\nu,\chi_{*}(s_{1})\rangle\langle\nu,\sigma_{*}(s_{2})\rangle-\langle\nu,\chi_{*}(s_{2})\rangle\langle\nu,\sigma_{*}(s_{1})\rangle\,\,\mbox{for
any $s_{1},s_{2}\in{\mathbb{R}}$}.$

We therefore apply Theorem 9.1 to deduce the strong convergence of sequence
$(\rho^{\varepsilon},v^{\varepsilon})$ a.e. and in $L^{r}_{{\text{\rm
loc}}}({\mathbb{R}}^{2}_{+})$ for all $r\in[1,\infty)$. Then the proof of
Theorem 1.3 concludes in the usual way.

The proof of the Newtonian limit of Theorem 1.3 is similar.

Acknowledgements. The research of Gui-Qiang G. Chen was supported in part by
the UK Engineering and Physical Sciences Research Council Award EP/L015811/1,
and the Royal Society–Wolfson Research Merit Award (UK). The research of
Matthew Schrecker was supported in part by the UK Engineering and Physical
Sciences Research Council Award EP/L015811/1. This paper is a continuation of
the joint program initiated by Gui-Qiang G. Chen, Philippe LeFloch, and Yachun
Li in [4, 5, 6, 7]. We think Professors Ph. LeFloch and Yachun Li for their
inputs and helpful discussions.

## References

  * [1] Ball, J. M., A version of the fundamental theorem for Young measures, PDEs and Continuum Models of Phase Transitions $($Nice, 1988$)$, 207–215, Lecture Notes in Physics, 344, Springer: Berlin (1989).
  * [2] Chen, G.-Q., Convergence of the Lax-Friedrichs scheme for isentropic gas dynamics (III), Acta Math. Sci. 8 (1988), 243–276.
  * [3] Chen, G.-Q., The compensated compactness method and the system of isentropic gas dynamics, Preprint MSRI-00527-91, Math. Sci. Res. Inst., Berkeley, (1990).
  * [4] Chen, G.-Q. and LeFloch, P. G., Compressible Euler equations with general pressure law, Arch. Ration. Mech. Anal. 153 (2000), 221–259.
  * [5] Chen, G.-Q. and LeFloch, P. G., Existence theory for the isentropic Euler equations, Arch. Ration. Mech. Anal. 166 (2003), 81–98.
  * [6] Chen, G.-Q. and Li, Y., Relativistic Euler equations for isentropic fluids: stability of Riemann solutions with large oscillation, Z. Angew. Math. Phys. 55 (2004), 903–926.
  * [7] Chen, G.-Q. and Li, Y., Stability of Riemann solutions with large oscillation for the relativistic Euler equations, J. Diff. Eqs. 202 (2004), 332–353.
  * [8] Chen, G.-Q. and Perepelitsa, M., Vanishing viscosity limit of the Navier-Stokes equations to the Euler equations for compressible fluid flow, Comm. Pure Appl. Math. 63 (2010), 1469–1504.
  * [9] Dafermos, C. M., Hyperbolic Conservation Laws in Continuum Physics, Springer-Verlag: Berlin, 2016.
  * [10] Dal Maso, G., LeFloch, P. G. and Murat, F., Definition and weak stability of nonconservative products, J. Math. Pures Appl. 74 (1995), 483–548.
  * [11] Ding, X., Chen, G.-Q. and Luo, P., Convergence of the Lax-Friedrichs scheme for isentropic gas dynamics (I)–(II), Acta Math. Sci. 5 (1985), 483–500, 501–540 (in English); and 7 (1987), 467–480; 8 (1989), 61–94 (in Chinese).
  * [12] Ding, M. and Li, Y., Global existence and non-relativistic global limits of entropy solutions to the 1D piston problem for the isentropic relativistic Euler equations, J. Math. Phys. 54 (2013), 031506.
  * [13] Ding, M. and Li, Y., Non-relativistic limits of rarefaction wave to the 1-D piston problem for the isentropic relativistic Euler equations, J. Math. Phys. 58 (2017), 081510.
  * [14] DiPerna, R. J., Convergence of the viscosity method for isentropic gas dynamics, Comm. Math. Phys. 91 (1983), 1–30.
  * [15] Gelfand, I. M. and Shilov, G. E., Generalised Functions, Vol I, Academic Press (1964).
  * [16] Hsu, C.-H., Lin, S.-S. and Makino, T., On the relativistic Euler equation, Methods Appl. Anal. 8 (2001), 159–207.
  * [17] Li, Y., Feng, D. and Wang, Z., Global entropy solutions to the relativistic Euler equations for a class of large initial data, Z. Angew. Math. Phys. 56 (2005), 239–253.
  * [18] Liang, E. P. T., Relativistic simple waves: shock damping and entropy production, Astrophys. J. 211 (1977), 361–376.
  * [19] Lions, P.-L., Perthame, B. and Tadmor, E., Kinetic formulation for the isentropic gas dynamics and p-system, Comm. Math. Phys. 163 (1994), 415–431.
  * [20] Lions, P.-L., Perthame, B. and Souganidis, P. E., Existence and stability of entropy solutions for the hyperbolic systems of isentropic gas dynamics in Eulerian and Lagrangian coordinates, Comm. Pure Appl. Math. 49 (1996), 599–638.
  * [21] Murat, F., Compacité par compensation, Ann. Scuola Norm. Sup. Pisa Sci. Fis. Mat. 5 (1978), 489–507.
  * [22] Murat, F., L’injection du cône positif de $H^{-1}$ dans $W^{-1,q}$ est compacte pour tout $q<2$, J. Math. Pures Appl. 60 (1981), 309–322.
  * [23] Olver, F. W. J., Lozier, D. W., Boisvert R. F. and Clark, C. W. (eds.), NIST Handbook of Mathematical Functions, Cambridge University Press (2010).
  * [24] Smoller, J. and Temple, B., Global solutions of the relativistic Euler equations, Comm. Math. Phys. 156 (1993), 67–99.
  * [25] Taub, A. H., Relativistic Rankine-Hugoniot equations, Physical Rev. (2), 74 (1948), 328–334.
  * [26] Tartar, L., Compensated compactness and applications to partial differential equations, In: Research Notes in Mathematics, Nonlinear Analysis and Mechanics, Herriot-Watt Symposium, Vol. IV, 136–212, Res. Notes in Math. 39, Pitman, Boston, Mass.-London (1979).
  * [27] Thorne, K. S., Gravitational collapse and the death of a star, Science 150 (1965), 1671–1679.


# Growing Steerable Neural Cellular Automata

Ettore Randazzo,1 Alexander Mordvintsev,1 Craig Fouts2

1 Google Research
{etr<EMAIL_ADDRESS>

2 Columbia University
<EMAIL_ADDRESS>
###### Abstract

Neural Cellular Automata (NCA) models have shown remarkable capacity for
pattern formation and complex global behaviors stemming from local
coordination. However, in the original implementation of NCA, cells are
incapable of adjusting their own orientation, and it is the responsibility of
the model designer to orient them externally. A recent isotropic variant of
NCA (Growing Isotropic Neural Cellular Automata) makes the model orientation-
independent - cells can no longer tell up from down, nor left from right - by
removing its dependency on perceiving the gradient of spatial states in its
neighborhood. In this work, we revisit NCA with a different approach: we make
each cell responsible for its own orientation by allowing it to "turn" as
determined by an adjustable internal state. The resulting Steerable NCA
contains cells of varying orientation embedded in the same pattern. We observe
how, while Isotropic NCA are orientation-agnostic, Steerable NCA have
chirality: they have a predetermined left-right symmetry. We therefore show
that we can train Steerable NCA in similar but simpler ways than their
Isotropic variant by: (1) breaking symmetries using only two seeds, or (2)
introducing a rotation-invariant training objective and relying on
asynchronous cell updates to break the up-down symmetry of the system.

## Introduction

Recently, Neural Cellular Automata (NCA) have shown remarkable capacity for
discovering complex rules for a variety of tasks, including digit
classification (Randazzo et al.,, 2020; Walker et al.,, 2022), texture
generation (Niklasson et al., 2021b, ; Mordvintsev et al.,, 2021),
morphogenesis (Mordvintsev et al.,, 2020, 2022; Sudhakaran et al.,, 2021;
Grattarola et al.,, 2021), adversarial attacks (Randazzo et al.,, 2021;
Cavuoti et al.,, 2022), and control systems (Variengien et al.,, 2021). For
the most part, implementations of such NCA are anisotropic: cells are
directionally-dependent on some external factors, such as a predetermined
orientation and coordinate system or an oriented perception kernel. Anisotropy
may be an undesirable property, as it makes it impossible to produce groups of
cells with inherently different orientations unless externally defined, and it
relies on non-local information (such as a grid- or axis-alignment) that may
or may not be available in practical applications.

Figure 1: Rotation-invariant Steerable NCA training pipeline. The target
pattern is augmented with extra channels (a) to break symmetries that may
interfere with training. Both the target- and NCA-grown patterns are sharpened
(b) to steer optimization into preserving fine details. A polar transformation
(c) is applied to turn the unknown rotation between two images into a
horizontal shift. Fourier-domain image matching then enables efficient
computation of pixel-matching losses across all orientations (d). The blue
plot shows losses with respect to all rotations of the original pattern. The
minimal loss value is selected for backpropagation.

### Isotropy in Cellular Automata

Isotropic cellular automata is not a novel concept: indeed, Conway’s Game of
Life is isotropic (Gardner,, 1970). Generalizations, such as Lenia (Chan,,
2019), are also inherently isotropic and rely on the initial configuration of
cells for the emergence of complex behaviours. Until very recently, however,
discovering complex isotropic rules has been achieved with limited success,
with procedures mostly dependent on manual search and evolutionary strategies.

With the advent of NCA, end-to-end differentiable optimization has provided a
powerful way of discovering more complex rules. An example of fully isotropic
NCA is presented in Mordvintsev et al., (2021), where the NCA can be seen as a
reaction-diffusion system. However, the loss used is a rotation-invariant
variation of the texture loss used in Niklasson et al., 2021b , which is not
invariant to reflections, and the resulting textures are not capable of
creating very complex shapes. In a more recent work (Mordvintsev et al.,,
2022), a fully isotropic NCA (IsoNCA) is introduced for complex pattern
formation. This is the first study on NCA that explicitly addresses issues
with training fully isotropic pattern growth, such as the inherent rotation
and reflection symmetry of any target pattern. The authors solved these issues
with two alternative approaches: (1) by breaking the up-down and left-right
symmetries of the system with the initial conditions by means of structured
seeds, and (2) by introducing a novel rotation-reflection invariant loss,
alongside some other enhancements to facilitate training.

We also want to highlight that, since NCA can be seen as a special case of
Graph Neural Networks, one can construct them to be isotropic (Grattarola et
al.,, 2021; Gala et al.,, 2023). The focus of this paper, however, is on the
more traditional grid-based NCA, such as in IsoNCA, where cells are immovable
and their perception is restricted by a local neighborhood that cannot be
arbitrarily distant.

### Chirality and Steerable NCA

While IsoNCA might be the right choice of architecture for a variety of tasks,
sometimes the property of reflection-invariance is unnecessary or even
undesirable. Many biological molecules and structures, such as proteins,
sugars, nucleic acids, and amino acids, exhibit chirality. It is possible that
chirality at this level helps break asymmetry at higher levels such as organs,
which are highly asymmetric in most animals (Inaki et al.,, 2016).

Inspired by biological cells, we develop and train NCA that have an inherent
cellular chirality, which we call Steerable NCA due to their capacity to self-
orient their perception field. Steerable NCA are capable of perceiving
gradients of states in their neighborhood, like in Mordvintsev et al., (2020),
but the orientation of the perception field is dependent on internal states
which may vary from cell to cell. Since the orientation is exclusively defined
by an angle for each cell, the x-y gradients have an inherent chirality.

Similarly to IsoNCA, we show we can also train Steerable NCA in two possible
ways: (1) Steerable NCA can be initialized by placing two different seeds (as
opposed to the minimum of 3 required for IsoNCA) to break the up-down
symmetry, and (2) a single seed can be used alongside a rotation-invariant
loss (as opposed to the rotation-reflection variant used for IsoNCA) and
stochastic updates are responsible for breaking the up-down symmetry of the
system. These simplifications are possible because Steerable NCA need only
break a single up-down symmetry, as opposed to both up-down and left-right
required by IsoNCA, resulting in simpler and faster training regimes.

We released a set of Google Colabs to easily reproduce the entirety of the
experiments presented in this paper, available at https://github.com/google-
research/self-organising-systems/tree/master/isotropic_nca.

## Steerable Neural CA model

The Steerable NCA model described here can be seen as a variant of the Growing
NCA model (Mordvintsev et al.,, 2020), differing only in its perception phase.
This section covers key features of the model design. Below, we adapt the
equivalent section in (Mordvintsev et al.,, 2022) to better showcase the
similarities and differences between the two approaches.

##### Grid

Cells exist on a regular Cartesian grid; the state of each cell is represented
by a vector

$\mathbf{s}=[s^{0}=R,s^{1}=G,s^{2}=B,s^{3}=A,...,s^{C-1}]$

where $C$ is the number of channels and the first four channels represent a
visible RGBA image. The whole grid is initialized with zero values, except for
the seed cell in which $A=1$. Cells then iteratively update their states using
only the information collected from their 3x3 Moore neighbourhood.

##### Stochastic updates

Cell updates happen stochastically: at each NCA step, cells are updated with
probability $p_{\text{upd}}$ (we use $p_{\text{upd}}=0.5$ in our experiments).
This stochasticity is meant to eliminate dependence on a global shared clock
that synchronizes the updates between cells. Previous work on NCA discusses
the impact of this strategy on NCA robustness (Niklasson et al., 2021a, ). In
Isotropic and Steerable NCA models, asynchronicity plays a critical role in
the symmetry-breaking process; this asynchronicity can be seen as a strategy
for generating noise, which has been documented to help biological systems
construct complex functions in simple ways (Samoilov et al.,, 2006).

##### “Alive” and “empty” cells

The alpha channel ($s^{3}=A$) plays a special role in determining whether a
cell is currently “alive” or “empty”; here, a cell is alive if $A>0.1$ or if
it has at least one alive cell in its 3x3 neighbourhood. The state of empty
cells is explicitly set to zeros after each NCA step.

##### Perception

Each cell collects information about the state of its neighborhood using three
filters: a per-channel discrete 3x3 Laplacian filter $K_{lap}$, and Sobel
filters $K_{x},K_{y}$ for the x- and y-axis respectively. The Laplacian filter
computes the difference between the state of the cell and the average state of
its neighbours, where:

$K_{lap}=\begin{bmatrix}1&2&1\\\ 2&-12&2\\\ 1&2&1\\\ \end{bmatrix},\quad
K_{x}=\begin{bmatrix}-1&0&1\\\ -2&0&2\\\ -1&0&1\\\ \end{bmatrix}$

and $K_{y}=K_{x}^{\top}$.

The Sobel filters accumulate information about state gradients along their
respective axis, resulting, for each cell, in a gradient vector for the
x-axis, $g_{x}$, and one for the y-axis, $g_{y}$. In Growing NCA, these state
gradients were rotated by a predetermined angle, $\theta$, to perceive
correspondingly rotated x- and y-axis gradients. In Steerable NCA, each cell
determines its own rotation angle, which it uses to rotate its perception
field accordingly. In this work, we show two variants of Steerable NCA that
accomplish the same goal in different ways: angle- and gradient-based
rotations.

We highlight that the addition of the Laplacian filter (missing in the
original Growing NCA paper) greatly improves training speed and performance
with respect to models whose perception is solely made up of Sobel filters.

##### Angle-based rotation

In this variant, each cell maintains its orientation angle,
$\mathbf{s}_{\theta}$, as one of its state channels. The cell cannot perceive
this angle in any way, but the update rule can still modify it, causing the
cell to rotate its perception over time. For each step, the final perception
state $\mathbf{p}$ is computed as follows, with input
$\mathbf{x}=(\mathbf{s},\mathbf{s}_{\theta})$:

$\displaystyle g_{x}\,$ $\displaystyle=K_{x}\ast\mathbf{s},\quad
g_{y}=K_{y}\ast\mathbf{s}$ $\displaystyle\mathbf{p}_{x}$
$\displaystyle=g_{x}cos(\mathbf{s}_{\theta})+g_{y}sin(\mathbf{s}_{\theta})$
$\displaystyle\mathbf{p}_{y}$
$\displaystyle=g_{y}cos(\mathbf{s}_{\theta})-g_{x}sin(\mathbf{s}_{\theta})$
$\displaystyle\mathbf{p}\;\>$
$\displaystyle=concat(\mathbf{s},\,K_{lap}\ast\mathbf{s},\,\mathbf{p}_{x},\,\mathbf{p}_{y})$

This variant is very intuitive, but it requires attention - rotating the
cells’ angle states accordingly - while performing external rotations on
groups of cells to prevent them from behaving unexpectedly. The gradient-based
variant is designed to circumvent this issue.

##### Gradient-based rotation

In this variant, the orientation angle for each cell is dynamically inferred
from the cell’s respective neighborhood. To accomplish this, a concentration
state, $\mathbf{s}_{conc}$, is passed through the Sobel filters to determine
its respective gradients in the x- and y-direction. This pair of numbers is
then normalized, and the resulting x and y normalized gradients are
respectively used as the cosine and sine of the cell’s implicit rotation
angle. To stabilize computations, we clip the minimum value of the dividend to
1 during normalization. This has the effect that when the gradient of
$\mathbf{s}_{conc}$ is minimal, the perceived state gradients are effectively
multiplied by zero. For each step, the final perception state $\mathbf{p}$ is
computed as follows, with input $\mathbf{x}=(\mathbf{s},\mathbf{s}_{conc})$:

$\displaystyle g_{conc}$
$\displaystyle=[K_{x}\ast\mathbf{s}_{conc},\,K_{y}\ast\mathbf{s}_{conc}]$
$\displaystyle[cos_{\theta},sin_{\theta}]$
$\displaystyle=\frac{g_{conc}}{\textit{Clip}(\textit{Norm}(g_{conc}),1,\infty)}$
$\displaystyle g_{x}\,$ $\displaystyle=K_{x}\ast\mathbf{s},\quad
g_{y}=K_{y}\ast\mathbf{s}$ $\displaystyle\mathbf{p}_{x}$
$\displaystyle=g_{x}cos_{\theta}+g_{y}sin_{\theta}$
$\displaystyle\mathbf{p}_{y}$
$\displaystyle=g_{y}cos_{\theta}-g_{x}sin_{\theta}$
$\displaystyle\mathbf{p}\;\>$
$\displaystyle=concat(\mathbf{x},\,K_{lap}\ast\mathbf{x},\,\mathbf{p}_{x},\,\mathbf{p}_{y})$

Note how in this variant, cells can perceive $\mathbf{s}_{conc}$ without
incurring any issues; as introduced in the previous paragraph, if one were to
rotate a group of cells externally, the resulting system would not abruptly
vary its behaviour.

##### Update rule

Cells stochastically update their states using a learned rule that is
represented by a two-layer neural network:

$\mathbf{s}_{t+1}=\mathbf{s}_{t}+relu(\mathbf{p}_{t}W_{0}+b_{0})W_{1}$

where parameters $W_{0}$, $b_{0}$, and $W_{1}$ have shapes $(32,192)$, $(192)$
and $(192,16)$, respectively, giving a total of 9408 learned parameters.

## Training Steerable NCA

As introduced before, Steerable NCA can be trained in two different ways: (1)
by breaking symmetries during initialization using two different seeds, or (2)
by starting with a single seed and using a rotation-invariant loss function to
train the system to break up-down symmetry itself.

### Two seeds strategy

By initializing the system with two different seeds - distinguished by their
RGB configuration and placed slightly apart from one another - we define a
more restrictive initial condition which implicitly resolves up-down symmetry.
Doing so enables us to choose an arbitrary "training alignment" for the target
pattern with which we can perform L2 loss, identically to the training
procedure in Growing NCA. Here, just like IsoNCA, each seed serves as a point
of origin for one of two sub-regions of the resultant shape. However, unlike
IsoNCA, the chirality of our model precludes the necessity of left-right
symmetry breaking, allowing us to train the model with two seeds instead of
three (see Figure 2 for one possible two seed mapping).

Figure 2: Two-seed mapping paradigm from which Steerable NCA can be trained to
grow with the traditional L2 loss and without any additional auxiliary losses.

To generate visually- and spatially-differentiated seeds, we map their RGB
encodings to distinct HSV hues and separate them by eight cells. While the
exact specification is arbitrary, this procedure guarantees that (1) the seeds
lie within the target pattern, and (2) they can be distinguished from one
another by the model. As we’ll see in the next section, this differentiation
enables us to implicitly specify the model’s orientation by imposing a
rotation on the seed structure.

### Single seed strategy

In the single seed strategy, we start from a single seed cell and train
Steerable NCA using a loss dependent on the target pattern. Similarly to
IsoNCA, we use pixel-wise differences to match the pattern produced by the
trained model to the target, but modify the loss function to make it rotation-
invariant only and not rotation-reflection invariant. Because we cannot know a
priori what rotation the resulting pattern will have (decided by the random
up-down symmetry breaking that could happen on any axis in the 2d plane), we
instead select an individual rotation of the target that minimizes the pixel-
wise loss value for each NCA-generated sample in the training batch. Note that
there is no need for reflection-invariance in the training loss because once
one symmetry is broken, the perception scheme of Steerable NCA automatically
determines an orthogonal axis always in the same way and, therefore,
deterministically breaks left-right symmetry. Figure 3 shows a visual
representation of how the minimal loss is selected.

Figure 3: How the rotation-invariant loss gets computed: to the left, as a
reference, we calculate the loss with respect to the target pattern for every
rotation and see the minimum value is at 0 degrees. To the right, we calculate
the loss of a partially-formed pattern during training. The minimum loss can
be seed at roughly 60 degrees, and it gets selected for backpropagation.

We use the same algorithm presented in IsoNCA to optimize computation, using a
polar coordinate transformation and Fast Fourier Transform (FFT) to
efficiently compute the discrepancy across different target rotations. We
refer to Mordvintsev et al., (2022) for more implementation details.

#### Auxiliary channels

Figure 4: Target patterns and corresponding auxiliary targets used for the
rotation-invariant training: (a) the spiderweb does not use auxiliary targets;
(b) the lizard uses a binary auxiliary target; (c) the heart uses both binary
and radial-encoding auxiliary targets.

This rotation-invariant loss suffers from the same issues described in IsoNCA,
where some patterns exhibit strong local minima within different rotations. We
solve this issue in the same way as IsoNCA, using auxiliary channels. Here, we
mirror the single-seed experiments in IsoNCA and showcase the same three
example patterns: the spiderweb, where no auxiliary loss is needed; the
lizard, where a binary auxiliary channel is needed; and the heart, where both
the binary and radial-encoding auxiliary channels are needed (see Figure 4).

## Results

### Two seeds experiments

By adjusting the rotation of our seed structure, we see in Figure 5 that
Steerable NCA are able to self-orient based on their initial conditions. In
this figure, we present the angle-based variant, but the gradient-based
variant demonstrates identical behavior.

Figure 5: Starting in the top-left corner and moving clockwise, 0 degree, 30
degree, 170 degree, and 280 degree patterns produced by angle-based Steerable
NCA trained with L2 loss to grow at 0 degrees. Nonzero degree rotations are
inferred by the model during growth based on the arrangement of the starting
condition’s two constituent seeds.

Modifying the diameter of the seed structure after training further highlights
the model’s dependence on the intercellular message passing that occurs during
the formative stage of model growth. We see this in Figure 6, where the cells
attempt to align with one another based on the level of separation between the
two seeds.

Figure 6: Starting in the top-left corner and moving clockwise, seeds with a
diameter of 0 cells, 2 cells, 6 cells, and 8 cells and their resultant
patterns produced by angle-based Steerable NCA trained with L2 loss to grow
from a seed with a 4 cell diameter at 0 degrees.

We further analyze this dependence in Figure 7, wherein we highlight the
directional state of angle-based Steerable NCA at snapshots before and after
the growths originating from the two seeds intersect. Here, the seed angles
are initialized randomly, and the emerging cells progressively self-orient
into a coherent vector field that describes the model’s internal
representation of its own orientation. Moreover, we see that distributing
growth between the seeds this way produces an inflection point at which the
model "decides" the orientation of the resulting pattern.

Figure 7: Snapshots of angle-based Steerable NCA angle states during the early
alignment phase in which the growths produced by each seed intersect and self-
orient.

It’s worth noting that the influence of the two seeds decays over time, in
contrast to the single-seed case presented in Figure 11, in which a "residue"
from the seed remains as a divergent point in the vector field. This suggests
some additional structure necessitated by the rotation-invariant training,
dependent upon an orientation selected by the model rather than one imposed on
it.

### Single seed experiments

Figure 8: Unfolding of three trained angle-based Steerable NCA rules. Figure
9: Comparison between angle-based Steerable NCA (top) and gradient-based
Steerable NCA (bottom). Figure 10: Step 5000 for different runs of three
trained angle-based Steerable NCA models.

Figure 8 shows the execution in time of three different trained angle-based
Steerable NCA models, starting from a single seed and running for up to 5000
steps. Here, we can observe how trained patterns tend to rotate very slowly.
This is an artifact of the rotation-invariant loss used, which does not
enforce patterns to stand still; every possible rotation of the pattern has
the same loss and, thus, a dynamically rotating final state will have the same
loss as a static, non-rotating, one. Note in Figure 10 how every run of the
same model results in a different pattern rotation, but always the same
reflection. Figure 8 and Figure 10 show the angle-based model, but gradient-
based models manifest similar behaviours. Figure 9 compares the lizard pattern
produced by the angle- and gradient-based models.

Training Steerable NCA appears to be faster than training IsoNCA. For
instance, training a stable and well-formed lizard takes only 10000 training
steps for the Steerable variant, as opposed to 30000 steps for the Isotropic
variant with the exact same training configurations. Moreover, the Steerable
models empirically grow their target shape faster. This is likely due to
Steerable NCA having to break only one symmetry as opposed to the two required
by IsoNCA.

Figure 11: Internal angles for each cell in the Steerable NCA lizard models,
snapshot at 5000 steps. Left: an angle-based model in which the angles form a
spiral originating at the center. Right: a gradient-based model in which
angles cannot have discontinuities, and instead are suppressed, forming a
manifold of adirectional cells within the lizard.

Figure 11 shows the internal angles of both Steerable NCA model variants for
the lizard pattern, after a time unfolding of 5000 steps. Here, it appears
that cells favour having diverging angles from one another. The angles in the
angle-based model (left) seem to generate a spiral originating at the center
of the pattern field. On the other hand, angles in gradient-based model
(right) cannot create discontinuities, as the angle gradients are continuous
by construction. However, there appear to be some uniform manifolds where
angles are suppressed.

## Conclusion and future work

In this work we demonstrated the capability of chiral NCA models to reliably
grow complex asymmetric patterns. We show that they are capable of both
inferring a specified orientation based on initial conditions as well as
learning to break symmetries independently using a rotation-invariant loss
function. Steerable NCA are shown to be simpler to train than their Isotropic
counterparts, while maintaining local coordination properties. We believe that
this added tool in the differentiable morphogenetic toolbox may be useful in
morphogenesis simulations and multi-agent-based modeling, where agents might
benefit from having an inner sense of orientation.

## References

  * Cavuoti et al., (2022) Cavuoti, L., Sacco, F., Randazzo, E., and Levin, M. (2022).  Adversarial takeover of neural cellular automata.  volume ALIFE 2022: The 2022 Conference on Artificial Life of ALIFE 2022: The 2022 Conference on Artificial Life.  38\.
  * Chan, (2019) Chan, B. W.-C. (2019).  Lenia: Biology of artificial life.  Complex Systems, 28(3):251–286.
  * Gala et al., (2023) Gala, G., Grattarola, D., and Quaeghebeur, E. (2023).  E(n)-equivariant graph neural cellular automata.
  * Gardner, (1970) Gardner, M. (1970).  Mathematical games.  Scientific American, 223(4):120–123.
  * Grattarola et al., (2021) Grattarola, D., Livi, L., and Alippi, C. (2021).  Learning graph cellular automata.  In Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W., editors, Advances in Neural Information Processing Systems, volume 34, pages 20983–20994. Curran Associates, Inc.
  * Inaki et al., (2016) Inaki, M., Liu, J., and Matsuno, K. (2016).  Cell chirality: its origin and roles in left-right asymmetric development.  Philos. Trans. R. Soc. Lond. B Biol. Sci., 371(1710).
  * Mordvintsev et al., (2022) Mordvintsev, A., Randazzo, E., and Fouts, C. (2022).  Growing isotropic neural cellular automata.  volume ALIFE 2022: The 2022 Conference on Artificial Life of ALIFE 2022: The 2022 Conference on Artificial Life.  65\.
  * Mordvintsev et al., (2021) Mordvintsev, A., Randazzo, E., and Niklasson, E. (2021).  Differentiable programming of reaction-diffusion patterns.  In ALIFE 2021: The 2021 Conference on Artificial Life. MIT Press.
  * Mordvintsev et al., (2020) Mordvintsev, A., Randazzo, E., Niklasson, E., and Levin, M. (2020).  Growing neural cellular automata.  Distill.  https://distill.pub/2020/growing-ca.
  * (10) Niklasson, E., Mordvintsev, A., and Randazzo, E. (2021a).  Asynchronicity in neural cellular automata.  volume ALIFE 2021: The 2021 Conference on Artificial Life of ALIFE 2021: The 2021 Conference on Artificial Life.  116\.
  * (11) Niklasson, E., Mordvintsev, A., Randazzo, E., and Levin, M. (2021b).  Self-organising textures.  Distill.  https://distill.pub/selforg/2021/textures.
  * Randazzo et al., (2021) Randazzo, E., Mordvintsev, A., Niklasson, E., and Levin, M. (2021).  Adversarial reprogramming of neural cellular automata.  Distill.  https://distill.pub/selforg/2021/adversarial.
  * Randazzo et al., (2020) Randazzo, E., Mordvintsev, A., Niklasson, E., Levin, M., and Greydanus, S. (2020).  Self-classifying mnist digits.  Distill.  https://distill.pub/2020/selforg/mnist.
  * Samoilov et al., (2006) Samoilov, M. S., Price, G., and Arkin, A. P. (2006).  From fluctuations to phenotypes: the physiology of noise.  Sci. STKE, 2006(366):re17.
  * Sudhakaran et al., (2021) Sudhakaran, S., Grbic, D., Li, S., Katona, A., Najarro, E., Glanois, C., and Risi, S. (2021).  Growing 3d artefacts and functional machines with neural cellular automata.  The 2021 Conference on Artificial Life.
  * Variengien et al., (2021) Variengien, A., Nichele, S., Glover, T., and Pontes-Filho, S. (2021).  Towards self-organized control: Using neural cellular automata to robustly control a cart-pole agent.
  * Walker et al., (2022) Walker, K., Palm, R. B., Garcia, R. M., Faina, A., Stoy, K., and Risi, S. (2022).  Physical neural cellular automata for 2d shape classification.