# Learning dynamic representations of the functional connectome in
neurobiological networks

Luciano Dyballa1, Samuel Lang2, Alexandra Haslund-Gourley1,
Eviatar Yemini2∗, Steven W. Zucker1
1Dept. Computer Science, Yale University; 2Dept. Neurobiology, UMass Chan
Medical School
Correspondence to<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The static synaptic connectivity of neuronal circuits stands in direct
contrast to the dynamics of their function. As in changing community
interactions, different neurons can participate actively in various
combinations to effect behaviors at different times. We introduce an
unsupervised approach to learn the dynamic affinities between neurons in live,
behaving animals, and to reveal which communities form among neurons at
different times. The inference occurs in two major steps.111 Code is available
at https://github.com/dyballa/dynamic-connectomes. First, pairwise non-linear
affinities between neuronal traces from brain-wide calcium activity are
organized by non-negative tensor factorization (NTF). Each factor specifies
which groups of neurons are most likely interacting for an inferred interval
in time, and for which animals. Finally, a generative model that allows for
weighted community detection is applied to the functional motifs produced by
NTF to reveal a dynamic functional connectome. Since time codes the different
experimental variables (e.g., application of chemical stimuli), this provides
an atlas of neural motifs active during separate stages of an experiment
(e.g., stimulus application or spontaneous behaviors). Results from our
analysis are experimentally validated, confirming that our method is able to
robustly predict causal interactions between neurons to generate behavior.

##  1 Introduction

The connectome in neurobiology might seem roughly analogous to the
architecture of artificial neural networks (ANNs) in artificial intelligence
(AI), in the sense that it specifies the relevant structure of connections. In
AI, engineering this architecture is an engaging part of problem formulation;
in neurobiology, the structure of the connectome complicates the process of
circuit dissection for several reasons. First, there are different types of
connections between neurons: some are due to chemical synapses or electrical
ones (gap junctions) (White et al., 1986; Jarrell et al., 2012; Cook et al.,
2019; Witvliet et al., 2021), and others induced by, e.g., extrasynaptic
signaling of monoamines or neuropeptides (Bentley et al., 2016; Ripoll-Sánchez
et al., 2022; Beets et al., 2023). Each one defines a distinct connectome
(Fig. 1a). Secondly, unlike ANNs, in biology different connections can be
engaged at any given moment (e.g., as a result of neuromodulatory signals that
alter system-wide states (Bargmann & Marder, 2013)). As in complex social
networks, the effective—or functional—connectome varies with task and animal
state; these dynamics render the predictive power of the biological
connectomes problematic. We agree with others working in different domains,
e.g. Bassett et al. (2011); Skarding et al. (2021), that a functional
connectome that is dynamic across time is required, given that a number of
phenomena have been observed wherein neurons flexibly change how they encode
behavior in a state-dependent manner (Atanas et al., 2023). Our goal is to
learn a representation of these time-varying communities of neurons organized
by behavioral responses in a “simple” organism, Caenorhabditis elegans.

We present a novel algorithm for learning the dynamic community organization
within a weighted connectome, based on brain-wide activity measurements of
each individual neuron. Our main contribution is to organize the
connections—or estimated similarities between neural activity—rather than the
neurons themselves. It differs from previous attempts by considering the full
system of similarities and animals across time, rather than taking a step-by-
step statistical approach (Varshney et al., 2011; Bentley et al., 2016; Cook
et al., 2019; Uzel et al., 2022). Our approach combines two previously
unrelated components: (i) the use of tensor factorization to reveal putative
groups of related neurons/animals in time, which are then (ii) passed to a
community detection algorithm. Although the emphasis here is on the algorithm,
our results are confirmed with experiments that silence specific neurons
predicted by our method to directly measure their impact on behavior.
Moreover, even though the dynamic functional connectomes we reveal are
specific to C. elegans, we believe our approach is widely applicable to
learning representations of dynamic communities of neurons in other organisms,
as well as in other domains, such as fMRI and social networks.

Figure 1: From static to dynamic connectomes. (a) Different modes of neural
communication are represented by substantially different connectomes. Shown
here are chemical synapses (Syn), gap junctions (GJ), monoamines (MA), and
neuropeptide (NP) connectomes for the nematode worm C. elegans (image adapted
from Bentley et al. (2016)). How these relate to behavior has been
problematic, however. The immediate neighborhoods for neuron AWB are shown (b)
from Syn and (c) from GJ; which governs AWB communication and when? (d, top)
The dataset consists of activity traces from each neuron in C. elegans across
time; dashed lines illustrate stimulus presentation of e.g. a repulsive
chemical that is directly sensed by the worm. Note how the traces differ
across time and across worms. (d, bottom) Different modules, or communities of
neurons, become active at different time periods, revealing how neurons
interact dynamically to encode behavior. Our goal is to infer these
communities comprising the dynamic functional connectome.

##  2 Methods

It is now possible, using the “NeuroPAL” method (Yemini et al., 2021), to
determine both the cell-type-identity of every neuron and also record its
activity in C. elegans across time in awake, and even behaving animals (Dag et
al., 2023). Previously, neural activity traces from multiple worms have been
analyzed using traditional statistics such as correlation (e.g., Kato et al.
(2015); Yemini et al. (2021); Susoy et al. (2021); Randi et al. (2023)), small
local motifs (Uzel et al., 2022), or more complex predictive models (Cecere et
al., 2021; Dag et al., 2023) which assume that most neurons will respond in
roughly the same way when exposed to the same stimulus. However, both past and
recent studies indicate different ways that neurons can change how they encode
behavior in a state-dependent manner (Bargmann & Marder, 2013; Gordus et al.,
2015; Atanas et al., 2023), despite identical repeated presentations of
chemosensory stimuli intended to produce stereotypical behavior. Multiple
worms are therefore needed to determine whether “typical” functional networks
might activate after the introduction of a particular stimulus.

The third dimension (individual worms) in the data suggests using a tensor
rather than a matrix for representation of neural activity (Williams et al.,
2018; Dyballa et al., 2024). A naive tensor formulation would be to build it
using the temporal activity directly: worms $\times$ neurons $\times$ temporal
traces (e.g., Williams et al. (2018)). However, this would restrict the
similarity computation between traces to a multi-linear form, which can be
ambiguous for imaging data. We propose instead a novel tensor formulation with
dimensions: time $\times$ worms $\times$ pairwise affinities, where
instantaneous affinities (that may vary over time) are precomputed from the
temporal traces of each pair of neurons. In this way, non-linearities can be
introduced by performing the affinity computations between neurons separately
from the factorization (see below). Moreover, given that neural affinities are
expected to be state-dependent and changing over time, the tensor components
will be able to provide a summarized description of which communities or sub-
networks are active and when. We demonstrate this in Results.

###  2.1 Interpretable, instantaneous affinities; Fig. 3a

Experimental recordings of C. elegans currently span several minutes to nearly
half an hour in duration, during which it is likely that different groups of
neurons will become active at different times. Affinities should thus be a
function of time, rather than the more traditional method of computing a
global-time measure such as correlation (e.g., Yemini et al. (2021); Randi et
al. (2023)). We compute, instead, local non-linear similarities between pairs
of neural trace curves, which we term differential affinity.

Building a time-varying similarity measure between traces is a subtle problem.
For example, two neurons that are mostly silent throughout an experiment will
have small Euclidean distance (and high cosine similarity) between them, even
though there may be no evidence that the neurons participate in the same
circuit. Or, if they happen to maintain high, constant levels of activity at
some point in time, then their distances will also be locally small, but if
they reached those levels several seconds apart, it is unlikely that the two
neurons are actually interacting. The instantaneous trace values are not
meaningful individually: the temporal history of how they arrived at that
level is highly relevant. We rely on the information present in the rates of
monotonic change in activity, and compute what we term local differential
affinity.

We compare two neurons’ derivatives during intervals in which both had a
constant sign, i.e., periods of monotonic increase or decrease in activity.
The intuition is that two neurons with coinciding changes in activity are
likely to be influencing one another, or being influenced by a common
phenomenon, and thereby participating in the same functional circuit. The
details of this computation are given in the Appendix; examples of the
instantaneous local differential affinity are shown in Fig. 2.

Figure 2: Computation of dynamic differential affinity. (a, top) Activity
traces of two neurons, AWCL and ASER, across time. Note that they share marked
periods of similar trends in activity. (a, middle) Smoothed derivatives of the
top traces. Note how the derivatives agree during periods in which both traces
vary together. (a, bottom) Affinity trace across time computed between AWCL
and ASER (see section 2.1 and Appendix A.1 for details). Note how the “bumps”
coincide with regions of potential interest. (b) Similar plots for a different
neuron pair. Again, our non-linear affinity measure has low value except for
periods when both neurons increase or decrease their activities, i.e., are
more likely to be interacting.

Affinities show how similar two curves are, locally, in terms of their
absolute derivatives, which is motivated by the fact that two neurons with
very similar but opposite sign derivatives are still likely to be interacting
by means of some inhibitory mechanism. As a result, they can be interpreted as
how likely it is for the two neurons to be interacting (either directly or
indirectly, perhaps via common input). The instantaneous pairwise affinities
$a_{ij}^{(t)}$ at a time $t=1,\ldots,T$ may be bundled into an affinity
matrix, $A^{(t)}$, which is used in the tensor incorporating information
across multiple worms.

###  2.2 Dynamic affinity tensor decomposition; Fig. 3b

We now have a set of affinity matrices $\\{A^{(t)}\\}$, each one encoding an
affinity network at a time $t$, on which a community detection algorithm could
be run. This would yield a collection of “instantaneous” network communities.
However, which instants should be selected? How long should they be? Are they
independent? Rather than addressing these questions individually, we exploit
the summarization power of non-negative tensor factorization to automatically
cluster affinity networks across time (and animals!): we seek contiguous
intervals during which a given affinity pattern is approximately preserved
across several worms. Recall that our tensor has dimensions: time $\times$
worms $\times$ pairwise affinities. (We note as an aside that this tensor can
also help with “completing” affinity matrices containing missing data from a
few neurons.) See Fig. 3 for a graphical illustration of the following
discussion.

If there are $n$ neurons, $A(t)$ has dimension $n\times n$. Note that there is
no expectation for parts of the connectivity structure in $A(t)$ to be
decomposable as the product of two matrices. Therefore, instead of stacking
the 2-D matrices in a 4-way tensor, we vectorize each $A(t)$ by reshaping it
as a vector of dimension $n(n-1)/2$, the total number of pairs (for
computational efficiency, since they are symmetric, only the upper-triangular
part is used).

Figure 3: The procedure for inferring dynamic functional connectomes in three
stages. (a) Stage 1: Building affinities. Neurons in behaving worms are imaged
to yield a matrix of activity (in time) for each neuron for each worm. The
differential affinity computation compares traces to yield a tensor of
affinities (in time) for each neuron pair for each worm. (b) Stage 2: Non-
negative tensor factorization yields components that reveal affinity patterns
found for certain worms (bar plot) over different time intervals (temporal
curve). For example, component 2 indicates that neuron pairs with high
affinity in worms 4 and 6 are active early in the experiment, and the last
factor applies mainly to worm 2. (c) Stage 3: From affinities to a functional
connectome. The affinity factor from a selected tensor component is rearranged
as an affinity matrix, from which an equivalent weighted graph is implicitly
built. A community detection algorithm then reveals groupings of neurons
behaving similarly over the indicated time interval.

In a Canonical Polyadic (CP) decomposition (Carroll & Chang, 1970; Harshman et
al., 1970; Kolda & Bader, 2009), the goal is to approximate a $n$-way tensor
$\bm{T}\in\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times I_{n}}$ by a sum of
rank-1 tensors. Letting $R$ be the number of components chosen, we may express
the specific case of a 3-way tensor as (Kolda, 2006):

$\tilde{\bm{T}}=\sum_{r=1}^{R}{v^{(1)}_{r}\circ v^{(2)}_{r}\circ
v^{(3)}_{r}},$ (1)

where $\circ$ stands for the vector outer product and $X^{(k)}$ is called a
_factor matrix_ containing the factors $v^{(k)}_{r}$ as its columns. Each
rank-1 tensor is thus formed by the outer product between each factor in the
same component (a _component_ refers to each set of associated factors). Most
tensor decomposition algorithms use squared reconstruction error as objective
function (Cichocki et al., 2009):

$\min_{X^{(1)},X^{(2)},X^{(3)}}\frac{1}{2}\|\bm{T}-\tilde{\bm{T}}\|^{2}.$ (2)

where $\|\cdot\|$ is the norm of the vectorized tensor.

Since the data consist of non-negative neural affinities, we adopt a non-
negative tensor decomposition, also known as non-negative tensor factorization
(NTF) (Bro & De Jong, 1997), which adds a non-negativity constraint
$X^{(k)}\geq 0,\forall k$ to eq. 2. Non-negative factorization is popularly
used for matrices (NMF), and achieves a parts-based representation with more
easily-interpretable components (Lee & Seung, 1999; Cichocki et al., 2009);
additionally, it has been shown to produce more stable factors, is less prone
to overfitting, and has high parameter efficiency (Williams et al., 2018). We
use a hierarchical alternating least squares (HALS) algorithm (Cichocki &
Phan, 2009), as implemented in the tensortools Python library (Williams,
2024). Its efficiency when applied to real-world datasets has been compared
against several other algorithms in Phan & Cichocki (2008).

Each component resulting from our tensor will be constituted by three factors:
one representing a pattern $f_{a}$ of (vectorized) pairwise affinities; a
vector $f_{t}$ representing periods of time where such pattern was observed;
and a vector $f_{w}$ whose loadings, or coefficients, indicate how much each
worm’s affinities can be described by $f_{a}$ at those times expressed in
$f_{t}$. These components, combined, can be seen as a soft multi-clustering of
worms in terms of their activity patterns and the times at which they occur.
See Fig. 3b and Appendix A.4.

Figure 4: Tensor components for a real worm during an experiment in which
repellent and attractive stimuli were applied. The application of these
stimuli is important, since it provides correspondence between the time axes.
(a) One tensor component that selects out the interval during which repulsive
NaCl (salt) was applied. All worms reacted, but worms 5, 6, and 7 showed
higher responses for certain pairs of neurons (from the affinity factor, shown
as an unsorted matrix). The weighted graph visualized with a traditional
force-directed embedding is difficult to interpret, but the nested community
structure revealed by the NWSBM algorithm (b) provides a clear view of the
neuronal modules involved in salt perception. (c) A different tensor component
that selects several time intervals when stimuli were applied, but is
especially involved with the attractant 2,3-pentanedione. Most worms are again
implicated, but the community of neurons (d) is qualitatively different.

###  2.3 Neural affinity communities; Fig. 3c

Notice that, because the neurons are sorted in an arbitrary way, looking at
the affinity factor matrices should not be particularly enlightening. However,
because affinities are non-negative and bounded, they can be readily treated
as adjacencies (weighted edges) in a graph where the nodes correspond to
individual neurons. This enables us to use a community-detection algorithm to
cluster similar neurons based on their activity over time.

Although several methods for dynamic community detection (DCD) exist, using
both traditional ML (Rossetti & Cazabet, 2018) and deep learning (Skarding et
al., 2021), the vast majority are focused on unweighted networks (Li et al.,
2021). That is because in many applications a dynamic network refers to a
graph that changes its number of nodes or edges over time, not their weights
(with few exceptions, e.g., Guo et al. (2014)). Of course, for our purposes
the edge weights (i.e., affinity factors from NTF) are essential, as they
convey the likelihood that two neurons interact over some period of time
(section 2.1).

Since our tensor decomposition organizes the temporal dimension—i.e., which
networks are active at which times—we seek a community detection method
appropriate for handling edge weights. Additionally, we wish to prevent
spurious aggregations of nodes due to random fluctuations in the affinity
values (from possible noise in the activity traces). Therefore it is
preferable to infer a latent community structure that carries explanatory
power, as opposed to adopting classical descriptive methods that attempt to
find communities according to some preconceived notion of a good division of
the network into groups.

A classical model for community structure is the stochastic block model (SBM)
(Holland et al., 1983), which groups nodes according to their probabilities of
connection to the rest of the network. Statistical inference can be used to
find the best-fitting model to the data. A principled approach for this task
is to formulate generative models that allow this modular decomposition to be
found via statistical inference (Aicher et al., 2015; Peixoto, 2019). In
particular, the nested weighted-SBM algorithm (NWSBM) (Peixoto, 2018; 2014)
extends the SBM to be applied to weighted networks by treating edge weights as
covariates, and employs a Bayesian approach capable of inferring the number of
communities from the data in an unsupervised fashion. Moreover, it requires no
hyperparameter tuning.

For example, let the partition of the network into $C$ communities be denoted
as $\bm{c}=\\{c_{i}\\}$, where $c_{i}\in[0,C-1]$ is the community membership
of node $i$. A model can be defined that generates a network $\bm{G}$ with
probability

$P\left(\bm{G}|\bm{\theta},\bm{c}\right),$ (3)

where $\bm{\theta}$ is a set of parameters controlling how the node partition
affects the structure of $\bm{G}$. Assuming there is only one choice of
$\bm{\theta}$ that is compatible with the generated network, we may write the
Bayesian posterior probability as

$P\left(\bm{c}|\bm{G}\right)=\frac{P\left(\bm{G}|\bm{\theta},\bm{c}\right)P\left(\bm{\theta},\bm{c}\right)}{P\left(\bm{G}\right)}.$
(4)

The selected model is the one found to yield the smallest description length
(Peixoto, 2015), namely

$\mathcal{L}=-\ln P(\bm{G}|\bm{\theta},\bm{c})-\ln P(\bm{\theta},\bm{c}),$ (5)

since eq. 4 above can be rewritten as

$P\left(\bm{c}|\bm{G}\right)=\frac{\exp(-\mathcal{L})}{P\left(\bm{G}\right)}.$
(6)

Such a problem is NP-hard in general, so the likelihood is optimized
stochastically. Hence, instead of assuming that there is a single “best fit”
to the model, an efficient multi-flip Markov-Chain Monte Carlo (MCMC)
algorithm performs model averaging to increase the robustness of the results;
additional details can be found in Peixoto (2018). This algorithm was chosen
after running benchmark tests against several other popular community
detection methods (section 3.3, Table 1).

##  3 Results

Our approach was applied to the dataset from Yemini et al. (2021), in which
calcium activity from 189 neurons in the head of C. elegans were recorded from
a total of 21 individual worms. Each worm performed a series of behavioral
tasks in response to stimuli: one gustatory repellent (NaCl) and two olfactory
attractants (2-butanone and 2,3-pentanedione) were flushed into the worm
surroundings over periods of 10 s. Different worms were exposed to the stimuli
in 3 possible sequences. We restricted the data to sensory and inter-neurons,
since those are more likely to be directly responsive to the chemosensory
stimuli used.

###  3.1 Dynamic functional connectomes

Fig. 4a,c shows informative factors resulting from running NTF on our tensor
(section 2.3). Each affinity factor yields a vectorized affinity matrix, which
can be reorganized as the upper-triangular part of a symmetric matrix (for
ease of visualization). This matrix in turn can be interpreted as a weighted
adjacent matrix of a network in which nodes represent individual neurons
(section 2.3). The communities inferred by the NWSBM algorithm for each factor
then apply to certain time windows and worms, weighted by their respective
factor loadings, and represent a transient functional connectome. Two examples
are shown in Fig. 4b,d, with their respective temporal factors revealing that
such neural interactions are mostly active during the presence of NaCl and
pentanedione, respectively. A richly-connected individual community from the
NaCl-sensing connectome is further analyzed in Fig. 5.

Figure 5: A transient community for salt sensing. (a) We here highlight one
community from the dendrogram shown in Fig. 4a (cf. inset below). By reverting
back to the original activity traces, we can confirm that the individual
neuronal responses agree with a functional circuit organization, and that the
affinity measure is meaningful. (b) The temporal factor of the corresponding
tensor component shows that this particular circuit is most prominently
involved in a response to NaCl (salt), plus minor but non-negligible responses
to the other two stimuli (2-butanone and 2,3-pentanedione). (c) Individual
activity traces for the most strongly connected neurons within the community,
highlighted in proportion to the instantaneous temporal loadings. Note that
the traces are most similar precisely during the NaCl interval for several
worms; traces are plotted from the three worms with highest loadings in the
corresponding worm factor (see Fig. 4a). (d) Note the complexity of the traces
outside this interval. In particular, even though in worm 7 all the selected
neurons were correlated throughout the experiment, this was not true for other
worms, which explains the NaCl specificity in the temporal factor. Notice
that, because the affinities are computed from absolute derivatives, neurons
will be strongly connected even when their traces have opposite signs.

###  3.2 Validation experiments

Our algorithm predicted the involvement of multiple neurons in the salt-
sensing (NaCl) circuit, as those neurons found clustered with the canonical
NaCl-sensing neuron ASE. Among these predictions were several neurons not
previously known to play a role in this salt-sensation circuit. For example,
our algorithm predicted that the worm’s sole aversive olfactory neuron AWB
(Troemel et al., 1997) is substantially involved in salt sensation. Therefore,
to experimentally test the accuracy of this surprising of a previously unknown
role for AWB in sensing NaCl.

To do so, we compared salt avoidance in worms where AWB was functional versus
worms where AWB was silenced. Briefly, we used a standard “drop test”
(Bazzicalupo & Hilliard, 2005) to measure avoidance of salt, wherein a drop of
160 mM NaCl buffer (as used in Yemini et al. (2021)) was placed in front of
worms, and their response was recorded as avoidance if they reversed, and non-
avoidance if they continued moving forward. Each of our samples consisted of
25–50 worms that were drop tested, and collectively scored for their mean
avoidance. We tested 12 such samples of functional AWB controls versus 12
samples with AWB silenced. To silence AWB, we used a published collection of
worms (Wang et al., 2017), and crossed an AWB neuron-specific driver (syIs666)
to an inhibitory histamine-gated chloride channel (HisCl1–syIs373) effector
(Pokala et al., 2014). These AWB$>$HisCl1 worms (strain YYY24, with genotype
syIs666;syIs373) have a functional AWB that is silenced upon treatment with
histamine.

We found that silencing AWB substantially increased salt avoidance ($p$-value
= 5.7e-11, mean effect size = +25%); see Fig. A2 (Appendix A.3). Our results
are particularly striking since AWB is an aversive neuron and thus silencing
it would be expected to decrease avoidance rather than increase it as we
observed. This suggests the potential for oppositional interactions between
salt sensation and olfactory aversion circuits. Our results not only validate
the predictive power of our algorithm, they further highlight its strength in
predicting results that expert scientific researchers may find contradict
assumptions they have based on published canonical neuronal roles.

###  3.3 Comparison with other community detection methods

The NWSBM method was compared to other popular algorithms for community
detection, including classical methods such as Louvain (Blondel et al., 2008;
Rossetti et al., 2019) and Combo (Sobolevsky et al., 2014), based on
modularity maximization; asymptotic surprise (AS) (Traag et al., 2015); NNSED
(Sun et al., 2017; Rozemberczki et al., 2020), a non-negative encoder-decoder
approach; and GNNS (Sobolevsky & Belyi, 2022; Belyi, 2022), based on recurrent
graph neural networks. The quality of their results was evaluated using
normalized mutual information (NMI) (Danon et al., 2005), a widely adopted
measure of agreement between clustering assignments (Fortunato, 2010). Because
it requires ground truth—not available for the vast majority of real network
datasets with weighted edges—, the algorithms above were evaluated on a
benchmark of synthetic networks generated using the popular Lancichinetti-
Fortunato-Radicchi (LFR) algorithm (Lancichinetti & Fortunato, 2009), a
standard tool for creating benchmarks with various types of networks (Gopalan
& Blei, 2013; Fortunato & Hric, 2016; Yang et al., 2016). It provides several
parameters that control the how the connections and their weights will be
distributed within and across communities, allowing for a rich variety of
block structures to be generated. Our benchmark consisted of 9 different types
of network (see Appendix A.2 and Fig. A1 for details). Table 1 shows the
average NMI scores obtained by each algorithm over ten instances of each type
of network. Although most of the algorithms tested was able to score highest
for at least one network type, NWSBM outperformed the others in the majority
of cases.

Table 1: Mean NMI scores of algorithms applied to the weighted-LFR benchmark

.

Algorithm | Net 1 | Net 2 | Net 3 | Net 4 | Net 5 | Net 6 | Net 7 | Net 8 | Net 9
---|---|---|---|---|---|---|---|---|---
NWSBM | 0.28 | 0.33 | 0.39 | 0.65 | 0.39 | 0.54 | 0.25 | 0.51 | 0.64
Louvain | 0.08 | 0.19 | 0.13 | 0.04 | 0.23 | 0.29 | 0.12 | 1.0 | 0.19
Combo | 0.13 | 0.18 | 0.14 | 0.04 | 0.23 | 0.28 | 0.12 | 1.0 | 0.16
AS | 0.18 | 0.16 | 0.16 | 0.47 | 0.20 | 0.12 | 0.28 | 0.20 | 0.20
NNSED | 0.0 | 0.18 | 0.38 | 0.32 | 0.39 | 0.50 | 0.39 | 0.43 | 0.45
GNNS100 | 0.14 | 0.19 | 0.15 | 0.04 | 0.22 | 0.28 | 0.12 | 1.0 | 0.18

##  4 Discussion and Conclusion

We present here an algorithm for dynamic connectome discovery in the nematode
worm C. elegans. The key insight was to use the differential affinity between
neurons in a tensor factorization approach, rather than (the more common) use
of individual neurons. The advantages were (i) that non-linear methods could
be used in calculating affinities before they were tensorized; (ii) the
temporal factor in each tensor component revealed experimental epochs during
which functional circuits appeared; and (iii) the affinity factors could be
remade into weighted graphs on which community detection algorithms could be
run. In the end, the algorithm was able to make (surprising) predictions about
individual neurons that were involved in unexpected functional roles, which
were experimentally confirmed.

While this last point illustrates a role for machine learning in biological
research, we also showed that reverting the community structure back to the
original affinities (Fig. 5) can explain how the original traces caused a
given pair of nodes to be ultimately grouped together.

Finally, we believe that our approach could be used more widely in
understanding community behavior not only in neuroscience, but also in social
and ethological situations. Working with the higher-order affinities directly
informed the biology in our case; this is most likely true for other dynamical
interactions as well.

#### Acknowledgments

Supported by NIH Grant 1R01EY031059, NSF Grant 1822598.

## References

  * Aicher et al. (2015) Christopher Aicher, Abigail Z Jacobs, and Aaron Clauset.  Learning latent block structure in weighted networks.  _Journal of Complex Networks_ , 3(2):221–248, 2015.
  * Atanas et al. (2023) Adam A Atanas et al.  Brain-wide representations of behavior spanning multiple timescales and states in C. elegans.  _Cell_ , 186(19):4134–4151.e31, 2023.  ISSN 0092-8674.
  * Bargmann & Marder (2013) Cornelia I Bargmann and Eve Marder.  From the connectome to brain function.  _Nature Methods_ , 10(6):483–490, 2013.
  * Bassett et al. (2011) Danielle S Bassett et al.  Dynamic reconfiguration of human brain networks during learning.  _Proc Natl Acad Sci USA_ , 108(18):7641–7646, 2011.
  * Bazzicalupo & Hilliard (2005) Paolo Bazzicalupo and Massimo Hilliard.  Behavioral assays – Drop assay.  In Anne C. Hart (ed.), _WormMethods_. WormBook, 2005.  doi: 10.1895/wormbook.1.87.1.  URL http://www.wormbook.org/chapters/www˙behavior/behavior.html#sec4˙3.
  * Beets et al. (2023) Isabel Beets et al.  System-wide mapping of peptide-GPCR interactions in C. elegans.  _Cell Reports_ , 42(9), 2023.
  * Belyi (2022) Alexander Belyi.  GNNS.  https://github.com/Alexander-Belyi/GNNS, 2022.
  * Bentley et al. (2016) Barry Bentley et al.  The multilayer connectome of Caenorhabditis elegans.  _PLoS Computational Biology_ , 12(12):e1005283, 2016.
  * Blondel et al. (2008) Vincent D Blondel, Jean-Loup Guillaume, Renaud Lambiotte, and Etienne Lefebvre.  Fast unfolding of communities in large networks.  _Journal of statistical mechanics: theory and experiment_ , 2008(10):P10008, 2008.
  * Bro & De Jong (1997) Rasmus Bro and Sijmen De Jong.  A fast non-negativity-constrained least squares algorithm.  _J Chemom_ , 11(5):393–401, 1997.
  * Carroll & Chang (1970) J Douglas Carroll and Jih-Jie Chang.  Analysis of individual differences in multidimensional scaling via an n-way generalization of “Eckart-Young” decomposition.  _Psychometrika_ , 35(3):283–319, 1970.
  * Cecere et al. (2021) Zachary T Cecere et al.  State-dependent network interactions differentially gate sensory input at the motor and command neuron level in Caenorhabditis elegans.  _bioRxiv_ , pp. 2021–04, 2021.
  * Cichocki & Phan (2009) Andrzej Cichocki and Anh-Huy Phan.  Fast local algorithms for large scale nonnegative matrix and tensor factorizations.  _IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences_ , 92(3):708–721, 2009.
  * Cichocki et al. (2009) Andrzej Cichocki, Rafal Zdunek, Anh Huy Phan, and Shun-ichi Amari.  _Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation_.  John Wiley & Sons, 2009.
  * Cook et al. (2019) Steven J Cook et al.  Whole-animal connectomes of both Caenorhabditis elegans sexes.  _Nature_ , 571(7763):63–71, 2019.
  * Dag et al. (2023) Ugur Dag et al.  Dissecting the functional organization of the C. elegans serotonergic system at whole-brain scale.  _Cell_ , 186(12):2574–2592.e20, 2023.  ISSN 0092-8674.  doi: https://doi.org/10.1016/j.cell.2023.04.023.
  * Danon et al. (2005) Leon Danon, Albert Diaz-Guilera, Jordi Duch, and Alex Arenas.  Comparing community structure identification.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2005(09):P09008, 2005.
  * Dyballa et al. (2024) Luciano Dyballa et al.  Population encoding of stimulus features along the visual hierarchy.  _Proc Natl Acad Sci USA_ , 121(4):e2317773121, 2024.
  * Fortunato (2010) Santo Fortunato.  Community detection in graphs.  _Physics Reports_ , 486(3-5):75–174, 2010.
  * Fortunato & Hric (2016) Santo Fortunato and Darko Hric.  Community detection in networks: a user guide.  _Physics Reports_ , 659:1–44, 2016.
  * Gopalan & Blei (2013) Prem K Gopalan and David M Blei.  Efficient discovery of overlapping communities in massive networks.  _Proc Natl Acad Sci USA_ , 110(36):14534–14539, 2013.
  * Gordus et al. (2015) Andrew Gordus, Navin Pokala, Sagi Levy, Steven W Flavell, and Cornelia I Bargmann.  Feedback from network states generates variability in a probabilistic olfactory circuit.  _Cell_ , 161(2):215–227, 2015.
  * Guo et al. (2014) Chonghui Guo, Jiajia Wang, and Zhen Zhang.  Evolutionary community structure discovery in dynamic weighted networks.  _Physica A: Statistical Mechanics and its Applications_ , 413:565–576, 2014.
  * Harshman et al. (1970) Richard A Harshman et al.  Foundations of the PARAFAC procedure: Models and conditions for an “explanatory” multimodal factor analysis.  _UCLA Working Papers in Phonetics_ , 16:1–84, 1970.
  * Holland et al. (1983) Paul W Holland, Kathryn Blackmond Laskey, and Samuel Leinhardt.  Stochastic blockmodels: First steps.  _Social Networks_ , 5(2):109–137, 1983.
  * Jarrell et al. (2012) Travis A Jarrell et al.  The connectome of a decision-making neural network.  _Science_ , 337(6093):437–444, 2012.
  * Kato et al. (2015) Saul Kato et al.  Global brain dynamics embed the motor command sequence of Caenorhabditis elegans.  _Cell_ , 163(3):656–669, 2015.
  * Kolda (2006) Tamara G Kolda.  Multilinear operators for higher-order decompositions.  Technical Report SAND2006-2081, Sandia National Laboratories, April 2006.  URL http://www.osti.gov/scitech/biblio/923081.
  * Kolda & Bader (2009) Tamara G Kolda and Brett W Bader.  Tensor decompositions and applications.  _SIAM Rev_ , 51(3):455–500, 2009.
  * Lancichinetti & Fortunato (2009) Andrea Lancichinetti and Santo Fortunato.  Benchmarks for testing community detection algorithms on directed and weighted graphs with overlapping communities.  _Phys Rev E_ , 80(1):016118, 2009.
  * Lee & Seung (1999) Daniel D Lee and H Sebastian Seung.  Learning the parts of objects by non-negative matrix factorization.  _Nature_ , 401(6755):788–791, 1999.
  * Li et al. (2021) Shudong Li, Laiyuan Jiang, Xiaobo Wu, Weihong Han, Dawei Zhao, and Zhen Wang.  A weighted network community detection algorithm based on deep learning.  _Applied Mathematics and Computation_ , 401:126012, 2021.
  * Peixoto (2014) Tiago P Peixoto.  The graph-tool Python library.  _figshare_ , 2014.  doi: 10.6084/m9.figshare.1164194.  URL http://figshare.com/articles/graph˙tool/1164194.
  * Peixoto (2015) Tiago P Peixoto.  Model selection and hypothesis testing for large-scale network models with overlapping groups.  _Phys Rev X_ , 5(1):011033, 2015.
  * Peixoto (2018) Tiago P Peixoto.  Nonparametric weighted stochastic block models.  _Phys Rev E_ , 97(1):012306, 2018.
  * Peixoto (2019) Tiago P Peixoto.  Bayesian stochastic blockmodeling.  In _Advances in Network Clustering and Blockmodeling_ , chapter 11, pp. 289–332. John Wiley & Sons, 2019.
  * Phan & Cichocki (2008) Anh Huy Phan and Andrzej Cichocki.  Multi-way nonnegative tensor factorization using fast hierarchical alternating least squares algorithm (hals).  In _Proceedings of The 2008 International Symposium on Nonlinear Theory and its Applications_ , 2008.
  * Pokala et al. (2014) Navin Pokala, Qiang Liu, Andrew Gordus, and Cornelia I Bargmann.  Inducible and titratable silencing of Caenorhabditis elegans neurons in vivo with histamine-gated chloride channels.  _Proc Natl Acad Sci USA_ , 111(7):2770–2775, 2014.
  * Randi et al. (2023) Francesco Randi, Anuj K. Sharma, Sophie Dvali, and Andrew M. Leifer.  Neural signal propagation atlas of C. elegans, 2023.
  * Ripoll-Sánchez et al. (2022) Lidia Ripoll-Sánchez et al.  The neuropeptidergic connectome of C. elegans.  _bioRxiv_ , pp. 2022–10, 2022.
  * Rossetti & Cazabet (2018) Giulio Rossetti and Rémy Cazabet.  Community discovery in dynamic networks: a survey.  _ACM Computing Surveys (CSUR)_ , 51(2):1–37, 2018.
  * Rossetti et al. (2019) Giulio Rossetti, Letizia Milli, and Rémy Cazabet.  Cdlib: a python library to extract, compare and evaluate communities from complex networks.  _Applied Network Science_ , 4(1):1–26, 2019.
  * Rozemberczki et al. (2020) Benedek Rozemberczki, Oliver Kiss, and Rik Sarkar.  Karate Club: An API Oriented Open-source Python Framework for Unsupervised Learning on Graphs.  In _Proceedings of the 29th ACM International Conference on Information and Knowledge Management (CIKM ’20)_ , pp. 3125–3132. ACM, 2020.
  * Skarding et al. (2021) Joakim Skarding, Bogdan Gabrys, and Katarzyna Musial.  Foundations and modeling of dynamic networks using dynamic graph neural networks: a survey.  _IEEE Access_ , 9:79143–79168, 2021.
  * Sobolevsky & Belyi (2022) Stanislav Sobolevsky and Alexander Belyi.  Graph neural network inspired algorithm for unsupervised network community detection.  _Applied Network Science_ , 7(1):1–19, 2022.
  * Sobolevsky et al. (2014) Stanislav Sobolevsky, Riccardo Campari, Alexander Belyi, and Carlo Ratti.  General optimization technique for high-quality community detection in complex networks.  _Phys Rev E_ , 90(1):012811, 2014.
  * Sun et al. (2017) Bing-Jie Sun, Huawei Shen, Jinhua Gao, Wentao Ouyang, and Xueqi Cheng.  A non-negative symmetric encoder-decoder approach for community detection.  In _Proceedings of the 2017 ACM on Conference on Information and Knowledge Management_ , pp. 597–606, 2017.
  * Susoy et al. (2021) Vladislav Susoy et al.  Natural sensory context drives diverse brain-wide activity during C. elegans mating.  _Cell_ , 184(20):5122–5137, 2021.
  * Traag et al. (2015) V A Traag, R Aldecoa, and J-C Delvenne.  Detecting communities using asymptotical surprise.  _Phys Rev E_ , 92:022816, 2015.
  * Troemel et al. (1997) E R Troemel, B E Kimmel, and C I Bargmann.  Reprogramming chemotaxis responses: sensory neurons define olfactory preferences in C. elegans.  _Cell_ , 91(2):161–169, 1997.
  * Uzel et al. (2022) Kerem Uzel, Saul Kato, and Manuel Zimmer.  A set of hub neurons and non-local connectivity features support global brain dynamics in C. elegans.  _Current Biology_ , 32(16):3443–3459, 2022.
  * Varshney et al. (2011) Lav R Varshney, Beth L Chen, Eric Paniagua, David H Hall, and Dmitri B Chklovskii.  Structural properties of the Caenorhabditis elegans neuronal network.  _PLoS Computational Biology_ , 7(2):e1001066, 2011.
  * Wang et al. (2017) Han Wang et al.  cGAL, a temperature-robust GAL4-UAS system for Caenorhabditis elegans.  _Nature Methods_ , 14(2):145–148, 2017.
  * White et al. (1986) John G White et al.  The structure of the nervous system of the nematode Caenorhabditis elegans.  _Philos Trans R Soc Lond B Biol Sci_ , 314(1165):1–340, 1986.
  * Williams (2024) Alex H Williams.  tensortools.  https://github.com/neurostatslab/tensortools, 2024.
  * Williams et al. (2018) Alex H Williams et al.  Unsupervised discovery of demixed, low-dimensional neural dynamics across multiple timescales through tensor component analysis.  _Neuron_ , 98(6):1099–1115, 2018.
  * Witvliet et al. (2021) Daniel Witvliet et al.  Connectomes across development reveal principles of brain maturation.  _Nature_ , 596(7871):257–261, 2021.
  * Yang et al. (2016) Zhao Yang, René Algesheimer, and Claudio J Tessone.  A comparative analysis of community detection algorithms on artificial networks.  _Scientific Reports_ , 6(1):30750, 2016.
  * Yemini et al. (2021) Eviatar Yemini et al.  Neuropal: a multicolor atlas for whole-brain neuronal identification in C. elegans.  _Cell_ , 184(1):272–288, 2021.

##  Appendix A Appendix

###  A.1 Local Differential Affinities

We begin by computing each neuron $i$’s derivative $g_{i}(t)$ over time, pre-
smoothing the raw trace $\tau_{i}(t)$ with a Gaussian filter of bandwidth 0.5
sec to avoid the amplification of any measurement noise (Fig. 2, top panels).
Notice that periods where the derivative sign is constant must necessarily
start and end at points where the derivative is either exactly zero or changes
sign, and are commonly termed “zero-crossings”. Because sampling is discrete
over time, under a sampling interval $T$ a zero-crossing is computed as those
time points $t$ such that

$g_{i}(t-T/2)g_{i}(t+T/2)\leq 0,$ (7)

i.e., there is a change of sign between two consecutive samples, and

$|g_{i}(t-T/2)|+|g_{i}(t+T/2)|>0,$ (8)

i.e., the derivatives are not both exactly zero. When both conditions are met,
we say $g_{i}(t)$ is a zero-crossing, and denote it as

$g_{i}(t)\stackrel{{\scriptstyle\times}}{{\approx}}0.$ (9)

We then partition $g_{i}(t)$ by those intervals where it has constant sign
(i.e., those lying between zero-crossings); we call these _bumps_. Thus the
support of $g_{i}(t)$ (subset of the time domain during which the neuron’s
activity is not constant), may be represented as a set of $m$ time intervals
$b_{i}^{(1)},b_{i}^{(2)},\ldots,b_{i}^{(m)}$, where $b_{i}^{(t)}$ denotes the
interval $(t_{0},t_{1})$, such that

$g_{i}(t_{0})\stackrel{{\scriptstyle\times}}{{\approx}}0\stackrel{{\scriptstyle\times}}{{\approx}}g_{i}(t_{1}).$
(10)

Thus, two bumps $b_{i}^{(t)}$ and $b_{j}^{(t^{\prime})}$ from neurons $i$ and
$j$, respectively, will overlap when

$b_{i}^{(t)}\cap b_{j}^{(t^{\prime})}>0.$ (11)

The affinity $a_{ij}$ between neurons $i$ and $j$ at time $t$ is then defined
as the fraction of overlap between the areas under their derivative bumps
occurring at $t$: $b_{i}$ and $b_{j}$ respectively. This is analogous to the
Jaccard index:

$a_{ij}^{(t)}\equiv\frac{\mathrm{area}(b_{i})\cap\mathrm{area}(b_{j})}{\mathrm{area}(b_{i})\cup\mathrm{area}(b_{j})},$
(12)

where $\mathrm{area}(b_{i})$ is the unsigned area under the bump $b_{i}$:

$\mathrm{area}(b_{i})\equiv\int_{t_{b_{i}}}^{t_{b_{i}}^{\prime}}b_{i}^{(t)}dt,$
(13)

with the integral taken over the time interval under each bump. Note that
these affinity values lie in the range $[0,1]$, and because they are computed
locally, will be sensitive to spurious fluctuations of activity. To prevent
this, we take into account global information about the traces and weight this
by the relative change in trace levels during those local bumps compared to
the global range of values of $\tau_{i}(t)$.

Notice that the areas are unsigned, meaning we ignore the signs of the bumps.
This is motivated by the fact that two neurons with very similar but opposite
sign derivatives are still likely to be interacting by means of some
inhibitory mechanism. This further contributes to making the affinities
interpretable as indicating the likelihood of interactions between pairs of
neurons, regardless of the specific physiological connectivity mechanisms
involved.

###  A.2 Weighted communities benchmark

A total of 9 weighted networks with $N$=125 nodes (compatible with the size of
our affinity networks) were generated using the weighted LFR benchmark
(Lancichinetti & Fortunato, 2009). The mean neighborhood size $k$ influences
the resulting number of communities produced. How such neighborhoods are
distributed between intra- and inter-community connections is determined by
the topology ($\mu_{t}$) and weight ($\mu_{w}$) mixing parameters,
respectively. The weights within a node’s neighborhood are sampled from a
power-law distribution with exponent $\beta$. Finally, the degree sequence and
community sizes are drawn from power-law distributions with exponents
$-\tau_{1}$ and $-\tau_{2}$, respectively. By choosing different sets of
parameters for each network type, we aimed to vary several properties: the
number of communities (determined by the $k/N$ ratio ); the density of inter-
block connections ($\mu_{w}/\mu_{t}$ ratio); as well as the overall weight
distribution ($\beta$); we thus cover a variety of scenarios. Example weighted
adjacency matrices from each style and their respective parameters are shown
in Fig. A1.

Figure A1: Family of synthetic networks generated using the weighted LFR
benchmark. Each panel shows the weighted adjacency matrix of one instance of
each type of network, sorted by its ground-truth community membership. Default
parameter values used: $N=125$, $k=25$, $maxk=100$, $\mu_{w}=0.5$,
$\mu_{t}=0.5$, $\beta=2$, $\tau_{1}=2$, $\tau_{2}=1$. Each network type varies
a single parameter from this list (value shown above the matrix). (1–3)
Varying $\beta$ while keeping other parameters fixed dramatically changes the
weight distribution. (4–6) Smaller values of $\mu_{t}$ increase the intra-
block connections (irrespective of their weights). (7) A smaller $k$ creates a
larger number of small communities. (8–9) The $\mu_{w}/\mu_{t}$ ratio controls
the ratio between inter- and intra-block weights.

###  A.3 Experimental validation results

Silencing AWB substantially increased salt avoidance ($p$-value=5.7e-11, mean
effect size = +25%), validating the predictive power of our algorithm; see
Fig. A2. Surprisingly, this suggests the potential for oppositional
interactions between salt sensation and olfactory aversion circuits.

Figure A2: Results for avoidance to 160 mM NaCl for AWB controls (+AWB) vs.
silenced (-AWB). Two-tailed unpaired t-test P=5.7574e-11, mean effect size is
+25% avoidance when silencing AWB, N=12 vs. 12. red=95% CI, blue=1 SD.

###  A.4 Selection of the number of tensor components

Our logic for selecting the number of tensor components, $R$, was to use as
many as possible to minimize reconstruction error, provided the results across
multiple random initializations remained stable (i.e., small variance). Small
reconstruction error suggests a faithful representation, and small variability
guarantees that the retained components are robust. Based on Fig. A3a, the
error variability (std. dev. across 15 runs) reaches a minimum when
$R\approx$14–15 (shaded area), then increases sharply for $R>15$. We therefore
selected R=15. Fig. A3b shows that the distributions of error for each choice
of $R$ are separate until $R=15$, after which they being to overlap; this
further strengthens our confidence in our choice for the number of tensor
components.

Figure A3: Selecting the number of tensor components, $R$. (a) Reconstruction
error as a function of $R$. Although the mean error (blue curve) always
decreases with increasing $R$, its variability (std. dev. across 15 runs,
dashed gray curve) reaches a minimum when $R\approx 14$–15 (shaded area), then
increases sharply for $R>15$. (b) The empirical distributions of error for
each choice of $R$ remain separate until $R$=15, but starting at $R=16$ they
begin to overlap. These results motivate our choice for the number of tensor
components to use, namely 15.


††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Nonlinear Thouless Pumping of Solitons Across an Impurity

Xuzhen Cao  State Key Laboratory of Quantum Optics and Quantum Optics Devices,
Institute of Laser Spectroscopy, Shanxi University, Taiyuan 030006, China
Collaborative Innovation Center of Extreme Optics, Shanxi University, Taiyuan,
Shanxi 030006, China  Chunyu Jia  Department of Physics, Zhejiang Normal
University, Jinhua 321004, China  Ying Hu  State Key Laboratory of Quantum
Optics and Quantum Optics Devices, Institute of Laser Spectroscopy, Shanxi
University, Taiyuan 030006, China  Collaborative Innovation Center of Extreme
Optics, Shanxi University, Taiyuan, Shanxi 030006, China  Zhaoxin Liang
Department of Physics, Zhejiang Normal University, Jinhua 321004, China

###### Abstract

The nonlinear Thouless pumping is an exciting frontier of topological physics.
While recent works have revealed the quantized motion of solitons in Thouless
pumps, the interplay between the topology, nonlinearity and disorder remains
largely unexplored. Here, we investigate the nonlinear Thouless pumping of
solitons in the presence of an impurity in the context of a Bose–Einstein
condensate. Using both the Gross–Pitaevskii equation and Lagrangian
variational approach, we analyze the interaction between a moving soliton and
an impurity. Without the pump, the soliton can pass through a light impurity,
but gets trapped by the impurity with large mass. In marked contrast, we find
the soliton in Thouless pumps can always transit through the impurity, and its
motion is topologically quantized. Our result explicitly showcases the
robustness of topological soliton pumping against microscopic imperfections,
and opens a new perspective in the information processing with solitons.

##  I Introduction

Central to understanding quantized topological transport is the concept of
Thouless pumping [1, 2, 3, 4]. There, a quantum particle acted on by a
periodic potential that varies adiabatically and periodically in time shows a
quantized motion; the quantized displacement, stable against disorder, is
dictated by the Chern number of the underlying band structure in the momentum-
time space [4]. The experimental implementation of Thouless pumps and the
observations of topologically quantized transport have been achieved in a wide
variety of systems, including ultracold atoms [5, 6, 7, 8, 9, 10, 11],
photonics [12, 13, 14, 15, 16], and spin systems [17, 18]. Interestingly,
these synthetic topological systems can operate even beyond the linear regime,
i.e., in the presence of nonlinearities that arise from the Kerr effect in
optical platforms [19, 20, 21, 22, 23] or the interparticle interactions in
atomic setups [24, 25, 26]. At present, exploring nonlinear topological
pumping has attracted significant interests and experimental efforts.

A paradigm of the nonlinear Thouless pumping concerns the topological
transport of nonlinear excitations, known as solitons [27, 28, 29, 30, 31, 32,
33]. It has been shown that for weak nonlinearity, the motion of the soliton
can be topologically quantized, where the quantized displacement is directly
related to the topology of the underlying band structures [34, 12, 6, 7, 13,
8, 17, 15, 10, 11, 35]; whereas, the quantization breaks down in the strong
nonlinearity limit [3, 36, 37, 38, 39, 40, 20, 41, 42, 43, 44, 45, 46, 47, 48,
49, 50, 29, 28]. More recent studies have shown that the nonlinearity can
fundamentally modify the adiabatic dynamics and thus the Berry connection
[33]. So far, however, exploring the interplay of nonlinear excitations,
topology and disorder remains largely unchartered territory.

In this work, we explore the soliton transport in nonlinear Thouless pumps in
the presence of an impurity, based on a quasi-one-dimensional (1D) Bose-
Einstein condensate (BEC). We first numerically simulate the motion of the
soliton through the impurity using the Gross–Pitaevskii (GP) equation [51, 52,
53]. Further, we analytically study the interaction between the soliton and
the impurity using the Lagrangian variational approach. Without the pumps, the
fate of the soliton relies crucially on the effective mass of the impurity:
For light impurities, the soliton can easily pass through, but for impurities
with large masses, the soliton can get trapped. However, in the presence of
Thouless pumps, we find the soliton can always transit through the impurity,
and its motion becomes topologically quantized. Such topology-enhanced soliton
transport through the impurity provides insights into the interplay between
nonlinearity, topology and disorder, and explicitly elucidates the stability
of the topological nonlinear pumping against local imperfections. Our result
may enable novel possibilities in the information processing with solitons.

The paper is organized as follows. In Sec. II, we present detailed
descriptions of our model system. In Sec. III, by numerically solving GP
equation, we investigate the nonlinear Thouless pumping of solitons in the
presence of an impurity. In Sec. IV, we analytically derive the equations of
motions for the soliton using Lagrangian variational approach, to gain more
understandings of its interaction with the soliton. Both the numerical and
analytical results explicitly demonstrate the robustness of topological
soliton pumping against imperfections. In Sec. V, we summarize our work.

##  II Model system

We consider a 3D BEC with an attractive interatomic interaction in the
following trap geometry [5, 6]: in the $x$\- direction, the BEC is trapped in
an optical superlattice $V_{\text{OSL}}(x,t)$ and an impurity potential
$V_{{\text{imp}}}$, while in the $y$\- and $z$\- directions, the BEC is
tightly confined in the harmonic trap with the large trap frequency
$\omega_{\perp}$ such that the atomic motions in these directions are
affectively frozen. At the mean-field level, the considered BEC can be
described by the condensate wavefunction $\Psi(x,y,z)=\psi(x)\phi(y,z)$ with
$\phi(y,z)=1/(\sqrt{\pi}a_{\perp})e^{-[y^{2}+z^{2}]/2a_{\perp}^{2}}$ and
$a_{\perp}=\sqrt{\hbar/(m\omega_{\perp})}$. Averaging out the transverse
degrees of freedom results in an effective qusi-1D BEC described by the
following GP equation [54, 55], i.e.,

$\\!\\!\\!i\hbar\frac{\partial\psi}{\partial
t}=-\frac{\hbar^{2}}{2m}\frac{\partial^{2}\psi}{\partial
x^{2}}+g|\psi|^{2}\psi+[V_{{\text{imp}}}(x)+V_{{\text{OSL}}}(x,t)]\psi.$ (1)

Here the $m$ is the atomic mass, and $g=2\hbar^{2}a_{s}/ma_{\perp}^{2}$ is the
confinement-modified coupling constant, where $a_{s}<0$ is the negative s-wave
scattering length.

Figure 1: Soliton transport in the presence of a light impurity without and
with nonlinear Thouless pumps. In (a) and (c), the pump is absent, where
$V_{s}=V_{l}=0$. In (b) and (d), the pump is present, with $V_{s}=V_{l}=25$
and the periodicity $T$. In (a) and (b), the solitonic motion in space (x) and
time (t/T) is simulated via the numerical solutions of the dimensionless GP
Eq. (2) with $\gamma=5$ and the initial state
$\psi(x,0)=Ae^{-(x-x_{0})^{2}/l^{2}}e^{ikx}$, where $x_{0}=-3$, $l=0.4$ and
$A^{2}=15$. The color bar denotes the magnitude of $|\psi|^{2}$. The blue line
in (a) denotes the center of the impurity potential. The green line in (b)
denotes the pumping potential at $t=0$. Panels (c) and (d) present the results
for the center-of-mass position of the soliton based on the Euler-Lagrangian
Eqs. (7)-(11) with scaled parameters $\tilde{\gamma}=0.0258$ and $\eta(0)=1$,
$\kappa(0)=0.0258$, $\alpha(0)=\beta(0)=a(0)=0$.

In Eq. (1), we consider the time-dependent superlattice potential
$V_{\text{OSL}}/E_{R}=-V_{\text{s}}\cos^{2}\left(\frac{2\pi
x}{d}\right)-V_{\text{l}}\cos^{2}\left(\frac{\pi x}{d}-\Omega t\right)$, with
the recoil energy $E_{R}=\pi^{2}/2md^{2}$. It consists of a primary lattice
with the spatial period $d$ and lattice strength $V_{\text{s}}$, and a
superimposed second lattice with the period $d/2$ and the strength $V_{l}$;
the relative phase between the two lattices is $\Omega t$. Such a time-
dependent potential has been experimentally realized in quantum gases to
implement Thouless pumps [6].

Moreover, we consider the impurity potential in Eq. (1) as
$V_{\text{imp}}=-V_{0}\delta(x)$, corresponding to an impurity centered at
$x=0$ with the impurity strength $V_{0}$. Such an impurity trap has been
experimentally realized using bichromatic fields [56, 57], or ‘dark-state’
optical potentials [58, 59].

For subsequent analysis, we shall rescale GP equation (1) into the
dimensionless form via $x\rightarrow x/d$, $t\rightarrow
2E_{R}t/\pi^{2}\hbar$, $\psi\rightarrow\psi\sqrt{|g|\pi^{2}/2E_{R}}$,
$\gamma=V_{0}\pi^{2}/2d$, $V_{s(l)}=V_{s(l)}\pi^{2}/2$ and
$\nu=\Omega\hbar\pi^{2}/2E_{R}$. The resulting dimensionless GP equation reads

$\displaystyle i\frac{\partial}{\partial
t}\psi=-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}\psi-|\psi|^{2}\psi-\gamma\delta(x)\psi+V_{\text{OSL}}(x,t)\psi,$ (2)

with the normalization condition $\int|\psi|^{2}dx=N$. The dimensionless
pumping potential is explicitly given by

$\displaystyle V_{\text{OSL}}(x,t)=-V_{s}\cos^{2}(2\pi x)-V_{l}\cos^{2}(\pi
x-\nu t)$ (3)

which is periodic in time with the periodicity $T=\pi/\nu$. In addition, the
$\gamma$ in Eq. (2) reflects the strength of the defect, and can be
interpreted as the effective mass of the impurity [55].

##  III Solitonic Thouless pumping through an impurity: GP theory

In this section, we present numerical solutions of the soliton motion based on
the GP equation (2). Before proceeding, let us first briefly recapitulate
several familiar cases studied before. In the absence of both impurity
($\gamma=0$) and Thouless pump ($V_{\text{OSL}}=0$), the quasi-1D BEC
described by GP Eq. (2) is known to support solitons as the nonlinear
excitations. Without impurity ($\gamma=0$), Equation (2) with
$V_{\text{OSL}}\neq 0$ has been used to study the nonlinear Thouless pumping
of solitons such as in Ref. [29, 30]. There, the soliton is found to undergo a
topologically quantized transport even for the nonlinearity below some
critical value, and becomes trapped if the nonlinearity strength exceeds the
threshold. At the meantime, Equation (2) with $\gamma\neq 0$ and
$V_{\text{OSL}}=0$ has been used to study the interaction between a moving
soliton with an impurity. As shown in Ref. [54, 55], the behavior of the
soliton crucially depends on the impurity mass $\gamma$. Specifically, the
soliton can only move through a very light impurity, but becomes to be trapped
when $\gamma$ becomes large.

In this work, we are interested in the motion of the soliton when both
Thouless pump and the impurity are present, i.e., $V_{\text{OSL}}\neq
0,\gamma\neq 0$. To this end, we start with the numerical solutions of the GP
Eq. (2). We assume at time $t=0$, the soliton initially takes the form
$\psi(x,0)=Ae^{-(x-x_{0})^{2}/l^{2}}e^{ikx}$ [29, 55]. Here, the $k$ and
$x_{0}$, respectively, correspond to the initial momentum and the center-of-
mass position of the soliton; the $l$ represents the initial width of the
envelope and $|A|^{2}$ measures the nonlinear strength. In order to elucidate
the role of the topological nonlinear pump, in Figs. 1 (a) and (b) as well as
Fig. 2 (a) and (b), we compare typical simulations of the GP equations under
$V_{\text{OSL}}=0$ and $V_{\text{OSL}}\neq 0$ for various $\gamma$.

Consider first the regime of light impurities, as shown in Fig. 1(a) for
$\gamma=5$ without the pump. We see that the soliton directly tunnels through
the light impurity potential centered at $x=0$, where the collision within the
impurity only causes a small deviation of the soliton from its original path.
In the presence of the topological pump with $V_{s}=V_{l}=25$, as shown in
Fig. 1 (b), the soliton not only transmits through the impurity, but also the
solitonic displacement becomes quantized to the underlying Chern number of the
pump, despite the presence of the impurity.

In comparison, Figure 2(a) and (b) illustrate the soliton motion for large
impurity mass. Without the pump, as illustrated in Fig. 2(a), we see that the
soliton gets trapped in the impurity potential, as opposed to the case of
light impurity in Fig. 1(a). Interestingly, after applying the Thouless pump
with $V_{s}=V_{l}=25$ in Eq. (3), the soliton is able to pass through the
impurity instead of being localized by it [see Fig. 2 (b)]. Moreover, the
solitonic displacement in a period becomes quantized to unity. This phenomena
contrast strongly with the pump-free counterpart, and showcase the stability
of the topological nonlinear pumping against the static disorder.

##  IV Solitonic Thouless pumping through an impurity: variational approach

To gain further understandings how the interaction between the soliton and the
impurity is affected by the Thouless pump, in this sectionn we develop an
analytical study of the soliton motion using the Lagrangian variational
approach.

Writing $\psi(x,t)=A\phi(X,\tau)e^{i(V_{s}+V_{l})t/2}$ with the new variables
$\tau=A^{2}t$ and $X=Ax$, we recast Eq. (2) into the following form

$i\phi_{\tau}=-\frac{1}{2}\phi_{XX}-|\phi|^{2}\phi-\tilde{\gamma}\delta(X)\phi+\tilde{V}_{\text{OSL}}(X,\tau)\phi.$
(4)

Here, we have used the notations $\phi_{\tau}={\partial\phi}/{\partial\tau}$,
$\phi_{XX}={\partial^{2}\phi}/{\partial X^{2}}$, and the scaled parameters
$\tilde{V}_{\text{OSL}}(X,\tau)=-\tilde{V_{s}}\cos\left({4\pi
X}/{A}\right)-\tilde{V_{l}}\cos\left({2\pi X}/{A}-{2\nu\tau}/{A^{2}}\right)$
with $\tilde{V}_{s}=V_{s}/(2A^{2})$ and $\tilde{V_{l}}=V_{l}/(2A^{2})$, and
$\tilde{\gamma}=\gamma/A$.

As the benchmark, we recall that without the pump
($\tilde{V}_{\text{OSL}}=0$), there exist two special cases where known
solutions exist. (i) Without both the nonlinearity and the pump, Equation (4)
is reduced to the familiar linear Schrödinger equation with the delta-
potential, whose solution is a bound state localized in the defect potential
(see e.g. Ref. [54, 55] and the references therein), i.e.,
$\phi_{\text{im}}(x)\sim\sqrt{\tilde{\gamma}}e^{-\tilde{\gamma}|X|}$ with the
localization length $1/\tilde{\gamma}$. (ii) Without the impurity
($\tilde{\gamma}=0$) and the pump, Equation (4) supports moving solitons of
the form
$\phi_{\text{soliton}}\sim\eta\text{sech}(\eta(X-c\tau))\exp{(i(\eta^{2}-c^{2})\tau/2+icX)}$
[54], with the amplitude $\eta$, width $\eta$, and wavevector $c$.

Figure 2: Soliton transport in the presence of a heavy impurity without and
with nonlinear Thouless pumps. In (a) and (c), the pump is absent, where
$V_{s}=V_{l}=0$. In (b) and (d), the pump is present, with $V_{s}=V_{l}=25$.
In (a) and (b), the solitonic motion is simulated via the numerical solutions
of the GP equation (2) with $\gamma=75$ and the initial state
$\psi(x,0)=Ae^{-(x-x_{0})^{2}/l^{2}}e^{ikx}$, where $x_{0}=-3$, $l=0.4$ and
$A^{2}=15$. Panels (c) and (d) present the results for the center-of-mass
position of the soliton based on the Euler-Lagrangian equation Eqs. (7)-(11)
with scaled parameters $\tilde{\gamma}=0.258$ and $\eta(0)=1$,
$\kappa(0)=0.0258$, $\alpha(0)=\beta(0)=a(0)=0$.

In the presence of the Thouless pump ($\tilde{V}_{\text{OSL}}(X,\tau)\neq 0$),
neither the bound state $\phi_{\text{im}}$ nor the soliton
$\phi_{\text{soliton}}$ stands as the general solutions of Eq. (4). To seek
the solutions in this case, we use the Lagrangian variational approach along
the line of Refs. [54, 55]. We assume a trial wave function which is the
combination of the moving soliton and the impurity-induced bound state, i.e.,

$\\!\\!\\!\\!\\!\phi(X,\tau)=\left[\eta\text{sech}(\eta X-z)e^{i\kappa
X}+a\sqrt{\tilde{\gamma}}e^{-\tilde{\gamma}|X|+i\alpha}\right]e^{i\beta},$ (5)

where $\eta(\tau)$, $z(\tau)$, $\kappa(\tau)$, $a(\tau)$, $\alpha(\tau)$ and
$\beta(\tau)$ are the variational parameters to be determined below.
Specifically, $\eta(\tau)$ and $z(\tau)$ are the amplitude and the center-of-
mass position of the soliton, respectively, $\kappa(\tau)$ is the wavenumber
of the soliton, $a(\tau)$ is the amplitude of the bound-state component,
$\alpha(\tau)$ is the relative phase between the soliton and the bound state,
and $\beta(\tau)$ is the global phase of the trial function. The key
assumption underlying the ansatz (5) is that the functional forms of the
soliton and the impurity-induced bound state are preserved in the presence of
the Thouless pumping, but the corresponding parameters become slowly time-
dependent.

With the help of the ansatz (5), the Lagrangian of $L$ corresponding to the GP
Eq. (4) can be derived directly following standard procedures [55, 54] as
follows,

$\displaystyle L$ $\displaystyle=$
$\displaystyle\\!-\\!2\eta\frac{d\beta}{d\tau}\\!-\\!2z\frac{d\kappa}{d\tau}\\!-\\!\kappa^{2}\eta\\!-\\!\frac{\eta^{3}}{3}\\!+\\!\frac{2}{3}\eta^{3}\\!+\\!\tilde{\gamma}\eta^{2}\text{sech}^{2}(z)\\!+\\!\frac{2\pi^{2}}{A}[2\tilde{V_{s}}\cos(\frac{4\pi
z}{\eta A})\text{csch}(\frac{2\pi^{2}}{\eta A})+\tilde{V}_{l}\cos(\frac{2\pi
Az-2\eta\nu\tau}{\eta A^{2}})\text{csch}(\frac{\pi^{2}}{\eta A})]$ (6)
$\displaystyle+$ $\displaystyle 2\eta
a\tilde{\gamma}^{3/2}\text{sech}(z)\cos(\alpha)-a^{2}\frac{\partial}{\partial\tau}\beta-a^{2}\frac{\partial}{\partial\tau}\alpha+\frac{a^{2}\tilde{\gamma}^{2}}{2}+\frac{a^{2}\tilde{\gamma}^{2}A^{2}\tilde{V}_{s}}{4\pi^{2}+A^{2}\tilde{\gamma}^{2}}+\frac{a^{2}\tilde{\gamma}^{2}A^{2}\tilde{V}_{l}}{\pi^{2}+A^{2}\tilde{\gamma}^{2}}\cos\left(\frac{2\nu\tau}{A^{2}}\right)+\frac{1}{4}a^{4}\tilde{\gamma}.$

In the Lagrangian (6), the terms in the first line correspond to contribution
from the moving soliton; the terms in order of $O(a^{1})$ come from the mutual
interaction between the bound state and the moving soliton; the terms in order
of $O(a^{2})$ describe the impurity’s energy in the presence of the pumping
potential; the terms in order of $O(a^{4})$ origin from the self-interaction
of the bound state.

Thus the Euler-Lagrangian equations for the variational parameters are given
by $\partial L/\partial q_{i}-d/dt(\partial L/\partial\dot{q_{i}})=0$ with
$q_{i}=\eta,z,\kappa,a,\alpha,\beta$. Crucial for our study is the center-of-
mass position $z(\tau)$ of the soliton, whose equation of motion are given by

$\displaystyle\frac{d^{2}z}{d\tau^{2}}=$ $\displaystyle-$ $\displaystyle
a\eta\kappa\tilde{\gamma}^{3/2}\text{sech}(z)\sin(\alpha)-\frac{8\pi^{3}\tilde{V_{s}}}{A^{2}}\text{csch}\left(\frac{2\pi^{2}}{A\eta}\right)\sin\left(\frac{4\pi
z}{A\eta}\right)-\frac{2\pi^{3}\tilde{V_{l}}}{A^{2}}\text{csch}\left(\frac{\pi^{2}}{A\eta}\right)\sin\left(\frac{2A\pi
z-2\nu\tau\eta}{A^{2}\eta}\right)$ (7) $\displaystyle-$
$\displaystyle\tilde{\gamma}\eta^{3}\text{sech}^{2}(z)\text{tanh}(z)-a\eta^{2}\tilde{\gamma}^{3/2}\cos(\alpha)\text{sech}(z)\text{tanh}(z).$

The dynamics of $z(\tau)$ depends on the time-dependent parameters $\kappa$,
$\eta$, $a$, $\alpha$ and $\beta$, whose equations of motion are given as
follows

$\displaystyle\frac{d}{d\tau}a$ $\displaystyle=$
$\displaystyle\eta\tilde{\gamma}^{3/2}\text{sech}(z)\sin(\alpha),\frac{d}{d\tau}\eta=-a\eta\tilde{\gamma}^{3/2}\text{sech}(z)\sin(\alpha).$
(8) $\displaystyle\frac{d}{d\tau}\kappa$ $\displaystyle=$
$\displaystyle-\frac{8\pi^{3}\tilde{V_{s}}}{A^{2}\eta}\text{csch}\left(\frac{2\pi^{2}}{A\eta}\right)\sin\left(\frac{4\pi
z}{A\eta}\right)-\frac{2\pi^{3}\tilde{V_{l}}}{A^{2}\eta}\text{csch}\left(\frac{\pi^{2}}{A\eta}\right)\sin\left(\frac{2A\pi
z-2\nu\tau\eta}{A^{2}\eta}\right)-\tilde{\gamma}\eta^{2}\text{sech}^{2}(z)\text{tanh}(z)$
(9) $\displaystyle-$ $\displaystyle
a\eta{\tilde{\gamma}}^{3/2}\cos(\alpha)\text{sech}(z)\text{tanh}(z),$
$\displaystyle\frac{d}{d\tau}\alpha$ $\displaystyle=$
$\displaystyle\frac{1}{a}\cos(\alpha)\text{sech}(z)\eta{\tilde{\gamma}}^{3/2}+\frac{1}{2}a^{2}\tilde{\gamma}-\frac{d}{d\tau}\beta+\frac{\tilde{\gamma}^{2}}{2}+\frac{A^{2}\tilde{V}_{l}\tilde{\gamma}^{2}}{\pi^{2}+A^{2}\tilde{\gamma}^{2}}\cos\left(\frac{2\nu\tau}{A^{2}}\right)+\frac{A^{2}\tilde{V}_{s}\tilde{\gamma}^{2}}{4\pi^{2}+A^{2}\tilde{\gamma}^{2}},$
(10) $\displaystyle\frac{d}{d\tau}\beta$ $\displaystyle=$
$\displaystyle\frac{4\pi^{4}\tilde{V}_{s}}{A^{2}\eta^{2}}\cos\left(\frac{4\pi
z}{A\eta}\right)\text{coth}\left(\frac{2\pi^{2}}{A\eta}\right)\text{csch}\left(\frac{2\pi^{2}}{A\eta}\right)+\frac{8\pi^{3}\tilde{V}_{s}}{A^{2}\eta^{2}}\text{csch}\left(\frac{2\pi^{2}}{A\eta}\right)\sin\left(\frac{4\pi
z}{A\eta}\right)z$ (11) $\displaystyle+$
$\displaystyle\frac{\pi^{4}\tilde{V_{l}}}{A^{2}\eta^{2}}\cos\left(\frac{2A\pi
z-2\tau\nu\eta}{A^{2}\eta}\right)\text{coth}\left(\frac{\pi^{2}}{A\eta}\right)\text{csch}\left(\frac{\pi^{2}}{A\eta}\right)+\frac{2\pi^{3}\tilde{V}_{l}}{A^{2}\eta^{2}}\text{csch}\left(\frac{\pi^{2}}{A\eta}\right)\sin\left(\frac{2A\pi
z-2\tau\nu\eta}{A^{2}\eta}\right)z$ $\displaystyle+$
$\displaystyle\tilde{\gamma}\text{sech}^{2}(z)\eta+\frac{1}{2}\eta^{2}-\frac{\kappa^{2}}{2}+a{\tilde{\gamma}}^{3/2}\cos(\alpha)\text{sech}(z).$

Equation (7) provides the key analytical result of this work; it shows how the
motion of a soliton against an impurity is affected by the Thouless pump. Note
that without the pump (i.e., $\tilde{V}_{s}=\tilde{V}_{l}=0$), the above
equations recover the corresponding results of Ref. [54]. As shown in Eq. (7),
the displacement of the soliton is strongly affected by the time-dependent
superposition amplitudes of the soliton component $\eta(\tau)$ and the bound-
state component $a(\tau)$ in Eq. (8).

To solve the set of Euler-Lagrangian equations in Eqs. (7)-(11) for the
variational parameters, without loss of generality, we consider the following
initial conditions. At $\tau=0$, the bright soliton is centered at $z(0)=-3$,
far away from the impurity; the initial amplitude and velocity of the soliton
are chosen as $\eta=1$ and $\kappa=0.0258$, respectively; for other
parameters, we choose $a(0)=\alpha(0)=\beta(0)=0$.

To double-check whether our derivations reproduce the known results in Ref.
[54] in the absence of Thouless pumps, at first, we take $\tilde{V}_{s}=0$ and
$\tilde{V}_{l}=0$ in Eqs. (7)-(11). Indeed, as illustrated in Fig. 1 (c) for
the light impurity with $\tilde{\gamma}=0.0258$, the soliton directly transits
through the impurity potential without being localized. This can be
intuitively understood by noting that for small $\tilde{\gamma}$, the bound
state has a much smaller energy than the soliton, resulting in a small
admixture of the localized bound-state in Eq. (5). Therefore, the moving
soliton is only slightly perturbed by the impurity. With the increase of
$\tilde{\gamma}$, however, the bound state gains more weight in the
wavefunction described by Eq. (5). For $\tilde{\gamma}=0.258$ as in Fig. 2
(c), we see that the initially moving soliton becomes localized around the
impurity, consistent with previous findings.

Next, we turn on the Thouless pump with $\tilde{V}_{s}=\tilde{V}_{l}=25$, and
again solve Eqs. (7)-(11). The result for the light impurity is shown in Fig.
1 (d). We see that the soliton in the Thouless pump pass through a light
impurity as if the impurity does not exist, and its mean displacement in a
period becomes quantized. Remarkably, as shown in Fig. 2 (d), the nonlinear
topological pump enables the soliton to pass through even a very heavy
impurity, contrary to the pump-free case, and the motion is quantized. This
provides a direct evidence that the quantized transport of a quantum particle
in nonlinear topological Thouless pump is robust against the local
imperfections.

##  V Conclusion

In this work, we theoretically investigate the nonlinear Thouless pumping of a
soliton in the presence of an impurity based on the system of a quasi-1D BEC.
Using both the GP equation and the Lagrangian variational method, we compare
the interaction between the impurity and the soliton in the absence and
presence of the nonlinear topological pumps. Without the pump, whether the
soliton can move through the impurity depends crucially on the impurity mass.
In contrast, the soliton in the nonlinear Thouless pump can always transmit
through the impurity regardless, and exhibits a quantized motion. Our work
sheds light on the interplay of nonlinearity, topology and disorder, and
promises potential applications in information processing with solitons.

We thank Qidong Fu, Fangwei Ye, Qi Zhang, and Biao Wu for stimulating
discussions. This work was supported by the National Natural Science
Foundation of China (Nos. 12074344,12374246), the Zhejiang Provincial Natural
Science Foundation (Grant Nos. LZ21A040001) and the key projects of the
Natural Science Foundation of China (Grant No. 11835011).

## References

  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Quantized hall conductance in a two-dimensional periodic potential, Phys. Rev. Lett. 49, 405 (1982).
  * Thouless [1983] D. J. Thouless, Quantization of particle transport, Phys. Rev. B 27, 6083 (1983).
  * Niu and Thouless [1984] Q. Niu and D. J. Thouless, Quantised adiabatic charge transport in the presence of substrate disorder and many-body interaction, Journal of Physics A: Mathematical and General 17, 2453 (1984).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Berry phase effects on electronic properties, Rev. Mod. Phys. 82, 1959 (2010).
  * Wang  _et al._ [2013a] L. Wang, M. Troyer, and X. Dai, Topological charge pumping in a one-dimensional optical lattice, Phys. Rev. Lett. 111, 026802 (2013a).
  * Lohse  _et al._ [2016] M. Lohse, C. Schweizer, O. Zilberberg, M. Aidelsburger, and I. Bloch, A thouless quantum pump with ultracold bosonic atoms in an optical superlattice, Nat. Phys. 12, 350 (2016).
  * Nakajima  _et al._ [2016] S. Nakajima, T. Tomita, S. Taie, T. Ichinose, H. Ozawa, L. Wang, M. Troyer, and Y. Takahashi, Topological thouless pumping of ultracold fermions, Nat. Phys. 12, 296 (2016).
  * Lu  _et al._ [2016] H.-I. Lu, M. Schemmer, L. M. Aycock, D. Genkina, S. Sugawa, and I. B. Spielman, Geometrical pumping with a bose-einstein condensate, Phys. Rev. Lett. 116, 200402 (2016).
  * Lohse  _et al._ [2018] M. Lohse, C. Schweizer, H. M. Price, O. Zilberberg, and I. Bloch, Exploring 4d quantum hall physics with a 2d topological charge pump, Nature 553, 55 (2018).
  * Nakajima  _et al._ [2021] S. Nakajima, N. Takei, K. Sakuma, Y. Kuno, P. Marra, and Y. Takahashi, Competition and interplay between topology and quasi-periodic disorder in thouless pumping of ultracold atoms, Nat. Phys. 17, 844 (2021).
  * Minguzzi  _et al._ [2022] J. Minguzzi, Z. Zhu, K. Sandholzer, A.-S. Walter, K. Viebahn, and T. Esslinger, Topological pumping in a floquet-bloch band, Phys. Rev. Lett. 129, 053201 (2022).
  * Kraus  _et al._ [2012] Y. E. Kraus, Y. Lahini, Z. Ringel, M. Verbin, and O. Zilberberg, Topological states and adiabatic pumping in quasicrystals, Phys. Rev. Lett. 109, 106402 (2012).
  * Ke  _et al._ [2016] Y. Ke, X. Qin, F. Mei, H. Zhong, Y. S. Kivshar, and C. Lee, Topological phase transitions and thouless pumping of light in photonic waveguide arrays, Laser & Photonics Reviews 10, 995 (2016).
  * Zilberberg  _et al._ [2018] O. Zilberberg, S. Huang, J. Guglielmon, M. Wang, K. P. Chen, Y. E. Kraus, and M. C. Rechtsman, Photonic topological boundary pumping as a probe of 4d quantum hall physics, Nature 553, 59 (2018).
  * Cerjan  _et al._ [2020] A. Cerjan, M. Wang, S. Huang, K. P. Chen, and M. C. Rechtsman, Thouless pumping in disordered photonic systems, Light: Science & Applications 9, 178 (2020).
  * Wang  _et al._ [2022] P. Wang, Q. Fu, R. Peng, Y. V. Kartashov, L. Torner, V. V. Konotop, and F. Ye, Two-dimensional thouless pumping of light in photonic moiré lattices, Nat. Commun. 13, 6738 (2022).
  * Ma  _et al._ [2018] W. Ma, L. Zhou, Q. Zhang, M. Li, C. Cheng, J. Geng, X. Rong, F. Shi, J. Gong, and J. Du, Experimental observation of a generalized thouless pump with a single spin, Phys. Rev. Lett. 120, 120501 (2018).
  * Bastidas [2022] V. M. Bastidas, Topological thouless pumping in arrays of coupled spin chains, Phys. Rev. B 106, L220308 (2022).
  * Leykam and Chong [2016] D. Leykam and Y. D. Chong, Edge solitons in nonlinear-photonic topological insulators, Phys. Rev. Lett. 117, 143901 (2016).
  * Tangpanitanon  _et al._ [2016] J. Tangpanitanon, V. M. Bastidas, S. Al-Assam, P. Roushan, D. Jaksch, and D. G. Angelakis, Topological pumping of photons in nonlinear resonator arrays, Phys. Rev. Lett. 117, 213603 (2016).
  * Bisianov  _et al._ [2019] A. Bisianov, M. Wimmer, U. Peschel, and O. A. Egorov, Stability of topologically protected edge states in nonlinear fiber loops, Phys. Rev. A 100, 063830 (2019).
  * Ivanov  _et al._ [2020] S. K. Ivanov, Y. V. Kartashov, A. Szameit, L. Torner, and V. V. Konotop, Vector topological edge solitons in floquet insulators, ACS Photonics 7, 735 (2020).
  * Dai  _et al._ [2024] T. Dai, Y. Ao, J. Mao, Y. Yang, Y. Zheng, C. Zhai, Y. Li, J. Yuan, B. Tang, Z. Li, J. Luo, W. Wang, X. Hu, Q. Gong, and J. Wang, Non-hermitian topological phase transitions controlled by nonlinearity, Nat. Phys. 20, 101 (2024).
  * Solnyshkov  _et al._ [2017] D. D. Solnyshkov, O. Bleu, B. Teklu, and G. Malpuech, Chirality of topological gap solitons in bosonic dimer chains, Phys. Rev. Lett. 118, 023901 (2017).
  * González-Cuadra  _et al._ [2020] D. González-Cuadra, A. Dauphin, P. R. Grzybowski, M. Lewenstein, and A. Bermudez, Dynamical solitons and boson fractionalization in cold-atom topological insulators, Phys. Rev. Lett. 125, 265301 (2020).
  * Lo  _et al._ [2021] P.-W. Lo, C. D. Santangelo, B. G.-g. Chen, C.-M. Jian, K. Roychowdhury, and M. J. Lawler, Topology in nonlinear mechanical systems, Phys. Rev. Lett. 127, 076802 (2021).
  * Jürgensen  _et al._ [2021] M. Jürgensen, S. Mukherjee, and M. C. Rechtsman, Quantized nonlinear thouless pumping, Nature 596, 63 (2021).
  * Mostaan  _et al._ [2022] N. Mostaan, F. Grusdt, and N. Goldman, Quantized topological pumping of solitons in nonlinear photonics and ultracold atomic mixtures, Nat. Commun. 13, 5997 (2022).
  * Fu  _et al._ [2022a] Q. Fu, P. Wang, Y. V. Kartashov, V. V. Konotop, and F. Ye, Nonlinear thouless pumping: Solitons and transport breakdown, Phys. Rev. Lett. 128, 154101 (2022a).
  * Fu  _et al._ [2022b] Q. Fu, P. Wang, Y. V. Kartashov, V. V. Konotop, and F. Ye, Two-dimensional nonlinear thouless pumping of matter waves, Phys. Rev. Lett. 129, 183901 (2022b).
  * Jürgensen and Rechtsman [2022] M. Jürgensen and M. C. Rechtsman, Chern number governs soliton motion in nonlinear thouless pumps, Phys. Rev. Lett. 128, 113901 (2022).
  * Jürgensen  _et al._ [2023] M. Jürgensen, S. Mukherjee, C. Jörg, and M. C. Rechtsman, Quantized fractional thouless pumping of solitons, Nat. Phys. 19, 420 (2023).
  * Tuloup  _et al._ [2023] T. Tuloup, R. W. Bomantara, and J. Gong, Breakdown of quantization in nonlinear thouless pumping, New Journal of Physics 25, 083048 (2023).
  * Citro and Aidelsburger [2023] R. Citro and M. Aidelsburger, Thouless pumping and topology, Nat. Rev. Phys. 5, 87 (2023).
  * Xiang  _et al._ [2023] Z.-C. Xiang, K. Huang, Y.-R. Zhang, and et al, Simulating chern insulators on a superconducting quantum processor, Nat. Commun. 14, 5433 (2023).
  * Berg  _et al._ [2011] E. Berg, M. Levin, and E. Altman, Quantized pumping and topology of the phase diagram for a system of interacting bosons, Phys. Rev. Lett. 106, 110405 (2011).
  * Qian  _et al._ [2011] Y. Qian, M. Gong, and C. Zhang, Quantum transport of bosonic cold atoms in double-well optical lattices, Phys. Rev. A 84, 013608 (2011).
  * Wang  _et al._ [2013b] L. Wang, M. Troyer, and X. Dai, Topological charge pumping in a one-dimensional optical lattice, Phys. Rev. Lett. 111, 026802 (2013b).
  * Grusdt and Höning [2014] F. Grusdt and M. Höning, Realization of fractional chern insulators in the thin-torus limit with ultracold bosons, Phys. Rev. A 90, 053623 (2014).
  * Zeng  _et al._ [2016] T.-S. Zeng, W. Zhu, and D. N. Sheng, Fractional charge pumping of interacting bosons in one-dimensional superlattice, Phys. Rev. B 94, 235139 (2016).
  * Li and Fleischhauer [2017] R. Li and M. Fleischhauer, Finite-size corrections to quantized particle transport in topological charge pumps, Phys. Rev. B 96, 085444 (2017).
  * Ke  _et al._ [2017] Y. Ke, X. Qin, Y. S. Kivshar, and C. Lee, Multiparticle wannier states and thouless pumping of interacting bosons, Phys. Rev. A 95, 063630 (2017).
  * Taddia  _et al._ [2017] L. Taddia, E. Cornfeld, D. Rossini, L. Mazza, E. Sela, and R. Fazio, Topological fractional pumping with alkaline-earth-like atoms in synthetic lattices, Phys. Rev. Lett. 118, 230402 (2017).
  * Hayward  _et al._ [2018] A. Hayward, C. Schweizer, M. Lohse, M. Aidelsburger, and F. Heidrich-Meisner, Topological charge pumping in the interacting bosonic rice-mele model, Phys. Rev. B 98, 245148 (2018).
  * Nakagawa  _et al._ [2018] M. Nakagawa, T. Yoshida, R. Peters, and N. Kawakami, Breakdown of topological thouless pumping in the strongly interacting regime, Phys. Rev. B 98, 115147 (2018).
  * Stenzel  _et al._ [2019] L. Stenzel, A. L. C. Hayward, C. Hubig, U. Schollwöck, and F. Heidrich-Meisner, Quantum phases and topological properties of interacting fermions in one-dimensional superlattices, Phys. Rev. A 99, 053614 (2019).
  * Haug  _et al._ [2019] T. Haug, R. Dumke, L.-C. Kwek, and L. Amico, Topological pumping in aharonov–bohm rings, Commun. Phys. 2, 127 (2019).
  * Unanyan  _et al._ [2020] R. Unanyan, M. Kiefer-Emmanouilidis, and M. Fleischhauer, Finite-temperature topological invariant for interacting systems, Phys. Rev. Lett. 125, 215701 (2020).
  * Greschner  _et al._ [2020] S. Greschner, S. Mondal, and T. Mishra, Topological charge pumping of bound bosonic pairs, Phys. Rev. A 101, 053630 (2020).
  * Chen  _et al._ [2020] Y.-L. Chen, G.-Q. Zhang, D.-W. Zhang, and S.-L. Zhu, Simulating bosonic chern insulators in one-dimensional optical superlattices, Phys. Rev. A 101, 013627 (2020).
  * Pérez-García  _et al._ [1996] V. M. Pérez-García, H. Michinel, J. I. Cirac, M. Lewenstein, and P. Zoller, Low energy excitations of a bose-einstein condensate: A time-dependent variational analysis, Phys. Rev. Lett. 77, 5320 (1996).
  * Dalfovo  _et al._ [1999] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Theory of bose-einstein condensation in trapped gases, Rev. Mod. Phys. 71, 463 (1999).
  * Al Khawaja  _et al._ [2002] U. Al Khawaja, H. T. C. Stoof, R. G. Hulet, K. E. Strecker, and G. B. Partridge, Bright soliton trains of trapped bose-einstein condensates, Phys. Rev. Lett. 89, 200404 (2002).
  * Forinash  _et al._ [1994] K. Forinash, M. Peyrard, and B. Malomed, Interaction of discrete breathers with impurity modes, Phys. Rev. E 49, 3400 (1994).
  * Jia and Liang [2022] C. Jia and Z. Liang, Interaction between an impurity and nonlinear excitations in a polariton condensate, Entropy 24, 1789 (2022).
  * Fallani  _et al._ [2007] L. Fallani, J. E. Lye, V. Guarrera, C. Fort, and M. Inguscio, Ultracold atoms in a disordered crystal of light: Towards a bose glass, Phys. Rev. Lett. 98, 130404 (2007).
  * Chomaz  _et al._ [2022] L. Chomaz, I. Ferrier-Barbut, F. Ferlaino, B. Laburthe-Tolra, B. L. Lev, and T. Pfau, Dipolar physics: a review of experiments with magnetic quantum gases, Reports on Progress in Physics 86, 026401 (2022).
  * Łacki  _et al._ [2016] M. Łacki, M. A. Baranov, H. Pichler, and P. Zoller, Nanoscale “dark state” optical potentials for cold atoms, Phys. Rev. Lett. 117, 233001 (2016).
  * Wang  _et al._ [2018] Y. Wang, S. Subhankar, P. Bienias, M. Łacki, T.-C. Tsui, M. A. Baranov, A. V. Gorshkov, P. Zoller, J. V. Porto, and S. L. Rolston, Dark state optical lattice with a subwavelength spatial structure, Phys. Rev. Lett. 120, 083601 (2018).