# Adapting Brain-Like Neural Networks for Modeling Cortical Visual Prostheses

Jacob Granley
Department of Computer Science
University of California Santa Barbara
<EMAIL_ADDRESS>
&Alexander Riedel
Ernst-Abbe-University, Jena, Germany
<EMAIL_ADDRESS>
&Michael Beyeler
Department of Computer Science
Department of Psychology and Brain Sciences
University of California Santa Barbara
<EMAIL_ADDRESS>
###### Abstract

Cortical prostheses are devices implanted in the visual cortex that attempt to
restore lost vision by electrically stimulating neurons. Currently, the vision
provided by these devices is limited, and accurately predicting the visual
percepts resulting from stimulation is an open challenge. We propose to
address this challenge by utilizing ‘brain-like’ convolutional neural networks
(CNNs), which have emerged as promising models of the visual system. To
investigate the feasibility of adapting brain-like CNNs for modeling visual
prostheses, we developed a proof-of-concept model to predict the perceptions
resulting from electrical stimulation. We show that a neurologically-inspired
decoding of CNN activations produces qualitatively accurate phosphenes,
comparable to phosphenes reported by real patients. Overall, this is an
essential first step towards building brain-like models of electrical
stimulation, which may not just improve the quality of vision provided by
cortical prostheses but could also further our understanding of the neural
code of vision.

##  1 Introduction

Visual neuroprostheses are emerging as a promising treatment option for
restoring visual function lost to injury or disease. Analogous to cochlear
implants, cortical prostheses electrically stimulate neurons in the early
visual system, typically in the primary (V1) or secondary visual cortex (V2),
to elicit neuronal responses that the brain interprets as visual percepts
(‘phosphenes’) [11]. Several devices are currently in development, such as the
Orion Visual Cortical Prosthesis System [1], which is based on relatively
large surface electrodes, and CORTIVIS [12, 13], which is based on the Utah
array of intracortical penetrating electrodes.

Despite recent technological advances, the vision restored by these devices
remains restricted to white-ish or yellow-ish phosphenes of simple geometric
shape [6, 13, 1]. The inability to generate more complex visual patterns may
at least be partially due to our limited understanding of how visual
prostheses interact with the human visual system to shape perception. Current
devices might indiscriminately stimulate diverse populations of cortical
neurons, each with their own complex neuronal response properties that are
often not fully understood [8]. Constructing accurate phosphene models, which
predict the appearance of percepts elicited from electrical stimuli, directly
from patient data remains one of the biggest challenges in the field,
primarily due to the difficulty of data acquisition, the amount of noise in
perceptual measures, and patient-to-patient variability.

Alternatively, convolutional neural networks (CNNs) have emerged as promising
models of the visual system, allowing the simulation of the natural single-
unit cortical response to ecological visual stimuli. To this end, “Brain-
Score” [21, 22] provides a comprehensive benchmark that measures how closely
activation of intermediate layers in CNNs trained for object classification
corresponds to primate cortical and behavioral responses across multiple
datasets. If adapted properly, these neural networks could be powerful tools
for understanding the neuronal and perceptual effects of electrical
stimulation on the visual cortex.

Here, we investigate the feasibility of adapting brain-like deep neural
networks (DNNs) for modeling visual cortical prostheses. We make the following
contributions:

  * •

We develop techniques facilitating the use of brain-like DNNs trained on
object classification to predict visual responses to electrical stimulation.
We highlight several key challenges associated with this process, such as
strategies for mapping stimulation intensities and anatomical locations to DNN
neurons.

  * •

We implement a proof-of-concept phosphene model and compare its predictions to
phosphenes reported by real patients. We show that many of the qualitative
features that characterize clinically observed phosphenes naturally emerge
from a neurologically-inspired decoding of predicted neuronal responses, such
as yellow-white, circular phosphenes, and ‘shimmering’ or ‘fireworks’-like
qualities of phosphenes. In addition, our model reproduces observed and
expected trends in phosphene size.

  * •

We discuss the unique opportunities that brain-like phosphene models offer,
which could have implications both in neuroscience (i.e., illuminating the
decoding methods the brain uses to interpret unusual activation patterns) and
in machine learning (i.e., benchmarks measuring DNN’s ability to reproduce the
cortical behaviors observed in prosthesis patients).

##  2 Methods

In this section, we discuss a number of the challenges associated with
adapting an object classification DNN to be used for modeling of electrical
stimulation. To the best of our knowledge, we are the first to attempt to
model cortical prostheses with brain-like DNNs. Therefore, we also describe
many potential solution techniques, and motivate the choices used in our
phosphene model.

### Cortical Prostheses

We simulated electrical stimulation with two cortical implants: Orion [1] and
CORTIVIS [13]. The Orion implant has 60 large surface electrodes
($1\text{\,}\mathrm{mm}$ diameter) spaced $4.2\text{\,}\mathrm{mm}$ apart
horizontally, and $3\text{\,}\mathrm{mm}$ apart diagonally. CORTIVIS is based
on the Utah array, which consists of 96 small penetrating electrodes on a
square grid spaced $0.4\text{\,}\mathrm{mm}$ apart.

### Choice of Brain-Like Model

We used the `effnetb1` network [20], which is currently the best performing
model on brain-score.org [21, 22]. We chose this over other networks that
directly predict neural activations due to the correspondence of multiple
cortical and artificial layers (V1, V2, V4, and IT) and its thorough
evaluation across multiple electrophysiological and behavioral datasets. The
network is an EfficientNet B1, consisting of 7.8M parameters, with
adversarially-robust training as described in [20]. We focused on the layers
corresponding to cortical areas V1 and V2 (`blocks.3.0` and `blocks.3`), since
these are the most relevant sites for cortical prostheses.

### Mapping Artificial Neurons to Anatomical Locations

Mapping electrical stimulation patterns to DNN activations is a primary
challenge associated with our approach. One obvious obstacle is that there is
no 1:1 correspondence between artificial and biological neurons. In other
words, the corresponding anatomical distribution of the artificial neurons
across the cortical surface is unknown.

Previous works have successfully used a wiring length constraint to enforce
topographic organization of IT neurons [4, 18], but it is unclear how this
might extend to other cortical layers. Further, a method that infers
anatomical location after training of the DNN is preferable for its
compatibility with existing models. One such method would be to effectively
‘reverse’ the Brain-Score metrics by training a regression model to associate
measured cortical responses with artificial neuron activations. While
promising, this approach is limited by the specific cortical dataset and
regression model used, which might not generalize well in the absence of
patient-specific cortical responses.

Figure 1: Mapping artificial neurons from a brain-like DNN to anatomical
locations in visual cortex. _A_ : Examples of stimuli that maximally activate
V1 neurons. Each neuron’s receptive field (RF)center is defined as the center
of mass (marked in red). _B_ : Inferred RF centers of all 51,840 V1 and V2
artificial neurons in visual field coordinates (degrees of visual angle; dva).
_C_ : The corresponding anatomical locations in V1 and V2, under the Wedge-
Dipole visuotopic mapping [19], vertical line marking the discontinuous
boundary between hemispheres. The 4 clusters (corresponding to channels) of
neurons boxed in _B_ correspond to the circled clusters in _C_.

We therefore propose a new approach suited for the general patient, which can
be applied after training of the DNN, and can be made patient-specific without
requiring neuronal responses, as illustrated in Fig. 1. In this approach,
artificial neurons are mapped to the visual field using an activation
maximization technique, and then projected to the cortical surface using well-
established visuotopic mappings. First, gradient descent was used to find the
stimuli that maximally activate each artificial neuron [24], similar to
feature visualization techniques [9]. We used the Adam optimizer with an
initial learning rate of 0.1, decreasing by a factor of 5 whenever the loss
failed to decrease for 10 iterations. We also enforced L2 regularization on
the stimuli, encouraging irrelevant areas in the stimulus to be dark. The vast
majority of resulting ‘maximal stimuli’ were found to be circular and well
localized (Fig. 1a). These receptive fields (RFs) were slightly larger than
primate V1 and V2 RFs, but with an increase in size (33%) from V1 and V2 that
is similar to primate cortex [23]. We also input each stimulus back into the
network to verify that it selectively activated only the corresponding
artificial neuron. Each stimulus’s center of mass was used as the neuron’s
location in visual field (Fig. 1B). To convert this to a cortical location
(Fig. 1C), we used a general Wedge-Dipole mapping [19]—but note that if fMRI
data were available, then recent patient-specific mapping techniques (e.g.,
[2]) would serve as suitable replacements. The RF centers are nearly uniformly
spread across the visual field, but are disproportionately sparse in the
central foveal cortex region due to cortical magnification.

### Cortical Activation Patterns

Previous work has demonstrated that cortical electrodes sparsely activate
neurons within a certain radius of electrodes, sometimes as far as millimeters
away [16]. We therefore modeled cortical activation patterns as Gaussians
centered on each stimulating electrode, with standard deviation ($\rho$) as an
adjustable parameter, allowing us to account for a range of possible implants.
The stimulation intensity $I$ at cortical location $\mathbf{x}$ in mm was
given by

$I(\mathbf{x};\rho)=\sum_{e\in
E}\,a_{e}\exp\left(\frac{-||\mathbf{x}-\mathbf{e}||_{2}^{2}}{2\rho^{2}}\right),$
(1)

where $E$ is the set of electrodes, $\mathbf{e}$ denotes the location of each
electrode, and $a_{e}$ is the amplitude of each electrode. We set $\rho$ equal
to the electrode spacing of the two cortical prostheses studied (i.e.,
$4.2\text{\,}\mathrm{mm}$ for Orion and $0.4\text{\,}\mathrm{mm}$ for
CORTIVIS), and otherwise treated stimulation as identical.

### Electrical Stimulation and Phosphene Model

The final challenge is mapping cortical activation patterns to the artificial
neurons and reconstructing the resulting phosphenes. One straightforward
solution motivated by previous literature [3, 5, 14] would be to model current
spread on the cortical surface, and ‘copy and paste’ the current intensity
onto the corresponding DNN activation maps. However, what it means for an
artificial neuron to be ‘activated’ differs from biological neurons (e.g.,
artificial neurons can be positively or negatively activated, each neuron is
scaled differently), and is dependent on network architecture and each layer’s
activation function. Further we empirically found that artificial activation
patterns from this technique did not produce realistic phosphenes.

We instead propose a neurologically-inspired linear decoding technique that is
not architecture-dependent (Fig. 2). To visualize the phosphene ($P$)
resulting from stimulation, we average all of the previously computed maximal
stimuli ($S$), weighted by the stimulation intensity ($I$) at each neuron’s
corresponding location in V1 or V2, normalized by activity in a local
neighborhood:

$P_{i}=\frac{\sum_{j}I_{j}*S_{j,i}}{\sum_{j\in\mathcal{N}(i)}I_{j}},$ (2)

where the index $i$ iterates over pixels, $j$ iterates over different DNN
neurons, and $\mathcal{N}(i)$ denotes all DNN neurons within the local
neighborhood of $i$. We found that limiting each maximal stimulus to the 99th
percentile of measured RF sizes produced high quality reconstructions.

This decoding strategy is consistent with the idea that each neuron’s
activation level is indicative of the likelihood that said neuron encountered
its preferred stimulus; the resulting percept is thus a weighted sum of
preferred stimuli [17]. Normalization by activity within a local neighborhood
bears some resemblance to divisive normalization [7]. A benefit of this
strategy is its uniform treatment of different cortical areas, allowing for
seamless modeling of implants on the border of V1 and V2.

Figure 2: Simulating electrical stimulation of early visual cortex to produce
phosphenes. To obtain the predicted phosphene for a electrode activation
pattern (_A_), the stimulation intensity at each artificial neuron is computed
using Eq. 1. Shown is the stimulation intensity across the neurons in V1 (_B_)
and the same neurons in visual field (_C_). The intensities at each neuron are
then used as the weight of the corresponding neuron’s maximally activating
stimulus. The weighted average of all stimuli is output as the final predicted
percept (_D_), which spans $16\times 16$ degrees of visual angle.

##  3 Results

Fig. 3 shows representative examples of single-electrode phosphenes for Orion
and CORTIVIS implanted parafoveally in V1. Prediction for a single phosphene
took approximately 3s on CPU. In general, phosphenes were round yellow-white
blobs as expected, but with some subtle high frequency details, such as
intermittent patches of bright colors. While predicted Orion phosphenes tended
to be highly saturated and large blobs, predicted CORTIVIS phosphenes appeared
sparser, smaller, and slightly more colorful. When we simulated the two
implants at the same visual field location in V2, phosphenes were
qualitatively similar, but 41% larger on average. Due to cortical
magnification, electrodes near the fovea may produce smaller phosphenes than
electrodes near the periphery. This relationship was linear and strong in
Orion for V1 and V2 (t-test, $n=90$, $p<.001$) matching previous studies [5],
but not for the CORTIVIS ($p=0.803$).

Figure 3: Predicted phosphenes for Orion (_A_) and CORTIVIS (_B_). _C_ :
Predicted phosphene size in V1 to V2 (implant location not shown) increases
with eccentricity. _D_ : Predicted phosphene size (measured as radius of the
phosphene’s convex hull) in V1 as a function of current spread ($\rho$).

While it is well known that phosphene size saturates with increasing current
spread [5], it is less clear whether phosphenes can become arbitrarily small.
In theory, the smallest phosphene elicited by stimulating a single neuron
should be limited by that neuron’s RF size (because the neuron cannot sense
which subregion of its RF was stimulated). Consistent with these theoretical
considerations, our model predicts phosphene size to saturate both with high
and low current spreads (Fig. 3D).

##  4 Discussion

In this work we adapt a brain-like DNN [20] as a decoding model for visual
cortical prostheses. We show that our model makes a number of predictions that
agree well with the existing neuroscience literature. First, predicted
phosphenes are yellow-ish white blobs. Notably, this is not a manually
designed feature such as in other models, but arises naturally from the
proposed decoding scheme. Second, predicted phosphenes have high-frequency
details, such as ‘shimmering’ or ‘fireworks’ (patches of bright colors) that
qualitatively match many patients’ descriptions of observed phosphenes [10].
Third, the model reproduces trends in phosphene size that have either been
observed in cortical prostheses [5] or are in line with our theoretical
understanding of neuronal population codes.

Our model generally produces similar phosphenes to existing, simpler models
(e.g. linear-nonlinear models [6]). However, rather than being hard-coded,
these features naturally emerge from our population decoding scheme. Our
model’s ability to reproduce some of the subtler details of phosphenes
inspires hope that this approach might serve as a reliable phosphene model for
a range of experiments, enabling improvements in cortical prostheses. For
example, if the phosphenes are accurate, then a similar brain-like phosphene
model could be used in an end-to-end stimulus optimization framework [15] to
discover superior stimulus encoding algorithms.

This study is inherently preliminary as it is limited to a single DNN and has
yet to be comprehensively validated against a large dataset of real phosphenes
(which currently does not exist). Further, it must be demonstrated that our
results are not just due to peculiarities of the network presented in [20],
though our preliminary investigations suggest that other brain-like models
make similar predictions.

The intersection of cortical prostheses and DNNs offers a unique set of
opportunities for advancing our knowledge of neuroscience. On the one hand,
visual prostheses offer a safe way to stimulate and record from visual cortex
of awake humans, which may serve as a meaningful benchmark dataset for ‘brain-
like‘ DNNs. On the other hand, building brain-like models of electrical
stimulation may not only improve the quality of vision provided by cortical
prostheses but also offer a unique testbed to further our understanding of the
neural code of vision.

## References

  * Beauchamp et al. [2020] Michael S. Beauchamp, Denise Oswalt, Ping Sun, Brett L. Foster, John F. Magnotti, Soroush Niketeghad, Nader Pouratian, William H. Bosking, and Daniel Yoshor.  Dynamic Stimulation of Visual Cortex Produces Form Vision in Sighted and Blind Humans.  _Cell_ , 181(4):774–783.e5, May 2020.  ISSN 00928674.  doi: 10.1016/j.cell.2020.04.033.  URL https://linkinghub.elsevier.com/retrieve/pii/S0092867420304967.
  * Benson et al. [2014] Noah C. Benson, Omar H. Butt, David H. Brainard, and Geoffrey K. Aguirre.  Correction of Distortion in Flattened Representations of the Cortical Surface Allows Prediction of V1-V3 Functional Organization from Anatomy.  _PLOS Computational Biology_ , 10(3):e1003538, March 2014.  ISSN 1553-7358.  doi: 10.1371/journal.pcbi.1003538.  URL https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003538.  Publisher: Public Library of Science.
  * Beyeler et al. [2019] Michael Beyeler, Devyani Nanduri, James D. Weiland, Ariel Rokem, Geoffrey M. Boynton, and Ione Fine.  A model of ganglion axon pathways accounts for percepts elicited by retinal implants.  _Scientific Reports_ , 9(1):1–16, June 2019.  ISSN 2045-2322.  doi: 10.1038/s41598-019-45416-4.  URL https://www.nature.com/articles/s41598-019-45416-4.
  * Blauch et al. [2022] Nicholas M. Blauch, Marlene Behrmann, and David C. Plaut.  A connectivity-constrained computational account of topographic organization in primate high-level visual cortex.  _Proceedings of the National Academy of Sciences_ , 119(3):e2112566119, January 2022.  ISSN 0027-8424, 1091-6490.  doi: 10.1073/pnas.2112566119.  URL https://pnas.org/doi/full/10.1073/pnas.2112566119.
  * Bosking et al. [2017] William H. Bosking, Ping Sun, Muge Ozker, Xiaomei Pei, Brett L. Foster, Michael S. Beauchamp, and Daniel Yoshor.  Saturation in Phosphene Size with Increasing Current Levels Delivered to Human Visual Cortex.  _Journal of Neuroscience_ , 37(30):7188–7197, July 2017.  ISSN 0270-6474, 1529-2401.  doi: 10.1523/JNEUROSCI.2896-16.2017.  URL https://www.jneurosci.org/content/37/30/7188.  Publisher: Society for Neuroscience Section: Research Articles.
  * Bosking et al. [2022] William H. Bosking, Denise N. Oswalt, Brett L. Foster, Ping Sun, Michael S. Beauchamp, and Daniel Yoshor.  Percepts evoked by multi-electrode stimulation of human visual cortex.  _Brain Stimulation_ , 15(5):1163–1177, September 2022.  ISSN 1935-861X.  doi: 10.1016/j.brs.2022.08.007.  URL https://www.sciencedirect.com/science/article/pii/S1935861X22001851.
  * Carandini & Heeger [2012] Matteo Carandini and David J. Heeger.  Normalization as a canonical neural computation.  _Nature Reviews Neuroscience_ , 13(1):51–62, January 2012.  ISSN 1471-0048.  doi: 10.1038/nrn3136.  URL https://www.nature.com/articles/nrn3136.  Number: 1 Publisher: Nature Publishing Group.
  * Carandini et al. [2005] Matteo Carandini, Jonathan B. Demb, Valerio Mante, David J. Tolhurst, Yang Dan, Bruno A. Olshausen, Jack L. Gallant, and Nicole C. Rust.  Do We Know What the Early Visual System Does?  _Journal of Neuroscience_ , 25(46):10577–10597, November 2005.  ISSN 0270-6474, 1529-2401.  doi: 10.1523/JNEUROSCI.3726-05.2005.  URL http://www.jneurosci.org/content/25/46/10577.
  * Erhan et al. [2009] Dumitru Erhan, Y. Bengio, Aaron Courville, and Pascal Vincent.  Visualizing Higher-Layer Features of a Deep Network.  _Technical Report, Univeristé de Montréal_ , January 2009.
  * Erickson-Davis & Korzybska [2021] Cordelia Erickson-Davis and Helma Korzybska.  What do blind people “see” with retinal prostheses? Observations and qualitative reports of epiretinal implant users.  _PLoS ONE_ , 16(2):e0229189, February 2021.  ISSN 1932-6203.  doi: 10.1371/journal.pone.0229189.  URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7875418/.
  * Fernandez [2018] Eduardo Fernandez.  Development of visual Neuroprostheses: trends and challenges.  _Bioelectronic Medicine_ , 4(1):12, August 2018\.  ISSN 2332-8886.  doi: 10.1186/s42234-018-0013-8.  URL https://doi.org/10.1186/s42234-018-0013-8.
  * Fernández & Normann [2017] Eduardo Fernández and Richard A. Normann.  CORTIVIS Approach for an Intracortical Visual Prostheses.  In Veit Peter Gabel (ed.), _Artificial Vision: A Clinical Guide_ , pp. 191–201. Springer International Publishing, Cham, 2017.  ISBN 978-3-319-41876-6.  doi: 10.1007/978-3-319-41876-6_15.  URL https://doi.org/10.1007/978-3-319-41876-6_15.
  * Fernández et al. [2021] Eduardo Fernández, Arantxa Alfaro, Cristina Soto-Sánchez, Pablo González-López, Antonio M. Lozano Ortega, Sebastian Peña, María Dolores Grima, Alfonso Rodil, Bernardeta Gómez, Xing Chen, Pieter R. Roelfsema, John D. Rolston, Tyler S. Davis, and Richard A. Normann.  Visual percepts evoked with an Intracortical 96-channel microelectrode array inserted in human occipital cortex.  _The Journal of Clinical Investigation_ , October 2021.  ISSN 0021-9738.  doi: 10.1172/JCI151331.  URL https://www.jci.org/articles/view/151331.  Publisher: American Society for Clinical Investigation.
  * Granley & Beyeler [2021] Jacob Granley and Michael Beyeler.  A Computational Model of Phosphene Appearance for Epiretinal Prostheses.  In _2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC)_, pp. 4477–4481, November 2021.  doi: 10.1109/EMBC46164.2021.9629663.  ISSN: 2694-0604.
  * Granley et al. [2022] Jacob Granley, Lucas Relic, and Michael Beyeler.  A Hybrid Neural Autoencoder for Sensory Neuroprostheses and Its Applications in Bionic Vision, May 2022.  URL http://arxiv.org/abs/2205.13623.  arXiv:2205.13623 [cs].
  * Histed et al. [2009] M. H. Histed, V. Bonin, and R. C. Reid.  Direct activation of sparse, distributed populations of cortical neurons by electrical microstimulation.  _Neuron_ , 63(4):508–22, August 2009.  ISSN 1097-4199 (Electronic) 0896-6273 (Linking).  doi: 10.1016/j.neuron.2009.07.016.
  * Jazayeri & Movshon [2006] M. Jazayeri and J. A. Movshon.  Optimal representation of sensory information by neural populations.  _Nat Neurosci_ , 9(5):690–6, May 2006.  ISSN 1097-6256 (Print) 1097-6256 (Linking).  doi: 10.1038/nn1691.
  * Lee et al. [2020] Hyodong Lee, Eshed Margalit, Kamila M. Jozwik, Michael A. Cohen, Nancy Kanwisher, Daniel L. K. Yamins, and James J. DiCarlo.  Topographic deep artificial neural networks reproduce the hallmarks of the primate inferior temporal cortex face processing network.  preprint, Neuroscience, July 2020.  URL http://biorxiv.org/lookup/doi/10.1101/2020.07.09.185116.
  * Polimeni et al. [2006] J. R. Polimeni, M. Balasubramanian, and E. L. Schwartz.  Multi-area visuotopic map complexes in macaque striate and extra-striate cortex.  _Vision Research_ , 46(20):3336–3359, October 2006.  ISSN 0042-6989.  doi: 10.1016/j.visres.2006.03.006.  URL https://www.sciencedirect.com/science/article/pii/S0042698906001428.
  * Riedel [2022] Alexander Riedel.  Bag of Tricks for Training Brain-Like Deep Neural Networks.  March 2022.  URL https://openreview.net/forum?id=SudzH-vWQ-c.
  * Schrimpf et al. [2018] Martin Schrimpf, Jonas Kubilius, Ha Hong, Najib J. Majaj, Rishi Rajalingham, Elias B. Issa, Kohitij Kar, Pouya Bashivan, Jonathan Prescott-Roy, Franziska Geiger, Kailyn Schmidt, Daniel L. K. Yamins, and James J. DiCarlo.  Brain-Score: Which Artificial Neural Network for Object Recognition is most Brain-Like?  preprint, Neuroscience, September 2018.  URL http://biorxiv.org/lookup/doi/10.1101/407007.
  * Schrimpf et al. [2020] Martin Schrimpf, Jonas Kubilius, Michael J. Lee, N. Apurva Ratan Murty, Robert Ajemian, and James J. DiCarlo.  Integrative Benchmarking to Advance Neurally Mechanistic Models of Human Intelligence.  _Neuron_ , 108(3):413–423, November 2020.  ISSN 0896-6273.  doi: 10.1016/j.neuron.2020.07.040.  URL https://www.cell.com/neuron/abstract/S0896-6273(20)30605-X.  Publisher: Elsevier.
  * Van den Bergh et al. [2010] Gert Van den Bergh, Bin Zhang, Lutgarde Arckens, and Yuzo M. Chino.  Receptive-field Properties of V1 and V2 Neurons in Mice and Macaque monkeys.  _The Journal of comparative neurology_ , 518(11):2051–2070, June 2010.  ISSN 0021-9967.  doi: 10.1002/cne.22321.  URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2881339/.
  * Walker et al. [2019] Edgar Y. Walker, Fabian H. Sinz, Erick Cobos, Taliah Muhammad, Emmanouil Froudarakis, Paul G. Fahey, Alexander S. Ecker, Jacob Reimer, Xaq Pitkow, and Andreas S. Tolias.  Inception loops discover what excites neurons most using deep predictive models.  _Nature Neuroscience_ , pp. 1–6, November 2019.  ISSN 1546-1726.  doi: 10.1038/s41593-019-0517-x.  URL https://www.nature.com/articles/s41593-019-0517-x.


# Structural Inference of Networked Dynamical Systems with Universal
Differential Equations

J. Koch<EMAIL_ADDRESS>Z. Chen  A. Tuor  J. Drgona  D. Vrabie  Pacific
Northwest National Laboratory, Richland, WA

###### Abstract

Networked dynamical systems are common throughout science in engineering;
e.g., biological networks, reaction networks, power systems, and the like. For
many such systems, nonlinearity drives populations of identical (or near-
identical) units to exhibit a wide range of nontrivial behaviors, such as the
emergence of coherent structures (e.g., waves and patterns) or otherwise
notable dynamics (e.g., synchrony and chaos). In this work, we seek to infer
(i) the intrinsic physics of a base unit of a population, (ii) the underlying
graphical structure shared between units, and (iii) the coupling physics of a
given networked dynamical system given observations of nodal states. These
tasks are formulated around the notion of the Universal Differential Equation,
whereby unknown dynamical systems can be approximated with neural networks,
mathematical terms known a priori (albeit with unknown parameterizations), or
combinations of the two. We demonstrate the value of these inference tasks by
investigating not only future state predictions but also the inference of
system behavior on varied network topologies. The effectiveness and utility of
these methods is shown with their application to canonical networked nonlinear
coupled oscillators.

††preprint: APS/123-QED

> Data-driven modeling of dynamical systems has experienced a surge in method
> development in concert with advances in machine learning and artificial
> intelligence. In this work, we restrict our scope specifically to the
> problem of data-driven modeling of networked systems; a challenging problem
> in which one needs to elicit not only the intrinsic physics of individual
> nodes of the network, but also what nodes talk to whom and how that
> communication influences nodal behaviors. The implication of such
> methodology is far-reaching – one can begin to infer properties of networked
> systems with respect to network topology and/or external perturbations,
> answering hypotheticals such as “what would happen if we remove this node?”

##  I Introduction

Networked dynamical systems, i.e., those described by the evolution of
differential equations on a graph, are ubiquitous in science and engineering.
When applied in the context of modeling dynamics of a population of units
(e.g., a collection of fireflies or neural field activity), such systems
exhibit a wide variety of population-wide complex behaviors including
synchrony, pattern formation, chaos, and characteristic transitions between
these states.

Presently, we are concerned with the task of inferring (i) the physics of the
individual units that make up the broader population, (ii) the underlying
graphical structure of the units, and (iii) the coupling physics of a given
networked dynamical system given nodal observations. Knowledge of these
physics can allow for inference of system properties beyond future state
predictions, including system-level stability and the emergence of coherent
structures. Furthermore, in disambiguating between unit physics, coupling
physics, and graphical interactions, one can project learned physics onto new,
previously unexplored settings. For example, such tasks might include
investigating system stability with respect to node deletions or other
topological changes to the underlying graph.

In Section II, a brief overview of networked dynamical systems is presented in
the context of state-of-the-art data-driven modeling efforts. In Section III,
the methodology for the present work is laid out. Numerical experiments are
presented in Section IV and is followed by discussion in Section V.

\begin{overpic}[width=433.62pt]{img/fig_01.png} \put(1.0,44.0){(a)}
\put(67.0,44.0){(b)} \end{overpic} Figure 1: In (a), a small network of
unknown oscillators interact according to an adjacency matrix shown in (b). In
this pictoral example, the nodes are assumed to be second-order oscillators
with a directed graph structure, $A_{i,j}$. $\operatorname{NN}_{1}$
approximates the intrinsic nodal physics and $\operatorname{NN}_{2}$
approximates the coupling physics. The adjacency matrix in (b) is used in
Section IV.2.

##  II Background

Networked dynamical systems is an area of intense research in engineering and
the sciences. Many physical systems feature a collection of individual units,
or nodes, possessing common intrinsic physics. The manner in which these nodes
communicate - both how nodes communicate and with whom information is shared
or directed - is of great interest, as these physics can inform system design
and control.

Examples in the physical sciences and engineering range from biological
networks and pattern formation in neuroscience Earn, Rohani, and Grenfell
(1998); Klemm and Bornholdt (2005) to chemical reaction networks Ji and Deng
(2021); Owoyele and Pal (2022) to thermoacoustics of combustors Dange  _et
al._ (2019); Premraj  _et al._ (2021); Guan  _et al._ (2022); Manoj  _et al._
(2019). Coupled dynamical systems have also been used successfully to model
social and epidemiological networks Bakhtiarnia, Fahim, and Miandoab (2020);
Keeling and Eames (2005). Despite existing in disparate domains with distinct
physics, these systems all can exhibit common qualitative behaviors like
synchrony, pattern formation, and characteristic bifurcations between these
states Xiao  _et al._ (2008); Strogatz (2001); Arenas  _et al._ (2008); Li and
Chen (2006).

The Kuramoto family of coupled oscillators Kuramoto  are the simplest
nonlinear oscillator exhibiting many of these qualitative behaviors of
physical systems. The generic Kuramoto oscillator possesses an assumed natural
frequency $\omega$ and an all-to-all sinusoidal coupling:

$\frac{d\theta_{i}}{dt}=\omega+\frac{K}{N}\sum_{j=1}^{N}{\sin\left(\theta_{i}-\theta_{j}\right)}.$
(1)

The utility of the Kuramoto model is enshrined in its simplicity: basic
modifications to the network topology (encoded through an adjacency matrix)
and coupling term can produce simple, interpretable models for chaos,
turbulence, chimeras, and bifurcations between behaviors on an application-by-
application basis. As a phase oscillator, the Kuramoto oscillator lacks the
ability to capture higher-order dynamics. Nevertheless, the general
prescription for a graph-coupled nonlinear dynamical system follows from that
of the Kuramoto model:

$\frac{d\mathbf{x}_{i}}{dt}=\mathbf{f}(\mathbf{x}_{i})+\sum_{j\neq
i}^{N}\mathbf{A}_{i,j}\mathbf{g}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right),$
(2)

where $\mathbf{x}$ are the node states, $\mathbf{f}$ is the node physics,
$\mathbf{g}$ is the coupling physics, and the adjacency matrix
$\mathbf{A}_{i,j}$ directs the communication between nodes.

###  II.1 Related Work

Many methods exist in data science and engineering to model nonlinear
dynamical systems directly from observations. These methods range from so-
called “black-box” approaches, where interpretability and/or physicality are
not of concern, to physics-informed methods, where one uses prior knowledge of
a system to construct, constrain, or otherwise regularize an optimization or
regression task. Specifically for modeling dynamical systems, such methods
include black-box neural networks for time-stepping (e.g., a multi-layer
perceptron or recurrent neural network) Skomski  _et al._ (2021), symbolic
regression via a sparse selection of candidate functions Brunton, Proctor, and
Kutz (2016), or Neural Ordinary Differential Equations (NODEs) Chen  _et al._
(2018). In the latter, unknown dynamical systems can be approximated with
neural networks, mathematical terms known a priori \- potentially
parameterized, or combinations of the two. Incorporating prior knowledge in
NODEs – e.g., a known term in the right-hand side (RHS) of the differential
equation – can bridge the gap between “black-box” modeling paradigms and from-
first-principles forward modeling Rackauckas  _et al._ (2020).

In the broader context of learning dynamics on graphs, Graph Neural Networks
(GNNs) Scarselli  _et al._ (2009); Poli  _et al._ (2019a); Rusch  _et al._
(2022); Asikis, Böttcher, and Antulov-Fantulin (2022) aim to relate node
features to those of their adjoining neighbors through diffusive processes
(e.g. message passing Duvenaud  _et al._ (2015); Gilmer  _et al._ (2017)) in
an end-to-end differentiable algorithm. GNNs can process graph-structured data
of arbitrary size and shape, making them very attractive for applications in
physical chemistry and related fields. Although GNNs typically operate over
static graphs, nodal hidden states are dynamic. These dynamics can be captured
with traditional time-stepping models leading to natural extension of
continuous graph neural networks (CGNNs) Xhonneux, Qu, and Tang (2020) also
referred to as graph neural ordinary differential equations (GDEs) Poli  _et
al._ (2019b). Additionally, inference of a system’s underlying graph can be
accomplished within the GNN framework through generative models, i.e., graph
variational autoencoders Kipf and Welling (2016).

The evolution of dynamics on graphs is intimately related to space-time
dynamics governed by Partial Differential Equations (PDEs). Considering the
discretization of space in a continuously varying spatial field, the states
associated with the discrete nodes evolve according to some intrinsic physics
subject to spatial interaction. In the case of a PDE, the spatial interaction
follows the application (and truncation) of the Taylor series approximation
for the differential operators on the discretized domain, resulting in sparse
(and often symmetric) adjacency matrices. The corresponding data-driven method
for approximating PDEs in a consistent manner is Neural PDEs – a method-of-
lines technique where the spatial interaction (stencil of the discretization)
is learned through a convolutional neural network Rackauckas  _et al._ (2020);
Eliasof, Haber, and Treister (2021); Chamberlain  _et al._ (2021).

Networked dynamical systems do not necessarily exist over an explicit spatial
domain and instead have a sparse node-to-node coupling through fundamentally
different physics (e.g., resistive coupling Wang, Tang, and Ma (2019), or
pulse coupling Abbott and van Vreeswijk (1993)). In such systems, learning the
underlying graph is difficult. In Song, Choi, and Kahng (2021), a neural
network is trained on a large number of samples ($10^{6}$) to output an
adjacency matrix directly from observations of standard Kuramoto oscillators.
Other network reconstruction tasks have been attempted with sparse regression
techniques Han  _et al._ (2015), or by leveraging community structure Yan and
Gregory (2012). Compressive sensing techniques have been successfully used to
identify and place missing nodes Su, Wang, and Lai (2012); Wu, Wang, and Zheng
(2012). More recently, methods using the Gumbel Softmax technique Jang, Gu,
and Poole (2016) have been used to propose adjacency matrices from data Zhang
_et al._ (2019), where jointly trained network proposal models and dynamics
models work together to produce a candidate model for observed dynamics on a
graph.

###  II.2 Contributions

In this work, we are motivated by modeling tasks in which network connectivity
strongly influences qualitative behavior (i.e., limit cycling and bifurcation
structure). To that end, we emphasize learning physically-consistent models
that are interpretable and extensible. Our main contribution to this field is
the construction and demonstration of a modeling framework in which one can
recover system attributes required for informed modeling and control, namely
(i) the intrinsic nodal physics, (ii) the coupling physics, and (iii) the
underlying adjacency matrix. Our method utilizes universal differential
equations to construct a unified approach for cooperatively learning the
adjacency and physics components. We also show this method’s generalization
value by extending learned models to never-before-seen networks.

\begin{overpic}[width=162.6075pt]{img/fig_02.png} \end{overpic} Figure 2: Adjacency matrix for $11$-node directed network of nonlinear oscillators as defined in Eq. 10. Dark elements correspond to $\mathbf{A}_{i,j}=1$ with $\mathbf{A}_{i,j}=0$ otherwise. \begin{overpic}[width=216.81pt]{img/fig_03.png} \end{overpic} Figure 3: Single node data ($i=1$) for $t\in[0,50)$ for position $x$ and velocity $v$. The vertical dashed lines indicate intervals for data splits into thirds for train-development-test sets. \begin{overpic}[width=483.48647pt]{img/fig_04.png} \end{overpic} Figure 4: Sweep of sparsity promotion penalty weight, $\alpha$, over several orders of magnitude. The top row is the visualization of the learned adjacency matrix ($\hat{\mathbf{A}}$) for the different values of $\alpha$. The middle row corresponds to the element-wise squared error between the actual and learned matrices. The last row shows the thresholding step, i.e., rounding to the nearest integer. \begin{overpic}[width=216.81pt]{img/fig_05.png} \end{overpic} Figure 5: L2 error between learned and actual adjacency matrix over sweep of sparsity promotion penalty weight, $\alpha$. The lowest error occurs at $\alpha=1\cdot 10^{-6}$, although the error is insensitive to error until several orders of magnitude higher $\alpha$. \begin{overpic}[width=483.48647pt]{img/fig_06.png} \end{overpic} Figure 6: The learned nodal and coupling physics from the coupled system with the 11-node directed network of Fig. 6 is deployed on a never-before-seen 3-node network. The adjacency matrix for this network is shown in Fig. 1b. The learned physics models qualitatively capture the transients in the time signals, including the limit cycle behavior.  Table 1: Model selection metrics. Boldface represents the best (lowest) value among the candidate models. |  | $\alpha$
---|---|---
|  | 0 | $1\cdot 10^{-7}$ | $1\cdot 10^{-6}$ | $1\cdot 10^{-5}$ | $1\cdot 10^{-4}$ | $1\cdot 10^{-2}$ | $1\cdot 10^{-1}$
${||A||}_{1}$ | 24.59 | 24.58 | 24.46 | 23.82 | 22.16 | 26.74 | 26.74
$N_{f}$-step MSE | Train | 0.00789 | 0.01108 | 0.00761 | 0.00643 | 0.00837 | 0.01392 | 0.01392
Dev. | 0.01019 | 0.01442 | 0.00968 | 0.00911 | 0.01153 | 0.01790 | 0.01790
Test | 0.01002 | 0.01415 | 0.00951 | 0.00905 | 0.01148 | 0.01757 | 0.01757

##  III Methodology

Our goal is to identify a networked nonlinear dynamical system from
measurements of nodal states over time. We begin by defining $N$ assumed-
identical nodes with states $\mathbf{x}_{i}(t)\in\mathbb{R}^{d}$ where
$i=(1,2,...,N)$, $d$ is the state dimension of the $i$-th node, and $N\cdot d$
is the system state dimension. The connectivity of the nodes is prescribed by
an assumed time-invariant adjacency matrix $\mathbf{A}\in\mathbb{R}^{N\times
N}$, where entry $\mathbf{A}_{i,j}$ prescribes the learnable (parameterized by
$\theta_{g}$) interaction
$\mathbf{g}(\mathbf{x}_{i},\mathbf{x}_{j};\theta_{g}):\mathbb{R}^{2d}\rightarrow\mathbb{R}^{d}$
between nodes $i$ and $j$. This interaction is assumed to be directed
($\mathbf{A}$ not necessarily symmetric), unweighted (meaning our final
adjacency during deployment should have binary elements), and sparse.

Thus, we define an expected dynamical system (Figure 1) of the form:

$\frac{d\mathbf{x}_{i}}{dt}=\mathbf{f}\left(\mathbf{x}_{i};\theta_{f}\right)+\sum_{j\neq
i}^{N}\mathbf{A}_{i,j}\mathbf{g}\left(\mathbf{x}_{i},\mathbf{x}_{j};\theta_{g}\right)=\mathbf{F}\left(\mathbf{x};\Theta\right),$
(3)

where
$\mathbf{f}(\mathbf{x}_{i},\theta_{f}):\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$
is the learnable intrinsic physics of the nodes parameterized by $\theta_{f}$
and $\Theta$ is the concatenated vector of model parameters. In general,
$\mathbf{f}$ and $\mathbf{g}$ can be any differentiable surrogate model (e.g.
a neural network) mapping to the full nodal state dimension $d$. However,
these can be constructed to satisfy a priori knowledge of a system; e.g. a
parameterized equation, composite function, or other closure-type model
Rackauckas  _et al._ (2020).

A Neural ODE evolve states $\mathbf{x}$ through time $t$ according to the
transition dynamics specified by the RHS of Eq. 3. The evolution of
$\mathbf{x}$ corresponds to the initial value problem:

$\mathbf{x}\left(t_{end}\right)=\mathbf{x}\left(t_{0}\right)+\int_{t_{0}}^{t_{end}}\mathbf{F}\left(\mathbf{x};\Theta\right)dt.$
(4)

Equation 4 can be solved with standard numerical ODE solvers,

$\mathbf{x}\left(t_{end}\right)=\operatorname{ODESolve}\left(\mathbf{F}\left(\mathbf{x};\Theta\right),\mathbf{x}_{0},t_{0},t_{end}\right).$
(5)

Equation 5 can be made fully differentiable one of two ways: (i) via the
adjoint sensitivity method, whereby a second differential equation is
formulated and solved for the evolution of an error metric backward in time
from $t_{end}$ to $t_{0}$ (see Chen et al. Chen  _et al._ (2018)), or (ii)
through the computational graph of the forward pass of the ODE solver (e.g.
composition of differentiable elementary operations of the fourth-order Runge-
Kutta integrator, for example). The latter method is employed in this work for
its ease of implementation.

A necessary output of the training routine is a adjacency matrix with
connections between nodes chosen to fit the data. Assumed is that the network
to be modeled is unweighted, i.e., nodes $i$ and $j$ are either “connected”
with $A_{i,j}=1$ or “disconnected” with $A_{i,j}=0$. Several methods exist for
making such a selection task differentiable (e.g. the Gumbel-Softmax trick).
In this work, to begin training, we instead assume all-to-all coupling
activated with a sigmoid:

$\mathbf{A}=\sigma\left(\left[0\right]-I/\epsilon\right),$ (6)

with small number $\epsilon$, such that the off-diagonal entries of
$\mathbf{A}$ default to a value of one-half and are bounded to
$\mathbf{A}_{i,j}\in(0,1)$. After training is complete, the entries of
$\mathbf{A}$ are thresholded to either $0$ or $1$ based on some criteria; for
the present study, this is simply a rounding operation on the elements of
$\mathbf{A}$ post-training.

The training data is formatted in a batched tuple with $N_{f}$-step moving-
window forecasting horizon

$\mathcal{D}=\\{\mathbf{X}^{(k)}(t),\mathbf{X}^{(k)}(t+\Delta),...,\mathbf{X}^{(k)}(t+N_{f}\Delta)\\},$
(7)

where $\mathbf{X}=\\{\mathbf{x}_{1},...,\mathbf{x}_{N}\\}$; $k$ is the $k$-th
batch ($k=1,...,N_{b}$); $\Delta$ is the size of discrete timesteps. During
training, the surrogate model with parameters $\bm{\Theta}$ is asked to
predict the following $N_{f}-1$ timesteps provided $\mathbf{X}^{(k)}(t)$ for
each batch $k$. Therefore, the loss function to be minimized could be
constructed as the mean squared error between predictions and measurements,
given as

$\mathcal{L}(\bm{\Theta},\mathbf{A})=\frac{1}{N_{b}\cdot N_{f}\cdot N\cdot
d}||\hat{\mathbf{X}}^{(k)}(t;\bm{\Theta})-\mathbf{X}^{(k)}(t)||^{2}_{2}+\alpha||\mathbf{A}||_{1},$
(8)

where $\hat{\mathbf{X}}^{(k)}(t;\bm{\Theta})$ is prediction for all system
states at time $t$ in batch $k$; $||\mathbf{A}||_{1}$ is a $\ell_{1}$ penalty
intended to promote sparse connections in the adjacency matrix; and $\alpha$
is the weighting coefficient for this $\ell_{1}$ penalty.

##  IV Numerical Experiments

Presented are two numerical experiments designed to show the utility of the
presented methods. First, in Section IV.1, we perform the data-driven modeling
of a coupled oscillator system with unknown adjacency, node physics, and
coupling physics. Once these physics are learned, we shift our attention to
inferring system behavior on a never-before-seen network in Section IV.2.

###  IV.1 Modeling Coupled Nonlinear Oscillators

Consider a canonical nonlinear oscillator:

$\displaystyle\frac{d^{2}x_{i}}{dt^{2}}=$ $\displaystyle-
x_{i}-a_{1}\frac{dx_{i}}{dt}\left(a_{2}x_{i}^{4}-a_{3}x_{i}+a_{4}\right)$ (9)
$\displaystyle+\sum_{j=1}^{N}\mathbf{A}_{i,j}\left(\frac{d^{2}x_{i}}{dt^{2}}-\frac{d^{2}x_{j}}{dt^{2}}\right),$

or, written as a system of first-order ordinary differential equations:

$\displaystyle\frac{dx_{i}}{dt}=$ $\displaystyle v_{i}$ (10)
$\displaystyle\frac{dv_{i}}{dt}=$ $\displaystyle-
x_{i}-a_{1}v_{i}\left(a_{2}x_{i}^{4}-a_{3}x_{i}+a_{4}\right)$
$\displaystyle+\sum_{j=1}^{N}\mathbf{A}_{i,j}\left(v_{i}-v_{j}\right),$

with $N=11$, $a_{1,...,4}=\\{0.2,11.0,11.0,1.0\\}$, $\mu=0.2$, and adjacency
matrix specified by that of Fig. 2.

Depending on parameter values ($a_{1,...,4}$ and $\mu$), the system detailed
in Eq. 10 exhibits unique behaviors; for the specific parameter values and
network detailed, this system shows limit cycle behavior, as seen in the data
generated for this study and shown in Fig. 3. For demonstration of our
methodology, the system in Eq. 10 was simulated from a random initial
condition for 500 timesteps with temporal resolution of $\Delta t=0.1$. This
constitutes our dataset for this task.

We aim to fit the notional coupled dynamical system (Eq. 3) to this dataset.
The trainable surrogate for nodal intrinsic physics is
$\operatorname{NN}_{1}(x_{i},v_{i};\theta_{1}):\mathbb{R}^{2}\rightarrow\mathbb{R}^{1}$,
and similarly for the interaction physics:
$\operatorname{NN}_{2}(\mathbf{x}_{i},\mathbf{x}_{j};\theta_{2}):\mathbb{R}^{4}\rightarrow\mathbb{R}^{1}$.
Thus, by assuming second-order physics, our model takes the form:

$\displaystyle\frac{dx_{i}}{dt}=$ $\displaystyle v_{i}$ (11)
$\displaystyle\frac{dv_{i}}{dt}=$
$\displaystyle\operatorname{NN}_{1}(x_{i},v_{i};\theta_{1})+\sum_{j\neq
i}^{N}A_{i,j}\operatorname{NN}_{2}(\mathbf{x}_{i},\mathbf{x}_{j};\theta_{2})$
$\displaystyle A_{i,j}(\theta_{3})=$
$\displaystyle\sigma\left(\operatorname{vec}^{-1}(\theta_{3})-I/\epsilon\right),$

where the $\operatorname{vec}^{-1}\left(\cdot\right)$ operator reshapes the
vector $\theta_{3}\in\mathbb{R}^{N^{2}}$ to $\mathbb{R^{N\times N}}$. The
network $\operatorname{NN}_{1}$ is chosen to be a feed-forward network with
two hidden layers of 50 nodes each with Leaky ReLU activation functions.
$\operatorname{NN}_{2}$ has a single hidden layer of size 4 with Leaky ReLU
activation applied.

The integration scheme for solving Eq. 11 is a straightforward fourth-order
Runge-Kutta integrator. The forecast horizon is set to $N_{f}=5$ during
training. The optimizer used is Adam with a fixed learning rate of 0.02.
First, the model parameters are optimized for 1000 epochs to minimize the loss
as specified by Eq. 8 with $\alpha=0$; that is, without penalizing adjacency
sparsity. After 1000 epochs, the sparsity penalty weight $\alpha$ is increased
to a fixed value and the optimization continues with another 1000 epochs.

Table 1 summarizes the error metrics across the training, development, and
test data splits (Figure 3) through a sweep of values for $\alpha$ (log-10
scale). For the purposes of this work, the model selection criteria is the
trained model with the lowest Mean-Squared-Error (MSE) on the development data
split. For this set of experiments, a value of $\alpha=1\cdot 10^{-5}$
resulted in the lowest error for this metric. As shown in Fig. 5, the
corresponding learned adjacency matrix identically matches the ground-truth
adjacency matrix after a thresholding operation. In general, as the sparsity
promotion penalty increases, several elements in the adjacency matrix do not
meet the threshold and are set to zero. Note, however, that the error between
the learned and true adjacency matrices is relatively insensitive to $\alpha$
for low values of $\alpha$ (less than $\alpha=1\cdot 10^{-5}$) for our chosen
training procedure.

###  IV.2 Inference

The trained model of Section IV contains neural network approximations for the
intrinsic nodal physics ($\operatorname{NN}_{1}$) and nodal coupling physics
($\operatorname{NN}_{2}$). By virtue of the construction of our model as a
system of ODEs, these approximations should hold regardless of the underlying
graph structure of the system. To test this, these physics approximations are
deployed on a small, never-before-seen three-node network whose adjacency
matrix is given by Fig. 1b.

Figure 6 shows phase portraits and time series for the three nodes of the new
network. In the top row, position $x_{i}$ and velocity $v_{i}$ show the clear
limit cycle behavior for both a ground-truth simulation and the roll-out of
the learned model. In the second row, the time series reconstruction
(displayed is $x_{i}$) from a random initial condition for each node is shown.
Although this network was not included in the training data, the learned nodal
and coupling physics approximations perform well in reproducing the major
qualitative behavior of the system, including limit cycling and following
transient behavior of the nodal states as they approach the limit cycle (e.g.
for time $0<t<20$ for each of the nodes).

##  V Discussion

Although an exhaustive set of numerical experiments was not presented, we have
shown the general application and utility of universal differential equations
for structural inference of networked dynamical systems. Specifically, we have
shown that one can approximate (i) nodal physics, (ii) coupling physics, and
(iii) the underlying graphical structure of such systems directly from data.
Furthermore, we have demonstrated that these learned physics can be
successfully applied to never-before-seen networks and capture the general
qualitative behavior of these systems.

Several modeling choices were made to deliver the results contained within
this work; examples include the neural network architecture(s), the network
characteristics (weighted vs. unweighted, directed vs. undirected, etc.), and
the neural ODE model structure and integration scheme, among may others. While
each of these aspects of the modeling process can have large impacts on model
quality and performance, in this section detailed are those deserving of extra
attention. Specifically we focus on (i) potential shortcomings of the
framework and avoidance strategies and (ii) applications and future
development work.

### Shortcomings

The model form of Eq. 11 is generic in that assumed is second-order dynamics
with additive coupling - no other structure or physical prior is enforced in
the model construction. This contrasts so-called “gray” or “white”-box models
in which known terms are inserted into the model as a priori knowledge. In
this context, such prior knowledge can take the form of known nodal physics,
coupling physics, and/or adjacency matrix. The class of models presented in
this work fall into the category of “black”-box models, which are prone to
overfitting, especially in the context of neural ODEs. Model performance will
likely suffer for systems with (i) weakly nonlinear dynamics (e.g. a Kuramoto
oscillator network with small coupling coefficient) or (ii) data existing
solely on the system’s attractor. In the case of the former, disambiguation of
distinct signals from time series data may be difficult or impossible,
especially in the situation of global synchrony (a population of homogeneous
Kuramoto oscillators). Likewise, should observations of a system exist
exclusively on that system’s attractor, it is likely that a neural ODE-based
model may learn exclusively the shape of the attractor in phase space -
learned models are unlikely to generalize well off of the attractor. Inclusion
of prior knowledge will strongly regularize the neural ODE such that learned
closure terms or parameters will generalize better than their black-box
counterparts. For example, reformulating the learning task to learn
exclusively the adjacency matrix of a networked dynamical system given known
nodal and coupling physics will generally produce a better model than the
black-box counterpart.

Similarly, the networks presented in this work have $N=11$ for Section IV.1
and $N=3$ for IV.2 with uniform nodal physics. These are small networks with
sparse connections, though for the presented training routine, assumed is all-
to-all coupling meaning that the number of trainable parameters (for
$\mathbf{A}$) exhibits quadratic growth with the number of nodes, underscoring
the need for regularization during training (e.g. $\ell_{1}$ penalty for
$\mathbf{A}$).

### Applications and Future Work

This work was motivated by the open problem of inferring the graphical
structure of networked dynamical systems from data. Following the layout of
the presented results, this modeling framework can be directly applied to (i)
separate nodal and coupling physics from nodal interactions and (ii) estimate
the stability properties of an observed system with respect to perturbations
in the underlying network.

The class of physical systems amenable to the presented techniques is wide-
ranging. In closed combustion systems, thermoacoustic instabilities can arise
from coupling of injectors with fluid dynamic processes of the combustor-
injector placement and control is an area of intense study, which amounts to
perturbations of an adjacency matrix of the coupled injectorsDange  _et al._
(2019); Manoj  _et al._ (2019); Koch  _et al._ (2020). In power system network
dynamics, nodes obey the so-called swing equations for synchronous operation
of equipment. Such physics can be approximated by a network of coupled
oscillators, from which one can estimate stability properties of the system
directly from data. Extending the method further, time-varying adjacency
matrices can be incorporated to model molecular dynamics, whereby the
underlying graphical structure evolves with the nodes. In general, broad
applicability can be had by encoding physical priors. For example, one could
specify a known energy input/output balance in a model for a dissipative
system (e.g. building thermal dynamics and other lumped capacitance-type
models) such that the model can be readily constructed and analyzed in the
context of a particular graph structure. Our modeling framework can inform
system design and control in these contexts from a purely data-driven
perspective, e.g., the modification or evolution of the graph structure and
system control, among other tasks.

##  VI Conclusion

We have demonstrated the use of universal differential equations for the data-
driven modeling of networked dynamical systems. Specifically, we have
demonstrated the capability to extract (i) intrinsic nodal physics, (ii)
coupling physics, and (iii) the underlying graphical structure of the system.
This is accomplished by assuming a generic differential equation form that can
be applied to a wide variety of networked dynamical systems. The modeling
framework was demonstrated on an 11-node-directed network of nonlinear
oscillators. The utility of this framework lies in the disambiguation of the
physics from network structure. In doing so, the physics can be applied to a
network of arbitrary topology, allowing one to investigate system properties
with respect to perturbations in the underlying network.

## Acknowledgements

This work was supported through the Data Model Converge (DMC) initiative at
the Pacific Northwest National Laboratory (PNNL). PNNL is a multi-program
national laboratory operated for the U.S. Department of Energy (DOE) by
Battelle Memorial Institute under Contract No. DE-AC05-76RL0-1830.

## Conflicts of Interest

The authors have no conflicts to disclose.

## Author’s Contributions

J. Koch: conceptualization, software, investigation, visualization, draft
preparation. Z. Chen: software, investigation, visualization, draft
preparation. A. Tuor: software, supervision. J. Drgona: software, draft
editing, supervision. D. Vrabie: funding acquisition, advisor.

## Data Availability

Data available from the authors upon reasonable request.

## References

## References

  * Earn, Rohani, and Grenfell (1998) D. J. D. Earn, P. Rohani, and B. T. Grenfell, “Persistence, chaos and synchrony in ecology and epidemiology,” Proceedings of the Royal Society of London. Series B: Biological Sciences 265, 7–10 (1998).
  * Klemm and Bornholdt (2005) K. Klemm and S. Bornholdt, “Topology of biological networks and reliability of information processing,” Proceedings of the National Academy of Sciences 102, 18414–18419 (2005).
  * Ji and Deng (2021) W. Ji and S. Deng, “Autonomous discovery of unknown reaction pathways from data by chemical reaction neural network,” The Journal of Physical Chemistry A 125, 1082–1092 (2021).
  * Owoyele and Pal (2022) O. Owoyele and P. Pal, “ChemNODE: A neural ordinary differential equations framework for efficient chemical kinetic solvers,” Energy and AI 7, 100118 (2022).
  * Dange  _et al._ (2019) S. Dange, K. Manoj, S. Banerjee, S. A. Pawar, S. Mondal, and R. I. Sujith, “Oscillation quenching and phase-flip bifurcation in coupled thermoacoustic systems,” Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 093135 (2019).
  * Premraj  _et al._ (2021) D. Premraj, K. Manoj, S. A. Pawar, and R. I. Sujith, “Effect of amplitude and frequency of limit cycle oscillators on their coupled and forced dynamics,” Nonlinear Dynamics (2021), 10.1007/s11071-020-06135-y.
  * Guan  _et al._ (2022) Y. Guan, K. Moon, K. T. Kim, and L. K. Li, “Synchronization and chimeras in a network of four ring-coupled thermoacoustic oscillators,” Journal of Fluid Mechanics 938 (2022), 10.1017/jfm.2022.130.
  * Manoj  _et al._ (2019) K. Manoj, S. A. Pawar, S. Dange, S. Mondal, R. I. Sujith, E. Surovyatkina, and J. Kurths, “Synchronization route to weak chimera in four candle-flame oscillators,” Physical Review E 100, 062204 (2019).
  * Bakhtiarnia, Fahim, and Miandoab (2020) A. Bakhtiarnia, A. Fahim, and E. M. Miandoab, “Data-driven discovery of social network dynamics,” in  _2020 4th International Conference on Smart City, Internet of Things and Applications (SCIOT)_ (IEEE, 2020).
  * Keeling and Eames (2005) M. J. Keeling and K. T. Eames, “Networks and epidemic models,” Journal of The Royal Society Interface 2, 295–307 (2005).
  * Xiao  _et al._ (2008) Y. Xiao, M. Xiong, W. Wang, and H. Wang, “Emergence of symmetry in complex networks,” Physical Review E 77, 066108 (2008).
  * Strogatz (2001) S. H. Strogatz, “Exploring complex networks,” Nature 410, 268–276 (2001).
  * Arenas  _et al._ (2008) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, “Synchronization in complex networks,” Physics Reports 469, 93–153 (2008).
  * Li and Chen (2006) Z. Li and G. Chen, “Global synchronization and asymptotic stability of complex dynamical networks,” IEEE Transactions on Circuits and Systems II: Express Briefs 53, 28–33 (2006).
  * (15) Y. Kuramoto, “Self-entrainment of a population of coupled non-linear oscillators,” in  _International Symposium on Mathematical Problems in Theoretical Physics_ (Springer-Verlag) pp. 420–422.
  * Skomski  _et al._ (2021) E. Skomski, S. Vasisht, C. Wight, A. Tuor, J. Drgoňa, and D. Vrabie, “Constrained block nonlinear neural dynamical models,” in  _2021 American Control Conference (ACC)_ (2021) pp. 3993–4000.
  * Brunton, Proctor, and Kutz (2016) S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” Proceedings of the National Academy of Sciences 113, 3932–3937 (2016).
  * Chen  _et al._ (2018) R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, “Neural ordinary differential equations,” (2018), arXiv:1806.07366 [cs.LG] .
  * Rackauckas  _et al._ (2020) C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, A. Ramadhan, and A. Edelman, “Universal differential equations for scientific machine learning,” (2020), arXiv:2001.04385 [cs.LG] .
  * Scarselli  _et al._ (2009) F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini, “The graph neural network model,” IEEE Transactions on Neural Networks 20, 61–80 (2009).
  * Poli  _et al._ (2019a) M. Poli, S. Massaroli, J. Park, A. Yamashita, H. Asama, and J. Park, “Graph neural ordinary differential equations,” (2019a), arXiv:1911.07532 [cs.LG] .
  * Rusch  _et al._ (2022) T. K. Rusch, B. P. Chamberlain, J. Rowbottom, S. Mishra, and M. M. Bronstein, “Graph-coupled oscillator networks,” (2022), arXiv:2202.02296 [cs.LG] .
  * Asikis, Böttcher, and Antulov-Fantulin (2022) T. Asikis, L. Böttcher, and N. Antulov-Fantulin, “Neural ordinary differential equation control of dynamics on graphs,” Physical Review Research 4, 013221 (2022).
  * Duvenaud  _et al._ (2015) D. Duvenaud, D. Maclaurin, J. Aguilera-Iparraguirre, R. Gómez-Bombarelli, T. Hirzel, A. Aspuru-Guzik, and R. P. Adams, “Convolutional networks on graphs for learning molecular fingerprints,” (2015), arXiv:1509.09292 [cs.LG] .
  * Gilmer  _et al._ (2017) J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl, “Neural message passing for quantum chemistry,” (2017), arXiv:1704.01212 [cs.LG] .
  * Xhonneux, Qu, and Tang (2020) L.-P. A. C. Xhonneux, M. Qu, and J. Tang, “Continuous graph neural networks,” in  _Proceedings of the 37th International Conference on Machine Learning_ , ICML’20 (JMLR.org, 2020).
  * Poli  _et al._ (2019b) M. Poli, S. Massaroli, J. Park, A. Yamashita, H. Asama, and J. Park, “Graph neural ordinary differential equations,” CoRR abs/1911.07532 (2019b), 1911.07532 .
  * Kipf and Welling (2016) T. N. Kipf and M. Welling, “Variational graph auto-encoders,” (2016), arXiv:1611.07308 [stat.ML] .
  * Eliasof, Haber, and Treister (2021) M. Eliasof, E. Haber, and E. Treister, “Pde-gcn: Novel architectures for graph neural networks motivated by partial differential equations,” (2021), arXiv:2108.01938 [cs.LG] .
  * Chamberlain  _et al._ (2021) B. P. Chamberlain, J. Rowbottom, M. Gorinova, S. Webb, E. Rossi, and M. M. Bronstein, “Grand: Graph neural diffusion,” (2021), arXiv:2106.10934 [cs.LG] .
  * Wang, Tang, and Ma (2019) C. Wang, J. Tang, and J. Ma, “Minireview on signal exchange between nonlinear circuits and neurons via field coupling,” The European Physical Journal Special Topics 228, 1907–1924 (2019).
  * Abbott and van Vreeswijk (1993) L. F. Abbott and C. van Vreeswijk, “Asynchronous states in networks of pulse-coupled oscillators,” Physical Review E 48, 1483–1490 (1993).
  * Song, Choi, and Kahng (2021) J. U. Song, K. Choi, and B. Kahng, “Machine learning approaches for kuramoto coupled oscillator systems,” (2021), arXiv:2109.08918 [cond-mat.stat-mech] .
  * Han  _et al._ (2015) X. Han, Z. Shen, W.-X. Wang, and Z. Di, “Robust reconstruction of complex networks from sparse data,” Physical Review Letters 114, 028701 (2015).
  * Yan and Gregory (2012) B. Yan and S. Gregory, “Finding missing edges in networks based on their community structure,” Physical Review E 85, 056112 (2012).
  * Su, Wang, and Lai (2012) R.-Q. Su, W.-X. Wang, and Y.-C. Lai, “Detecting hidden nodes in complex networks from time series,” Physical Review E 85, 065201 (2012).
  * Wu, Wang, and Zheng (2012) X. Wu, W. Wang, and W. X. Zheng, “Inferring topologies of complex networks with hidden variables,” Physical Review E 86, 046106 (2012).
  * Jang, Gu, and Poole (2016) E. Jang, S. Gu, and B. Poole, “Categorical reparameterization with gumbel-softmax,” (2016), arXiv:1611.01144 [stat.ML] .
  * Zhang  _et al._ (2019) Z. Zhang, Y. Zhao, J. Liu, S. Wang, R. Tao, R. Xin, and J. Zhang, “A general deep learning framework for network reconstruction and dynamics learning,” Applied Network Science 4 (2019), 10.1007/s41109-019-0194-4.
  * Koch  _et al._ (2020) J. Koch, M. Kurosaka, C. Knowlen, and J. N. Kutz, “Mode-locked rotating detonation waves: Experiments and a model equation,” Physical Review E 101, 013106 (2020).


* [39] G.N. Candlish, R. Smith, and M. Fellhauer, 2015, RAyMOND: an N-body and hydrodynamics code for MOND. Mon. Not. R. Astron. Soc. 446, 1060
  * [40] G.F. Thomas, B. Famaey, R. Ibata, F. Lueghausen, and P. Kroupa, 2017, Stellar streams as gravitational experiments. I. The case of Sagittarius. Astron. Astrophys. 603, A65
  * [41] M. Bílek, I. Thies, P. Kroupa, and B. Famaey, 2018, MOND simulation suggests an origin for some peculiarities in the Local Group. Astron. Astrophys. 614, A59
  * [42] I. Banik, I. Thies, R. Truelove, G. Candlish, B. Famaey, M.S. Pawlowski, R. Ibata, and P. Kroupa, 2022, 3D hydrodynamic simulations for the formation of the Local Group satellite planes. Mon. Not. R. Astron. Soc. 513, 129
  * [43] I. Banik, S.T. Nagesh, H. Haghi, P. Kroupa, and HS. Zhao, 2022, Overestimated inclinations of Milgromian disc galaxies: The case of the ultradiffuse galaxy AGC 114905. Mon. Not. R. Astron. Soc. 513, 3541
  * [44] M. Milgrom, 1986, Solutions for the Modified Newtonian Dynamics Field Equation, Astrophys. J. 302, 617
  * [45] R. Brada and M. Milgrom, 2000, The Modified Dynamics is Conducive to Galactic Warp Formation. Astrophys. J. Lett. 531, L21
  * [46] R. Brada and M. Milgrom, 2000, Dwarf Satellite Galaxies in the Modified Dynamics. Astrophys. J. 541, 556
  * [47] M. Milgrom, 2015, MOND theory. Canad. J. Phys. 93, 107
  * [48] A. Yu. Ignatiev, 2007, Is Violation of Newton’s Second Law Possible? Phys. Rev. Lett. 98, 101101
  * [49] J. Bekenstein and J. Magueijo, 2006, Modified Newtonian dynamics habitats within the Solar System. Phys. Rev. D 73, 103513
  * [50] J. Magueijo and A. Mozaffari, 2012, Case for testing modified Newtonian dynamics using LISA pathfinder. Phys. Rev. D 85, 043527
  * [51] A.R. Penner, 2020, A proposed experiment to test gravitational anti-screening and MOND using Sun-Gas giant saddle points. Astrophys. Sp. Sc. 365, 154
  * [52] X. Wu and P. Kroupa, 2015, Galactic rotation curves, the baryon-to-dark-halo-mass relation and space-time scale invariance. Mon. Not. R. Astron. Soc. 446, 330
  * [53] A. Hees, B. Famaey, G.W. Angus, and G. Gentile, 2016, Combined Solar system and rotation curve constraints on MOND. Mon. Not. R. Astron. Soc. 455, 449
  * [54] H. Haghi, A.E. Bazkiaei, A.H. Zonoozi, and P. Kroupa, 2016, Declining rotation curves of galaxies as a test of gravitational theory. Mon. Not. R. Astron. Soc. 458, 4172
  * [55] K.H. Chae et al., 2020, Testing the Strong Equivalence Principle: Detection of the External Field Effect in Rotationally Supported Galaxies. Astrophys. J. 904, 51
  * [56] K.H. Chae et al., 2021, Testing the Strong Equivalence Principle. II. Relating the External Field Effect in Galaxy Rotation Curves to the Large-scale Structure of the Universe. Astrophys. J. 921, 104
  * [57] R. Stiskalek and H. Desmond, 2023, On the fundamentality of the radial acceleration relation for late-type galaxy dynamics, Mon. Not. R. Astron. Soc. 525, 6130
  * [58] H. Desmond, 2023, The underlying radial acceleration relation. Mon. Not. R. Astron. Soc. 526, 3342
  * [59] X. Hernandez, M.A. Jiménez, and C. Allen, 2012, Wide binaries as a critical test of classical gravity. Europ. Phys. J. 72, 1884
  * [60] I. Banik, Indranil, C. Pittordis, and W. Sutherland, 2021, Detailed numerical implementation of the wide binary test. arXiv:2109.03827
  * [61] K.H. Chae, 2023, Breakdown of the Newton-Einstein Standard Gravity at Low Acceleration in Internal Dynamics of Wide Binary Stars. Astrophys. J. 952, 128
  * [62] X. Hernandez, 2023, Internal kinematics of GAIA DR3 wide binaries: anomalous behaviour in the low acceleration regime. Mon. Not. R. Astron. Soc. 525, 1401
  * [63] K.H. Chae, 2023, Robust Evidence for the Breakdown of Standard Gravity at Low Acceleration from Statistically Pure Binaries Free of Hidden Companions. arXiv:2309.10404
  * [64] X. Hernandez, V. Verteletskyi, L. Nasser, and A. Aguayo-Ortiz, 2023, Statistical analysis of the gravitational anomaly in Gaia wide binaries. arXiv:2309.10995
  * [65] G.F. Thomas et al., 2018, Stellar streams as gravitational experiments. II. Asymmetric tails of globular cluster streams. Astron. Astrophys. 609, A44
  * [66] P. Kroupa et al., 2023, Asymmetrical tidal tails of open star clusters: stars crossing their cluster’s práh challenge Newtonian gravitation. Mon. Not. R. Astron. Soc. 517, 3613
  * [67] M. Milgrom, 2009, MOND effects in the inner Solar system. Mon. Not. R. Astron. Soc. 399, 474
  * [68] L. Blanchet and J. Novak, 2011, External field effect of modified Newtonian dynamics in the Solar system. Mon. Not. R. Astron. Soc. 412, 2530
  * [69] O. Bienaymé et al., 2014, Weighing the local dark matter with RAVE red clump stars. Astron. & Astrophys. 571, A92.
  * [70] S.R. Loebman et al., 2014, The Milky Way Tomography with Sloan Digital Sky Survey. V. Mapping the Dark Matter Halo. Astrophys. J. 794, 151.
  * [71] G.W. Angus et al., 2015, Mass models of disc galaxies from the DiskMass Survey in modified Newtonian dynamics. Mon. Not. R. Astron. Soc. 451, 3551
  * [72] M. Milgrom, 2015, Critical take on ”Mass models of disk galaxies from the DiskMass Survey in MOND”. arXiv:1511.08087.
  * [73] G.W. Angus, 2016, The dynamics of face-on galaxies in MOND. J. Phys. 718, 032001
  * [74] M. Milgrom, 2018, MOND and the dynamics of NGC 628, arXiv:1802.05140
  * [75] M. Lisanti et al., 2019, Testing dark matter and modifications to gravity using local Milky Way observables. Phys. Re. D 100, 083009.
  * [76] M. Milgrom, 1994, Dynamics with a Nonstandard Inertia-Acceleration Relation: An Alternative to Dark Matter in Galactic Systems. Ann. Phys. 229, 384
  * [77] M. Milgrom, 1999, The modified dynamics as a vacuum effect. Phys. Lett. A 253, 273
  * [78] M. Milgrom, 2002, MOND—theoretical aspects. New Astron. Rev. 46, 741
  * [79] M. Milgrom, 2006, MOND as modified inertia. arXiv:astro-ph/0510117
  * [80] M. Milgrom, 2011, MOND–particularly as modified inertia. Act. Phys. Polon. B 42, 2175; arXiv:1111.1611
  * [81] M. Milgrom, 2022, Models of modified-inertia formulation of MOND. Phys. Rev. D 106, 064060
  * [82] L.J. Wang, 2010, Recovering modified Newtonian dynamics by changing inertia. arXiv:1011.3618
  * [83] F. Namouni, 2015, Towards an interpretation of MOND as a modification of inertia. Mon. Not. R. Astron. Soc. 452, 210
  * [84] M. Alzain, 2017, Modified Newtonian Dynamics (MOND) as a Modification of Newtonian Inertia. J. Astrophys. Astron. 38, 59
  * [85] R. Costa, G. Franzmann, and J.P. Pereira, 2019, A local Lagrangian for MOND as modified inertia. arXiv:1904.07321
  * [86] M. Milgrom, 2014, General virial theorem for modified-gravity MOND. Phys. Rev. D 89, 024016
  * [87] M. Milgrom, 2006, Massive particles in acoustic space-times: Emergent inertia and passive gravity. Phys. Rev. D 73, 084005
  * [88] J.T. Cushing, 1981, Electromagnetic mass, relativity, and the Kaufmann experiments. Am. J. Phys.e 49, 1133
  * [89] O. Tiret, F. Combes, G.W. Angus, G. W. , B. Famaey, and H.S. Zhao, 2007, Velocity dispersion around ellipticals in MOND. Astron. Astrophys. 476, L1
  * [90] O. Müller, B. Famaey, and HS. Zhao, 2019, Predicted MOND velocity dispersions for a catalog of ultra-diffuse galaxies in group environments. Astron. Astrophys. 623, A36,
  * [91] H. Haghi, P. Kroupa, I. Banik, X.n Wu, A. Hasani Zonoozi, B. Javanmardi, A. Ghari, O. Müller, J. Dabringhausen, and HS. Zhao, 2019, A new formulation of the external field effect in MOND and numerical simulations of ultra-diffuse dwarf galaxies – application to NGC 1052-DF2 and NGC 1052-DF4. Mon. Not. R. Astron. Soc. 487, 2441
  * [92] E. Asencio et al., 2022, The distribution and morphologies of Fornax Cluster dwarf galaxies suggest they lack dark matter. Mon. Not. R. Astron. Soc. 515, 2981
  * [93] I, Banik and H.S. Zhao, 2015, The External Field Dominated Solution In QUMOND & AQUAL: Application To Tidal Streams. arXiv:1509.08457
  * [94] K.H. Chae and M. Milgrom 2022, Numerical Solutions of the External Field Effect on the Radial Acceleration in Disk Galaxies. Astrophys. J. 928, 24
  * [95] R. Pauc̆o and J. Klac̆ka, 2016, Sedna and the cloud of comets surrounding the Solar System in Milgromian dynamics. Astron. Astrophys. 589, A63
  * [96] R. Pauc̆o, 2017, Towards an explanation of orbits in the extreme trans-Neptunian region: The effect of Milgromian dynamics. Astron. Astrophys. 603, A11
  * [97] C. Migaszewski, 2023, On the origin of extreme trans-Neptunian objects within Modified Newtonian Dynamics. Mon. Not. Roy. Astron. Soc. 525, 805
  * [98] K. Brown and H. Mathur, 2023, Modified Newtonian Dynamics as an Alternative to the Planet Nine Hypothesis. Astron. J. 166, 168
  * [99] F.J. Sánchez-Salcedo, J. Reyes-Iturbide, and X. Hernandez, 2006, An extensive study of dynamical friction in dwarf galaxies: the role of stars, dark matter, halo profiles and MOND. Mon. Not. Roy. Astron. Soc. 370, 1829
  * [100] O. Tiret and F. Combes, 2007, Evolution of spiral galaxies in modified gravity. Astron. Astrophys. 464, 517
  * [101] C. Nipoti, L. Ciotti, J. Binney, and P. Londrillo, 2008, Dynamical friction in modified Newtonian dynamics. Mon. Not. Roy. Astron. Soc. 386, 2194
  * [102] F.J. Sánchez-Salcedo, 2009, Gaseous drag on a gravitational perturber in Modified Newtonian Dynamics and the structure of the wake. Mon. Not. Roy. Astron. Soc. 392, 1573
  * [103] G.W. Angus and A. Diaferio, 2009, Resolving the timing problem of the globular clusters orbiting the Fornax dwarf galaxy. Mon. Not. Roy. Astron. Soc. 396, 887.
  * [104] H. Vakili, P. Kroupa, and S. Rahvar, 2017, Type I Shell Galaxies as a Test of Gravity Models. Astrophys. J. 848, 55
  * [105] M. Bílek, HS. Zhao, B. Famaey, O. Müller, P. Kroupa, and R. Ibata, 2021, Evolution of globular-cluster systems of ultra-diffuse galaxies due to dynamical friction in MOND gravity. Astron. Astrophys 653, A170
  * [106] M. Bílek, HS. Zhao, B. Famaey, S.T. Nagesh, F. Combes, O. Müller, M. Hilker, P. Kroupa, and R. Ibata, 2023, Do old globular clusters in low mass galaxies disprove modified gravity? arXiv:2307.03202
  * [107] K. Vieira, et al. 2022, Milky Way Thin and Thick Disk Kinematics with Gaia EDR3 and RAVE DR5. Astrophys. J. 932, 28
  * [108] M. Sereno,and Ph. Jetzer, 2006, Dark matter versus modifications of the gravitational inverse-square law: results from planetary motion in the Solar system. Mon. Not. R. Astron. Soc. 371, 626
  * [109] M. Milgrom, 2020, The $a_{0}$ – cosmology connection in MOND. arXiv:2001.09729