# Mitigation of the onset of hosing in the linear regime through plasma
frequency detuning

Mariana Moreira<EMAIL_ADDRESS>GoLP/Instituto de
Plasmas e Fusão Nuclear, Instituto Superior Técnico, Universidade de Lisboa,
1049-001 Lisboa, Portugal  Patric Muggli  Max Planck Institute for Physics,
D-80805 Munich, Germany  CERN, CH-1211 Geneva, Switzerland  Jorge Vieira
GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior Técnico,
Universidade de Lisboa, 1049-001 Lisboa, Portugal

###### Abstract

The hosing instability poses a feasibility risk for plasma-based accelerator
concepts. We show that the growth rate for beam hosing in the linear regime
(which is relevant for concepts that use a long driver) is a function of the
centroid perturbation wavelength. We demonstrate how this property can be used
to damp centroid oscillations by detuning the plasma response sufficiently
early in the development of the instability. We also develop a new theoretical
model for the early evolution of hosing. These findings have implications for
the general control of an instability’s growth rate.

In the last decades, plasma-based acceleration has been explored as a possible
avenue for the next generation of particle accelerators. Plasma can sustain
electric fields about a thousand times higher than conventional accelerator
elements. This estimate is based on a plasma electron density
$n_{0}=10^{18}~{}\mathrm{cm}^{-3}$ and the cold nonrelativistic wavebreaking
field [1, 2]
$E_{0}=m_{e}c\,\omega_{p}/e\,\approx\,96\sqrt{n_{0}\left[\mathrm{cm}^{-3}\right]}\,\left[\mathrm{V/m}\right]$,
where $c$ is the speed of light, $m_{e}$ is the electron mass, $e$ is the
elementary charge, $\omega_{p}=\sqrt{e^{2}n_{0}/\varepsilon_{0}\,m_{e}}$ is
the electron plasma frequency, and $\varepsilon_{0}$ is the vacuum
permittivity.

One approach to leveraging plasma-based acceleration for high-energy physics
applications is to drive a plasma wave (and the associated wakefields) with a
high-energy particle bunch and accelerate a trailing bunch in a single, many-
meter-long plasma, thus avoiding the complexity of staging schemes [3, 4].
This concept has been demonstrated recently in the AWAKE experiment at CERN
using a 400-GeV proton bunch as the driver through a ten-meter-long plasma,
leading to the acceleration of injected electrons from 18 MeV up to 2 GeV [5].

A plasma wave is excited effectively by an impulse with a length of the order
of a plasma wavelength $\lambda_{p}=2\pi c/\omega_{p}$. It is currently
impractical to compress a 400-GeV proton bunch into this longitudinal size,
even for lower plasma densities on the order of $10^{14}~{}\mathrm{cm}^{-3}$.
Nevertheless, high-amplitude wakefields can be reached by letting the bunch
undergo the self-modulation instability (SMI) [6, 7] and controlling its onset
with seeding [8, 9]. The same transverse wakefields that modulate the beam
envelope at $\sim\lambda_{p}$ (self-modulation) can also modulate the centroid
of the bunch when it is not propagating axisymmetrically, leading to the
hosing instability (HI) [10]. This oscillation of the drive bunch centroid can
alter the structure of the wakefields significantly, thus impacting their
quality and that of an accelerated witness bunch.

The long-bunch HI and the SMI have comparable growth rates and can couple to
each other when growing from similar seed levels [11, 12]. Nevertheless,
hosing can be avoided in a fully self-modulated bunch. This requires strongly
seeding the SMI [13], which has been accomplished in experiments by letting an
ionizing laser pulse propagate with the proton bunch and create the plasma
[8]. In the future, SMI seeding may be achieved with a preceding short
electron bunch, such that the entire proton bunch self-modulates [9, 14, 15].
When misaligned, this arrangement may seed the growth of hosing. In the short-
bunch, nonlinear wakefield regime, it has been established that hosing
eventually reaches saturation [16, 17], leading to an increased emittance.
However, saturation in the long-bunch case is not yet well understood, both in
the presence and absence of the SMI. Since this could represent a major
impediment for high-quality acceleration, it would be useful to develop
further mitigation methods for hosing in the linear wakefield regime.

In this Letter we show that the growth rate at the onset of hosing in the
linear regime of plasma wakefield acceleration is highly sensitive to the
centroid perturbation wavelength, and we exploit this property to suppress the
growth of the HI. As a proof of concept of this mitigation approach, we
consider the front portion of a long, cold electron bunch propagating through
several plasma density steps. Nevertheless, these HI properties are general
and their application as a mitigation strategy can be scaled for any
relativistic particle driver through the betatron period. We also develop a
theoretical model that can describe the early evolution of hosing and
supplement existing asymptotic models. Our work opens up an unprecedented
avenue for the control of hosing, which is vital for the future of plasma-
based acceleration.

We consider a particle bunch propagating along $z$ with a density profile
$n_{b}(\zeta,y)=n_{b0}\cdot f(\zeta)\cdot g(y,y_{c})$, where $n_{b0}$ is the
peak density and $f$ and $g$ are the normalized longitudinal and transverse
profiles, respectively. The longitudinal coordinate $\zeta=z-c\>t$ moves at
the beam velocity $v_{b}\approx c$, while $y$ is the transverse coordinate and
$y_{c}(\zeta,z)=\left<y\right>$ is the bunch centroid, where the angle
brackets denote a bunch-weighted transverse average defined as
$\left<\cdot\right>=\int\cdot~{}n_{b}(\zeta,y)~{}dy\>/\int
n_{b}(\zeta,y)~{}dy$. An equation to describe hosing can be derived from the
transverse averaging of the motion of a single relativistic bunch particle:

$\frac{d^{2}y_{c}}{dz^{2}}=\frac{m_{e}}{\gamma
M_{b}}\frac{\left<F_{y}\right>}{eE_{0}}\quad,$ (1)

where $M_{b}$ is the bunch particle mass and $F_{y}$ is the transverse
component of the force associated with the plasma wakefields driven by the
bunch. All distances are normalized to the plasma skin depth
$k_{p}^{-1}=c/\omega_{p}$ in Eq. (1) and in the rest of this Letter. Note
that, for small centroid displacements ($y_{c}\ll 1$), the Green’s function
solution to the plasma operator contained in $\left<F_{y}\right>$ of Eq. (1)
can be Taylor-expanded into a linear function of $y_{c}$. This means that, in
this approximation, Eq. (1) can always be written in the form of a harmonic
oscillator (HO) equation:

$\left(\frac{d^{2}}{dz^{2}}+k_{\mathrm{HO}}^{2}(\zeta,z)\right)y_{c}(\zeta,z)=F(\zeta,z,y_{c})\quad,$
(2)

where the specific forms of the natural wavenumber
$k_{\mathrm{HO}}^{2}(\zeta,z)$ and the driving force $F(\zeta,z,y_{c})$ depend
on the choice of geometry and transverse profile of the bunch. We also note
that a centroid velocity $v_{c}=c\cdot dy_{c}/dz$ may be defined in
association with the development of hosing.

For a Gaussian transverse profile, and assuming a constant beam envelope
$\sigma_{y}$, $\left<F_{y}\right>$ in 2D Cartesian geometry is given by [18]:

$\displaystyle\frac{\left<F_{y}\right>}{eE_{0}}=\sqrt{\frac{\pi}{8}}\frac{n_{b0}}{n_{0}}\left(\frac{q_{b}}{e}\right)^{2}\sigma_{y}\>\exp(\sigma_{y}^{2})\int_{\zeta}^{\infty}d\zeta^{\prime}\sin(\zeta-\zeta^{\prime})\>f(\zeta^{\prime})$
$\displaystyle\left\\{\exp\left[y_{c}(\zeta^{\prime})-y_{c}(\zeta)\right]\>\mathrm{erfc}\left[\frac{y_{c}(\zeta^{\prime})-y_{c}(\zeta)+2\>\sigma_{y}^{2}}{2\>\sigma_{y}}\right]\right.$
$\displaystyle\left.-\exp\left[y_{c}(\zeta)-y_{c}(\zeta^{\prime})\right]\>\mathrm{erfc}\left[\frac{y_{c}(\zeta)-y_{c}(\zeta^{\prime})+2\>\sigma_{y}^{2}}{2\>\sigma_{y}}\right]\right\\}\quad,$
(3)

where $q_{b}$ is the bunch particle charge, and $\mathrm{erfc}(x)$ is the
complementary error function.

For the following, we consider an initial perturbation of the bunch centroid
of the form $y_{c0}(\zeta)\propto\sin(\hat{k}\>\zeta)$, where
$\hat{k}=k/k_{p}$ is an arbitrary wavenumber. In order to measure the HI
growth rate, we define the growth factor

$\Gamma(z)=\frac{\int_{L}|y_{c}(z,\zeta)|\>d\zeta}{\int_{L}|y_{c}(0,\zeta)|\>d\zeta}\quad,$
(4)

where $L$ is the length of some region of interest along the bunch. Here we
will consider a window measuring $L=140~{}k_{p}^{-1}$, which represents around
22 $\lambda_{p}$, and a bunch front at $\zeta_{s}=135~{}k_{p}^{-1}$.

During the initial propagation in plasma, we can assume a constant
($z$-independent) average transverse force $\left<F_{y0}\right>$ caused by
$y_{c0}$ [replacing $y_{c}$ with $y_{c0}$ in Eq. (Mitigation of the onset of
hosing in the linear regime through plasma frequency detuning)]. This
assumption is valid within the typical timescale for bunch evolution, given by
the betatron period $k_{\beta}^{-1}=c/\omega_{\beta}$ (defined by
$\omega_{\beta}=\omega_{b}/\sqrt{2\gamma}$, where
$\omega_{b}=\sqrt{n_{b0}q_{b}^{2}/\varepsilon_{0}M_{b}}$). In this case the
solution to Eq. (1) is simply
$y_{c}(\zeta,z)=y_{c0}(\zeta)+\tfrac{1}{2}~{}z^{2}\cdot\left(\frac{m_{e}}{\gamma
M_{b}}\>\frac{\left<F_{y0}\right>}{eE_{0}}\right)$. To evaluate $\Gamma(z)$,
we consider a relativistic ($\gamma=480$) electron bunch with
$n_{b0}=0.001~{}n_{0}$, a Gaussian transverse profile with the rms transverse
size $\sigma_{y}\approx 0.27\>k_{p}^{-1}$, and a longitudinal profile
$f(\zeta)=\frac{1}{2}~{}\left(1+\cos(\sqrt{\frac{\pi}{2}}\frac{\zeta-\zeta_{c}}{\sigma_{z}})\right)$,
where $\sigma_{z}\approx 160\>k_{p}^{-1}$,
$\zeta_{c}=\zeta_{s}-\sqrt{2\pi}\sigma_{z}$, and with the limits
$\zeta=\pm\sqrt{2\pi}\sigma_{z}+\zeta_{c}$. The resulting curve (calculated
numerically) is shown in Fig. 1a) (blue) for $z=k_{\beta}^{-1}$, which is
approximately $0.74~{}\mathrm{m}$ for the above electron bunch parameters and
$n_{0}=0.5\cdot 10^{14}~{}\mathrm{cm}^{-3}$. For a proton bunch with a
population $N_{b}=3\cdot 10^{11}$, rms sizes $\sigma_{z}=12~{}\mathrm{cm}$ and
$\sigma_{r}=200~{}\mathrm{\mu m}$, and $\gamma=426$, $k_{\beta}^{-1}\approx
3.3~{}\mathrm{m}$.

The theoretical curve is confirmed by 2D Cartesian particle-in-cell
simulations [18] with OSIRIS [19] [red cross symbols in Fig. 1a), where each
symbol represents a simulation initialized with a centroid perturbation at
$k$]. An identical theoretical curve can be obtained in cylindrical
coordinates (assuming $y_{c}\ll 1$), and in fact we use 3D simulations later
in this Letter to explore the consequences of this result.

\begin{overpic}[width=433.62pt]{growthrate_prl.pdf} \put(13.0,68.0){a)}
\put(13.0,36.0){b)} \end{overpic} Figure 1: a) Hosing growth factor at
$z=k_{\beta}^{-1}$ versus wavenumber of an initial sinusoidal perturbation,
according to theory (blue curve) and simulations (red cross symbols). b)
Theoretical hosing growth spectrum along $z$. The horizontal dashed line
corresponds to the blue line-out in a).

The hosing frequency response in Fig. 1a) (note that $k$ may be interpreted as
a frequency via $f=\tfrac{k\>c}{2\pi}$) implies that there is a finite growth
rate for initial perturbations at frequencies other than the plasma frequency.
The tail in the growth rate spectrum for $k<k_{p}$ has been discussed before
in the context of long laser pulses under the term “long-wavelength hosing”
[20, 21]. For short wavelengths ($k>k_{p}$), however, the growth factor is
smaller than 1, i.e., there is a decrease of the initial perturbation
amplitude instead of growth. This property has not been previously predicted
and underpins a new HI mitigation process that we explore in this Letter. A
further observation is that the maximum growth is not attained for $k=k_{p}$
(as generally assumed for the HI) during this initial phase of propagation,
but at a slightly lower value ($k\approx 0.98~{}k_{p}$).

After some propagation in plasma, growth around $k_{p}$ becomes dominant, as
illustrated in Fig. 1b), based on the theoretical $\left<F_{y0}\right>$. Not
only does the amplitude of resonant growth become several times that of
negative growth, but the locations of both extrema, $k_{\mathrm{max}}$ and
$k_{\mathrm{min}}$, shift with increasing $z$ ($k_{\mathrm{max}}\rightarrow
k_{p}$ and $k_{\mathrm{min}}\rightarrow\infty$). The typical duration of the
transient phase where the maximum positive and negative growth amplitudes are
similar [Fig. 1a)] depends on the bunch properties, though in general this
corresponds roughly to a propagation distance $z\lesssim k_{\beta}^{-1}$.

\begin{overpic}[width=433.62pt]{phase_prl.pdf} \put(0.0,88.0){a)}
\put(0.0,33.0){b)} \end{overpic} Figure 2:  a) Initial centroid (black) and
average transverse force (blue) for three different seed wavenumbers, obtained
from 2D OSIRIS simulations at $z=0$. b) Phase shift between the initial
$y_{c}$ and $\left<F_{y}\right>$ as a function of the perturbation wavenumber,
obtained from theory (blue curve) and simulations (red cross symbols).

We find that the phase shift between the bunch centroid perturbation and the
transverse plasma response, represented by $\left<F_{y}\right>$, determines
the early HI growth dynamics in Fig. 1a). Figure 2a) displays the initial
centroid and $\left<F_{y}\right>$ from three simulations at different
perturbation wavenumbers, thereby illustrating three regimes. For
$k=0.9~{}k_{p}$ (slow growth), $\left<F_{y}\right>$ is almost in phase with
$y_{c}$, while the damping regime results from the fully out-of-phase
transverse force, which acts in the opposite direction of $y_{c}$ at every
$\zeta$ [see Fig. 2a), $k=1.1~{}k_{p}$]. For $k=k_{p}$ (resonant growth), the
wakefield response lags the centroid perturbation by $\pi/2$.

Since $\left<F_{y}\right>$ initially oscillates at $k$ (along $\zeta$), a
phase shift $\Delta\phi$ between both periodic curves can be measured
straightforwardly. The relationship between $\Delta\phi$ and $k$ is shown in
Fig. 2b), as obtained from the 2D simulations and the theoretical
$\left<F_{y0}\right>$, using a cross-correlation method [18]. There is
excellent agreement between theory and simulations for $k>0.5~{}k_{p}$. The
simulation data points for $k\leq 0.5~{}k_{p}$ [grey cross symbols in Fig.
2b)] are not valid, since the simulation window length becomes comparable to
the perturbation wavelength ($L/\lambda\lesssim 10$). The theoretical curve
was obtained by scaling the nominal window and bunch lengths by $1/k$, thus
encompassing several periods in the analysis. The three growth regimes are
again evident: two asymptotes at $\Delta\phi=0$ and $\Delta\phi=\pi$, and a
transition region where $\Delta\phi$ crosses $\pi/2$.

We note that the physics described in this section is characteristic of a
sinusoidally driven damped HO: consider the phase jump when the driving
frequency crosses the resonance, or how the curve for $\Gamma(z)$ at large $z$
[see Fig. 1b)] resembles the HO amplitude curve as a function of driving
frequency.

We wish to predict the development of the centroid displacement and velocity
during the initial propagation in plasma and for an arbitrary seed wavenumber.
Previous asymptotic models are inappropriate for our purposes, since they only
apply far enough behind the front ($k_{p}\zeta\gg k_{\beta}z$) and assume
initial centroid fluctuations at $k_{p}$ [11].

The early evolution of $y_{c}$ can be approximated by truncating a series
solution to Eq. (1) [22], which may be written as
$\partial_{z}^{2}y_{c}=\mathcal{R}\\{y_{c}\\}$, with the operator
$\mathcal{R}$ in the $\zeta$ variable. For this method to be tractable,
$\mathcal{R}$ should be a linear function of $y_{c}$, which is always the case
assuming $y_{c}\ll 1$. Writing $y_{c}$ as a power series,
$y_{c}=\sum_{n=0}^{\infty}a_{n}\>z^{n}$, we arrive at the following general
solutions for $y_{c}$ and $v_{c}$ (normalized to $k_{p}^{-1}$ and $c$,
respectively):

$\displaystyle
y_{c}(\zeta,z)=\sum_{n=0}^{\infty}\frac{z^{n}}{n!}\begin{cases}\mathcal{R}^{(\frac{n}{2})}\\{y_{c0}\\}\;,&n\text{
is even}\\\ \mathcal{R}^{(\frac{n-1}{2})}\\{v_{c0}\\}\;,&n\text{ is
odd}\end{cases}$ (5) $\displaystyle
v_{c}(\zeta,z)=\sum_{n=0}^{\infty}\frac{z^{n}}{n!}\begin{cases}\mathcal{R}^{(\frac{n}{2})}\\{v_{c0}\\}\;,&n\text{
is even}\\\ \mathcal{R}^{(\frac{n+1}{2})}\\{y_{c0}\\}\;,&n\text{ is
odd}\end{cases}$ (6)

where $y_{c0}(\zeta)$ and $v_{c0}(\zeta)$ are arbitrary initial conditions for
the centroid and centroid velocity, respectively, and the superscript in the
operator $\mathcal{R}$ signifies consecutive applications (e.g.
$\mathcal{R}^{(2)}\\{X\\}=\mathcal{R}\\{\mathcal{R}\\{X\\}\\}$ and
$\mathcal{R}^{(0)}\\{X\\}=X$). In cylindrical coordinates and for a transverse
flat-top profile $g(r)=\Theta[r_{b0}-r]$, where $r_{b0}$ is a constant bunch
radius,
$\mathcal{R}\\{X(\zeta,z)\\}=2\hat{k}_{\beta}^{2}I_{1}(r_{b0})K_{1}(r_{b0})\int_{\zeta}^{\infty}d\zeta^{\prime}\sin(\zeta-\zeta^{\prime})f(\zeta^{\prime})\left[X(\zeta,z)-X(\zeta^{\prime},z)\right]$
[11].

\begin{overpic}[width=433.62pt]{evolmodel_prl.pdf} \put(92.0,45.5){a)}
\put(92.0,21.0){b)} \end{overpic} Figure 3:  Comparison of the series model
truncated at $n$ terms (colored curves) with a 3D OSIRIS simulation (black
curves), a) along $\zeta$ at $z=2~{}k_{\beta}^{-1}$ and b) along $z$ at
$\zeta=30~{}k_{p}^{-1}$.

Figure 3 compares the model to a 3D simulation of a transversely flat-top
electron bunch with $r_{b0}\approx 0.27~{}k_{p}^{-1}$ and an initial centroid
$y_{c0}(\zeta)=0.027\sin(\zeta)$ propagating in constant-density plasma [18].
For the range in $\zeta$ we are considering, the truncated solutions Eq. (5)
and (6) are valid up to propagation distances $z\sim k_{\beta}^{-1}$ [see Fig.
3b)]. A larger number of series terms improves the agreement between the model
and simulation, as demonstrated in Fig. 3.

One could use the frequency response shown in Fig. 1a) to, for example, damp
the HI. In reality we may not control $k$, but we can control the local plasma
density and therefore the ratio of $k$ (fixed in the centroid perturbation) to
the local plasma wavenumber $k_{p}$, thereby operating in different growth
regimes (e.g. at $\Gamma<1$). Since $y_{c}$ and $v_{c}$ react differently to
each regime, the disruption of the bunch can only be minimised by alternating
between regimes. A mitigation set-up would therefore require propagating the
hosed or seeded bunch through plasma sections with different densities. The
series model can be applied to this set-up as long as we assume that the
plasma density does not change at the length scale of $k_{\beta}^{-1}$. In
that case we can normalize the entire model to the local plasma density.

Mitigation of a sinusoidal hosing seed at $k_{p}$ can be achieved by a series
of plasma density steps with densities above and below $n_{0}$. In fact, the
use of a density step to control an instability has been proposed before (e.g.
[23]). Note that, for an arbitrarily shaped initial centroid, the fastest
growing mode is sinusoidal, which means that this strategy can suppress the
growing component of any seed. The two configurations chosen to demonstrate
the concept are illustrated in Fig. 4 (inset). The parameters for these
density steps were found by experimenting with the series model.

We performed three 3D simulations (same parameters as before [18]): two where
a bunch seeded with $y_{c0}(\zeta)$ propagates through each of the two density
profiles shown in Fig. 4 (inset) and in resonant plasma ($n_{0}$) afterwards,
and one where the same bunch propagates exclusively in resonant plasma.

Figure 4:  Transverse energy of the bunch relative to the initial value along
the plasma (integrated over the simulation domain): plasma at $n_{0}$ (grey),
with two density steps (red), and with four density steps (blue). Results from
3D OSIRIS. Inset: plasma density profiles of the two density step
configurations. The vertical dash-dotted lines indicate the boundaries of each
density step.

To quantify the effect of the density steps on the bunch and its hosing seed,
we define a transverse energy
$\mathcal{E}(\zeta,z)=\frac{1}{2}\>v_{c}^{2}(\zeta,z)+\frac{1}{2}\>k_{\mathrm{HO}}^{2}(\zeta,z)\>y_{c}^{2}(\zeta,z)$,
in analogy with the HO. We note that the energy conservation equation for
$\mathcal{E}(\zeta,z)$ can be obtained by multiplying Eq. (2) by $v_{c}$. The
wavenumber $k_{\mathrm{HO}}^{2}(\zeta,z)$ is computed directly from the
$z$-evolving simulation data. In this geometry,
$k_{\mathrm{HO}}^{2}(\zeta,z)=2\>\hat{k}_{\beta}^{2}\;I_{1}[r_{b}(\zeta,z)]/r_{b}(\zeta,z)\int_{\zeta}^{\infty}d\zeta^{\prime}\sin(\zeta-\zeta^{\prime})\>f(\zeta^{\prime},z)\;K_{1}[r_{b}(\zeta^{\prime},z)]\;r_{b0}^{2}(\zeta^{\prime},z)/r_{b}(\zeta^{\prime},z)$,
where $f(\zeta^{\prime},z)$ is measured on axis.

Figure 4 compares $\mathcal{E}(\zeta,z)$ (integrated along $\zeta$ for a
complete picture of the bunch at every time step and normalised to the initial
value $\mathcal{E}_{0}$) for the three simulation cases. After two density
steps (red curves), the energy is smaller than in the resonant case by almost
two orders of magnitude (and lower than $\mathcal{E}_{0}$). It is possible to
stack further pairs of density steps with tuned parameters so as to extend the
suppressive effect on the HI, as demonstrated by the blue curves in Fig. 4.
After the steps, however, the growth rate tends to be exacerbated and
saturation is reached at similar energy values. Further work is necessary to
explore different density profile configurations and their ability to control
the growth and saturation of the HI.

In conclusion, we show that the growth rate of hosing in the linear regime
depends on the frequency of the centroid perturbation, and that the transient
growth spectrum at the onset of the instability can be harnessed to suppress
its growth. We demonstrate how this could be achieved by detuning the
background plasma frequency, i.e., by introducing a series of plasma density
steps. These findings apply both to long particle bunches and laser pulses,
and can be generalized for the control of both the hosing and self-modulation
instability growth rates.

###### Acknowledgements.

M.M. thanks John Farmer for fruitful discussions and acknowledges the
supervision of Bernhard Holzer. This work was partially supported by the
Portuguese Foundation for Science and Technology (FCT) through grants
CERN/FIS-TEC/0032/2017, CERN/FIS-TEC/0017/2019, IF/01635/2015 and
UID/FIS/50010/2023; and by the European Union’s Framework Programme for
Research and Innovation Horizon 2020 (2014–2020) under Grant Agreement No.
730871. We acknowledge PRACE for awarding us access to Piz Daint at CSCS,
Switzerland and MareNostrum at Barcelona Supercomputing Center (BSC), Spain,
where the simulations were performed.

## References

  * Akhiezer and Polovin [1956] A. Akhiezer and R. Polovin, Theory of wave motion of an electron plasma, J. Exp. Theor. Phys. 3, 696 (1956).
  * Dawson [1959] J. M. Dawson, Nonlinear electron oscillations in a cold plasma, Phys. Rev. 113, 383 (1959).
  * Steinke  _et al._ [2016] S. Steinke, J. van Tilborg, C. Benedetti, C. G. R. Geddes, J. Daniels, K. K. Swanson, A. J. Gonsalves, K. Nakamura, B. H. Shaw, C. B. Schroeder, E. Esarey, and W. P. Leemans, Staging of laser-plasma accelerators, Phys. Plasmas 23, 056705 (2016).
  * Lindstrøm [2021] C. A. Lindstrøm, Staging of plasma-wakefield accelerators, Phys. Rev. Accel. Beams 24, 014801 (2021).
  * AWAKE Collaboration [2018] AWAKE Collaboration, Acceleration of electrons in the plasma wakefield of a proton bunch, Nature 561, 363 (2018).
  * Kumar  _et al._ [2010] N. Kumar, A. Pukhov, and K. Lotov, Self-modulation instability of a long proton bunch in plasmas, Phys. Rev. Lett. 104, 255003 (2010).
  * Schroeder  _et al._ [2011] C. B. Schroeder, C. Benedetti, E. Esarey, F. J. Grüner, and W. P. Leemans, Growth and phase velocity of self-modulated beam-driven plasma waves, Phys. Rev. Lett. 107, 145002 (2011).
  * Batsch  _et al._ [2021] F. Batsch, P. Muggli, _et al._ (AWAKE Collaboration), Transition between instability and seeded self-modulation of a relativistic particle bunch in plasma, Phys. Rev. Lett. 126, 164802 (2021).
  * Verra  _et al._ [2022] L. Verra, G. Zevi Della Porta, J. Pucek, T. Nechaeva, S. Wyler, M. Bergamaschi, E. Senes, E. Guran, J. T. Moody, M. A. Kedves, E. Gschwendtner, P. Muggli, _et al._ (AWAKE Collaboration), Controlled growth of the self-modulation of a relativistic proton bunch in plasma, Phys. Rev. Lett. 129, 024802 (2022).
  * Whittum  _et al._ [1991] D. H. Whittum, W. M. Sharp, S. S. Yu, M. Lampe, and G. Joyce, Electron-hose instability in the ion-focused regime, Phys. Rev. Lett. 67, 991 (1991).
  * Schroeder  _et al._ [2012] C. B. Schroeder, C. Benedetti, E. Esarey, F. J. Grüner, and W. P. Leemans, Coupled beam hose and self-modulation instabilities in overdense plasma, Phys. Rev. E 86, 026402 (2012).
  * Schroeder  _et al._ [2013] C. B. Schroeder, C. Benedetti, E. Esarey, F. J. Grüner, and W. P. Leemans, Coherent seeding of self-modulated plasma wakefield accelerators, Phys. Plasmas 20, 056704 (2013).
  * Vieira  _et al._ [2014] J. Vieira, W. B. Mori, and P. Muggli, Hosing instability suppression in self-modulated plasma wakefields, Phys. Rev. Lett. 112, 205001 (2014).
  * P. Muggli for the AWAKE Collaboration [2020] P. Muggli for the AWAKE Collaboration, Physics to plan AWAKE Run 2, in  _J. Phys.: Conf. Ser._, Vol. 1596 (IOP Publishing, 2020) p. 012008.
  * Gschwendtner [2021] E. Gschwendtner, Awake Run 2 at CERN, in  _Proc. IPAC’21_, 12th Int. Particle Acc. Conf (JACoW Publishing, 2021) pp. 1757–1760.
  * Mehrling  _et al._ [2017] T. J. Mehrling, R. A. Fonseca, A. Martinez de la Ossa, and J. Vieira, Mitigation of the hose instability in plasma-wakefield accelerators, Phys. Rev. Lett. 118, 174801 (2017).
  * Lehe  _et al._ [2017] R. Lehe, C. B. Schroeder, J.-L. Vay, E. Esarey, and W. P. Leemans, Saturation of the hosing instability in quasilinear plasma accelerators, Phys. Rev. Lett. 119, 244801 (2017).
  * [18] See Supplemental Material for a detailed derivation of the average transverse force for a Gaussian transverse profile in 2D Cartesian geometry, the parameters of the particle-in-cell simulations, and a description of the phase shift measurement method.
  * Fonseca  _et al._ [2002] R. A. Fonseca, L. O. Silva, F. S. Tsung, V. K. Decyk, W. Lu, C. Ren, W. B. Mori, S. Deng, S. Lee, T. Katsouleas, and J. C. Adam, Osiris: A three-dimensional, fully relativistic particle in cell code for modeling plasma based accelerators, in  _Computational Science — ICCS 2002_ , edited by P. M. A. Sloot, A. G. Hoekstra, C. J. K. Tan, and J. J. Dongarra (Springer Berlin Heidelberg, 2002) pp. 342–351.
  * Duda  _et al._ [1999] B. J. Duda, R. G. Hemker, K. C. Tzeng, and W. B. Mori, A long-wavelength hosing instability in laser-plasma interactions, Phys. Rev. Lett. 83, 1978 (1999).
  * Kaluza  _et al._ [2010] M. C. Kaluza, S. P. D. Mangles, A. G. R. Thomas, Z. Najmudin, A. E. Dangor, C. D. Murphy, J. L. Collier, E. J. Divall, P. S. Foster, C. J. Hooker, A. J. Langley, J. Smith, and K. Krushelnick, Observation of a long-wavelength hosing modulation of a high-intensity laser pulse in underdense plasma, Phys. Rev. Lett. 105, 095003 (2010).
  * Chao  _et al._ [1980] A. W. Chao, B. Richter, and C.-Y. Yao, Beam emittance growth caused by transverse deflecting fields in a linear accelerator, Nuclear Instruments and Methods 178, 1 (1980).
  * Lotov and Tuev [2021] K. V. Lotov and P. V. Tuev, Plasma wakefield acceleration beyond the dephasing limit with 400 GeV proton driver, Plasma Phys. Control. Fusion 63, 125027 (2021).


# Group SELFIES: A Robust Fragment-Based Molecular String Representation

Austin H. Cheng
University of Toronto Andy Cai11footnotemark: 1
University of Toronto
Santiago Miret
Intel Labs
Gustavo Malkomes
SigOpt
Mariano Phielipp
Intel Labs
Alán Aspuru-Guzik
University of Toronto
Equal contribution.Correspondence to<EMAIL_ADDRESS>
###### Abstract

We introduce Group SELFIES, a molecular string representation that leverages
group tokens to represent functional groups or entire substructures while
maintaining chemical robustness guarantees. Molecular string representations,
such as SMILES and SELFIES, serve as the basis for molecular generation and
optimization in chemical language models, deep generative models, and
evolutionary methods. While SMILES and SELFIES leverage atomic
representations, Group SELFIES builds on top of the chemical robustness
guarantees of SELFIES by enabling group tokens, thereby creating additional
flexibility to the representation. Moreover, the group tokens in Group SELFIES
can take advantage of inductive biases of molecular fragments that capture
meaningful chemical motifs. The advantages of capturing chemical motifs and
flexibility are demonstrated in our experiments, which show that Group SELFIES
improves distribution learning of common molecular datasets. Further
experiments also show that random sampling of Group SELFIES strings improves
the quality of generated molecules compared to regular SELFIES strings. Our
open-source implementation of Group SELFIES is available at
https://github.com/aspuru-guzik-group/group-selfies, which we hope will aid
future research in molecular generation and optimization.

##  1 Introduction

The discovery of functional molecules for drugs and energy materials is
crucial to tackling global challenges in public health and climate change.
Different types of generative models can suggest potential molecules to
synthesize and test, but the performance of the models and molecules heavily
relies on the underlying molecular representation. Several models generate
molecules represented as SMILES strings [1, 2, 3, 4, 5], but their generated
output can be invalid due to syntax errors or incorrect valency. SELFIES [6]
is a molecular string representation that overcomes chemical invalidity
challenges by ensuring that any string of SELFIES characters can be decoded to
a molecule with valid valency. This not only makes it a natural representation
for chemical language models that output molecular strings, but also for
genetic algorithms such as GA+D, STONED, and JANUS [7, 8, 9] for molecular
optimization.

SELFIES improves string-based molecular generation by encoding prior knowledge
of valency constraints into the representation independently of the
optimization method. The representation has been shown to improve distribution
learning by language models [10], as well as image2string and string2string
translation [11, 12] and molecular generation in data-sparse regimes [13].
Additionally, simple add/replace/delete edits to SELFIES strings can generate
new but similar molecules, enabling genetic algorithms that directly
manipulate strings to generate molecules [7]. Alternatively, guiding these
simple string edits with Tanimoto similarity can interpolate between molecules
as performed in STONED by Nigam et al. [8], which can then be applied as
crossover operations in genetic algorithms such as JANUS [9]. Molecular
interpolation has also been used to find counterfactual decision boundaries
that explain a molecular classifier’s decisions [14].

While SMILES and SELFIES represent molecules at the individual atom and bond
level, human chemists typically think about molecules in terms of the
substructures that they contain. Human chemists can distinguish molecular
substructures based solely on the image of a molecule and induce the molecular
properties those substructures usually imply. Many fragment-based generative
models take advantage of this inductive bias [15, 16, 17, 18, 19, 20, 21, 22,
23, 24] by constructing custom representations amenable to fragment-based
molecular design. However, these approaches are not string-based, thereby
losing desirable properties of string representations: easy manipulation, and
direct input into established language models.

Similar to how SELFIES incorporates prior knowledge of valency constraints, we
can also incorporate prior knowledge in the form of functional groups and
molecular substructures into the representation. In this work, we combine the
flexibility of string representations with the chemical robustness of SELFIES
and the interpretability and inductive bias of fragment-based approaches into
a novel string representation: Group SELFIES, a robust string representation
that extends SELFIES to include tokens which represent functional groups or
entire substructures.

Figure 1: Visual overview of SELFIES and Group SELFIES. SELFIES is robust, so
shuffling tokens around will yield new molecules with correct valency. Group
SELFIES maintains robustness while adding group tokens, highlighted in color.
When Group SELFIES tokens are shuffled, structures like benzene rings are more
often preserved, while shuffled SELFIES strings rarely ever preserve
structures. Incidentally, Group SELFIES also improves the readability of
molecular string representations since chemists can see what substructures are
present.

In Section 2, we discuss how Group SELFIES fits into related research and then
formally introduce the representation in Section 3. Specifically, we outline
how molecules are encoded into and decoded from Group SELFIES, and we show
that arbitrary Group SELFIES strings can be decoded to molecules with valid
valency. The representation enables users to easily specify their own groups
or extract fragments from a dataset, leveraging the wide area of
cheminformatics research available there. In Section 4, we find that Group
SELFIES is more compact than SMILES or SELFIES and improves distribution
learning. Additionally, we compare molecules generated via randomly sampling
SELFIES and Group SELFIES strings and find that Group SELFIES improves the
quality of generated molecules. Molecular generation via random sampling
provides greater emphasis on the representation itself by abstracting away the
complexities of the type of generative method used, which we leave to future
work as described in Section 5.

Overall, Group SELFIES provides the flexibility of _group representation_ ,
the ability to represent _extended chirality_ via chiral group tokens and
_chemical robustness_ as summarized in Table 1.

Representation | Robustness |  | Substructure
---
Control
| Extended
---
Chirality
| Distribution
---
Learning
SMILES | no | no | no | ~
SELFIES | yes | no | no | ~
Group SELFIES | yes | yes | yes | improved
Table 1:  Comparison of the capabilities of SMILES, SELFIES, and Group
SELFIES. Group SELFIES provides _group representation_ , representation of
_extended chirality_ and _chemical robustness_. Additionally as shown in
Section 4, Group SELFIES improves distribution learning compared to other
representations.

##  2 Related Work

Fragment-Based String Representations: Group SELFIES is not the first
fragment-based string representation that has been proposed. Historical string
representations, such as Wiswesser Line Notation (WLN) [25, 26, 27], Hayward
Notation [28], and Skolnik Notation [29], all predate SMILES and represent
molecules non-atomically. They use tokens that represent functional groups,
such as carboxyls or phenyls, as well as ring systems. WLN strings are usually
shorter and sometimes easier for trained humans to understand than SMILES, as
it is easier to recognize functional groups encoded as single characters than
functional groups encoded atomically. SYBYL Line Notation (SLN) [30] allows
for “macro atoms” which specify multiple atoms in a substructure. The
Hierarchical Editing Language for Macromolecules (HELM) [31] represents
complex biomolecules by declaring monomers and then connecting them in a
polymer line notation. Human-Readable SMILES [32] applies common abbreviations
for chemical substituents to process and compress SMILES strings into a more
human-readable format. SMILES Pair Encoding [33] breaks down SMILES strings by
tokenizing them in a data-driven way that recognizes common substructures.

Genetic Programming: A string representation of molecules such as SELFIES can
be thought of as a programming language where programs specify how to
construct molecules. Genetic programming [34] uses genetic algorithms to
design programs that fulfill desired constraints. In particular for linear
genetic programming [35], programs are represented as linear sequences of
atomic instructions, rather than as a tree of expressions. Linear sequences of
atomic instructions are easily mutated by changing any instruction in the
sequence, since any sequence of atomic instructions will still be a runnable
program. In this way, SELFIES and Group SELFIES can be thought of as domain-
specific languages for linear genetic programming.

Learned Grammars: Data-Efficient Graph Grammar Learning (DEG) [22] is a recent
approach for extracting useful formal graph grammars from small datasets of
molecules. In this context, a “useful grammar” means that molecules generated
by applying random applicable production rules usually have high scores. The
learned production rules of the graph grammar can be thought of as similar to
functional groups applied in Group SELFIES. Group SELFIES allows for
flexibility and fine control of substructures, which can extracted from any
procedure including DEG.

##  3 Representation

###  3.1 SELFIES Framework

Figure 2: Visual explanation of Group SELFIES encoding/decoding of celecoxib.
Top-left: molecular structure of celecoxib. Top-middle: the structure of
celecoxib colored by its groups and atoms, with arrows and attachment indices
indicating how the encoder and decoder navigate around the groups. Top-right:
index overload table in Group SELFIES, indicating how tokens are interpreted
as numbers. Bottom: Celecoxib represented in Group SELFIES. Tokens are colored
by the groups and atoms they act on. Index overloads are shown where
interpreted. Arrows indicate how branches return to their original
branchpoints and how rings form bonds with previously placed atoms (light-
blue).

Before introducing in Group SELFIES in greater detail, we summarize the
primary features of SELFIES and the reasons underlying its chemical
robustness. SELFIES is equipped with an encoder and a decoder. The encoder
takes in a molecule and converts it to a SELFIES string, and the decoder takes
in a SELFIES string and converts it to a molecule. To encode a molecule in
SELFIES, one traverses its molecular graph and outputs the processed traversal
as a string of SELFIES tokens. To decode a SELFIES string, one reads through
the string token-by-token, building the molecular graph along the way until
arriving at the finished graph. Since the encoding and decoding process alone
does not guarantee chemical robustness, the SELFIES decoder further includes
two important features:

  1. 1.

Each token in SELFIES is overloaded to ensure that it can be interpreted
sensibly in all contexts. For instance, all tokens in SELFIES can also be
interpreted as numbers, which is useful when expressing branch and ring
lengths.

  2. 2.

SELFIES keeps track of the available valency at each step in the decoding
process; if a bond would be formed that would exceed this valency, it changes
the bond order or ignores the bond. For instance, when decoding [C][O][=C],
adding [=C] would exceed the valency of [O], so SELFIES changes the bond order
and adds [C] instead.

By preserving these properties in the Group SELFIES decoder, we ensure
chemical robustness is preserved.

###  3.2 Basic Tokens in Group SELFIES

Group SELFIES strings consist of the following fundamental tokens:

  * •

[X] adds an atom with the atomic symbol X.

  * •

[Branch] creates a new branch off the current atom and saves the current atom
as a branchpoint to return to later, and is analogous to an opening
parenthesis ( in SMILES. [pop] exits the current branch, returning to the most
recent branchpoint, and is analogous to a closing parenthesis ) in SMILES.
Unlike in SMILES, however, [Branch] and [pop] tokens need not come in pairs,
which helps maintain robustness. This token is also different from the
[BranchX] tokens in SELFIES. Experiments in Appendix A.5 indicate this change
does not substantially affect the performance of Group SELFIES.

  * •

[RingX] indicates that a ring bond will be formed from the current atom. The
next X tokens immediately following [RingX] will be interpreted as a number
$N$, and we will count backwards $N$ atoms in placement order to determine the
target of the ring bond. For example, [Ring2] indicates that the next 2 tokens
will be interpreted as a 2-digit base-16 number $N$. Ring bonds are stored
until after all tokens have been read by the decoder; only then are ring bonds
placed, and only if it would not violate valency. Due to the addition of
groups, it is sometimes necessary to form ring bonds to atoms that are added
after the current atom (e.g. ring bonds within groups). To indicate this, we
use the [->] token before the number token to specify that we will count
forwards instead of backwards.

All tokens can be modified by adding =, #, \ or / to change the bond order or
stereochemistry of their parent bond (e.g. [#Branch] or [/C]). The parent bond
is the bond to the previous atom.

###  3.3 Groups

The primary difference between SELFIES and Group SELFIES is the addition of
groups. Each group is defined as a set of atoms and bonds representing the
molecular group with its attachment points, indicating how the group can
participate in bonding. Each attachment point has a specified maximum valency,
which allows us to continue tracking available valency while decoding. These
attachment points are labeled by attachment indices, and the encoder and
decoder will navigate around these attachment indices as described in Section
3.4.

Users must specify the groups they want to use using a dictionary that maps
group names to groups. This tells the encoder what groups to recognize, and
tells the decoder how to map group tokens to groups. We call this dictionary a
“group set”, and every group set defines its own distinct instance of Group
SELFIES. In particular, the decoder will not recognize a Group SELFIES string
that contains group tokens not present in the current group set.

To distinguish group tokens from other tokens, we include a : character at the
front of the token (e.g. [:4benzene]). All group tokens are of the form
[:S<group-name>], where S is the starting attachment index of the group, and
<group-name> is any alphanumeric string that does not start with a number.

Optionally, a priority value can be specified for each group, indicating the
priority with which the group should be recognized when encoding into Group
SELFIES. Priority affects the Group SELFIES encoder as described in Section
3.4. For each group, one can also specify its index overload value, which is
the value the group token takes when the decoder must interpret the token as a
number.

Figure 3: Representation of a possible “benzene” group and its use in a
corresponding group token. The *N notation represents an attachment point with
valency N. Each attachment point is labeled with its attachment index
0th,1st... The starting attachment point represented in the token is also
highlighted.

###  3.4 Encoding and Decoding

Encoding: To encode a molecule in Group SELFIES, the encoder first recognizes
and replaces substructure matches of groups from the molecule. By default, the
encoder iterates through the group set and recognizes the largest groups
first, but users can override this by specifying a priority for each group.
Setting a high priority value for a group indicates that it will be recognized
first when encoding into Group SELFIES, ensuring that other group replacements
will not overlap with this group. This encoding strategy iterates implies that
increasing the size of the group set will increase the running time of the
encoder linearly. We then traverse the graph similar to the encoding process
for SMILES and SELFIES, while also placing the correct tokens for tracking the
attachment indices of where the encoder entered and exited a group.

Decoding: When decoding Group SELFIES, the process is essentially the same as
regular SELFIES except when reading group tokens. When a group token is read
by the decoder, the group set dictionary determines the corresponding group.
Subsequently, all atoms of the group are placed and the main chain is
connected to the starting attachment point. The decoder selects the next
attachment point to branch off from by reading in the next token as a relative
index. By adding the current attachment index to a relative index modulo the
total number of attachment points in the group, the decoder selects the next
attachment point. From the specified attachment point, the decoder implicitly
branches off of the group and continues traversing until a [pop] token is
read. Once the branch is “popped”, the decoder returns to the group and can
navigate to the next attachment point using another relative index. If the
selected attachment point is occupied, then the next available attachment
point is used. If all attachment points have been used up, then the group
itself is immediately “popped”, returning to the most recent branchpoint
before the group was placed.

We verified the robustness of Group SELFIES by encoding and decoding 25M
molecules from the eMolecules database [36]. We provide a detailed example of
encoding/decoding the molecule celecoxib in Appendix A.1.

Group SELFIES manages chirality differently than SMILES and SELFIES. Rather
than use @-notation to specify tetrahedral chirality, all chiral centers must
be specified as groups. We provide an “essential set” of 23 groups which
encode all relevant chiral centers in the eMolecules database. Equipped with
this essential set, every molecule can be encoded-decoded while maintaining
chirality. It is also an option to not use the essential set, or only use a
subset of it, depending on what chiral centers are relevant to the problem at
hand. If a molecule has a chiral center not specified in the group set, then
encode-decode will not preserve chirality.

###  3.5 Determining Fragments

Group SELFIES has a built-in flexibility for assigning the set of fragments
that make up a group set. Hence, the construction of a useful group set often
remains an open design choice. Users can specify groups using a SMILES-like
syntax, which could be useful if one knows what groups are synthetically
available or are expected to be useful for their particular design task.
Fragments can also be obtained from several fragment libraries found in the
literature [37, 38, 39]. Generally, a useful set of groups will appear in many
molecules in the dataset and replace many atoms, with similar fragments merged
together to reduce redundancy.

In our experiments, we also tested various fragmentation algorithms that
extract fragments from a dataset, including a naïve technique that cleaves
side chains from rings and a method based on matched molecular pair analysis
[40]. Several other fragmentation algorithms from cheminformatics can be
readily applied [41, 42, 43, 44, 45, 22] and the essential features of Group
SELFIES outlined in Section 1 remain robust to different fragment discovery
methods.

##  4 Experiments

The experiments in the subsequent sections outline some of the advantages of
Group SELFIES compared to SMILES and regular SELFIES representations.
Concretely, we show that: (1) Group SELFIES is shorter and more compressible
than SMILES and SELFIES; (2) Group SELFIES preserves useful properties during
generation; (3) Group SELFIES improves distribution learning.

###  4.1 Compactness

Group SELFIES strings are typically shorter than their SMILES and SELFIES
equivalents when using a generic set of groups. In Figure 4, this generic set
was generated by taking a random selection of 10,000 molecules from ZINC-250k
[46] and fragmenting them into 30 useful groups using various algorithms (see
Determining Fragments). We then combined these 30 groups with the 23 groups of
the essential set. Figure 4 shows histograms of the lengths of
SMILES/SELFIES/Group SELFIES strings of the entire ZINC-250k dataset. Length
is the number of characters in SMILES strings, and the number of tokens in
(Group) SELFIES strings. Group SELFIES strings are usually shorter than their
SELFIES and SMILES counterparts because group tokens can represent multiple
atoms in a molecule.

Since Group SELFIES has a larger alphabet than SMILES or SELFIES, we estimate
the complexity of each representation with the compressed filesize of
ZINC-250k. We find that out of all representations, Group SELFIES can be
compressed the most (see Appendix A.2).

Figure 4: Histogram of lengths of SMILES, SELFIES, and Group SELFIES strings
of the ZINC-250k dataset. Here, Group SELFIES uses a group set of 53 groups.

###  4.2 Molecular Generation

Algorithm 1 Sampling random (Group) SELFIES strings

1: given molecule dataset $D$, and encoder, decoder from SELFIES or Group
SELFIES

2: $S\leftarrow\\{$encoder$($molecule$)$ for molecule $\in D\\}$ # (Group)
SELFIES strings

3: $B\leftarrow$ concatenate$\\{$tokenize$(s)$ for $s\in S\\}$ # bag of tokens

4: repeat

5: select random string $s\in S$

6: length $l\leftarrow$ length$(s)$

7: tokens $\leftarrow$ sample $l$ tokens from $B$

8: string $\leftarrow$ concatenate$($tokens$)$

9: output decoder$($string$)$ # generated molecule

10: until enough molecules are generated

To specifically compare the suitability of the SELFIES and Group SELFIES
representations for molecular generation, we use a primitive generative model
which samples random strings. First, we convert a subset of $N=$ 100,000
molecules from ZINC-250k into (Group) SELFIES strings. Then, we tokenize all
strings and combine them into a single bag of tokens. To generate a new
string, we first pick a random (Group) SELFIES string from our chosen subset
and take its length $l$. We then randomly sample $l$ tokens from the bag, and
concatenate into a generated string. We generate $N$ random strings for each
representation. For Group SELFIES, we use the same 53 groups used for the
length histogram in Section 4.1.

We show histograms of the SAScore [47] and QED [48] of molecules generated
from ZINC in Figure 5. The distributions of generated Group SELFIES more
closely overlap with the original ZINC dataset than the generated SELFIES,
showing that even with an extremely simplistic generative model, Group SELFIES
can preserve important structural information. We perform a similar analysis
for a dataset of nonfullerene acceptors (NFA) [49] in Appendix 10 and find
that Group SELFIES preserves many aromatic rings in contrast to SELFIES, which
rarely ever preserves aromatic rings.

Figure 5: Molecules generated by our primitive generative model are binned by
SAScore and QED. For both properties, generated Group SELFIES have greater
overlap with the original ZINC distribution. Bracketed values indicate the
Wasserstein distance (a measure of overlap) to the ZINC distribution. Dashed
lines indicate the means.

###  4.3 Distribution Learning

To further quantify the effectiveness of Group SELFIES in generative models,
we use the MOSES benchmarking framework [50] to evaluate variational
autoencoders (VAEs) trained with both Group SELFIES and SELFIES strings.
Models were trained for 125 epochs. The group set for the Group SELFIES VAE
was created by fragmenting the training set provided by MOSES and selecting
the 300 most diverse groups. A set of 100,000 molecules was then generated
from each model and evaluated on the metrics provided by MOSES.

Model | Valid ($\uparrow$) | Unique@1k ($\uparrow$) | Unique@10k ($\uparrow$) | FCD ($\downarrow$) | SNN ($\uparrow$)
---|---|---|---|---|---
Test | TestSF | Test | TestSF
Train | 1.0 | 1.0 | 1.0 | 0.008 | 0.4755 | 0.6419 | 0.5859
Group-VAE-125 | 1.0(0) | 1.0(0) | 0.9985(4) | 0.1787(29) | 0.734(109) | 0.6051(4) | 0.5599(3)
SELFIES-VAE-125 | 1.0(0) | 0.9996(5) | 0.9986(4) | 0.6351(43) | 1.3136(128) | 0.6014(3) | 0.5566(2)

Model | Frag ($\uparrow$) | Scaf ($\uparrow$) | IntDiv ($\uparrow$) | IntDiv2 ($\uparrow$) | Filters ($\uparrow$) | Novelty ($\uparrow$)
---|---|---|---|---|---|---
Test | TestSF | Test | TestSF
Train | 1.0 | 0.9986 | 0.9907 | 0.0 | 0.8567 | 0.8508 | 1.0 | 1.0
Group-VAE-125 | 0.9995(0) | 0.9977(1) | 0.9649(21) | 0.0608(65) | 0.8587(1) | 0.8528(1) | 0.9623(7) | 0.7187(11)
SELFIES-VAE-125 | 0.9989(0) | 0.9965(1) | 0.9588(15) | 0.0675(37) | 0.8579(1) | 0.8519(1) | 0.96(4) | 0.7345(16)

Table 2: Group SELFIES VAE and SELFIES VAE evaluated on MOSES metrics. The
Group SELFIES VAE mostly matches or outperforms the SELFIES VAE.

For most metrics, Group SELFIES performs approximately the same as SELFIES.
Validity is the percentage of generated molecules that are accepted by RDKit’s
parser. Uniqueness is the percentage of generated molecules that are not
identical to any other generated molecule. Similarity to nearest neighbor
(SNN) is the average Tanimoto similarity between generated molecules and the
nearest neighbor in the reference set. Fragment similarity (Frag) is a cosine
similarity based on the distribution of BRICS fragments [45] of generated and
reference molecules. Scaffold similarity (Scaf) is a cosine similarity based
on the distribution of Bemis-Murcko scaffolds [51] of generated and reference
molecules. Internal diversity (IntDiv) measures the chemical diversity of the
generated molecules using Tanimoto similarity. Filters is the fraction of
generated molecules that pass filters for unwanted fragments. Novelty is the
fraction of generated molecules not in the training set.

The Group SELFIES model performs especially well on the Fréchet ChemNet
Distance (FCD) metric [52] when compared to SELFIES. FCD measures the
difference between the activations of the penultimate layer of ChemNet (a
model trained to predict the bioactivity of molecules) in the validation set
and in the generated set. Due to how ChemNet was trained, the activations are
likely to encode a mixture of biological and chemical properties important to
drug likelihood. This makes comparing these activations more informative than
comparing standard properties like logP or molecular weight, where the
correlation to bioactivity is weaker and less deliberate. To visualize FCD,
some indices of the penultimate activations of ChemNet are graphed in Figure
6. Generated Group SELFIES match these distributions more closely than
generated SELFIES.

Figure 6: Distribution of some values from the second-to-last layer of ChemNet
for molecules generated by Group SELFIES and SELFIES compared to the
validation set. The difference in distributions is used to calculate FCD.
Bracketed values in the legend represent the Wasserstein distance to the
original MOSES distribution.

##  5 Discussion

Our experiments show that Group SELFIES has noticeable advantages compared to
SMILES and SELFIES representations, including greater readability provided by
the group tokens. With regards to SMILES, the primary advantage is chemical
robustness. The comparison with SELFIES is more nuanced, as discussed in the
section below.

###  5.1 Group SELFIES vs SELFIES

Substructure Control: Group SELFIES provides more fine-grained control of
substructures, which creates the following advantages: (1) An important
scaffold can be preserved during optimization; (2) Chiral and charged groups
can be preserved during optimization, ensuring that charged tokens do not
proliferate and create radicals; (3) Synthetically accessible building blocks
can be chosen as groups to improve synthesizability.

Substructure Control with SELFIES: Various techniques applied to SELFIES can
mitigate the challenges of preserving structure. One such example is to simply
combine substrings of SELFIES strings together. Indeed, further experiments in
Appendix A.4 show that simply replacing all group tokens by their SELFIES
substrings shows similar performance to Group SELFIES. Within the SELFIES
framework, however, an inserted substring will not necessarily have that exact
substructure when decoded because the first token of the inserted substring
may need to be interpreted as a number, which can have cascading effects for
the rest of the substring. It is also likely that upon further insertions,
that substructure will not be preserved. Additionally, it is also not clear
how an insertion based approach can create groups with 3 or more branches,
since creating a third branch requires insertion in the middle of the group
substring.

Extended Chirality: Group SELFIES is theoretically capable of representing
molecules with extended chirality which are traditionally not able to be
represented with SMILES and SELFIES. These representations can only handle
local chirality - that is, chirality with a single atom as the chiral center.
This is in contrast to global chirality, where there may be an axis or plane
of chirality. Group SELFIES can handle global chirality by taking an entire
complex or chiral substructure and abstracting it into a group, leaving
attachment points on the outside for varying functionalization. Figure 7 shows
examples of groups with local and global chirality that Group SELFIES can
handle. We leave the proper implementation of representing global chirality to
future work.

Figure 7: Examples of chiral groups that can be represented in group SELFIES.
Tetrahedral carbon, hexoses, and octahedral metal centers have local
chirality, whereas allenes, metallocenes, helicenes, and substituted biphenyls
have global chirality.

Computational Speed: One tradeoff of Group SELFIES is that encoding and
decoding is usually slower than with SELFIES, likely due to overhead of RDKit
operations. The encoder is particularly slow as it relies on performing a
substructure match for every group in the group set. The decoder is faster
than the encoder, though still slower than the SELFIES decoder. See Appendix
A.6 for timing. To improve computational performance in future work, one could
exploit substructure control of Group SELFIES to reduce the number of
encode/decode calls needed to obtain high performers. Additionally, the speed
of encoding and decoding operations can be improved with distributed
computing, since Group SELFIES is trivially parallelizable for a fixed group
set.

###  5.2 Future Work

One promising extension of Group SELFIES is to incorporate more flexibility
into the representation. In such a case, a group token can represent an entire
scaffold, except without the atom identities. Other tokens can then identify
the atom types on the scaffolds. This would allow optimization of the atom
types while maintaining the topological structure of the scaffold. Another
current limitation of Group SELFIES is that groups cannot overlap; more work
is needed to develop a representation that acknowledges how groups might
overlap, particularly for generating polycyclic compounds. A sequence-based
representation of cellular complexes [53] or hypergraphs might suggest a
promising direction.

Finally, given this paper’s focus on the molecular representation, we only
applied rather simple generative modeling methods. We hope that future work
can leverage Group SELFIES to perform molecular generation with more advanced
generative methods, including chemical language models, deep generative
models, and evolutionary methods.

## Acknowledgments and Disclosure of Funding

We thank Luca Thiede, Akshat Nigam, Robert Pollice, Kjell Jorner, Gary Tom,
Nathanael Kusanda, Edwin Yu, Naruki Yoshikawa, and Felix Strieth-Kalthoff for
useful discussions.

Calculations testing the robustness of Group SELFIES were performed on the
Niagara supercomputer at the SciNet HPC Consortium. SciNet is funded by: the
Canada Foundation for Innovation; the Government of Ontario; Ontario Research
Fund - Research Excellence; and the University of Toronto.

This research was partially supported by funds from Intel Corporation.

A.A.-G. thanks Anders G. Frøseth for his generous support. A.A.-G. also
acknowledges the generous support of Natural Resources Canada and the Canada
150 Research Chairs program.

## References

  * Chithrananda et al. [2020] Seyone Chithrananda, Gabriel Grand, and Bharath Ramsundar.  Chemberta: large-scale self-supervised pretraining for molecular property prediction.  _arXiv preprint arXiv:2010.09885_ , 2020.
  * Gómez-Bombarelli et al. [2018] Rafael Gómez-Bombarelli, Jennifer N Wei, David Duvenaud, José Miguel Hernández-Lobato, Benjamín Sánchez-Lengeling, Dennis Sheberla, Jorge Aguilera-Iparraguirre, Timothy D Hirzel, Ryan P Adams, and Alán Aspuru-Guzik.  Automatic chemical design using a data-driven continuous representation of molecules.  _ACS central science_ , 4(2):268–276, 2018.
  * Blaschke et al. [2020] Thomas Blaschke, Josep Arús-Pous, Hongming Chen, Christian Margreitter, Christian Tyrchan, Ola Engkvist, Kostas Papadopoulos, and Atanas Patronov.  Reinvent 2.0: an ai tool for de novo drug design.  _Journal of chemical information and modeling_ , 60(12):5918–5922, 2020.
  * Moss et al. [2020] Henry Moss, David Leslie, Daniel Beck, Javier Gonzalez, and Paul Rayson.  Boss: Bayesian optimization over string spaces.  _Advances in neural information processing systems_ , 33:15476–15486, 2020.
  * Sanchez-Lengeling et al. [2017] Benjamin Sanchez-Lengeling, Carlos Outeiral, Gabriel L Guimaraes, and Alan Aspuru-Guzik.  Optimizing distributions over molecular space. an objective-reinforced generative adversarial network for inverse-design chemistry (organic).  2017\.
  * Krenn et al. [2020] Mario Krenn, Florian Häse, AkshatKumar Nigam, Pascal Friederich, and Alan Aspuru-Guzik.  Self-referencing embedded strings (selfies): A 100% robust molecular string representation.  _Machine Learning: Science and Technology_ , 1(4):045024, 2020.
  * Nigam et al. [2019] AkshatKumar Nigam, Pascal Friederich, Mario Krenn, and Alán Aspuru-Guzik.  Augmenting genetic algorithms with deep neural networks for exploring the chemical space.  _arXiv preprint arXiv:1909.11655_ , 2019.
  * Nigam et al. [2021] AkshatKumar Nigam, Robert Pollice, Mario Krenn, Gabriel dos Passos Gomes, and Alan Aspuru-Guzik.  Beyond generative models: superfast traversal, optimization, novelty, exploration and discovery (stoned) algorithm for molecules using selfies.  _Chemical science_ , 12(20):7079–7090, 2021.
  * Nigam et al. [2022] AkshatKumar Nigam, Robert Pollice, and Alán Aspuru-Guzik.  Parallel tempered genetic algorithm guided by deep neural networks for inverse molecular design.  _Digital Discovery_ , 2022.
  * Flam-Shepherd et al. [2022a] Daniel Flam-Shepherd, Kevin Zhu, and Alán Aspuru-Guzik.  Language models can learn complex molecular distributions.  _Nature Communications_ , 13(1):1–10, 2022a.
  * Rajan et al. [2020] Kohulan Rajan, Achim Zielesny, and Christoph Steinbeck.  Decimer: towards deep learning for chemical image recognition.  _Journal of Cheminformatics_ , 12(1):1–9, 2020\.
  * Rajan et al. [2021] Kohulan Rajan, Achim Zielesny, and Christoph Steinbeck.  Stout: Smiles to iupac names using neural machine translation.  _Journal of Cheminformatics_ , 13(1):1–14, 2021\.
  * Frey et al. [2022] Nathan C Frey, Vijay Gadepally, and Bharath Ramsundar.  Fastflows: Flow-based models for molecular graph generation.  _arXiv preprint arXiv:2201.12419_ , 2022.
  * Wellawatte et al. [2022] Geemi P Wellawatte, Aditi Seshadri, and Andrew D White.  Model agnostic generation of counterfactual explanations for molecules.  _Chemical science_ , 13(13):3697–3705, 2022.
  * Jensen [2019] Jan H Jensen.  A graph-based genetic algorithm and generative model/monte carlo tree search for the exploration of chemical space.  _Chemical science_ , 10(12):3567–3572, 2019.
  * Jin et al. [2018] Wengong Jin, Regina Barzilay, and Tommi Jaakkola.  Junction tree variational autoencoder for molecular graph generation.  In _International conference on machine learning_ , pages 2323–2332. PMLR, 2018.
  * Xie et al. [2021] Yutong Xie, Chence Shi, Hao Zhou, Yuwei Yang, Weinan Zhang, Yong Yu, and Lei Li.  Mars: Markov molecular sampling for multi-objective drug discovery.  _arXiv preprint arXiv:2103.10432_ , 2021.
  * Bengio et al. [2021] Emmanuel Bengio, Moksh Jain, Maksym Korablyov, Doina Precup, and Yoshua Bengio.  Flow network based generative models for non-iterative diverse candidate generation.  _Advances in Neural Information Processing Systems_ , 34:27381–27394, 2021.
  * Jin et al. [2020] Wengong Jin, Regina Barzilay, and Tommi Jaakkola.  Multi-objective molecule generation using interpretable substructures.  In _International conference on machine learning_ , pages 4849–4859. PMLR, 2020.
  * Yang et al. [2021] Soojung Yang, Doyeong Hwang, Seul Lee, Seongok Ryu, and Sung Ju Hwang.  Hit and lead discovery with explorative rl and fragment-based molecule generation.  _Advances in Neural Information Processing Systems_ , 34:7924–7936, 2021.
  * Flam-Shepherd et al. [2022b] Daniel Flam-Shepherd, Alexander Zhigalin, and Alán Aspuru-Guzik.  Scalable fragment-based 3d molecular design with reinforcement learning.  _arXiv preprint arXiv:2202.00658_ , 2022b.
  * Guo et al. [2022] Minghao Guo, Veronika Thost, Beichen Li, Payel Das, Jie Chen, and Wojciech Matusik.  Data-efficient graph grammar learning for molecular generation.  _arXiv preprint arXiv:2203.08031_ , 2022.
  * Liu et al. [2018] Qi Liu, Miltiadis Allamanis, Marc Brockschmidt, and Alexander Gaunt.  Constrained graph variational autoencoders for molecule design.  _Advances in neural information processing systems_ , 31, 2018.
  * Polishchuk [2020] Pavel Polishchuk.  Crem: chemically reasonable mutations framework for structure generation.  _Journal of Cheminformatics_ , 12(1):1–18, 2020\.
  * Wiswesser [1951] WJ Wiswesser.  Simplified chemical coding for automatic sorting and printing machinery.  _Reading, Pa., Willson Products Inc_ , 1951.
  * Wiswesser [1968] William J Wiswesser.  107 years of line-formula notations (1861-1968).  _Journal of Chemical Documentation_ , 8(3):146–150, 1968.
  * Vollmer [1983] John J Vollmer.  Wiswesser line notation: an introduction.  _Journal of Chemical Education_ , 60(3):192, 1983\.
  * Hayward [1961] H Winston Hayward.  _A new sequential enumeration and line formula notation system for organic compounds_.  Number 21. Office of Research and Development, Patent Office, 1961.
  * Skolnik and Clow [1964] Herman Skolnik and Avis Clow.  A notation system for indexing pesticides.  _Journal of Chemical Documentation_ , 4(4):221–227, 1964.
  * Homer et al. [2008] R Webster Homer, Jon Swanson, Robert J Jilek, Tad Hurst, and Robert D Clark.  Sybyl line notation (sln): a single notation to represent chemical structures, queries, reactions, and virtual libraries.  _Journal of chemical information and modeling_ , 48(12):2294–2307, 2008.
  * Zhang et al. [2012] Tianhong Zhang, Hongli Li, Hualin Xi, Robert V Stanton, and Sergio H Rotstein.  Helm: a hierarchical notation language for complex biomolecule structure representation, 2012.
  * Garay-Ruiz et al. [2021] Diego Garay-Ruiz, Carles Bo, and Diego Garay Ruiz.  Human-readable smiles: Translating cheminformatics to chemistry.  2021\.
  * Li and Fourches [2021] Xinhao Li and Denis Fourches.  Smiles pair encoding: a data-driven substructure tokenization algorithm for deep learning.  _Journal of Chemical Information and Modeling_ , 61(4):1560–1569, 2021.
  * Koza and Poli [2005] John R Koza and Riccardo Poli.  Genetic programming.  In _Search methodologies_ , pages 127–164. Springer, 2005.
  * Brameier et al. [2007] Markus Brameier, Wolfgang Banzhaf, and Wolfgang Banzhaf.  _Linear genetic programming_ , volume 1.  Springer, 2007.
  * [36] Products: Structure data downloads.  URL https://search.emolecules.com/info/products-data-downloads.html.
  * Zdrazil and Guha [2017] Barbara Zdrazil and Rajarshi Guha.  The rise and fall of a scaffold: a trend analysis of scaffolds in the medicinal chemistry literature.  _Journal of Medicinal Chemistry_ , 61(11):4688–4703, 2017.
  * Ertl and Schuhmann [2019] Peter Ertl and Tim Schuhmann.  A systematic cheminformatics analysis of functional groups occurring in natural products.  _Journal of natural products_ , 82(5):1258–1263, 2019.
  * Sharif et al. [2022] Suliman Sharif, Ruibin Liu, Asuka A. Orr, Daniel Khavrutskii, Sunhwan Jo, Bettina Lier, Anastasia Croitoru, Chris Burke, Aziza Frank, Jacob Weiner, Nathaniel McClean, Aarion Romany, Mingtian Zhao, Takayuki Serizawa, Jared Deacon, Ian Jones, Shaoqi Zhan, Anmol Kumar, Mike Woster, Rebecca Pinette-Dorin, Elena Y. Chow, Sevien Schulhoff, and Alexander D. MacKerell Jr.  Global-Chem: A Chemical Knowledge Graph of common small molecules and their IUPAC/SMILES/SMARTS for selection of compounds relevant to diverse chemical communities, 8 2022.
  * Hussain and Rea [2010] Jameed Hussain and Ceara Rea.  Computationally efficient algorithm to identify matched molecular pairs (mmps) in large data sets.  _Journal of chemical information and modeling_ , 50(3):339–348, 2010.
  * Sheng and Zhang [2013] Chunquan Sheng and Wannian Zhang.  Fragment informatics and computational fragment-based drug design: an overview and update.  _Medicinal Research Reviews_ , 33(3):554–598, 2013.
  * Liu et al. [2017] Tairan Liu, Misagh Naderi, Chris Alvin, Supratik Mukhopadhyay, and Michal Brylinski.  Break down in order to build up: decomposing small molecules for fragment-based drug design with e molfrag.  _Journal of chemical information and modeling_ , 57(4):627–631, 2017.
  * Kutchukian et al. [2015] Peter S Kutchukian, Sung-Sau So, Christian Fischer, and Chris L Waller.  Fragment library design: Using cheminformatics and expert chemists to fill gaps in existing fragment libraries.  In _Fragment-Based Methods in Drug Discovery_ , pages 43–53. Springer, 2015.
  * Müller [2019] Simon Müller.  Flexible heuristic algorithm for automatic molecule fragmentation: application to the unifac group contribution model.  _Journal of cheminformatics_ , 11(1):1–12, 2019\.
  * Degen et al. [2008] Jörg Degen, Christof Wegscheid-Gerlach, Andrea Zaliani, and Matthias Rarey.  On the art of compiling and using’drug-like’chemical fragment spaces.  _ChemMedChem: Chemistry Enabling Drug Discovery_ , 3(10):1503–1507, 2008.
  * Irwin et al. [2020] John J Irwin, Khanh G Tang, Jennifer Young, Chinzorig Dandarchuluun, Benjamin R Wong, Munkhzul Khurelbaatar, Yurii S Moroz, John Mayfield, and Roger A Sayle.  Zinc20—a free ultralarge-scale chemical database for ligand discovery.  _Journal of chemical information and modeling_ , 60(12):6065–6073, 2020.
  * Ertl and Schuffenhauer [2009] Peter Ertl and Ansgar Schuffenhauer.  Estimation of synthetic accessibility score of drug-like molecules based on molecular complexity and fragment contributions.  _Journal of cheminformatics_ , 1(1):1–11, 2009\.
  * Bickerton et al. [2012] G Richard Bickerton, Gaia V Paolini, Jérémy Besnard, Sorel Muresan, and Andrew L Hopkins.  Quantifying the chemical beauty of drugs.  _Nature chemistry_ , 4(2):90–98, 2012.
  * Lopez et al. [2017] Steven A Lopez, Benjamin Sanchez-Lengeling, Julio de Goes Soares, and Alán Aspuru-Guzik.  Design principles and top non-fullerene acceptor candidates for organic photovoltaics.  _Joule_ , 1(4):857–870, 2017.
  * Polykovskiy et al. [2020] Daniil Polykovskiy, Alexander Zhebrak, Benjamin Sanchez-Lengeling, Sergey Golovanov, Oktai Tatanov, Stanislav Belyaev, Rauf Kurbanov, Aleksey Artamonov, Vladimir Aladinskiy, Mark Veselov, Artur Kadurin, Simon Johansson, Hongming Chen, Sergey Nikolenko, Alan Aspuru-Guzik, and Alex Zhavoronkov.  Molecular Sets (MOSES): A Benchmarking Platform for Molecular Generation Models.  _Frontiers in Pharmacology_ , 2020.
  * Murcko and Bemis [1996] Mark A Murcko and Guy W Bemis.  The properties of known drugs. 1.  _Molecular Frameworks. J. Med. Chem_ , 39:2887–2893, 1996\.
  * [52] Kristina Preuer, Philipp Renz, Thomas Unterthiner, Sepp Hochreiter, and Günter Klambauer.  Fréchet ChemNet distance: A metric for generative models for molecules in drug discovery.  58(9):1736–1741.  ISSN 1549-9596.  doi: 10.1021/acs.jcim.8b00234.  URL https://doi.org/10.1021/acs.jcim.8b00234.  Publisher: American Chemical Society.
  * Hajij et al. [2020] Mustafa Hajij, Kyle Istvan, and Ghada Zamzmi.  Cell complex neural networks.  _arXiv preprint arXiv:2010.00743_ , 2020.

##  Appendix A Appendix

###  A.1 Example

As an example, we will go through encoding and decoding the molecule
celecoxib.

Figure 8: Celecoxib

We will define our group set as consisting of 4 main groups: trifluoromethane,
toluene, benzene, and sulfonamide.

Figure 9: Group set used in this example.

We will now extract occurrences of the groups in the molecule. Since
priorities were not specified, we match groups by size in descending order.

We will now traverse the graph. The following diagrams demonstrate this
process.

We then start traversing at an arbitrary atom, in this case the methyl carbon
of the toluene.

We find that this atom belongs to a group, and so place the appropriate group
token. We arbitrarily select 0 to be the starting attachment point, since
there is no previous connection into this group.

We then place an index token of [+2] to go from the 0th attachment point to
the 2nd attachment point. For readability, index tokens are shown as numbers,
but in the actual encoder output a token with the appropriate overload value
would be used (e.g. in this case [Ring2] is overloaded to [+2]).

We again find that our current atom is in a group. This time, we place the
group token with the starting attachment point at the 2nd attachment index,
since that is where the last bond entered from.

We then use the [pop] token to exit the current group. In this case, we go
from the trifluoromethane group back to the pyrazole group.

The next few steps build out the carbon chain of the benzene using atomic
tokens, including placing a branch and a ring bond.

After the benzene is completed, the molecule is finished by placing a
sulfonamide group.

The decoding process for this Group SELFIES string will read in token-by-
token, building the molecular graph in the same order as shown in this
encoding example.

###  A.2 Complexity and compression

To take into account the alphabet size of SMILES, SELFIES, and Group SELFIES,
we compress each representation using zlib to estimate the complexity of each
representation. We first do index encoding to convert each string into a list
of ints by replacing each unique token with a unique int. We then compress
this representation using zlib.

Representation |  $\\#$ unique tokens | strings | +index encoding |  | +index encoding
---
+zlib compression
SMILES | 34 | 11.80 MB | 23.11 MB | 3.70 MB
SELFIES | 107 | 46.44 MB | 19.68 MB | 3.97 MB
Group SELFIES | 247 | 42.05 MB | 15.90 MB | 3.61 MB
Table 3: Filesize of ZINC-250k when represented in SMILES, SELFIES, and Group
SELFIES. When using strings, SMILES takes up the least space. When using index
encoding and then compressing using zlib, Group SELFIES takes up the least
space.

###  A.3 NFA experiments

Figure 10: Generated Group SELFIES and regular SELFIES binned by SAScore and
number of aromatic atoms. Molecules were generated from the NFA dataset.

We repeat Experiment 4.2 but for a dataset of nonfullerene acceptors (NFA)
[49]. This dataset contains many conjugated aromatic systems. In this case, we
histogram by SAScore and the number of aromatic atoms, as shown in Figure 10.
Generated SELFIES are rarely able to preserve the aromatic systems found in
NFA, but generated Group SELFIES preserve much more of the aromatic systems,
and can produce molecules with lower and higher SAScore. When binning by
number of aromatic atoms, molecules with 0 aromatic atoms are omitted for
clarity. About 49% / 98% of molecules generated via Group SELFIES / SELFIES
have 0 aromatic atoms. Group SELFIES uses a group set containing 74 groups
obtained via naïve fragmentation. $N=$ 51281 molecules are generated via Group
SELFIES / SELFIES, equal to the size of the NFA dataset.

###  A.4 Substring SELFIES

Substring SELFIES is generated by taking a Group SELFIES string and replacing
every group token with a SELFIES substring corresponding to it. We compare
Substring SELFIES and Group SELFIES in the context of Experiment 4.2 in Figure
11 and find that generated Substring SELFIES are on par with generated Group
SELFIES.

Figure 11: Generated Group SELFIES and Substring SELFIES binned by SAScore and
QED. Molecules were generated from ZINC-250k.

###  A.5 No-Group SELFIES

No-Group SELFIES is generated by using the Group SELFIES encoder with an empty
group set. We compare regular SELFIES to No-Group SELFIES to determine the
effect of small changes in the base encoding process, such as branch tokens.
In SELFIES, all [BranchX] expect a number token after to determine the length
of the branch, whereas in Group SELFIES, all [Branch] tokens do not need a
number token, instead using a [pop] token to exit branches. In Figure 12,
generated No-Group SELFIES are on par with or worse than generated SELFIES,
indicating that the modified encoding process does not contribute much to the
overall performance of Group SELFIES.

Figure 12: Generated regular SELFIES and No-Group SELFIES binned by SAScore
and QED. Molecules were generated from ZINC-250k.

###  A.6 Timing

Representation | Encode | Decode
---|---|---
SELFIES | 0.199 ms | 0.133 ms
Group SELFIES | 12.860 ms | 2.494 ms
Table 4: Running time of encoder/decoder for SELFIES and Group SELFIES. The
SELFIES encoder is 65x faster than the Group SELFIES encoder, and the SELFIES
decoder is 19x faster than the Group SELFIES decoder. Timing is averaged per
molecule over the ZINC-250k dataset, using the same group set of 53 groups as
the length histogram experiments. Computations were done on a 2021 Apple
MacBook Pro with a M1 Pro chip and 32 GB RAM.