# Information Bottleneck in Peptide Conformation Determination by X-ray
Absorption Spectroscopy

E. A. Eronen<EMAIL_ADDRESS>University of Turku, Department of
Physics and Astronomy, FI-20014 Turun yliopisto, Finland  A. Vladyka
University of Turku, Department of Physics and Astronomy, FI-20014 Turun
yliopisto, Finland  F. Gerbon  ESRF, The European Synchrotron, 71 Avenue des
Martyrs, CS40220, 38043 Grenoble Cedex 9, France  Ch. J. Sahle  ESRF, The
European Synchrotron, 71 Avenue des Martyrs, CS40220, 38043 Grenoble Cedex 9,
France  J. Niskanen<EMAIL_ADDRESS>University of Turku, Department
of Physics and Astronomy, FI-20014 Turun yliopisto, Finland

###### Abstract

We apply a recently developed technique utilizing machine learning for
statistical analysis of computational nitrogen K-edge spectra of aqueous
triglycine. This method, the emulator-based component analysis, identifies
spectrally relevant structural degrees of freedom from a data set filtering
irrelevant ones out. Thus tremendous reduction in the dimensionality of the
ill-posed nonlinear inverse problem of spectrum interpretation is achieved.
Structural and spectral variation across the sampled phase space is notable.
Using these data, we train a neural network to predict the intensities of
spectral regions of interest from the structure. These regions are defined by
the temperature-difference profile of the simulated spectra, and the analysis
yields a structural interpretation for their behavior. Even though the
utilized local many-body tensor representation implicitly encodes the
secondary structure of the peptide, our approach proves that this information
is irrecoverable from the spectra. A hard X-ray Raman scattering experiment
confirms the overall sensibility of the simulated spectra, but the predicted
temperature-dependent effects therein remain beyond the achieved statistical
confidence level.

##  I Introduction

Proteins are formed of amino acids, often cited as the building blocks of
life. The backbone of the amino acid chain consists of a -CO-NH-Cα\- sequence,
where adjacent residues are bound by the peptide bond Liljas2009 . The
geometry of the backbone contains pairs of flexible dihedral angles, the
Ramachandran angles, which together define the macroscopic structure and the
function of the protein. The related folding of proteins is a complex and
complicated question Dill2012 , affected by both intramolecular and
intermolecular interactions. The process is ultimately dependent on the
primary order of the amino acids Liljas2009  and the surrounding solvent
Ball2017 ; Canchi2013 .

While diffraction experiments allow for structure determination of proteins,
they require a crystalline sample. Owing to its localized mechanism, X-ray
spectroscopy maintains its sensitivity to atomistic structure also in the soft
condensed phase, although statistical variation of the individual spectra in
such environments has been found to be huge Wernet2004 ; Ottosson2011 ;
Niskanen2016a ; Niskanen2016b ; Niskanen2017 ; Sahle2018 ; Niskanen2019 ;
VazDaCruz2019 ; Pietzsch2022 . This variation can be accounted for by
evaluation of the spectra at numerous structures, sampled from the respective
statistical ensemble, and averaging over them. Furthermore, X-ray spectroscopy
is sensitive to the local structure of amino acids and peptides
Kaznacheyev2002 ; Gordon2003 ; Zubavichus2004 ; Zubavichus2005 ; Messer2005 ;
Ottosson2011 ; Blum2012 ; Weinhardt2019 . For example, different protonation
forms of glycine yield different resonant inelastic X-ray scattering spectra
Blum2012  and X-ray photoelectron spectra Ottosson2011 ; Niskanen2013 . It has
also been found that near-edge X-ray absorption fine structure is sensitive to
a few neighbouring amino acids in poly-Gly peptides Gordon2003 ;
Zubavichus2004 ; Weinhardt2019 .

To predict the native form of proteins in their biological environment,
classical molecular dynamics (MD) is often used Hollingsworth2018 . In this
technique, the atomistic composition of the system is accounted for, but the
forces between the atoms are computed from a prior parametrization rather than
from the electronic-nuclear system that ultimately defines them. Being
computationally much lighter, and allowing for time and size scales not
otherwise accessible, this makes classical MD subject to the quality of the
parametrization – the force field – used in the calculations
LindorffLarsen2012 ; MartinGarcia2015 . For the development of these models in
the solution environment, their performance assessment in the atomistic level
would be valuable.

Interpretation of X-ray spectra is typically not straightforward. This is due
to the complicated relation between spectra and structures, originating from
quantum nature of the electronic system. In addition, the statistical aspects
of the liquid state cause an experiment to probe only the ensemble average. To
tackle these problems, machine learning and related emulator-based component
analysis (ECA) Niskanen2022  may prove useful. The ECA algorithm carries out
dimensionality reduction in structural space to identify structural variations
having the strongest influence on the spectra. Thus the method identifies
recoverable (and irrecoverable) structural information without a prior
hypothesis.

In this work we report a theoretical study of the smallest tripeptide
triglycine in aqueous solution (see Figure 1) at temperatures of 300 K and 350
K, expecting a change in the distribution of the Ramachandran angles. We
calculate an extensive set of N K-edge X-ray absorption spectra from MD
trajectories with three different force fields, observing notable structural
and spectral variation. Next, we study the dependency between intensities of
spectral regions of interest (ROI) and the corresponding structures using a
neural network (NN) and ECA. The results show that while X-ray absorption ROIs
are sensitive to the nearest neighboring atoms of the absorption site, they
are not sensitive to the secondary structure i.e. the Ramachandran angles of
the system, in agreement with Schwartz et al. Schwartz2010 . The simulated
spectra are in good agreement with an X-ray Raman scattering (XRS) experiment
we performed. However, the predicted spectral difference profile as a function
of temperature remains unconfirmed with the achieved statistical uncertainty.

Figure 1: A schematic illustration of aqueous triglycine. The absorption sites
are denoted as N1, N2 and N3. The backbone dihedrals known as the Ramachandran
angles $\phi$ and $\psi$ are shown for the central residue. The molecular plot
was prepared using the VMD software HUMP96 .

##  II Methods

###  II.1 Simulations

We simulated the NPT ensemble by classical MD at 300 K and at 350 K and 1 bar
using three different force fields: AMBER-03 Duan2003 , Charmm27 MacKerell1998
; MacKerell2004  and OPLS-AA Kaminski2001 . For the first two force fields we
applied the TIP3P Jorgenssen1983  water model and for the latter the TIP4P
Jorgenssen1983  model. The cubic simulation cell had an edge length of $\sim
50$ Å. The input files were prepared using a desktop installation of GROMACS
Abraham2015  package (version 2020.1), and the simulations were performed
using the version 2020.5 gromacs2020.5  on a computing cluster. Rigid water
molecules were used whereas no constraints were applied on the triglycine
molecule. After 1 ns of initial thermalization with Berendsen thermostat and
barostat Berendsen1984 , we ran MD for 51 ns (timestep 0.5 fs) using the Nosé-
Hoover thermostatting Nose1984 ; Hoover1985  (time constant $\tau_{T}$=0.5 ps)
and Isotropic Parrinello-Rahman barostatting Andersen1980 ; Parrinello1980
(time constant $\tau_{P}$=5.0 ps, compressibility $4.5\times 10^{-5}$ bar-1).
To avoid the “hot-solvent/cold-solute” problem, we used separate thermostats
for the solute and for the solvent in the latter run, with conserved energy
drift of the order 1 kJ mol-1ns-1atom-1. The last 50 ns of the run were
sampled with 10 ps spacing for spectrum calculations.

From the MD snapshots we calculated the N K-edge X-ray absorption spectra of
30006 structures using the projector-augmented-wave (PAW) method Enkovaara2010
with plane wave basis and density functional theory (DFT), as implemented in
GPAW version 22.1.0 Mortensen2005 ; Enkovaara2010 ; Larsen2017 . As the
simulation cell of MD is too large to be used in the quantum mechanical
electron structure calculations, we used cut structures including only the
water molecules within 3.0 Å of the solute. A vacuum with a radius of 3.0 Å
was added around each structure. Excitations to the lowest 1500 valence
single-electron states were evaluated in the transition potential half-hole
(TP-HH) approximation Triguero1998  for each nitrogen site (Nex). The spectrum
onset was corrected using $\Delta$-DFT method for the lowest core-excited
state. The calculations utilized plane wave basis with the energy the cutoff
of 350 eV and the Perdew–Burke–Ernzerhof (PBE) functional pbe . To aid
convergence, occupation smearing by the Fermi-Dirac distribution (width 0.25
eV) was used.

The simulations resulted in energy–intensity pairs for the transitions of each
absorption site (denoted as N1, N2 and N3; see Figure 1), from which the
spectra were obtained by convolution with Gaussian functions of increasing
width analogously to the procedure presented in Leetmaa2010 . The full width
at half maximum of the function was obtained by a numerical grid search for
the best match with the experiment; we used 1.4 eV for the lowest state and
increased the value linearly in energy to 4.5 eV for states 5.5 eV above it,
or higher. An alternative set of convolution parameters was tested (from 0.2
eV to 4.25 eV for states 10 eV above the lowest state, or higher) without
qualitative changes in the results (see Supplementary Information). The
spectrum of a single snapshot was evaluated as the sum of the spectra from the
three nitrogen atoms and the resulting ensemble mean spectra of all the
systems matched the experiment well as seen in Figure 2. In this set, 11
spectra were omitted in the subsequent analysis due to obviously nonsensical
results.

To validate the computational results, sets of a few hundred to a few thousand
spectra were calculated while varying one of the following simulation
parameters: the number of water molecules included (both with a 6.0 Å cutoff,
and also without any water), the plane wave energy cutoff (600 eV) or the
functional (RPBE rpbe ). Some of these calculations were troubled with
convergence issues, but the successful ones show none of the parameter changes
to have a significant effect on the temperature difference profile (see
Supplementary Information). In total, the spectrum calculations took
approximately 300 000 CPU hours on Intel Xeon Gold 6148 processors.

###  II.2 Data analysis

We divided the spectra into three ROIs, the selection of which was based on
mean zero-passing location in the T-difference profiles $\Delta$ of the three
models. These areas roughly correspond to pre peak (I), main edge (II), and
continuum (III) in the spectrum, and each of them has a consistent
T-dependence in the simulations. Using the full spectrum instead of ROIs would
probably allow more information to be captured, but the risk of
overinterpretation also grows with the grid tightness, as simulations are
always erroneous in reproducing the experiment Niskanen2022a . To conclude, we
analyze the dependency between the ROI intensities $\mathbf{S}$ and the
corresponding structure $\mathbf{R}$.

We apply emulator-based component analysis Niskanen2022  to this data for its
interpretation. The algorithm searches for a few orthogonal basis vectors of
structural space by maximizing the generalized covered variance (R2-score) of
the prediction for projections of the data onto the spanned subspace. This
procedure utilizes the known target values of the original data points in the
evaluation of the score. The components of the ECA basis vectors indicate
dominant structural features in terms of the variation of the resulting output
(in this work the spectral ROIs). The method needs a suitable machine learning
(ML) based emulator $\mathbf{S}_{\mathrm{emu}}$ to predict the intensities of
the three ROIs of any given structure $\mathbf{R}$. We used an emulator
consisting of nine separate NNs, one for each nitrogen atom and ROI,
implemented by scikit-learn scikit-learn .

The individual structures $\mathbf{R}$ were encoded using the local many-body
tensor representation (LMBTR) Huo2022  implemented in the DScribe package
Himanen2020 . The system is represented as distances between pairs and angles
between triples of each element combination that includes the central atom –
the absorption site. The internal hyperparameters of the LMBTR descriptor and
the emulator were chosen by an alternating search, for best average ROI
prediction (see Supplementary Information). In the search for the best
descriptor, some LMBTR features were zero depending on the according set of
hyperparameters. These features were manually removed before training the
final model, without a significant effect on performance. The total
dimensionality of the final LMBTR vector $\mathbf{D}(\mathbf{R})$ describing
the local neighborhoods of the three nitrogen atoms was $1140$. The ECA
algorithm aims at dimensionality reduction, an optimization task complicated
by the notably large number of dimensions of this space. We applied an
alternating partial optimization with respect to 10% of the ECA vector
components at a time. The routine was completed with a full optimization of
all components at once. For these tasks we used the optimization toolbox of
the SciPy Scipy  Python library and the trust-region interior point method
Byrd1999  therein.

For the analysis of general spectrum–structure relationships, we combined all
the data from the three force fields and the two temperatures into one data
set, which was randomly divided for model selection and training (80%) and for
testing and application (20%) of the emulator.

###  II.3 Experiment

We measured XRS spectra schulke2007book  at the N K-edge of aqueous triglycine
using the multi-element XRS endstation Huotari2017  at beamline ID20 at the
European Synchrotron Radiation Facility (ESRF). In the experiment scattering
signal proportional to the double differential scattering cross section
(DDSCS) Sahle2014a
$\mathrm{d}^{2}\sigma/\mathrm{d}\Omega\mathrm{d}\omega_{2}$ at finite angle
element $\mathrm{d}\Omega$ and outgoing photon energy element
$\mathrm{d}\omega_{2}$ is recorded with a monochromatic and collimated
incident beam of hard X-rays. When carried out at small scattering angles
(forward direction) the XRS core-level signal is reduced to the dipole
spectrum, and thus yields a spectrum equivalent to XAS theory1967 . Four
analyzer modules, 12 spherically bent Si(660) crystals (bending radius
$R_{b}=1$ m) in each, were used to detect scattering in near-forward
direction. We used the spacial imaging property of the bent analyzer crystals
to include signal only from the interaction region by selecting the
corresponding pixels of the CCD detector in the analysis. We observed
different signal from the bulk liquid region and from the capillary walls and
therefore excluded the latter pixels from the analysis.

To minimize the effect from radiation damage, a heatable liquid flow cell
(modified version of a cell described elsewhere sahle2015miniature ) was used.
In the system, a magnetically driven pump rotor is used to circulate $\sim$ 5
ml of sample liquid through a capillary (2 mm outer diameter, 0.01 mm wall
thickness), where inelastic X-ray scattering takes place. Moreover, the sample
was replaced regularly to further reduce possible degradation. The sample
solution of molality 0.3 mol/kg of aqueous triglycine was prepared by
dissolving the powder sample (Sigma–Aldrich, purity $\geq$ 99.0%, lot #
BCCB1590) into de-ionized water ($\rho\approx
18.2\,\mathrm{M}\Omega\,\mathrm{cm}$, $\mathrm{TOC}\approx 2$ ppb).

##  III Results

Figure 2: The calculated and experimental results for aqueous triglycine. a–f:
Scatter plot of the Ramachandran angles of the central residues from the
simulated trajectories and the definitions Lovell2003  of the allowed (99.95%)
and favored (99.8%) regions of glycine residue according to the Top8000 data
set Hintze2016 . g–i: Computational ensemble mean and experimental spectra
(background removed) with the respective temperature difference profiles
$\Delta$ and spectral regions of interest. For the experiment, the difference
profile has been 4-fold binned from the spectra. The error bars and shading of
the simulated curve indicate the statistical uncertainty $\sigma$ (confidence
level 68%). The computational spectra have been shifted by $-$2.9 eV in all
cases for the pre peak to match with the experiment. The experiment is
presented scaled for the same main peak height as the respective 300 K
simulation.

The Ramachandran angles of the central residue, $\phi$ and $\psi$ (see Figure
1), can be visualized in a corresponding scatter plot. Depending on these two
dihedrals, the residue can either exhibit an $\alpha$-helix or $\beta$-sheet
secondary structure. Our MD simulations show triglycine appears in both of
these structural classes (Figure 2a-f), but the distributions are different
between the different force fields. As the molecule has no restrictive side
chains, the plot is symmetric for the left-handed and right-handed isomers.
The conformational variation of the liquid system gives rise to significant
spectral variation. Using simulations allows studying the dependency between
the shape of a spectrum and the underlying structure.

Figures 2g-i present the experimental spectra as well as the simulated
ensemble-mean spectra for the corresponding force fields and the two
temperatures. The temperature difference profile $\Delta$ is also shown. The
simulated curves were shifted by $-2.9$ eV to match the pre peak of the
experiment, and the experiment was scaled for the same main peak height with
the respective 300 K simulation. Due to the huge Compton scattering
background, approximately subtracted from the spectra, our experiment suffers
from a rather large statistical uncertainty. This prohibits further
conclusions as the predicted difference profiles lie within. However, an
impressive match between the simulated and experimental spectra is obtained.

We did not observe a significant linear correlation between any of the
structural features (internal coordinates and water related parameters) and
the ROI intensities (see Supplementary Information). However, a well-
functioning emulator will enable the use of ECA decomposition as a more
advanced analysis tool. In this case, the apparent collaborative action of the
structural degrees of freedom is, indeed, captured by an NN-based emulator as
depicted in Figure 3. The emulator achieved covered variances (R2 score) of
0.714, 0.945, and 0.953 for each region I–III with the test data. The
prediction accuracy is systematically worse for region I, the origin of which
is the N1s$\rightarrow\pi^{*}$ resonance of the N atoms of the peptide bonds.
Combining the regions together, the total covered variance stood at 0.874 for
z-score standardized ROIs and 0.942 for the absolute intensity ROIs. The
difference can be attributed to much larger total areas of ROIs II and III
dominating the absolute score. Going forward, a score of 0.874 is the upper
limit for covered spectral variance for the standardized ROI intensities by
ECA decomposition, as even with complete structural coverage the respective
ML-emulator-induced error will remain.

Figure 3: Emulator evaluation with standardized data. a-c: Predicted
intensities for ROIs I–III, respectively, plotted against the known ones.
Covered variance R2, Pearson’s r, and the mean squared error MSE are given for
each panel. Overall covered variance of the three ROIs is 0.874.

We split the test data further into two parts of the same size for the ECA.
The first part was used for the optimization of the standardized structural
descriptor space basis vectors (fit) and the latter was used for an
independent test (validation). The covered variance as a function of the rank
of the ECA expansion are given in Table 1 for both data parts. While the
covered ROI intensity variance improves up to rank five (emulator limit is
achieved), convergence for the validation data is practically reached with
three components.

Table 1: Covered ROI intensity variances as a function of the rank of the ECA decomposition. Rank | ECA | ECA
---|---|---
| fit | validation
1 | 0.536 | 0.519
2 | 0.753 | 0.712
3 | 0.858 | 0.813
4 | 0.865 | 0.817
5 | 0.877 | 0.817
Figure 4: The interatomic distance distributions for each of the nitrogen
atoms deduced from the first component vector transformed into the descriptor
space. a–d: Distances from the absorption site (denoted as Nex) to the
neighboring hydrogen, carbon, nitrogen and oxygen atoms, respectively. The
ROIs of the total N K-edge spectrum is sensitive to the hydrogen and oxygen
atoms at around 3 Å, which can mainly be contributed to the water molecules.
In addition, the ROIs are sensitive to the nearest carbon and nitrogen atoms
of each absorption site. Finally, the ROIs are also sensitive to the three
nearest hydrogen atoms of N1.

The ECA components can be used for structural interpretation Niskanen2022 ;
Vladyka2023 . We analyze the first component (Figure 4) in terms of convoluted
distributions of interatomic distances, a subset of the entire LMBTR
descriptor. This component shows the predominant structural changes associated
with spectral ROI variation: increase in ROIs I and III, and decrease in ROI
II (see Supplementary Information). There is a striking shift in N1–H, Nex–C
and Nex–O interatomic distances, and Nex–N distances, which however remain
somewhat inconclusive due to the limited distance range of the descriptor. The
surrounding water affects the spectrum as for both Nex–H and Nex–O there is a
shift at around 3 Å near the first H2O solvation shell. Interpretation of the
angular structural features from the LMBTR is more complicated because
exponential distance weighting was found to be required for best emulation
performance. The weighting mode was one of the numerous hyperparameters varied
in the model selection phase.

We next turn our focus on the features representing the secondary structure,
the Ramachandran angles. We assigned each simulated structure in the test set
(includes ECA fit and ECA validation data) to one of three classes:
$\alpha$-helix, $\beta$-sheet or ‘other’ based on the contours Lovell2003
from the Top8000 data set Hintze2016 . Borderline cases were assigned manually
as shown in Figure 5a. The fractions of $\alpha$-helices and $\beta$-sheets
are given in Table 2. We then used principal component analysis (PCA) to
reduce the dimensionality of the LMBTR-encoded structures to two (spectral
R${}^{2}=0.002$) as shown in Figure 5b. This plot shows clear clusters of
$\alpha$-helices and $\beta$-sheets with ‘others’ found in between the two.
The LMBTR can thus encode this information, even though the dihedrals $\phi$
and $\psi$ are not directly included in it. Correspondingly, we used the ECA
approach as a 2D dimensionality reduction tool for the same data (spectral
R${}^{2}=0.733$), as depicted in Figure 5c, without clearly observable
clusters. This difference can be explained by PCA focusing only on covered
structural variance, whereas ECA is guided by the spectral one. The result
therefore indicates that, for reasonable structures, XAS (or XRS) is
insensitive to the Ramachandran angles, which in turn cannot be reconstructed
from the spectral ROIs.

Table 2: The relative occurence (in %) of $\alpha$-helices and $\beta$-sheets at 300 K and 350 K with the three different force fields. | 300 K | 350 K
---|---|---
Force field | $\quad\alpha\quad$ | $\quad\beta\quad$ | $\quad\alpha\quad$ | $\quad\beta\quad$
AMBER03 | 31 | 61 | 27 | 62
Charmm27 | 66 | 30 | 58 | 34
OPLS-AA | 27 | 64 | 24 | 63
Figure 5: Analysis of the structural classes of the test data with PCA of
LMBTR features and ECA reconstruction. a: Data points classified to structural
classes. The allowed and favored regions according to the Top8000 data set
Hintze2016  are shown as contours. b: A 2-component PCA decomposition of the
LMBTR features of the data identifies the structural classes. c: A 2-component
ECA-coordinate reconstruction, based on spectral ROIs, covers drastically more
spectral variation but does not distinguish the structural classes. For
details, see text.

Feature importance is a metric used for the significance of an input feature
with respect to the output of an ML model Geron2019 . We define an importance
score for each LMBTR feature as its absolute value in the first ECA vector.
Next, we trained a new emulator with the same architecture as the original one
(except for the input layer) using a given number of LMBTR features of the
highest importance score. We tested the according models using the ECA
validation set and find that most LMBTR groups (e.g. a pair-wise distance
distribution) contain highly meaningful features in terms of performance
which, in turn, indicates strong interplay of the structural features of all
kinds. Improvement of covered spectral ROI variance is monotonous with respect
to the importance score and strongly saturates already with 300 features (see
Supplementary Information). This shows that magnitude of a feature in the ECA
vector directly measures its spectral significance.

##  IV Discussion

Although the different force fields give different structural distributions
(including the those of the Ramachandran angles) the spectral effect predicted
upon change of temperature is similar for all of them. This complicates the
validation of the force fields by core-level excitation spectra. Moreover, the
predicted difference profiles are small and cannot be confirmed or rejected by
the current experiment, as all predicted changes are within the error bars
that also include the possibility of zero effect. However, the experiment
supports the validity of the spectrum simulations, together with the
respective convergence checks.

Schwartz and co-workers used a liquid jet to measure the nitrogen K-edge X-ray
absorption spectrum of triglycine(aq) in ambient temperature Schwartz2010  by
total electron yield. They obtained a spectrum with roughly similar features,
but notably different pre peak to main edge ratio. Interestingly, our study is
in agreement with X-ray spectroscopy of solid triglycine Gordon2003 ;
Zubavichus2004 ; Weinhardt2019 . The bump at the post edge seen only in our
350 K spectrum seems to be present in the other studies Schwartz2010 ;
Gordon2003 ; Zubavichus2004 ; Weinhardt2019  and is also present in the N
K-edge spectrum of the similar molecule diglycine Gordon2003 . The cause of
this feature is not explained by our simulations or the measured 300 K
spectrum, but based on these references the possibility of some solid
triglycine in our 350 K heated sample cannot be excluded. The aforementioned
results have been obtained with numerous yield techniques not necessarily
equivalent to XRS used here, or the definition of XAS. Although the XRS has a
rather low count rate, it is known to be bulk sensitive and extremely stable.

The descriptor used in an ML work sets the ‘language’ for the resulting
scientific discussion. Although several descriptors for the atomistic
structure of molecules have been developed Himanen2020 ; Huo2022 ; Behler2011
; Bartok2013 ; Rupp2012 ; Hansen2015 ; Chandrasekaran2019 , we see three
general requirements for the spectrum analysis by ECA:

  1. 1.

The descriptor must allow for accurate ML with the available data.

  2. 2.

The descriptor must allow for ECA decomposition of high spectral variance
coverage with only a few dominant components.

  3. 3.

The descriptor must allow for interpretation in terms of structure; preferably
it is translatable into simple structural information.

The architecture of a well-performing descriptor gives insight into the
physical system, as it can encode relevant structural information. With LMBTR,
the distance-based weighting of angles improved the ML prediction accuracy,
which is understandable as the nearby atoms should, by intuition, affect the
spectrum more than the distant ones. Unfortunately, the weighting also makes
the angular features of the descriptor harder to interpret and (3) is not
reached. Thus, effective encoding of physical information for condition (1)
may render some of it unrecoverable, and moreover, a trade-off between all
three conditions can be expected. Previously, we have applied a similar method
for glassy GeO2 Vladyka2023 , where we found a variant of a Coulomb matrix
Rupp2012  (similar to the Bag of Bonds Hansen2015 ) a suitable descriptor for
ML and ECA. Arguably the well-defined covalent bond topology renders LMBTR
well-suited for aqueous triglycine.

The ECA algorithm applies projection of the structural descriptor vector onto
a limited number of basis vectors and relies on an emulator to predict the
spectral information for these new points much faster than the corresponding
electron structure calculation would do. This allows for using iterative
algorithms for the search of the the optimal basis vectors. The emulator must
be complex enough of a function (likely non-linear) yet generalizable to
sufficiently capture the relation between any relevant structure $\mathbf{R}$
and its spectrum $\mathbf{S}$. The selection of this mapping
$\mathbf{S}_{\mathrm{emu}}(\mathbf{D}(\mathbf{R}))$ is a complicated problem
with plenty of tunable hyperparameters, including those of the descriptor
$\mathbf{D}(\mathbf{R})$. In this work we settled for nine independent neural
networks, one for each atomic site and ROI with a combined output, as they
worked best for the LMBTR-encoded triglycine system with points (1)–(3) in
mind. As often the case in machine learning, the choice cannot be shown to be
the global optimum, but only the best in the particular hyperparameter search.
Up to date there are no generally accepted performance criteria in the X-ray
spectroscopic community.

The results of the importance score analysis show that almost all groups of
distance and angular distributions are necessary to maximise the covered
spectral ROI variance. On the other hand, a significant number of features
from each group could be ignored, as only 300 of them in total (originally
1140) are necessary to cover nearly all the accessible variance, which itself
rises quickly and monotonously with respect to the number of selected
features. This shows that the first ECA vector can indeed find the most
relevant features of the descriptor in the order of their spectral
significance and could serve as a tool for feature selection. Using the ECA,
the dependence of ROI intensities, as captured by the descriptor–emulator–ROI
mapping $\mathbf{S}_{\mathrm{emu}}(\mathbf{D}(\mathbf{R}))$, can practically
be condensed into three structural degrees of freedom when a 1140-dimensional
LMBTR is used.

The quantitative analysis of our simulations proves that while information
about the Ramachandran angles is present in the LMBTR-encoded structural
space, it is largely lost in the subspace covering a significant portion of
the spectral ROI variance. Therefore, our results support the conclusion of
Schwartz and co-workers Schwartz2010  in that XAS is insensitive to the
secondary structure of triglycine (within the reasonably expected structural
space). Furthermore, our results are at least partly in contradiction with
that of Gordon and co-workers Gordon2003  as we see a small dependency between
the most relevant part of the backbone conformation and its spectrum (Figure
5c) at most. This conclusion is manifested by also taking the mean spectra of
all $\alpha$-helices and all $\beta$-sheets in the data (see Supplementary
Information), which show a deviating difference profile from those presented
in Figure 2. The advantage of K-edge spectra, locality, seems to be a
limitation when it comes to secondary structure of proteins, for which there
is an information bottleneck for reasonably expected structures. The ECA shows
the ROIs of the total spectrum to mainly be sensitive to the nearest C and N
atoms from the solute together with the O and H atoms mostly from the water;
and additionally, to the nearest H atoms in the case of N-terminus of the
peptide.

##  V Conclusions

Experimental N K-edge X-ray Raman scattering spectra of aqueous triglycine can
be modelled by classical MD and TP-DFT calculations for a good match. However,
the experiment is not able to confirm or reject the predicted temperature
difference effect. A machine learning emulator can predict the intensities of
spectral regions of interest from the corresponding local many-body tensor
representation encoded structures. This enabled the application of the
emulator-based component analysis, which was able to condense the
structure–spectrum dependency practically into three degrees of freedom.
Moreover, the spectral significance of structural features was found to be
ordered along the magnitude of the according component in the first obtained
basis vector. From the analysis we conclude that the details of the secondary
structure of aqueous triglycine, expressed by the Ramachandran angles, are
lost in X-ray absorption spectral regions of interest due to an information
bottleneck. Instead, in the structural information available for analysis from
the used descriptor, the distance distributions between the absorption site
and its nearest neighbors significantly account for the spectral region
variance. Each of the tried force fields results in a similar temperature-
difference profile, and therefore distinguishing between them by X-ray
absorption seems improbable. On the other hand, the spectrum can be predicted
equally well with all of them.

In this work we have pushed structural decomposition by the ECA algorithm to
maximal size scales that X-ray spectra can be hoped to have sensitivity to.
The method identifies spectrally relevant structural subspace in a complicated
system without prior knowledge. Future prospects of this approach include
reconstructing the maximum obtainable structural information from the spectra,
the first steps of which have already been taken Vladyka2023 . We also note
that the method is not limited to X-ray spectroscopy, but can be applied to a
wide variety of inverse problems for which an emulator for the forward problem
is known.

## Author contributions

E.A.E. Data analysis, machine learning, experiment, simulations and writing
manuscript. A.V. Data analysis, experiment and writing manuscript. F.G.
experiment. C.J.S Experiment, data analysis and writing manuscript. J.N.
Research design, experiment, simulations, writing manuscript and funding.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

E.A.E. acknowledges Jenny and Antti Wihuri Foundation for funding. E.A.E.,
A.V. and J.N. acknowledge Academy of Finland for funding via project 331234.
The authors acknowledge CSC – IT Center for Science, Finland, and the FGCI –
Finnish Grid and Cloud Infrastructure for computational resources. We
acknowledge the European Synchrotron Radiation Facility (ESRF) for provision
of synchrotron radiation facilities. Prof. M. Metsä-Ketelä is thanked for
sharing his insight in sample environment materials used in biochemical
experiments. Dr. H. Müller at the ESRF is thanked for help regarding chemistry
lab and sample preparation during the experiment. Dr. Y. Watier at the ESRF is
thanked for support with the sample environment.

## Data availability

The data and relevant scripts are available in Zenodo: 10.5281/zenodo.8239300.

## References

  * [1] Anders Liljas, Lars Liljas, Jure Piskur, Göran Lindblom, Poul Nissen, and Morten Kjeldgaard.  Textbook of structural biology.  World Scientific Publishing, 2009.
  * [2] Ken A. Dill and Justin L. MacCallum.  The protein-folding problem, 50 years on.  Science, 338(6110):1042–1046, 2012.
  * [3] Philip Ball.  Water is an active matrix of life for cell and molecular biology.  Proceedings of the National Academy of Sciences, 114(51):13327–13335, 2017.
  * [4] Deepak R. Canchi and Angel E. García.  Cosolvent effects on protein stability.  Annual Review of Physical Chemistry, 64(1):273–293, 2013.
  * [5] Ph. Wernet, D. Nordlund, U. Bergmann, M. Cavalleri, M. Odelius, H. Ogasawara, L. Å. Näslund, T. K. Hirsch, L. Ojamäe, P. Glatzel, L. G. M. Pettersson, and A. Nilsson.  The structure of the first coordination shell in liquid water.  Science, 304(5673):995–999, 2004.
  * [6] Niklas Ottosson, Knut J. Børve, Daniel Spångberg, Henrik Bergersen, Leif J. Sæthre, Manfred Faubel, Wandared Pokapanich, Gunnar Öhrwall, Olle Björneholm, and Bernd Winter.  On the Origins of Core-Electron Chemical Shifts of Small Biomolecules in Aqueous Solution: Insights from Photoemission and ab Initio Calculations of Glycine (aq).  Journal of the American Chemical Society, 133(9):3120–3130, 2011\.
  * [7] J. Niskanen, Ch. J. Sahle, K. O. Ruotsalainen, H. Müller, M. Kavčič, M. Žitnik, K. Bučar, M. Petric, M. Hakala, and S. Huotari.  Sulphur k$\beta$ emission spectra reveal protonation states of aqueous sulfuric acid.  Scientific Reports, 6:21012, 2016.
  * [8] J. Niskanen, K. Kooser, J. Koskelo, T. Käämbre, K. Kunnus, A. Pietzsch, W. Quevedo, M. Hakala, A. Föhlisch, S. Huotari, and E. Kukk.  Density functional simulation of resonant inelastic x-ray scattering experiments in liquids: acetonitrile.  Physical Chemistry Chemical Physics, 18:26026–26032, 2016.
  * [9] Johannes Niskanen, Christoph J. Sahle, Keith Gilmore, Frank Uhlig, Jens Smiatek, and Alexander Föhlisch.  Disentangling structural information from core-level excitation spectra.  Physical Review E, 96:013319, 2017.
  * [10] Ch. J. Sahle, J. Niskanen, K. Gilmore, and S. Jahn.  Exchange-correlation functional dependence of the o 1s excitation spectrum of water.  Journal of Electron Spectroscopy and Related Phenomena, 222:57–62, 2018.
  * [11] J. Niskanen, M. Fondell, Ch. J. Sahle, S. Eckert, R. M. Jay, K. Gilmore, A. Pietzsch, M. Dantz, X. Lu, D. McNally, T. Schmitt, V. Vaz da Cruz, V. Kimberg, F. Gel’mukhanov, and A. Föhlisch.  Compatibility of quantitative x-ray spectroscopy with continuous distribution models of water at ambient conditions.  Proceedings of the National Academy of Science of USA, 116:4058–4063, 2019.
  * [12] V. Vaz da Cruz, F. Gel’mukhanov, S. Eckert, M. Iannuzzi, E. Ertan, A. Pietzsch, R. C. Couto, J. Niskanen, M. Fondell, M. Dantz, T. Schmitt, X Lu, D. McNally, R. M. Jay, V. Kimberg, A. Föhlisch, and M. Odelius.  Probing hydrogen bond strength in liquid water by resonant inelastic x-ray scattering.  Nature Communications, 10:1013, 2019.
  * [13] Annette Pietzsch, Johannes Niskanen, Vinicius Vaz da Cruz, Robby Büchner, Sebastian Eckert, Mattis Fondell, Raphael M. Jay, Xingye Lu, Daniel McNally, Thorsten Schmitt, and Alexander Föhlisch.  Cuts through the manifold of molecular H2O potential energy surfaces in liquid water at ambient conditions.  Proceedings of the National Academy of Sciences, 119(28):e2118101119, 2022.
  * [14] K. Kaznacheyev, A. Osanna, C. Jacobsen, O. Plashkevych, O. Vahtras, Ågren, V. Carravetta, and A. P. Hitchcock.  Innershell absorption spectroscopy of amino acids.  The Journal of Physical Chemistry A, 106(13):3153–3168, 2002.
  * [15] Michelle L. Gordon, Glyn Cooper, Cynthia Morin, Tohru Araki, Cássia C. Turci, Konstantin Kaznatcheev, and Adam P. Hitchcock.  Inner-shell excitation spectroscopy of the peptide bond: comparison of the C 1s, N 1s, and O 1s spectra of glycine, glycyl-glycine, and glycyl-glycyl-glycine.  The Journal of Physical Chemistry A, 107(32):6144–6159, 2003.
  * [16] Y. Zubavichus, M. Zharnikov, A. Schaporenko, and M. Grunze.  Nexafs study of glycine and glycine-based oligopeptides.  Journal of Electron Spectroscopy and Related Phenomena, 134(1):25–33, 2004.
  * [17] Yan Zubavichus, Andrey Shaporenko, Michael Grunze, and Michael Zharnikov.  Innershell absorption spectroscopy of amino acids at all relevant absorption edges.  The Journal of Physical Chemistry A, 109(32):6998–7000, 2005.
  * [18] B. M. Messer, C. D. Cappa, J. D. Smith, W. S. Drisdell, C. P. Schwartz, R. C. Cohen, and R. J. Saykally.  Local hydration environments of amino acids and dipeptides studied by x-ray spectroscopy of liquid microjets.  The Journal of Physical Chemistry B, 109(46):21640–21646, 2005\.
  * [19] M. Blum, M. Odelius, L. Weinhardt, S. Pookpanratana, M. Bär, Y. Zhang, O. Fuchs, W. Yang, E. Umbach, and C. Heske.  Ultrafast proton dynamics in aqueous amino acid solutions studied by resonant inelastic soft x-ray scattering.  The Journal of Physical Chemistry B, 116(46):13757–13764, 2012\.
  * [20] L. Weinhardt, A. Benkert, F. Meyer, M. Blum, D. Hauschild, R. G. Wilks, M. Bär, W. Yang, M. Zharnikov, F. Reinert, and C. Heske.  Local electronic structure of the peptide bond probed by resonant inelastic soft x-ray scattering.  Physical Chemistry Chemical Physics, 21(24):13207–13214, 2019.
  * [21] Johannes Niskanen, N. Arul Murugan, Zilvinas Rinkevicius, Olav Vahtras, Cui Li, Susanna Monti, Vincenzo Carravetta, and Hans Ågren.  Hybrid density functional–molecular mechanics calculations for core-electron binding energies of glycine in water solution.  Physical Chemistry Chemical Physics, 15:244–254, 2013.
  * [22] Scott A. Hollingsworth and Ron O. Dror.  Molecular dynamics simulation for all.  Neuron, 99(6):1129–1143, 2018.
  * [23] Kresten Lindorff-Larsen, Paul Maragakis, Stefano Piana, Michael P. Eastwood, Ron O. Dror, and David E. Shaw.  Systematic validation of protein force fields against experimental data.  PLOS ONE, 7(2):1–6, 2012.
  * [24] Fernando Martín-García, Elena Papaleo, Paulino Gomez-Puertas, Wouter Boomsma, and Kresten Lindorff-Larsen.  Comparing molecular dynamics force fields in the essential subspace.  PLOS ONE, 10(3):1–16, 2015.
  * [25] J. Niskanen, A. Vladyka, J. Niemi, and C.J. Sahle.  Emulator-based decomposition for structural sensitivity of core-level spectra.  Royal Society Open Science, 9(6):220093, 2022.
  * [26] Craig P. Schwartz, Janel S. Uejio, Andrew M. Duffin, Alice H. England, Daniel N. Kelly, David Prendergast, and Richard J. Saykally.  Investigation of protein conformation and interactions with salts via x-ray absorption spectroscopy.  Proceedings of the National Academy of Sciences, 107(32):14008–14013, 2010.
  * [27] William Humphrey, Andrew Dalke, and Klaus Schulten.  VMD – Visual Molecular Dynamics.  Journal of Molecular Graphics, 14:33–38, 1996.
  * [28] Yong Duan, Chun Wu, Shibasish Chowdhury, Mathew C. Lee, Guoming Xiong, Wei Zhang, Rong Yang, Piotr Cieplak, Ray Luo, Taisung Lee, James Caldwell, Junmei Wang, and Peter Kollman.  A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations.  Journal of Computational Chemistry, 24(16):1999–2012, 2003.
  * [29] A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin, and M. Karplus.  All-atom empirical potential for molecular modeling and dynamics studies of proteins.  The Journal of Physical Chemistry B, 102(18):3586–3616, 1998.
  * [30] Alexander D. Mackerell Jr., Michael Feig, and Charles L. Brooks III.  Extending the treatment of backbone energetics in protein force fields: Limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations.  Journal of Computational Chemistry, 25(11):1400–1415, 2004.
  * [31] George A. Kaminski, Richard A. Friesner, Julian Tirado-Rives, and William L. Jorgensen.  Evaluation and reparametrization of the opls-aa force field for proteins via comparison with accurate quantum chemical calculations on peptides.  The Journal of Physical Chemistry B, 105(28):6474–6487, 2001.
  * [32] William L. Jorgensen, Jayaraman Chandrasekhar, Jeffry D. Madura, Roger W. Impey, and Michael L. Klein.  Comparison of simple potential functions for simulating liquid water.  The Journal of Chemical Physics, 79(2):926–935, 1983.
  * [33] Mark James Abraham, Teemu Murtola, Roland Schulz, Szilárd Páll, Jeremy C. Smith, Berk Hess, and Erik Lindahl.  Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers.  SoftwareX, 1-2:19–25, 2015.
  * [34] Lindahl, Abraham, Hess, and van der Spoel.  Gromacs 2020.5 manual, 2021.
  * [35] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak.  Molecular dynamics with coupling to an external bath.  The Journal of Chemical Physics, 81(8):3684–3690, 1984.
  * [36] Shuichi Nosé.  A molecular dynamics method for simulations in the canonical ensemble.  Molecular Physics, 52(2):255–268, 1984.
  * [37] William G. Hoover.  Canonical dynamics: Equilibrium phase-space distributions.  Physical Review A, 31:1695–1697, 1985.
  * [38] Hans C. Andersen.  Molecular dynamics simulations at constant pressure and/or temperature.  The Journal of Chemical Physics, 72(4):2384–2393, 1980.
  * [39] M. Parrinello and A. Rahman.  Crystal structure and pair potentials: A molecular-dynamics study.  Physical Review Letters, 45:1196–1199, 1980.
  * [40] J Enkovaara, C Rostgaard, J J Mortensen, J Chen, M Dułak, L Ferrighi, J Gavnholt, C Glinsvad, V Haikola, H A Hansen, H H Kristoffersen, M Kuisma, A H Larsen, L Lehtovaara, M Ljungberg, O Lopez-Acevedo, P G Moses, J Ojanen, T Olsen, V Petzold, N A Romero, J Stausholm-Møller, M Strange, G A Tritsaris, M Vanin, M Walter, B Hammer, H Häkkinen, G K H Madsen, R M Nieminen, J K Nørskov, M Puska, T T Rantala, J Schiøtz, K S Thygesen, and K W Jacobsen.  Electronic structure calculations with gpaw: a real-space implementation of the projector augmented-wave method.  Journal of Physics: Condensed Matter, 22(25):253202, 2010.
  * [41] J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen.  Real-space grid implementation of the projector augmented wave method.  Physical Review B, 71(3):035109, 2005.
  * [42] Ask Hjorth Larsen, Jens Jørgen Mortensen, Jakob Blomqvist, Ivano E Castelli, Rune Christensen, Marcin Dułak, Jesper Friis, Michael N Groves, Bjørk Hammer, Cory Hargus, Eric D Hermes, Paul C Jennings, Peter Bjerre Jensen, James Kermode, John R Kitchin, Esben Leonhard Kolsbjerg, Joseph Kubal, Kristen Kaasbjerg, Steen Lysgaard, Jón Bergmann Maronsson, Tristan Maxson, Thomas Olsen, Lars Pastewka, Andrew Peterson, Carsten Rostgaard, Jakob Schiøtz, Ole Schütt, Mikkel Strange, Kristian S Thygesen, Tejs Vegge, Lasse Vilhelmsen, Michael Walter, Zhenhua Zeng, and Karsten W Jacobsen.  The atomic simulation environment—a python library for working with atoms.  Journal of Physics: Condensed Matter, 29(27):273002, 2017.
  * [43] L. Triguero, L. G. M. Pettersson, and H. Ågren.  Calculations of near-edge x-ray-absorption spectra of gas-phase and chemisorbed molecules by means of density-functional and transition-potential theory.  Physical Review B, 58:8097–8110, 1998.
  * [44] J. P. Perdew, K. Burke, and M. Ernzerhof.  Generalized gradient approximation made simple.  Physical Review Letters, 77:3865–3868, 1996.
  * [45] M. Leetmaa, M. P. Ljungberg, A. Lyubartsev, A. Nilsson, and L. G.M. Pettersson.  Theoretical approximations to X-ray absorption spectroscopy of liquid water and ice.  Journal of Electron Spectroscopy and Related Phenomena, 177(2-3):135–157, 2010.
  * [46] B. Hammer, L. B. Hansen, and J. K. Nørskov.  Improved adsorption energetics within density-functional theory using revised perdew-burke-ernzerhof functionals.  Physical Review B, 59(11):7413–7421, 1999.
  * [47] Johannes Niskanen, Anton Vladyka, J. Antti Kettunen, and Christoph J. Sahle.  Machine learning in interpretation of electronic core-level spectra.  Journal of Electron Spectroscopy and Related Phenomena, 260:147243, 2022.
  * [48] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [49] Haoyan Huo and Matthias Rupp.  Unified representation of molecules and crystals for machine learning.  Machine Learning: Science and Technology, 3(4):045017, 2022.
  * [50] Lauri Himanen, Marc O. J. Jäger, Eiaki V. Morooka, Filippo Federici Canova, Yashasvi S. Ranawat, David Z. Gao, Patrick Rinke, and Adam S. Foster.  DScribe: Library of descriptors for machine learning in materials science.  Computer Physics Communications, 247:106949, 2020.
  * [51] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [52] Richard H. Byrd, Mary E. Hribar, and Jorge Nocedal.  An interior point algorithm for large-scale nonlinear programming.  SIAM Journal on Optimization, 9(4):877–900, 1999.
  * [53] Winfried Schülke.  Electron dynamics by inelastic X-ray scattering, volume 7.  Oxford University Press, 2007.
  * [54] S. Huotari, Ch. J. Sahle, Ch. Henriquet, A. Al-Zein, K. Martel, L. Simonelli, R. Verbeni, H. Gonzalez, M.-C. Lagier, C. Ponchut, M. Moretti Sala, M. Krisch, and G. Monaco.  A large-solid-angle X-ray Raman scattering spectrometer at ID20 of the European Synchrotron Radiation Facility.  Journal of Synchrotron Radiation, 24(2):521–530, 2017.
  * [55] Ch. J. Sahle, A. Mirone, J. Niskanen, J. Inkinen, M. Krisch, and S. Huotari.  Planning, performing and analyzing X-ray Raman scattering experiments.  Journal of Synchrotron Radiation, 22(2):400–409, 2015.
  * [56] Yukio Mizuno and Yoshihiro Ohmura.  Theory of X-Ray Raman Scattering.  Journal of the Physical Society of Japan, 22(2):445–449, 1967.
  * [57] Ch. J. Sahle, Christian Henriquet, MA Schroer, Ina Juurinen, J Niskanen, and Michael Krisch.  A miniature closed-circle flow cell for high photon flux X-ray scattering experiments.  Journal of synchrotron radiation, 22(6):1555–1558, 2015.
  * [58] Simon C. Lovell, Ian W. Davis, W. Bryan Arendall III, Paul I. W. de Bakker, J. Michael Word, Michael G. Prisant, Jane S. Richardson, and David C. Richardson.  Structure validation by C$\alpha$ geometry: $\phi$,$\psi$ and C$\beta$ deviation.  Proteins: Structure, Function, and Bioinformatics, 50(3):437–450, 2003.
  * [59] Bradley J. Hintze, Steven M. Lewis, Jane S. Richardson, and David C. Richardson.  Molprobity’s ultimate rotamer-library distributions for model validation.  Proteins: Structure, Function, and Bioinformatics, 84(9):1177–1189, 2016.
  * [60] Anton Vladyka, Christoph J. Sahle, and Johannes Niskanen.  Towards structural reconstruction from x-ray spectra.  Physical Chemistry Chemical Physics, 25(9):6707–6713, 2023.
  * [61] Aurélien Géron.  Hands-on machine learning with Scikit-Learn, Keras and TensorFlow: concepts, tools, and techniques to build intelligent systems.  O’Reilly Media, Inc., 2 edition, 2019.
  * [62] Jörg Behler.  Atom-centered symmetry functions for constructing high-dimensional neural network potentials.  The Journal of Chemical Physics, 134(7):074106, 2011.
  * [63] Albert P. Bartók, Risi Kondor, and Gábor Csányi.  On representing chemical environments.  Physical Review B, 87:184115, 2013.
  * [64] Matthias Rupp, Alexandre Tkatchenko, Klaus-Robert Müller, and O. Anatole von Lilienfeld.  Fast and accurate modeling of molecular atomization energies with machine learning.  Physical Review Letters, 108:058301, 2012.
  * [65] Katja Hansen, Franziska Biegler, Raghunathan Ramakrishnan, Wiktor Pronobis, O. Anatole von Lilienfeld, Klaus-Robert Müller, and Alexandre Tkatchenko.  Machine learning predictions of molecular properties: Accurate many-body potentials and nonlocality in chemical space.  The Journal of Physical Chemistry Letters, 6(12):2326–2331, 2015\.
  * [66] Anand Chandrasekaran, Deepak Kamal, Rohit Batra, Chiho Kim, Lihua Chen, and Rampi Ramprasad.  Solving the electronic structure problem with machine learning.  npj Computational Materials, 5(1):22, 2019.

## Supplementary information

### Consistency of the spectrum simulations

Consistency of the spectrum simulations was studied with respect to geometric
water cutoff, the plane-wave energy cutoff, the exchange-correlation
potential, and the presence of the solvent. The results are presented in
Figure 6, which indicates similar shape for the obtained difference profile in
all cases. The statistical uncertainty $\sigma$ of the difference profile
$\Delta$ was calculated using a 10000-fold bootstrap procedure. The
calculations suffered from convergence issues: one third of evaluations with
6.0 Å solvent cutoff did not converge creating a potential selection bias.
Furthermore, calculations failed for a few snapshots in the other cases.

Figure 6: The difference profile remains qualitatively the same with a: 6.0 Å
water radial cutoff (one third of the snapshots did not converge), b: RPBE
functional (4 snapshots did not converge) and c: larger plane-wave (PW) energy
cutoff (2 snapshots did not converge), apart from an energy shift associated
with a more complete basis set. d: With the water completely removed the shape
of the simulated spectra differ from the experimental and other computational
results (1 snapshot did not converge). e: The spectra with the selected
simulation parameters.

### The best LMBTR hyperparameters within the DScribe package

A full list of the selected LMBTR hyperparameters within the implementation by
DScribe [50] is presented in Table 3.

Table 3: The selected LMBTR parameters within the DScribe package [50]. Parameter | Distances (k2) | Angles (k3)
---|---|---
geometry function | distance | angle
grid min | $0.7~{}\text{\AA}$ | $0^{\circ}$
grid max | $4.0~{}\text{\AA}$ | $180^{\circ}$
grid n | 30 | 10
grid $\sigma$ | $0.2~{}\text{\AA}$ | $12.5^{\circ}$
weighting function | unity | exponential
weighting scale | — | 0.8
weighting threshold | — | 1e–7

### Correlation between structural features and ROI intensities

We evaluated Pearson’s r coefficients between the ROIs and internal
coordinates i.e. bond lengths $r$, angles $\theta$, and dihedrals $\phi$. In
addition, the following water parameters were investigated: number of donated
$D$ and accepted $A$ hydrogen bonds, and the number of water molecules in each
of the two solvation shell $SS1$ and $SS2$ with respect to the absorbing
nitrogen. The definitions of the water parameters are described in [9]. Figure
7 shows the naming convention of the atoms in this analysis.

The results reveal that there are only a few features whose correlation is
consistent with respect to the temperature difference profile for all ROIs
(i.e. the sign of the correlation with a ROI intensity is always either equal
or opposite to the sign of the difference profile) within 2$\sigma$ error
obtained from a 10000-fold bootstrap resampling (Table 4). There are also many
features whose correlation is not strictly inconsistent with respect to the
difference profile (Table 5). Overall, the correlation analysis of the
internal coordinates and ROIs reveals notable complexity of the problem and
necessitates a more sophisticated analysis, that accounts for collaborative
action of several structural features.

Figure 7: The naming of the atoms. Molecular plot made by the Jmol software. Table 4: The features whose correlation is consistent with respect to the temperature difference profile with 2$\sigma$ uncertainty. Absolute Pearson’s r coefficients greater than or equal to 0.1 are given in bold. | I | II | III
---|---|---|---
$r_{\mathrm{C1,C2}}$ | $\hphantom{-}0.03$ | $-0.04$ | $\mathbf{\hphantom{-}0.11}$
$r_{\mathrm{N1,C1}}$ | $-0.01$ | $\mathbf{\hphantom{-}0.49}$ | $\mathbf{-0.36}$
$r_{\mathrm{C4,O2}}$ | $\mathbf{\hphantom{-}0.26}$ | $-0.06$ | $\hphantom{-}0.05$
$r_{\mathrm{C2,O1}}$ | $\mathbf{\hphantom{-}0.13}$ | $-0.06$ | $\hphantom{-}0.05$
$\theta_{\mathrm{C6,O2,O3}}$ | $\hphantom{-}0.04$ | $-0.02$ | $\hphantom{-}0.02$
$\phi_{\mathrm{N3,C5,H11,C6}}$ | $-0.02$ | $\hphantom{-}0.03$ | $-0.02$
$D_{\mathrm{N2}}$ | $\hphantom{-}0.04$ | $-0.07$ | $\mathbf{\hphantom{-}0.10}$
$D_{\mathrm{N3}}$ | $\hphantom{-}0.06$ | $-0.07$ | $\hphantom{-}0.09$
Table 5: The features whose correlation is not inconsistent with respect to the temperature difference profile with 2$\sigma$ uncertainty. Absolute Pearson’s r coefficients greater than or equal to 0.1 are given in bold. | I | II | III
---|---|---|---
$r_{\mathrm{C6,O4}}$ | $\hphantom{-}0.02$ | $\hphantom{-}0.00$ | $-0.01$
$r_{\mathrm{C4,N3}}$ | $\mathbf{\hphantom{-}0.27}$ | $-0.01$ | $\mathbf{\hphantom{-}0.40}$
$r_{\mathrm{C2,N2}}$ | $\mathbf{\hphantom{-}0.10}$ | $\hphantom{-}0.00$ | $\mathbf{\hphantom{-}0.38}$
$r_{\mathrm{C1,C2}}$ | $\hphantom{-}0.03$ | $-0.04$ | $\mathbf{\hphantom{-}0.11}$
$r_{\mathrm{N1,C1}}$ | $-0.01$ | $\mathbf{\hphantom{-}0.49}$ | $\mathbf{-0.36}$
$r_{\mathrm{O3,C6}}$ | $\hphantom{-}0.02$ | $\hphantom{-}0.00$ | $-0.01$
$r_{\mathrm{C4,O2}}$ | $\mathbf{\hphantom{-}0.26}$ | $-0.06$ | $\hphantom{-}0.05$
$r_{\mathrm{C3,H8}}$ | $-0.01$ | $-0.01$ | $\hphantom{-}0.02$
$r_{\mathrm{C2,O1}}$ | $\mathbf{\hphantom{-}0.13}$ | $-0.06$ | $\hphantom{-}0.05$
$r_{\mathrm{C1,H4}}$ | $-0.01$ | $-0.02$ | $\hphantom{-}0.03$
$r_{\mathrm{C1,H5}}$ | $-0.01$ | $-0.01$ | $\hphantom{-}0.02$
$r_{\mathrm{N1,H1}}$ | $-0.00$ | $\mathbf{\hphantom{-}0.17}$ | $-0.03$
$r_{\mathrm{N1,H2}}$ | $\hphantom{-}0.00$ | $\mathbf{\hphantom{-}0.16}$ | $-0.03$
$r_{\mathrm{N1,H3}}$ | $-0.01$ | $\mathbf{\hphantom{-}0.16}$ | $-0.04$
$\theta_{\mathrm{N2,C3,C4}}$ | $\mathbf{-0.19}$ | $\hphantom{-}0.03$ | $\hphantom{-}0.01$
$\theta_{\mathrm{C1,C2,N2}}$ | $-0.03$ | $\hphantom{-}0.01$ | $\hphantom{-}0.00$
$\theta_{\mathrm{C6,O2,O3}}$ | $\hphantom{-}0.04$ | $-0.02$ | $\hphantom{-}0.02$
$\theta_{\mathrm{C3,H8,C4}}$ | $\mathbf{\hphantom{-}0.14}$ | $-0.02$ | $-0.01$
$\theta_{\mathrm{N2,H6,C3}}$ | $-0.01$ | $-0.01$ | $\hphantom{-}0.07$
$\theta_{\mathrm{N1,H1,C1}}$ | $\hphantom{-}0.01$ | $-0.03$ | $\hphantom{-}0.03$
$\theta_{\mathrm{N1,H2,C1}}$ | $-0.00$ | $-0.03$ | $\hphantom{-}0.03$
$\theta_{\mathrm{N1,H3,C1}}$ | $-0.01$ | $-0.03$ | $\hphantom{-}0.02$
$\phi_{\mathrm{N3,C5,H10,C6}}$ | $\hphantom{-}0.01$ | $-0.02$ | $\hphantom{-}0.01$
$\phi_{\mathrm{N3,C5,H11,C6}}$ | $-0.02$ | $\hphantom{-}0.03$ | $-0.02$
$\phi_{\mathrm{N2,C3,H7,C4}}$ | $-0.02$ | $-0.01$ | $-0.00$
$\phi_{\mathrm{N2,C3,H8,C4}}$ | $\hphantom{-}0.01$ | $\hphantom{-}0.01$ | $\hphantom{-}0.00$
$\phi_{\mathrm{N1,C1,H4,C2}}$ | $-0.01$ | $-0.02$ | $\hphantom{-}0.00$
$\phi_{\mathrm{N1,C1,H5,C2}}$ | $\hphantom{-}0.02$ | $\hphantom{-}0.01$ | $\hphantom{-}0.00$
$SS1_{\mathrm{N1}}$ | $-0.01$ | $-0.01$ | $-0.03$
$SS1_{\mathrm{N2}}$ | $\hphantom{-}0.01$ | $\hphantom{-}0.01$ | $-0.02$
$SS1_{\mathrm{N3}}$ | $\hphantom{-}0.06$ | $-0.02$ | $-0.00$
$D_{\mathrm{N1}}$ | $\hphantom{-}0.00$ | $\mathbf{-0.10}$ | $\mathbf{\hphantom{-}0.17}$
$D_{\mathrm{N2}}$ | $\hphantom{-}0.04$ | $-0.07$ | $\mathbf{\hphantom{-}0.10}$
$A_{\mathrm{N2}}$ | $\hphantom{-}0.01$ | $-0.01$ | $\hphantom{-}0.01$
$D_{\mathrm{N3}}$ | $\hphantom{-}0.06$ | $-0.07$ | $\hphantom{-}0.09$
$A_{\mathrm{N3}}$ | $\hphantom{-}0.02$ | $-0.01$ | $\hphantom{-}0.01$

### Analysis of the spectral ROI intensity changes along the first ECA
component

Figure 8 shows an increasing trend for the ROI I and III intensities, and a
decreasing one for the ROI II intensity, along the first ECA component vector.

Figure 8: Intensity behaviour of the ROIs along the first ECA component a–c:
An increase is observed for ROI I and ROI III, and and a decrease is observed
for ROI II.

### Spectral significance of the components in the ECA vector

We carried out analysis of features by ordering them along the importance
score (component magnitude in the first ECA vector). We trained an emulator
multiple times iteratively removing 30 of the least important features at each
cycle. Only the input layer of the emulator architecture was modified between
the iterations.

The results show that most and almost all spectral ROI variance is covered by
150 (Figure 9) and by 300 (Figure 10) features of the highest importance
score, respectively. Almost every group of LMBTR features is present in both
cases (Figure 9a and Figure 10a). The covered variance rises quickly and
monotonously with respect to the number of selected features (Figure 9b and
Figure 10b), which shows that ECA can find and order the most relevant
structural features with respect to the spectral ROI variance. The PCA (Figure
9c and Figure 10c) of the selected features shows that the information about
the Ramachandran angles is still present in the reduced feature space (150 or
300 features, respectively).

Figure 9: a: The frequency of occurrence of a feature from different LMBTR
groups when 150 features of the highest importance score are used. b: The
covered spectral ROI intensity variation as a function of the set size. Rapid
and monotonous saturation is observed and the set of 150 features covers most
of the spectral variance. c: 2-component PCA of the 150-dimensional feature
vectors from feature selection by the importance score. The information about
the Ramachandran angles is still present in this structural descriptor. Figure
10: a: The frequency of occurrence of a feature from different LMBTR groups
when 300 features of the highest importance score are used. b: The covered
spectral ROI intensity variation as a function of the set size. Rapid and
monotonous saturation is observed and the set of 300 features covers most of
the spectral variance. c: 2-component PCA of the 300-dimensional feature
vectors from feature selection by the importance score. The information about
the Ramachandran angles is still present in this structural descriptor.

### Spectral difference between $\beta$-sheets and $\alpha$-helices

The mean spectra of all the structures representing $\beta$-sheets or
$\alpha$-helices is shown in Figure 11. The difference profile between the two
is also shown. The $\alpha$/$\beta$-difference profile differs drastically
from the temperature difference profile.

Figure 11: The mean spectra of all the structures classified as $\beta$-sheets
or $\alpha$-helices presented together with their difference profile.

### Evaluating the main results with an alternative convolution

We evaluated the results with an alternative convolution: Gaussians with FWHM
of 0.2 eV for the lowest state increasing linearly in energy to 4.25 eV for
states 10 eV above it, or higher. While the shape of the computational spectra
differ more drastically from the experimental ones, the difference profile is
qualitatively the same (Fig. 12). Compared to the main text, the same
conclusions can be made from the interatomic distance distributions from the
first ECA vector (Fig. 13) and the analysis of the $\alpha$/$\beta$ structural
classes (Fig. 14).

Figure 12: The calculated (and experimental) spectra for aqueous triglycine
with alternative convolution parameters. a–c: Computational ensemble mean and
experimental spectra (background removed) with the respective temperature
difference profiles $\Delta$ and spectral regions of interest. For the
experiment, the difference profile has been 4-fold binned from the spectra.
The error bars and shading of the simulated curve indicate the statistical
uncertainty $\sigma$ (confidence level 68%). The computational spectra have
been shifted by $-$2.9 eV in all cases for the pre peak to match with the
experiment. The experiment is presented scaled for the same main peak height
as the respective 300 K simulation. Figure 13: The interatomic distance
distributions for each of the nitrogen atoms deduced from the first component
vector obtained with an alternative convolution compared to the main text and
transformed into the descriptor space. a–d: Distances from the absorption site
(denoted as Nex) to the neighboring hydrogen, carbon, nitrogen and oxygen
atoms, respectively. The ROIs of the total N K-edge spectrum is sensitive to
the hydrogen and oxygen atoms at around 3 Å, which can mainly be contributed
to the water molecules. In addition, the ROIs are sensitive to the nearest
carbon and nitrogen atoms of each absorption site. Finally, the ROIs are also
sensitive to the three nearest hydrogen atoms of N1.  Figure 14: Analysis of
the structural classes of the test data with PCA of LMBTR features and ECA
reconstruction with an alternative convolution. a: Data points classified to
structural classes. The allowed and favored regions according to the Top8000
data set [59] are shown as contours. b: A 2-component PCA decomposition of the
LMBTR features of the data identifies the structural classes. c: A 2-component
ECA-coordinate reconstruction, based on spectral ROIs, covers drastically more
spectral variation but does not distinguish the structural classes.


# Learning deterministic hydrodynamic equations from stochastic active
particle dynamics

Suryanarayana Maddu1,2,3,4, Quentin Vagne2,3,6,7, Ivo F. Sbalzarini1,2,3,4,5 1
Technische Universität Dresden, Faculty of Computer Science, 01069 Dresden,
Germany  2 Max Planck Institute of Molecular Cell Biology and Genetics, 01307
Dresden, Germany  3 Center for Systems Biology Dresden, 01307 Dresden, Germany
4 Center for Scalable Data Analytics and Artificial Intelligence ScaDS.AI,
Dresden/Leipzig, Germany  5 Cluster of Excellence Physics of Life, TU Dresden,
Germany  6 Max Planck Institute for the Physics of Complex Systems, 01187
Dresden, Germany  7 Department of Genetics and Evolution, University of
Geneva, Switzerland

###### Abstract

We present a principled data-driven strategy for learning deterministic
hydrodynamic models directly from stochastic non-equilibrium active particle
trajectories. We apply our method to learning a hydrodynamic model for the
propagating density lanes observed in self-propelled particle systems and to
learning a continuum description of cell dynamics in epithelial tissues. We
also infer from stochastic particle trajectories the latent phoretic fields
driving chemotaxis. This demonstrates that statistical learning theory
combined with physical priors can enable discovery of multi-scale models of
non-equilibrium stochastic processes characteristic of collective movement in
living systems.

††preprint: APS/123-QED

Collective movement is ubiquitous in living systems and is observed at all
spatial scales from groups of animals [1] and shoal of fish [2] to the motion
of cells within tissues [3] and cytoskeletal molecules within cells [4, 5].
Despite their diverse nature, these systems exhibit common emerging
properties, including density-dependent transitions to ordered phases [6],
persistent trajectories [7], large density fluctuations [2], and
spatiotemporal patterning [8]. Naturally then, a fundamental question of
collective self-organized motion is how interactions between constituents at
the microscopic scale lead to the emergent dynamics at the macroscopic scale.

This question has been addressed by physical theories of self-organized, non-
equilibrium active systems both at the microscopic scale of individual
constituents [9] and at the macroscopic scale in the form of nonlinear
hydrodynamic equations [10]. Explaining the mechanism of self-organization,
however, requires linking the two levels of description in order to, e.g.,
predict how the microscopic interaction parameters define the hydrodynamic
transport coefficients [11] in a multi-scale model [12]. This has been
addressed by constructing continuum theories based on coarse-graining a
microscopic model [6, 13], by symmetry arguments [14], and by non-equilibrium
thermodynamics close to equilibrium [15]. These approaches, however, require
closure assumptions, are restricted to microscopic models that are
analytically tractable, and tend to produce over-complete continuum models
that are difficult to understand.

Figure 1: Principle of data-driven coarse-graining: A: Microscopic data from
stochastic active particle (agent) systems. B: Mean-field quantities of
interest computed from particle trajectory data. C: Statistical learning by
solving a penalized regression problem using robustness to data perturbations
for model selection. D: Validation of the learned hydrodynamic models using
criteria like data fidelity, interpretability, universality, or (statistical)
recovery probability.

Analytical coarse-graining is therefore increasingly complemented with data-
driven approaches using machine learning [16]. Statistical learning frameworks
have been used to infer effective dynamics from data both in space and time
[17, 18, 19], and it has recently been shown how to enforce consistency of the
learned models with first principles [20]. So far, however, most data-driven
approaches are applied on one scale, for example to learn particle interaction
potentials from particle distributions [21], to learn force fields and non-
equilibrium currents from stochastic trajectories [22] and Brownian movies
[23], or to learn hydrodynamic equations from hydrodynamic fields [24]. In
contrast to these single-scale applications, data-driven coarse-graining can
be used to learn hydrodynamic transport coefficients from microscopic data, to
verify closures in kinetic approaches, and to reveal novel principles of self-
organization and collective behavior. This has been impressively demonstrated
using specialized model designs [16, 25], but a generic framework guaranteeing
physical consistency and statistical robustness of the learned models has so
far been lacking.

Here, we provide a principled statistical learning framework for data-driven
coarse-graining under a wide class of physical priors. In doing so, we extend
the concept of group sparsity [20] to cross-scale problems. This enables us to
learn coarse-grained equations with spatiotemporally varying coefficients,
extract the latent stochastic drift fields from particle trajectories, and
infer hidden dependences intrinsic to the microscopic system.

Problem statement: Given observed trajectories $\bm{r}_{p}(t_{i})$ of active
particles (or agents) $p=1,\ldots,N_{p}$ along with properties like particle
velocity $\bm{v}_{p}(t_{i})$ or orientation $\theta_{p}(t_{i})$ at T discrete
times $t_{i}=\\{t_{0},t_{0}+\Delta t,\ldots,t_{0}+\text{T}\Delta t\\}$, infer
the functional form of a partial differential equation for hydrodynamic
variables of interest (Fig. 1A,B). Assuming that slowly varying hydrodynamic
quantities can be extracted from particle trajectories by averaging, the first
step of equation inference is to construct an over-complete dictionary of all
possible right-hand-side terms and numerically evaluate their values on the
data [17, 19]. The canonical form of a one-dimensional model with a single
scalar state variable $u$ and a dictionary
$\bm{\Theta}\in\mathbb{R}^{\text{N}\times\text{P}}$ of $\text{P}\in\mathbb{N}$
possible terms is:

$\underbrace{\begin{bmatrix}|\\\ u_{t}\\\ |\end{bmatrix}}_{\text{{\hbox{
\bm{U}_t \in\mathbb{R}^{ \text{N}\times
1}}}}}=\underbrace{\begin{bmatrix}\>\>|\qquad|\qquad|\qquad|\qquad|\qquad|\quad\\\
u\quad uu_{x}\>\>\>\>\ldots u_{xx}\>\ldots\>\>\>\ldots\\\
\>\>|\qquad|\qquad|\qquad|\qquad|\qquad|\quad\end{bmatrix}}_{\text{{\hbox{\bm{\Theta}\in\mathbb{R}^{
\text{N}\times\text{P}}}}}}\underbrace{\bm{\xi}}_{\in\mathbb{R}^{\text{P}\times
1}}.$ (1)

The left-hand-side vector $\bm{U}_{t}$ contains the discrete approximations of
the time derivatives of $u$ at N selected space and time points, and each
column of $\bm{\Theta}$ contains the discrete approximations of one possible
term of the right-hand side at the same data points. The problem then amounts
to finding a sparse (ideally the sparsest) vector $\bm{\xi}$ such that this
model explains the data (Fig. 1C). The sparsest $\bm{\xi}$ defines the
simplest sufficient model with the fewest right-hand-side terms. Physical
priors, e.g., about symmetries or conservation laws are incorporated by
grouping columns of $\bm{\Theta}$ [20]. Thus, the goal is to solve the
optimization problem:

$\displaystyle\hat{\bm{\xi}}^{\lambda}=\arg\min_{\bm{\xi}}$
$\displaystyle\frac{1}{2}\|\bm{U}_{t}-\sum_{j=1}^{m}\mathbf{\Theta}_{g_{j}}\bm{\xi}_{g_{j}}\|_{2}^{2}+$
$\displaystyle\lambda\sum_{j=1}^{m}\sqrt{p_{g_{j}}}\>\mathbb{1}\left(\|\bm{\xi}_{g_{j}}\|_{2}\neq
0\right),$ (2)

where $m$ is the number of groups,
$\bm{\Theta}_{g_{j}}\in\mathbb{R}^{\text{N}\times p_{g_{j}}}$ is the sub-
matrix of $\bm{\Theta}$ formed by all columns belonging to group
$g_{j}\subseteq\\{1,\ldots,\text{P}\\}$,
$\bm{\xi}_{g_{j}}=\\{\bm{\xi}_{i}:i\in g_{j}\\}$ is the vector $\bm{\xi}$
restricted to the index set $g_{j}$ of size $p_{g_{j}}$, i.e.,
$|g_{j}|=p_{g_{j}}$, and $\mathbb{1}(\cdot)$ is the indicator function. We
solve this optimization problem using the gIHT algorithm [20] based on
approximate proximal operators.

Stability selection chooses the regularization: Meaningful models are only
found if the regularization coefficient $\lambda$ is well chosen [17, 18, 19].
We choose it automatically so as to maximize the stability of the inference,
i.e., to favor models that are robust to random variations in the data [20].
The probability that group $g_{j}$ is selected (i.e., the corresponding
entries in $\bm{\xi}$ are non-zero) at a given value of $\lambda$ is
$\widehat{\Pi}_{j}^{\lambda}=B^{-1}\sum_{b=1}^{B}\mathbb{1}(g_{j}\cap\hat{S}^{\lambda}[I_{b}^{*}]\neq\emptyset)$
with $I^{*}$ an independent random sub-sample of the data of size
$|I^{*}|=\lfloor\text{N}/2\rfloor$ and
$\hat{S}^{\lambda}[I_{b}^{*}]=\\{g_{j}:\|\bm{\xi}_{g_{j}}^{\lambda}\|_{2}\neq
0\\}$ the inferred group coefficients for this sub-sample. This process is
independently repeated for $B\in\mathbb{N}>1$ $(\approx 100)$ different sub-
samples of the data. Using the statistical theory of stability selection [26]
the stable groups can be extracted as the set
$\mathbf{S}(\lambda)=\\{{g_{j}}:\widehat{\Pi}_{j}^{\lambda}\geq\pi_{\text{th}}\\}$
with the threshold probability $\pi_{\text{th}}$ chosen to control the number
of expected false positives, $E_{\text{fp}}$, corresponding to spurious right-
hand-side terms [26, 20]; see Supplement for details. For every identified
model $\mathbf{S}^{*}$, we also compute its recovery probability
$\mathrm{P}(\exists\lambda:\pm\mathbf{S}(\lambda)=\pm\mathbf{S}^{\ast})$ for
increasing sample size N.

Figure 2: Learning hydrodynamic equations for the SPP model. A: Stability plot
for learning polarity dynamics with dictionary ($\Theta$) size
$\text{N}=3000$, $\text{P}=14$. The blue and green regions indicate the
regularization ranges where the model $\mathbf{S}_{1,\mathbf{w}}$ is recovered
with $E_{\text{fp}}\leq 1$ and $E_{\text{fp}}\leq 2$, respectively. In the
yellow region, the model $\mathbf{S}_{2,\mathbf{w}}$ is recovered with
$E_{\text{fp}}\leq 3$ and in the violet region with $E_{\text{fp}}>3$. In all
cases, the threshold is $\pi_{\text{th}}=0.95$ (red dashed line). Line markers
identify model terms as given in the inset legend. B: Stability plot for
learning density dynamics with $\text{N}=3000$, $\text{P}=10$. In the colored
regularization range, the model $\mathbf{S}_{1,\rho}$ is recovered with
$E_{\text{fp}}\leq 1$ for $\pi_{\text{th}}=0.8$, and the correct average lane
speed ($\approx 2.1$) is estimated. C: Achievability plot to gauge the
identifiability of the models $\mathbf{S}_{1,\rho}$,
$\mathbf{S}_{1,\mathbf{w}}$,
$\mathbf{S}_{2,\mathbf{w}},\mathbf{S}_{A,\mathbf{w}}$ for increasing sample
size N. Colored bands are Bernoulli standard deviations for 20 independent
trials. D: Inferred latent density dependence of the particle speed $v_{p}$
for varying $\gamma$ and initial density $\rho_{0}=30$. The dotted lines show
the true relation $v_{p}(\rho)=v_{0}\mathrm{e}^{-\gamma\rho\pi}+v_{1}$.

Application to the self-propelled particles (SPP) model: We consider the self-
propelled stochastic particle system with simple alignment interactions [6]:

$\dot{\bm{r}}_{p}=v_{p}\bm{e}_{\theta_{p}};\>\dot{\theta}_{p}=\beta\sum_{q=1}^{N_{p}}F(\theta_{q}-\theta_{p},r_{qp})+\sqrt{2\epsilon}\,\eta_{p}(t),\vspace{-1em}$
(3)

where $\bm{e}_{\theta_{p}}=\\{\cos\theta,\sin\theta\\}$ is the alignment
vector of particle $p$ moving at speed $v_{p}=\|\bm{v}_{p}\|_{2}$, and
$r_{pq}=\|\bm{r}_{q}-\bm{r}_{p}\|_{2}$. The parameters $\beta$ and $\epsilon$
describe the alignment and fluctuation strengths, respectively; $\eta_{p}(t)$
is a Gaussian white noise with zero mean and unit variance. The particle
alignment function has the form $F(\theta,r_{pq})=\sin(\theta)/\pi R^{2}$ if
$r_{pq}<R$ and $0$ otherwise. Density dependent motility is introduced by
making the particle speed $v_{p}$ depend on the local density $\rho$ as
$v_{p}(\rho)=v_{0}\mathrm{e}^{-\gamma\rho\pi}+v_{1}$, where $v_{0}$ and
$v_{1}$ are the speeds in the dilute and crowded limits, respectively. The
strength of this dependence is set by the parameter $\gamma$. This simple
microscopic model is able to generate a variety of patterns, ranging from
propagating density lanes to moving particle clumps and aster formation. None
of those are observed in the standard Vicsek model [9] nor its hydrodynamic
limit [2, 27].

We start by exploring the SPP model close to the Vicsek limit (here
$\gamma=10^{-5}$), where existing hydrodynamic theories are available to
compare with [2]. Figure 2A shows $\widehat{\Pi}(\lambda)$ for the
hydrodynamic polarization density $\mathbf{w}(\bm{r})=\int
f(\bm{r},\theta)\mathbf{e}_{\theta}\,\mathrm{d}\theta$ and the density
$\rho(\bm{r})=\int f(\bm{r},\theta)\,\mathrm{d}\theta$ with particle
distribution function
$f(\bm{r},\theta)=\sum_{p=1}^{N_{p}}\delta(\bm{r}-\bm{r}_{p})\delta(\theta-\theta_{p})$.
Above a model probability threshold of $\pi_{\text{th}}=0.95$, we observe two
models for polarization density
$\mathbf{S}_{1,\mathbf{w}},\mathbf{S}_{2,\mathbf{w}}$ (see Table 1) that are
consistent and robust with varying levels of fidelity (shown by different
color shades) across different $\lambda$ values.For the density field (see
Fig. 2B), we recover the continuity equation
$\partial_{t}\rho=d_{4}\nabla\cdot\mathbf{w}$ as model $\mathbf{S}_{1,\rho}$
with the coefficient $d_{4}$ approximating the constant particle speed, i.e.,
$d_{4}\approx v$.

$\mathbf{S}_{1,\rho}$ | $\partial_{t}\rho=$ | $c_{4}\nabla\cdot\mathbf{w}$ |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---
$\mathbf{S}_{1,\mathbf{w}}$ | $\partial_{t}\mathbf{w}=$ | $d_{3}\nabla\rho$ | $+d_{12}\mathbf{w}\cdot\nabla\mathbf{w}$ | $+d_{7}\nabla\textrm{w}^{2}$ |  |  |  |  |
$\mathbf{S}_{2,\mathbf{w}}$ | $\partial_{t}\mathbf{w}=$ | $d_{3}\nabla\rho$ | $+d_{12}\mathbf{w}\cdot\nabla\mathbf{w}$ | $+d_{7}\nabla\textrm{w}^{2}$ | $+d_{2}\rho\mathbf{w}$ | $+d_{1}\mathbf{w}$ | $+d_{4}\textrm{w}^{2}\mathbf{w}$ |  |
$\mathbf{S}_{A,\mathbf{w}}$ | $\partial_{t}\mathbf{w}=$ | $d_{3}\nabla\rho$ | $+d_{12}\mathbf{w}\cdot\nabla\mathbf{w}$ | $+d_{7}\nabla\textrm{w}^{2}$ | $+d_{2}\rho\mathbf{w}$ | $+d_{1}\mathbf{w}$ | $+d_{4}\textrm{w}^{2}\mathbf{w}$ | $+d_{5}\mathbf{w}\nabla\cdot\mathbf{w}$ |

Table 1: Comparison of hydrodynamic models derived from kinetic theory
($\mathbf{S}_{A,\mathbf{w}},\mathbf{S}_{1,\rho}$) with models
($\mathbf{S}_{1,\mathbf{w}},\mathbf{S}_{2,\mathbf{w}},\mathbf{S}_{1,\rho}$)
learned from microscopic simulation data of the self-propelled particle system
(Eq. 3) in the Vicsek limit $\gamma\ll 1$.

We check the statistical consistency of the learned models by computing their
recovery probability across all sparse models. We find that the three models
$\mathbf{S}_{1,\mathbf{w}},\mathbf{S}_{2,\mathbf{w}},\mathbf{S}_{1,\rho}$ can
be learned with recovery probability increasing with sample size N, see Fig.
2C. We compare this with the over-complete hydrodynamic model
$\mathbf{S}_{A,\mathbf{w}}$ (see Table 1) derived from kinetic theory [6],
which cannot be recovered without lowering $\pi_{\text{th}}$ or including more
false positives into the model set. This is because the convective term
$\mathbf{w}\nabla\cdot\mathbf{w}$ is statistically exchangeable with other
terms in the dictionary, since its dynamics are effectively captured by terms
like $\mathbf{w}\cdot\nabla\mathbf{w}$ and $\nabla|\mathbf{w}|^{2}$.
Therefore, our approach finds a simpler model that is consistent with physical
priors and is robustly recoverable from the data.

Numerical simulation of the model
$\mathbf{S}_{1,\rho}\\!+\mathbf{S}_{2,\mathbf{w}}$ reveals propagating density
stripes or “lanes” that accurately capture the speed observed in microscopic
SPP simulations (see Supplementary Video). The sparser model
$\mathbf{S}_{1,\mathbf{w}}$ does not produce stripe patterns in simulations,
suggesting a necessity for alignment interactions through $\mathbf{w}$ and
$\mathbf{w}|\mathbf{w}|^{2}$.

Finally, we look at higher values of $\gamma$, where crowding effects are not
negligible. There, learning hydrodynamic models using generic dictionaries
without additional structural priors results in an ill-posed problem, as the
coefficients of the underlying hydrodynamic model can be a function of local
particle density. In our framework, however, we can explicitly accommodate for
this by grouping dictionary columns based on local density and using a block-
diagonal dictionary design. The problem then becomes well posed, and we
recover a stable model
$\mathbf{S}_{1,\rho}=d_{4}(\rho)\nabla\cdot\mathbf{w}|_{\rho_{k}},\>\forall
k=1,\ldots,|g|$, with a coefficient $d_{4}(\rho)$ that depends on density (see
Supplement). From dimensional analysis, we can interpret $d_{4}(\rho)$ as the
local particle speed depending (nonlinearly) on the density and the crowding
parameter $\gamma$ (see Fig. 2D). Therefore, this constitutive relation can be
learned automatically from data.

Figure 3: Learning hydrodynamic equations from stochastic cellular dynamics.
A: Stability plot for $\text{N}=5000$, $\text{P}=14$. The orange and yellow
regions indicate the regularization ranges where the model
$\mathbf{S}_{1,\bm{u}}$ is recovered with $E_{\text{fp}}\leq 1$ and
$E_{\text{fp}}\leq 2$, respectively. In the blue region, the model
$\mathbf{S}_{2,\bm{u}}$ recovered with $E_{\text{fp}}\leq 2$. Line markers
identify model terms as given in the inset legend. B: Achievability plot for
the models $\mathbf{S}_{1,\bm{u}},\mathbf{S}_{2,\bm{u}}$ for increasing sample
size N. Colored bands are Bernoulli standard deviations for 20 independent
trials. C: Comparison of the velocity magnitude $|\bm{u}|$ distribution from
both learned hydrodynamic models and the stochastic cellular simulation. D, E:
Same comparison for the velocity power spectrum and the velocity correlation
function (VCF). The vertical dashed line indicates the wavenumber
(${k}^{\prime}=2\pi/h$) associated with the bandwidth $h$ of the Gaussian
coarse-graining kernel.

Learning continuum models of collective cell dynamics: Collective cell
migration is a hallmark of morphogenetic events in living tissues, for example
during embryogenesis, wound healing, and cancer invasion [28]. We consider a
stochastic model that describes cells as particles moving with velocity
$\bm{v}_{p}$ [29, 30] in a thin epithelial tissue sheet, as governed by the
Langevin-like process:

$\frac{d\bm{v}_{p}}{dt}=-\alpha\bm{v}_{p}+\sum_{q\in\mathcal{N}(p)}\Big{[}\frac{\beta}{n_{p}}(\bm{v}_{q}-\bm{v}_{p})+\bm{f}_{pq}\Big{]}+\sigma(\rho_{p})\bm{\eta}_{p}$
(4)

for $p=1,\ldots,N_{p}$. Each particle $p$ has $n_{p}$ neighbors within a
neighborhood $\mathcal{N}(p)$. The pair-wise forces between particles are
modeled through $\bm{f}_{pq}$, and $\alpha,\beta$ control the substrate
friction and velocity alignment, respectively. Cell motion is driven by a
noise term of magnitude $\sigma(\rho)$, where $\rho$ is the particle density,
generated from an Ornstein-Uhlenbeck process $\bm{\eta}_{p}$ with correlation
time $\tau$; see Supplement for details. Since we are interested in learning a
mean-field hydrodynamic model of the process, we consider the regime of low
density fluctuations. We want to find a model for the coarse-grained velocity
field $\bm{u}(\bm{r})=\frac{\overline{v}(\bm{r})}{\rho(\bm{r})}\int
f(\bm{r},\theta)\mathbf{e}_{\theta}\,\mathrm{d}\theta$, where
$\overline{v}(\bm{r})$ is the speed at position $\bm{r}$ interpolated from the
closest particles.

$\mathbf{S}_{1,\bm{u}}$ | $\partial_{t}\bm{u}=$ | $d_{1}\bm{u}$ | $+d_{8}\Delta\bm{u}$ | $+d_{6}\nabla(\nabla\cdot\bm{u})$ |
---|---|---|---|---|---
$\mathbf{S}_{2,\bm{u}}$ | $\partial_{t}\bm{u}=$ | $d_{1}\bm{u}$ | $+d_{8}\Delta\bm{u}$ | $+d_{6}\nabla(\nabla\cdot\bm{u})$ | $+d_{4}\bm{u}|\bm{u}|^{2}$
Table 2: Minimal hydrodynamic models learned for the coarse-grained velocity
field $\bm{u}$ from simulation data of stochastic cellular dynamics (Eq. 4) in
epithelial tissue sheets.

The inference results in Fig. 3A suggest the existence of two stable models
$\mathbf{S}_{1,\bm{u}},\mathbf{S}_{2,\bm{u}}$ (see Table 2) that can be
recovered with high statistical robustness. The model $\mathbf{S}_{2,\bm{u}}$
has better recovery probability with increasing sample size N, as shown in
Fig. 3B. This model lacks the damping term $\bm{u}|\bm{u}|^{2}$. Analyzing the
models, the terms $\bm{u}$ and $\bm{u}|\bm{u}|^{2}$ control the order–disorder
transition of the system and set the magnitude of the velocity $\bm{u}$. The
terms $\Delta\bm{u}$ and $\nabla(\nabla\cdot\bm{u})$ can be interpreted as the
active counterparts of the bend and splay moduli from nemato-hydrodynamics
[31]. They microscopically originate from relative alignment interactions in
the particle system of Eq. 4.

We numerically validate the sparse models
$\mathbf{S}_{1,\bm{u}},\mathbf{S}_{2,\bm{u}}$ and find quantitative agreement
between the hydrodynamic models and the stochastic microscopic simulations.
The algorithmically inferred hydrodynamic equations are able to capture the
velocity magnitude distribution (Fig. 3C), power spectrum (Fig. 3D), and the
velocity correlation function
$\textrm{VCF}(|\bm{r}-\bm{r}^{\prime}|)=\big{\langle}\bm{u}(\bm{r},t)\bm{u}(\bm{r}^{\prime},t)\big{\rangle}_{t}$
(Fig. 3E) of the stochastic particle dynamics. We found the model
$\mathbf{S}_{2,\bm{u}}$ to be more numerical stable than
$\mathbf{S}_{1,\bm{u}}$ due to the presence of damping term
$\bm{u}|\bm{u}|^{2}$ that prevents velocity magnitude from exploding.

Learning hydrodynamic equations of non-equilibrium particle dynamics driven by
phoretic fields: We consider the Phoretic Brownian Particle (PBP) model [32,
33] with self-propelled particles driven by an imposed chemotactic
concentration field $c(\bm{r})$, i.e.

$\dot{\bm{r}}_{p}=\mathbf{e}_{\theta_{p}};\quad\dot{\theta}_{p}=\beta\mathbf{e}_{\theta_{p}}\times\nabla
c(\bm{r}_{p})+\sqrt{2}\,\eta_{p}(t),$ (5)

for $p=1,\ldots,N_{p}$. For $\beta>0$, the active particles turn towards the
phoretic gradients, for $\beta<0$ they turn against the gradient. This
alignment interaction is countered by rotational Brownian motion with Gaussian
white noise $\eta_{p}$ with zero mean and unit variance. Contrary to the two
previous examples, activity in the PBP model is governed by the imposed
phoretic field and is devoid of direct particle interactions.

We challenge our framework to recover a statistically consistent sparse
hydrodynamic model given only the positions $\bm{r}_{p}$ and orientations
$\theta_{p}$ of the particles, but withholding the underlying phoretic field
$c(\bm{r})$. Using symmetry arguments, we construct block diagonal
dictionaries that can account for a missing latent scalar or vector field (see
Supplement). For the density field $\rho$, we recover the term
$\nabla\cdot\mathbf{w}$ with an estimated coefficient $(c_{4}\approx 1)$
approximating the right self-advection speed of the particles. For the
polarization density $\mathbf{w}$, we consistently identify the model
$\mathbf{S}_{1,\mathbf{w}}$ given in Table 3. The estimated phoretic gradient
field closely approximates the true one (see Fig. 4A). The estimate is
obtained by smoothness-constrained least-squares regression; see Supplement.

In Fig. 4B–D, we compare the learned model $\mathbf{S}_{1,\mathbf{w}}$ with
the model $\mathbf{S}_{A,\mathbf{w}}$ (see Table 3) derived from kinetic
theory [32], and with direct PBP simulations. The agreement suggests that for
the particular parameter values chosen, the simpler model
$\mathbf{S}_{1,\mathbf{w}}$ is sufficient to describe the hydrodynamics of the
PBP system. As shown in the Supplement, the model $\mathbf{S}_{1,\mathbf{w}}$
can also be derived from kinetic theory under the assumption of low
anisotropy. Thus, at low anisotropy (i.e., small drift, small phoretic
gradient), it is not surprising that the model $\mathbf{S}_{1,\mathbf{w}}$ is
sufficient to capture the coarse-grained dynamics.

$\mathbf{S}_{1,\rho}$ | $\partial_{t}\rho=$ | $c_{4}\nabla\cdot\mathbf{w}$ |  |  |  |  |  |
---|---|---|---|---|---|---|---|---
$\mathbf{S}_{1,\mathbf{w}}$ | $\partial_{t}\mathbf{w}=$ | $d_{1}\mathbf{w}$ | $+d_{3}\nabla\rho$ | $+\rho\nabla c$ |  |  |  |
$\mathbf{S}_{A,\mathbf{w}}$ | $\partial_{t}\mathbf{w}=$ | $d_{1}\mathbf{w}$ | $+d_{3}\nabla\rho$ | $+\rho\nabla c$ | $+d_{7}\nabla\textrm{w}^{2}$ | $+(\nabla\mathbf{w})^{\top}\cdot\nabla c$ | $+\nabla c\cdot\nabla\mathbf{w}$ | $+(\nabla\cdot\mathbf{w})\nabla c$

Table 3: Hydrodynamic models of the PBP system (Eq. 5) learned from
microscopic simulation data ($\mathbf{S}_{1,\rho}$,
$\mathbf{S}_{1,\mathbf{w}}$) and derived from kinetic theory
($\mathbf{S}_{A,\mathbf{w}}$). Figure 4: Learning hydrodynamic equations for
the PBP model with $\mathrm{Pe}=1,\beta=0.3$. A: Smoothness-constrained
regression estimate of the latent phoretic field $\nabla c$ in the simulation
domain $[0,2\pi]\times[0,2\pi]$. The ground truth $\nabla c$ (red) compared
with the inferred phoretic field $\widetilde{\nabla c}$ (blue). The Relative
Mean Squared Error $\left(\frac{\|\nabla c-\widetilde{\nabla
c}\|_{2}^{2}}{\|\nabla c\|_{2}^{2}}\right)$ is $\text{RMSE}=0.2654$. B–D: Time
trajectories from the learned ($\mathbf{S}_{1,\mathbf{w}}$) and analytical
($\mathbf{S}_{A,\mathbf{w}}$) models of the density $\rho$ and both components
of the polarization density $\mathbf{w}$ at two random points in the spatial
domain with coordinates $(0.6,2.7)$ (dashed) and $(4.4,2.5)$ (solid),
highlighted by circle symbols in (A), compared with those from direct
microscopic PBP simulations.

Conclusions: We have presented a principled data-driven strategy for automated
learning of hydrodynamic equations from stochastic simulations of non-
equilibrium active particle dynamics. We inferred sparse hydrodynamic models
from simulation data of well-studied SPP and PBP models and found agreement
with analytical models derived by kinetic approaches. Numerical simulations of
the learned models reproduced relevant metrics of the stochastic particle
dynamics, providing quantitative validation. Importantly, our strategy allows
for encoding prior knowledge about the system via group sparsity. We have
shown that this enables recovery of hydrodynamic equations with spatially or
temporally varying coefficients and inference of spatially varying latent
fields.

Beyond benchmarking with existing analytical models, we applied our framework
to a stochastic non-equilibrium model of cell motility in living tissues,
where coarse-graining via kinetic approaches is not feasible due to the
analytically intractable particle interactions. Our data-driven framework
found two simple models with interpretable terms and long length- and time-
scale spectral characteristics in agreement with the stochastic microscopic
dynamics. For the PBP system, we used the learned models to identify a
sufficient closure for the kinetic approach to quantitatively reproduce the
particle dynamics.

Future work could incorporate stochastic force inference techniques [22, 23]
into our framework in order to learn hydrodynamic equations from experimental
observations of non-equilibrium particle dynamics.

###### Acknowledgements.

This work was supported by the German Research Foundation (Deutsche
Forschungsgemeinschaft, DFG) under Germany’s Excellence Strategy –
EXC-2068-390729961 – Cluster of Excellence “Physics of Life” of TU Dresden,
and by the Center for Scalable Data Analytics and Artificial Intelligence
(ScaDS.AI) Dresden/Leipzig, funded by the Federal Ministry of Education and
Research (Bundesministerium für Bildung und Forschung, BMBF).

## References

  * Couzin  _et al._ [2005] I. D. Couzin, J. Krause, N. R. Franks, and S. A. Levin, Effective leadership and decision-making in animal groups on the move, Nature 433, 513 (2005).
  * Toner  _et al._ [2005] J. Toner, Y. Tu, and S. Ramaswamy, Hydrodynamics and phases of flocks, Annals of Physics 318, 170 (2005).
  * Riedel  _et al._ [2005] I. H. Riedel, K. Kruse, and J. Howard, A self-organized vortex array of hydrodynamically entrained sperm cells, Science 309, 300 (2005).
  * Nédélec  _et al._ [1997] F. Nédélec, T. Surrey, A. Maggs, and S. Leibler, Self-organization of microtubules and motors, Nature 389, 305 (1997).
  * Karsenti [2008] E. Karsenti, Self-organization in cell biology: a brief history, Nature reviews Molecular cell biology 9, 255 (2008).
  * Farrell  _et al._ [2012] F. Farrell, M. Marchetti, D. Marenduzzo, and J. Tailleur, Pattern formation in self-propelled particles with density-dependent motility, Physical review letters 108, 248101 (2012).
  * Wensink  _et al._ [2012] H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, Meso-scale turbulence in living fluids, Proceedings of the national academy of sciences 109, 14308 (2012).
  * Schaller  _et al._ [2010] V. Schaller, C. Weber, C. Semmrich, E. Frey, and A. R. Bausch, Polar patterns of driven filaments, Nature 467, 73 (2010).
  * Vicsek  _et al._ [1995] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel type of phase transition in a system of self-driven particles, Physical review letters 75, 1226 (1995).
  * Bär  _et al._ [2020] M. Bär, R. Großmann, S. Heidenreich, and F. Peruani, Self-propelled rods: Insights and perspectives for active matter, Annual Review of Condensed Matter Physics 11, 441 (2020).
  * Großmann  _et al._ [2014] R. Großmann, P. Romanczuk, M. Bär, and L. Schimansky-Geier, Vortex arrays and mesoscale turbulence of self-propelled particles, Physical review letters 113, 258104 (2014).
  * Gao  _et al._ [2015] T. Gao, R. Blackwell, M. A. Glaser, M. D. Betterton, and M. J. Shelley, Multiscale polar theory of microtubule and motor-protein assemblies, Physical review letters 114, 048101 (2015).
  * Bertin  _et al._ [2009] E. Bertin, M. Droz, and G. Grégoire, Hydrodynamic equations for self-propelled particles: microscopic derivation and stability analysis, Journal of Physics A: Mathematical and Theoretical 42, 445001 (2009).
  * Kruse  _et al._ [2005] K. Kruse, J.-F. Joanny, F. Jülicher, J. Prost, and K. Sekimoto, Generic theory of active polar gels: a paradigm for cytoskeletal dynamics, The European Physical Journal E 16, 5 (2005).
  * Jülicher  _et al._ [2018] F. Jülicher, S. W. Grill, and G. Salbreux, Hydrodynamic theory of active matter, Reports on Progress in Physics 81, 076601 (2018).
  * Supekar  _et al._ [2021] R. Supekar, B. Song, A. Hastewell, A. Mietke, and J. Dunkel, Learning hydrodynamic equations for active matter from particle simulations and experiments, arXiv preprint arXiv:2101.06568 (2021).
  * Rudy  _et al._ [2017] S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Data-driven discovery of partial differential equations, Science Advances 3, e1602614 (2017).
  * Both  _et al._ [2021] G.-J. Both, S. Choudhury, P. Sens, and R. Kusters, Deepmod: Deep learning for model discovery in noisy data, Journal of Computational Physics 428, 109985 (2021).
  * Brunton  _et al._ [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 , 201517384 (2016).
  * Maddu  _et al._ [2021] S. Maddu, B. L. Cheeseman, C. L. Müller, and I. F. Sbalzarini, Learning physically consistent differential equation models from data using group sparsity, Physical Review E 103, 042310 (2021).
  * Helmuth  _et al._ [2010] J. A. Helmuth, G. Paul, and I. F. Sbalzarini, Beyond co-localization: inferring spatial interactions between sub-cellular structures from microscopy images, BMC bioinformatics 11, 1 (2010).
  * Frishman and Ronceray [2020] A. Frishman and P. Ronceray, Learning force fields from stochastic trajectories, Physical Review X 10, 021009 (2020).
  * Gnesotto  _et al._ [2020] F. S. Gnesotto, G. Gradziuk, P. Ronceray, and C. P. Broedersz, Learning the non-equilibrium dynamics of brownian movies, Nature communications 11, 1 (2020).
  * Reinbold  _et al._ [2021] P. A. Reinbold, L. M. Kageorge, M. F. Schatz, and R. O. Grigoriev, Robust learning from noisy, incomplete, high-dimensional experimental data via physically constrained symbolic regression, Nature communications 12, 1 (2021).
  * Romeo  _et al._ [2021] N. Romeo, A. Hastewell, A. Mietke, and J. Dunkel, Learning developmental mode dynamics from single-cell trajectories, eLife 10, e68679 (2021).
  * Meinshausen and Bühlmann [2010] N. Meinshausen and P. Bühlmann, Stability selection, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72, 417 (2010).
  * Mishra  _et al._ [2010] S. Mishra, A. Baskaran, and M. C. Marchetti, Fluctuations and pattern formation in self-propelled particles, Physical Review E 81, 061916 (2010).
  * Alert and Trepat [2020] R. Alert and X. Trepat, Physical models of collective cell migration, Annual Review of Condensed Matter Physics 11, 77 (2020).
  * Deforet  _et al._ [2014] M. Deforet, V. Hakim, H. Yevick, G. Duclos, and P. Silberzan, Emergence of collective modes and tri-dimensional structures from epithelial confinement, Nature communications 5, 1 (2014).
  * Sepúlveda  _et al._ [2013] N. Sepúlveda, L. Petitjean, O. Cochet, E. Grasland-Mongrain, P. Silberzan, and V. Hakim, Collective cell motion in an epithelial sheet can be quantitatively described by a stochastic interacting particle model, PLoS Comput Biol 9, e1002944 (2013).
  * De Gennes and Prost [1993] P.-G. De Gennes and J. Prost,  _The physics of liquid crystals_ , Vol. 83 (Oxford university press, 1993).
  * Liebchen  _et al._ [2017] B. Liebchen, D. Marenduzzo, and M. E. Cates, Phoretic interactions generically induce dynamic clusters and wave patterns in active colloids, Physical review letters 118, 268001 (2017).
  * Liebchen and Löwen [2018] B. Liebchen and H. Löwen, Synthetic chemotaxis and collective behavior in active matter, Accounts of chemical research 51, 2982 (2018).