Article https://doi.org/10.

1038/s41467-024-46598-w

Predicting multiple observations in

complex systems through low-dimensional
embeddings
1
Received: 11 November 2023 Tao Wu , Xiangyun Gao 2,3 , Feng An 4 , Xiaotian Sun2, Haizhong An2,3,
5,6 5,7
Zhen Su , Shraddha Gupta , Jianxi Gao 8,9 & Jürgen Kurths 5,7
Accepted: 4 March 2024

Forecasting all components in complex systems is an open and challenging

Check for updates task, possibly due to high dimensionality and undesirable predictors. We
1234567890():,;
1234567890():,;

bridge this gap by proposing a data-driven and model-free framework, namely,

feature-and-reconstructed manifold mapping (FRMM), which is a combination
of feature embedding and delay embedding. For a high-dimensional dynami-
cal system, FRMM finds its topologically equivalent manifolds with low
dimensions from feature embedding and delay embedding and then sets the
low-dimensional feature manifold as a generalized predictor to achieve pre-
dictions of all components. The substantial potential of FRMM is shown for
both representative models and real-world data involving Indian monsoon,
electroencephalogram (EEG) signals, foreign exchange market, and traffic
speed in Los Angeles Country. FRMM overcomes the curse of dimensionality
and finds a generalized predictor, and thus has potential for applications in
many other real-world systems.

Prediction of future states of a complex dynamical system is a chal- one-step and multistep ahead predictions of a target time series
lenging task across various disciplines1–3. System details are often variable16.
unknown, and only their time series data are accessible. Therefore, a Despite considerable efforts in the study of prediction tasks in
variety of data-driven techniques are designed for the prediction complex systems, it is still unsolved to design a generalized framework
task4,5, including traditional statistical models (e.g., autoregressive for the predictions of all components in a complex system. Since real-
integrated moving average (ARIMA))6, state space-based methods world systems often consist of many interconnected units, e.g., mul-
(e.g., sequential locally weighted global linear maps (S-maps)7 and tiple spatiotemporal observations in climate systems17 and thousands
multiview embedding (MVE))8, machine learning algorithms (e.g., of functionally connected neurons in the brain18, they therefore output
support vector machine (SVM)9, long short-term memory (LSTM)10, a large number of time series variables, and the interactions between
and reservoir computing (RC)11,12, and state-of-the-art combination these variables intrinsically contribute to the dynamical evolution of a
frameworks (e.g., multitask learning-based Gaussian process regres- complex system. A practical way to predict complex systems (espe-
sion machine (MT-GPRM)13, randomly distribution embedding cially for high-dimensional systems), as an approximation, is to study
(RDE)14 and autoreservoir neural network (ARNN)15). These advanced the dynamics of partial units, e.g., representative observations19.
approaches have shown potential for several significant tasks, e.g., However, identifying such representative variables remains a

1
College of Management Science, Chengdu University of Technology, Chengdu 610059, China. 2School of Economics and Management, China University of
Geosciences, Beijing 100083, China. 3Key Laboratory of Carrying Capacity Assessment for Resource and Environment, Ministry of Land and Resources,
Beijing 100083, China. 4School of Economics and Management, Beijing University of Chemical Technology, Beijing 100029, China. 5Potsdam Institute for
Climate Impact Research (PIK)–Member of the Leibniz Association, Potsdam 14473, Germany. 6Department of Computer Science, Humboldt University at
Berlin, Berlin 12489, Germany. 7Department of Physics, Humboldt University at Berlin, Berlin 12489, Germany. 8Department of Computer Science, Rensselaer
Polytechnic Institute, Troy, NY 12180, USA. 9Network Science and Technology Center, Rensselaer Polytechnic Institute, Troy, NY 12180, USA.
e-mail: gxy5669777@126.com; af15910602135@126.com; gaoj8@rpi.edu; kurths@pik-potsdam.de

Nature Communications | (2024)15:2242 1

Article https://doi.org/10.1038/s41467-024-46598-w

challenging task. Moreover, one should be cautious to ignore those original system from a single time series29. Given that low-dimensional
‘unimportant variables’, in which small perturbations may be amplified representations (in different coordinates) from two approaches show
and propagated to all components, resulting in heavy changes in sys- isomorphic structures with the original system. This enables predic-
tem behaviors (known as cascading effects)20,21. Instead, the capacity to tion tasks by one-to-one mapping between two low-dimensional
predict the future states of all components can help to better estimate representations. Additionally, in a dynamical system, each time series
the future behavior of a complex system. However, many existing variable can reconstruct a low-dimensional representation via delay
approaches present typical limitations for this task. (a) The uncertainty embedding30. Therefore, the low-dimensional representation from
of the predictor, which means that for a target variable, the predictors feature embedding can be practically selected as a generalized pre-
are often selected empirically22, e.g., several target-related observa- dictor to potentially identify the future dynamics of all components in
tions. If regarding all the remaining variables as predictors, some complex systems.
redundant information may negatively affect the performance (e.g.,
noise and irrelevant variables to the target variable15), especially for Results
high-dimensional real-world systems8. (b) The uncertainty of the pre- Low-dimensional representation from delay embedding
dictive model, which means that for different targets, some approa- According to Takens’ embedding theory, it is possible to reconstruct a
ches may train different models (e.g., the predictive models for y1 and low-dimensional attractor by a single time series from a high-
y2 are completely independent)23,24. There may need N models for an N- dimensional dynamical system30. Particularly, for an N-dimensional
dimensional system, leading to a high computational cost. (c) The system M, one can reconstruct a topologically isomorphic manifold
challenge in forecasting multiple observations typically results in ver- M xi (namely, reconstructed manifold) from every time series x i ðtÞ
ifying methods over only a single or possibly few observations13,14,23. within the system (i = 1,2,    ,N, t = 1,2,    ,L, L is the length of the ser-
Therefore, designing a unified and reliable framework to forecast all ies), and each state point on M xi is represented as
components in complex systems is still an open and challenging issue. X~ i ðtÞ = ðx i ðtÞ,x i ðt + τÞ,    ,x i ðt + ðE  1ÞτÞÞ, where E is the embedding
In this work, we develop a data-driven and model-free framework dimension and τ is the time lag. For example, the attractors of the
by a combination of manifold learning and delay embedding, namely, 3-dimensional Lorenz system and Rössler system are reconstructed in
feature-to-reconstructed manifold mapping (FRMM). The FRMM fra- 2-dimensional space from individual time series (Fig. 1a, b, d, e).
mework yields reliable predictions for all components via a general- M xi has an isomorphic topological structure with the original
ized and practical predictor, i.e., the system’s low-dimensional system M. It indicates that for every state point X ðtÞ on M, one can find
representation from manifold learning (feature embedding). The the- a corresponding state point X~ i ðtÞ on M xi through a smooth mapping
oretical foundation of the FRMM is based on the ground truth that φi . According to Takens30, φi is a one-to-one mapping, we therefore
high-dimensional systems often contain redundant information and identify a corresponding state point X ðtÞ on M for every X~ i ðtÞ on M xi via
that their essential dynamics or structures can be characterized by low- the inverse mapping φ1 ðX~ i ðtÞÞ. These processes can be represented as
dimensional representations25–27, e.g., the meaningful structure of a (1).
4096-dimensional image (64 pixels by 64 pixels) can be characterized
in a three-dimensional manifold with two pose variables and an azi- φi : M ! M xi ,φi ðX ðtÞÞ = X~ i ðtÞ,φi 1 ðX~ i ðtÞÞ = X ðtÞ, ð1Þ
muthal lighting angle28. These low-dimensional representations can be
sufficiently identified from two powerful techniques: feature embed- where X ðtÞ = ðx 1 ðtÞ,x 2 ðtÞ,    ,x N ðtÞÞ and X ðtÞ 2 M,X~ i ðtÞ 2 M xi .
ding and delay embedding. (a) Feature embedding finds a low-
dimensional representation by preserving the geometric features (e.g., Low-dimensional representation from feature embedding
nearest neighbor) of the original system as much as possible25. (b) Delay embedding can reconstruct low-dimensional representations
Delay embedding reconstructs an isomorphic structure with the of the original systems. Additionally, such low-dimensional

Fig. 1 | Low-dimensional embeddings of complex systems. The dynamical 2-dimensional representations (e and f) of the 3-dimensional Rössler system (d).
structure of the 3-dimensional Lorenz system (a) is represented in 2-dimensional Even with additive noise, one can also find their low-dimensional embeddings (cf. SI
space via delay embedding (b) (from time series x, whereE = 2,τ = 10) and feature Appendix Fig. S1).
embedding (c) (i.e., diffusion map algorithm). Analogously, one can find

Nature Communications | (2024)15:2242 2

Article https://doi.org/10.1038/s41467-024-46598-w

representations can be obtained from manifold learning algorithms. data as training samples, and the others are test samples. Since the
For example, based on the diffusion map algorithm31,32, we find embedding dimension and time lag are E = 2 and τ = 10, FRMM yields
2-dimensional representations that show equivalent structures with reliable 10-step (T = ðE  1Þτ) ahead predictions for all units, where the
the 3-dimensional Lorenz and Rössler systems (Fig. 1c, f). These tech- average ρ reaches 0.91 and the average error remains at a low level
niques embed a high-dimensional system in a low-dimensional space (RMSE = 0:38). Still, FRMM achieves accurate 10-step-ahead predic-
by retaining the essential geometric features (e.g., the neighboring tions for all units in the Rössler system (SI Appendix Fig. S2).
points in high-dimensional space are also adjacent in the low- To further evaluate the FRMM framework in high-dimensional
dimensional representations) (feature embedding) of the original systems, we select the 90-dimensional coupled Lorenz system as a
system. Since the embedding is a one-to-one mapping31, it can be benchmark (see Methods)14. Traditional regression-based predictions
written as (2) encounter the “curse of dimensionality”. Some neural network-based
frameworks set all the observations as input, leading to relatively high
ϕ : M ! M 0 ,ϕðX ðtÞÞ = Y ðtÞ,ϕ1 ðY ðtÞÞ = X ðtÞ, ð2Þ computational costs. Due to the sensitivity of initial states as well as
complex nonlinear dynamics between units, predicting all compo-
where M 0 represents an E-dimensional manifold (namely, feature nents is indeed a challenging task. Our FRMM framework embeds the
manifold), and Y ðtÞ 2 M 0 ,X ðtÞ 2 M. Real-world systems often show 90-dimensional Lorenz system into a relatively lower space (E = 11) via
diverse dynamical structures and geometry features, and five feature embedding and delay embedding and sets the feature mani-
algorithms are selected alternatively to identify feature embedding, fold as a generalized predictor to find the future states of all compo-
i.e., isometric feature mapping (ISOMAP)28, locally linear embedding nents (Fig. 3d–f). FRMM performs reliably in that all the ρ values are
(LLE)33, Laplacian34, diffusion map31, and local tangent space alignment higher than 0.6 and all the errors are lower than 0.8 (Fig. 3f). The
(LTSA)35. More details are provided in Methods. average ρ reaches 0.74, and the average RMSE is 0.6.
Real-world systems are often influenced by various external fac-
Prediction via mapping between low-dimensional tors, and they may behave with time-dependent dynamics, e.g., the
representations couplings among components are not constant but time-varying. For
Through delay embedding and feature embedding, a high-dimensional this case, we evaluate the FRMM by setting the coupling in the 90-
system is represented by its two low-dimensional manifolds: the dimensional Lorenz system to be increased by 0.2 after ten time
reconstructed manifold (M x ) and the feature manifold (M 0 ). This intervals (see Methods)15. FRMM remains reliable with an average ρ of
indicates a one-to-one mapping between the feature manifold and the 0.73 and RMSE of 0.63 (Fig. 3g–i).
reconstructed manifold. Then, for every state point Y ðtÞ on M 0 , we can
find its corresponding state point X~ i ðtÞ on M xi by a smooth mapping Performance of FRMM on real-world systems
(e.g., ψi ). To illustrate the FRMM in real-world systems, we use several bench-
mark samples across different disciplines, including climate systems,
ψi : M 0 ! M xi ,ψi ðY ðtÞÞ = X~ i ðtÞ,i = 1,2,    ,N, ð3Þ neuroscience, financial systems, and traffic systems (details of all
datasets are provided in SI Appendix Chapter 1.2). a) For the climate
where ψi ðxÞ = φi ϕ1 ðxÞ(see Eqs. (1) and (2)). system, we consider the Indian monsoon, which is a typical phenom-
Note that X~ i ðtÞ = ðx i ðtÞ,x i ðt + τÞ,    ,x i ðt + ðE  1ÞτÞÞ, we deduce a enon that generates dramatic influences on India’s agriculture and
spatiotemporal transformation from state points on M 0 to a temporal economy37. Skillful ahead prediction of this phenomenon is of
series (final component in X~ i ðtÞ): importance. However, it remains a challenging task due to the complex
spatial-temporal interactions among multiple observations. To this
_
ψi ðy1 ðtÞ,y2 ðtÞ,    ,yE ðtÞÞ = x i ðt + ðE  1ÞτÞ,i = 1,2,    ,N,t = 1,2,    ,L, ð4Þ end, we select the lower-level (850 hPa) zonal daily wind component
from region IMI2 (70E-90E, 20N-30N), provided on a spatial grid with a
where ðy1 ðtÞ,y2 ðtÞ,    ,yE ðtÞÞ 2 M 0 and x i ðt + ðE  1ÞτÞ 2 X~ i ðtÞ 2 M xi . resolution of 10 × 10 38. Wind speeds interact spatially and form a 231-
When t = L,_it yields at mostðE  1Þτ-forward _
dynamics of each variable dimensional subsystem. The FRMM performs reliable 20-day forward
x i ðtÞ once ψi is identified (more details of ψi are provided in Methods). predictions for all observations, and the average ρ and RMSE are 0.86
In this work, _
we employ the classical Gaussian process regression to and 0.41, respectively (Fig. 4a–c and Supplementary Fig. S3). Addi-
train every ψi (cf. SI Appendix Chapter 1.3)36. To guarantee robustness, tionally, the FRMM is certified by monthly observations in the same
we validate the performance by randomly dividing the observed series region (SI Appendix Fig. S4).
into a training set and a test set (i.e., cross-validation). Two widely used b) In neuroscience, electroencephalogram (EEG) signals have
metrics are employed to measure the performance, i.e., the Pearson been extensively used to study the underlying mechanisms of the
correlation between observed values and predicted values (ρ) and the human brain as well as some typical diseases39. Therefore, ahead pre-
normalized root mean square error (RMSE) (normalized by the stan- dictions of EEG signals are expected to deliver efficient early warnings
dard deviation of input series)15. The main architecture of FRMM is for related diseases. EEG signals are often captured from different
given in Fig. 2. regions in the brain and show spatial-temporal dynamics. To test
FRMM, we utilize a 64-dimensional subsystem that consists of EEG
Performance of FRMM on model systems signal series from 64 channels from a healthy participant40. By setting
To illustrate the mechanism of the FRMM framework, we start with the the low-dimensional (E = 5) feature manifold of the system as a gen-
benchmark Lorenz system15. For the 3-dimensional ordinary Lorenz eralized predictor, we achieve accurate 20-second ahead predictions
system (see Methods), the FRMM first identifies its 2-dimensional for all signals, where the average ρ and RMSE are 0.92 and 0.42
manifolds (i.e., feature manifold and reconstructed manifold) via fea- (Fig. 4d–f and Supplementary Fig. S5).
ture embedding and delay embedding (Fig. 1a–c). Both feature mani- c) Financial systems are typically complex systems influenced
folds and reconstructed manifolds show isomorphic structures with by numerous internal and external factors via various channels,
the original system, which indicates an isomorphism between the resulting in high uncertainty and instability, which in turn makes
feature manifold and reconstructed manifold. Then, the feature prediction a difficult and challenging task41. We start with a 70-
manifold can be utilized as a generalized predictor for the predictions dimensional subsystem from the foreign exchange market, which
of three units. According to the generic cross-validation (see Meth- includes the daily closing prices of 70 currencies against the US
ods), we validate the performance by randomly selecting 50% of the dollar. FRMM performs accurate 20-day-ahead predictions for all

Nature Communications | (2024)15:2242 3

Article https://doi.org/10.1038/s41467-024-46598-w

Fig. 2 | Sketch of FRMM framework. To forecast all components in an N-dimen- one-to-one mapping between feature manifold M 0 and reconstructed manifold
_
sional system (a), we find its E-dimensional representations from delay embedding M x i . Then, it is possible to find a mapping ψi (i = 1,2,    ,N) from the feature mani-
(b) and feature embedding (c). Thus, an N-dimensional dynamical system M is fold M 0 to the final coordinate of the reconstructed manifold M x i . Therefore, the
represented by two isomorphic low-dimensional manifolds (i.e., feature manifold feature manifold M 0 can be utilized as a generalized predictor to find the future
M 0 and reconstructed manifold M xi ). The foundation of an isomorphism suggests a dynamics (purple elements) of all components in complex systems (d).

observations, where the average ρ and RMSE are 0.89 and 0.41 temporal dependencies among the elements of the system42. Our
(Fig. 4g–i and Supplementary Fig. S6). In addition, the FRMM out- FRMM achieves 10-step ahead predictions (ρ ≥ 0:6) for 86% (179) of the
puts reliable predictions of 46 stock indices from global stock components, where the average ρ and RMSE are 0.78 and 0.68 (the
markets; see SI Appendix Fig. S7. average accuracies for all components are 0.73 and 0.7) (Fig. 4j–l and
d) Finally, we demonstrate the FRMM in a 207-dimensional traffic Supplementary Fig. S8). As reported in Fig. 4l, the FRMM exhibits
system, which consists of the traffic speeds collected from 207 loop relatively poor performance for a few components, whose time series
detectors in Los Angeles County. Forecasting the time evolution of this involve many abrupt changes, like tipping points, possibly caused by
traffic system is a challenging task due to the complex spatial and rush hours or accidents.

Nature Communications | (2024)15:2242 4

Article https://doi.org/10.1038/s41467-024-46598-w

Fig. 3 | Model systems. a–c The performances of the 3-dimensional Lorenz system. demonstrate that the FRMM framework can output reliable multistep ahead pre-
d–f The performances of the 90-dimensional coupled Lorenz system, where the dictions for all components in the Lorenz system. Note: We validate the accuracy by
prediction accuracies of all 90 components are distributed in f, in which the red randomly selecting 50% of the series as a training sample (green shaded area), and
dotted line shows the average values of ρ and RMSE. g–i The predictions of the 90- the others are test samples. k represents the randomly selected k-th data. STD
dimensional coupled Lorenz system with time-varying dynamics. The results represents the standard deviation of the observed series.

Discussion long-term predictions. Despite this, our framework achieves at

Robustness tests most 30-step-ahead predictions for all components in the 90-
Generally, a predictive model performs better with longer train- dimensional ordinary Lorenz system (Fig. 5d).
ing samples and shorter test samples, and the performance
decreases sharply when training short samples but testing long Remark on feature embedding
samples. Our FRMM performs robustly even when inputting a Identifying the low-dimensional feature manifold is critical for
short training sample (10% (140)) and verifying on a longer test our prediction. However, embedding a false feature manifold
sample (90% (1260)), where the average ρ and RMSE are 0.69 and (which has an inequivalent topology to the original system) may
0.75 (Fig. 5a, Supplementary Figs. S10 and S11). In addition, the result in poor prediction due to the one-to-one mapping not
FRMM is robust with deteriorating noise (Fig. 5b and Supple- being maintained between the feature manifold and the recon-
mentary Fig. S12). For the length of input data, the FRMM also structed manifold. For example, LLE and Laplacian fail to identify
outputs reliable predictions with short input data (Fig. 5c). For the 2-dimensional feature manifold of the 3-dimensional Lorenz
the predicted step, the FRMM can find at most ðE  1Þτ steps. system, resulting in poor predictions (variable z) by adding them
Theoretically, each τ can be employed to reconstruct an iso- to the FRMM framework (Fig. 6a, b). Conversely, despite the
morphic attractor if the observed series is long enough30. Often, diffusion map algorithm, ISOMAP and LTSA also output a reliable
only limited data are accessible, leading to poor reconstruction feature manifold of the original attractor and could be used to
with an even larger τ. Therefore, our framework is unable to make perform accurate predictions for all variables by substituting the

Nature Communications | (2024)15:2242 5

Article https://doi.org/10.1038/s41467-024-46598-w

Fig. 4 | Performance in real-world datasets. The predictions of daily wind speed FRMM is shown reliable for accurate multistep predictions in representative real-
(m/s) (a–c), per second EEG signal (d–f), daily exchange rate (g–i), and traffic speed world systems, where the predicted horizons are T = 20 (a–i) and T = 10(j–l).
(5-minute interval) (j–l). By randomly selecting 50% of the data as a test sample, the

diffusion map in the FRMM framework (Fig. 6c, d). However, the Nevertheless, the five powerful techniques in this work make
diffusion map outperforms ISOMAP and LTSA. Due to the diver- sense in many high-dimensional systems.
sity of dynamic structures in various real-world systems and the
lack of sufficient details, only from their time series, there is no Comparison with traditional methods
golden rule to select an optimal feature embedding algorithm Many existing predictive models perform well when training on long
(More discussions are given in SI Appendix Chapter 1.7). samples and verifying on short samples, but the performances often

Nature Communications | (2024)15:2242 6

Article https://doi.org/10.1038/s41467-024-46598-w

Fig. 5 | Robustness tests. We conduct tests of length of training sample (a), ahead predictions, but remains challenging for even longer horizons (d). Overall,
additive noise (b), length of input series (c), and predicted step (d). The 90- the results demonstrate that the FRMM framework is robust for several funda-
dimensional ordinary Lorenz system is used here, and the performance of all mental factors. η represents the proportion of randomly selected training samples.
components is distributed as violins. FRMM performs better with longer training As with cross-validation, the longer the training sample is, the shorter the test
samples and shorter test samples, it remains reliable when training short samples sample. σ represents the strength of additive white noise. L is the length of the
but testing long samples (a). Given a target variable, FRMM can achieve multistep observed series. T denotes the predicted step.

decrease sharply with short training samples and long test samples. We models directly train and fit relations between time series from a
first compare the robustness of the length of training and test samples system, and they may perform poorly due to the inconstant corre-
with some traditional methods (e.g., ARIMA6, MVE8, SVM9, LSTM10, lations estimated from time series43 (the fitted parameters in a
RC11) and several advanced delay embedding-based frameworks (e.g., model are often time-varying). Instead, FRMM finds the mapping
ARNN15, RDE14, MT-GPRM13. As depicted in Fig. 7, ARIMA fails to predict between low-dimensional representations of a system, and this
variable z in the Lorenz system even when training on long samples, mapping is inherently supported. In summary, FRMM overcomes
whereas the other methods predict it accurately. Although the per- the curse of dimensionality, has higher interpretability, and shows
formance decreases as the training sample length decreases and the potential to be applied in various fields.
test sample length increases, FRMM always remains relatively robust Gaussian process regression is applied to find the mapping
compared to other methods. Second, we compare the performance of between the feature manifold and the reconstructed manifold. We
our FRMM with the aforementioned methods across different data- need to note that this mapping can be also trained by some neural
sets. The results indicate that FRMM yields relatively better predictions network algorithms, e.g., ARNN utilizes reservoir computing to train a
(Table 1). In summary, FRMM is shown more reliable for the predic- mapping from the original attractor to the delay attractor. Despite the
tions of all components in complex systems. satisfactory performance of neural networks, they often rely on suffi-
cient and rather large training samples. Besides, there remain unknown
Final remark on FRMM framework hidden details as black-box characters inside of some artificial neural
Our data-driven and model-free framework (FRMM) has been illu- networks. More importantly, the trade-offs between accuracy, cost,
strated by both representative models and real-world systems, and and interpretability are needed to be balanced in practical applica-
it has several advantages. First, FRMM performs predictions by tions. On this basis, it seems more satisfactory to integrate Gaussian
mapping between low-dimensional representations, which is well- process regression in our FRMM framework.
grounded in theory that the topological structure of a high- FRMM is developed based on a popular framework, namely spa-
dimensional dynamical system can be theoretically characterized tiotemporal information (STI) transformation44. Several advanced STI-
in low-dimensional space from delay embedding and feature based methods (e.g., MT-GPRM13, RDE14, and ARNN15) have been pro-
embedding. Second, for the uncertainties of the predictor and posed to predict various complex systems. FRMM shows individual
predictive model, FRMM sets the feature manifold as a generalized characteristics and meaningful improvements comparing with many
predictor to find future states of all components, and Gaussian existing STI-based methods. We clarify them from three aspects,
process regression is utilized as a fixed tool to train all mappings including the prediction task, the architecture, and the theoretical
between embedded manifolds. Third, many existing predictive foundation.

Nature Communications | (2024)15:2242 7

Article https://doi.org/10.1038/s41467-024-46598-w

Fig. 6 | Performance with other feature embedding techniques. Several repre- geometry of the original attractor, FRMM yields reliable predictions for all com-
sentative algorithms are considered, including LLE (a), Laplacian (b), ISOMAP (c), ponents. Several algorithms can be utilized for feature embedding in the
and LTSA (d). LLE and Laplacian algorithms fail to find faithful low-dimensional 3-dimensional Lorenz system, and an integration of Diffusion map performs the
representations of Lorenz attractor, resulting in poor performance for some best predictions among them (Fig. 3a).
components (e.g., variable z). While LTSA and ISOMAP preserve the fundamental

For the prediction task, it is still unsolved for the predictions of all highly related components as predictors. FRMM focuses on system’s
components in complex systems. Though some existing STI-based fundamental dynamics and sets the system’s low-dimensional feature
frameworks have the potential to address this issue, their abilities are manifold as a fixed and generalized predictor, which gives an efficient
often certified on partial components and not fully tested by all units in predictor when predicting different components in complex systems.
complex systems. Note that verifying a predictive model on fewer For the theoretical foundation, many existing STI-based fra-
observations of complex systems may be risky. Take the generic meworks create the STI equation by non-delay embedding and
3-dimensional Lorenz system as an example (Eq. (16)), it is possible to delay embedding, which originates from that a complex system
predict variables x and y through a linear regression model, but this can be approximately represented by different coordinates.
model fails to predict variable z45. It is uncritical to conclude that a Generally, the non-delay embedding of complex systems can be
linear regression model can predict the Lorenz system. In this direc- approximated in a space with either low or high dimension. (e.g.,
tion, FRMM is faithful and exhibits higher potential for the predictions MT-GPRM sets all selected observations as a representation of
of all components in complex systems. original systems, RDE finds non-delay embedding by randomly
For the architecture, the main difference between FRMM and selecting several observations). The theoretical foundation of
other STI-based frameworks is the selection of predictor. Some STI- FRMM is based on a well-accepted report that a high-dimensional
based frameworks set the original system as the fixed predictor, e.g., system often has redundant information, and the system’s fun-
MT-GPRM, while some frameworks may use different predictors for damental dynamics (e.g., the topology of complex systems) are
different targets, e.g., for each target variable, ARNN finds several restored in low-dimensional manifolds31–35. FRMM framework

Nature Communications | (2024)15:2242 8

Article https://doi.org/10.1038/s41467-024-46598-w

focuses on low-dimensional dynamics of complex systems, and abrupt, rapid, and even irreversible transitions (known as tipping
these low-dimensional dynamics are identified by feature points)46,47. The behavior of a system shifts between contrasting states,
embedding and delay embedding. The feature embedding is and the historical rules are often not held when a system crosses the
conducted by powerful manifold learning algorithms, and these threshold, leading to poor predictions of our framework. The phe-
methods can automatically extract and restore the fundamental nomenon of critical transitions, often caused by diverse external fac-
topology of the original system in a low-dimensional space. Thus, tors, is reported in numerous real-world systems. Despite this shortfall,
FRMM shows different theoretical foundations with existing STI- the FRMM framework also inspires to identify the tipping points of a
based frameworks. Additionally, identifying the fundamental system, e.g., the occurrence of poor performance may indicate an
dynamics of a high-dimensional system theoretically helps to underlying shift in the system.
reduce the negative impacts of redundant information in a high-
dimensional system, and will be beneficial for better predictions. Methods
These are also supported by the relatively higher performance FRMM framework
and robustness of FRMM in many real-world datasets (Table 1 Given an N-dimensional system with time series
and Fig. 7). x i ðtÞ(i = 1,2,    ,N,t = 1,2,    ,L), we aim to predict all units within the
However, like other STI-based frameworks, FRMM fails to predict system. For this task, the main structure of our FRMM framework is
different components synchronously. In other words, for each target listed as follows:
variable, one needs to train a suitable mapping. Additionally, FRMM (1) For each target variable x i ðtÞ, we estimate the embedding
also has limitations for situations in which a system experiences dimension E (cf. SI Appendix Chapter 3.3) and time lag τ based on the
false nearest neighbor algorithm and mutual information function,
respectively48,49. Then, this approach allows for reconstructing an
isomorphic manifold M xi in an E-dimensional space (E is usually much
smaller than N). The reconstructed manifold M xi is given as

2 3
x i ð1Þ x i ð1 + τÞ  x i ð1 + ðE  1ÞτÞ
6 . .. .. 7
6 . 7
6 . .  . 7
6 7
M xi = 6
6 x i ðhÞ x i ðh + τÞ  x i ðLÞ 7,
7 ð5Þ
6 . .. .. 7
6 .  7
4 . . . 5
x i ðLÞ x i ðL + τÞ  x i ðL + ðE  1ÞτÞ

where L is the length of the series and h = L  ðE  1Þτ. Note that the
elements from x i ð1Þ to x i ðLÞ are observed values from the original
system, others (x i ðL + τÞ,    ,x i ðL + ðE  1ÞτÞ) are unknown, and our goal
is to predict them. According to Takens, each time series variable can
be used to reconstruct an E-dimensional manifold30, which gives the
fundamental basis for predicting all components in complex systems.
(2) Moreover, the low-dimensional manifolds for the system can
also be identified by preserving their fundamental geometric features
(feature embedding). In this work, we provide several techniques to
find low-dimensional representations for the systems, e.g., isometric
Fig. 7 | Comparison of robustness with classic predictive models concerning
feature mapping (ISOMAP), locally linear embedding (LLE), Laplacian,
the length of the training sample, including ARIMA, SVM, MVE, LSTM, RC,
ARNN, RDE, and MT-GPRM. We compare the robustness of the one-step ahead diffusion map, and local tangent space alignment (LTSA), since real-
prediction of variable z from the 3-dimensional Lorenz system. Many methods world systems often behave with different dynamical structures and
show reliable predictions when inputting long training samples, while our FRMM have various geometric features.
remains robust when training on short samples and verifying on long samples. η We select LLE as an example to clarify the main idea of low-
gives the proportion of the training sample. dimensional embedding. Given an N-dimensional dimensional system
with observed vectors X i = ðx 1i ,x 2i ,    ,x Ni Þ, we approximate each point

Table 1 | Comparison of performance with several classic approaches

Real-world dataset Metric Method
FRMM ARNN RDE MVE SVM LTSM RC MT-GPRM ARIMA
Wind speed ρ 0.86 0.88 0.61 0.52 0.45 0.61 0.67 0.83 0.12
RMSE 0.41 0.51 0.79 0.97 0.88 0.85 0.7 0.52 0.89
EEG signal ρ 0.89 0.68 0.51 0.45 0.34 0.38 0.56 0.61 0.05
RMSE 0.6 0.74 0.79 1.03 0.91 0.82 0.89 0.87 1.07
Exchange rate ρ 0.89 0.81 0.66 0.61 0.42 0.51 0.67 0.84 0.23
RMSE 0.41 0.63 0.85 0.89 0.8 0.92 0.84 0.44 0.98
Traffic speed ρ 0.73 0.67 0.56 0.42 0.17 0.31 0.53 0.7 -0.18
RMSE 0.7 0.82 0.91 1.01 0.94 0.97 0.96 0.68 1.02
Note: We compare the average performance for all components in these systems, where 50% of the observed series are used as test samples (cross-validation). Two accuracy metrics are employed,
including ρ (Pearson correlation between predicted values and observed values, see the number in the first row for each dataset) and RMSE (Normalized by the standard deviation of the input series,
see the number in the second row for each dataset). All the datasets are the same to those in the main experiments in Results. All simulations are operated in Matlab 2018a, with the exception of MVE
prediction, which is conducted using package “rEDM’”, in R.

Nature Communications | (2024)15:2242 9

Article https://doi.org/10.1038/s41467-024-46598-w

_
by a linear function of its K nearest neighbors (e.g.,K = 8). In Eq. (12), ϕi can be easily obtained by the transform (13)
0 1 0 1
X
N X x i ð1Þ x i ð1 + τÞ    x i ð1 + ðE  1ÞτÞ 0 x i ð1 + ðE  1ÞτÞ
X~ i = wij X j , wij = 1, ð6Þ B . .. .. C . B .. C
B . C . B C
i=1 j B . .  . C . B . C
B C B C
B x i ðhÞ x i ðh + τÞ  x i ðLÞ C½ 0  = B x ðLÞ C:
where wij measures the weight_between the ith point and jth point. To B C B i C
_ B . .. .. C . B .. C
find the optimal set of weights W = ðwij Þ, we minimize the loss function B . C .. B C
@ . .  . A @ . A
x i ðLÞ x i ðL + τÞ    x i ðL + ðE  1ÞτÞ 1 x i ðL + ðE  1ÞτÞ
 2
_ N 
X X
N 
  ð13Þ
W = arg min X i  wij X j  : ð7Þ
 
i=1 j=1
According to Eqs. (11) and (12), we deduce a form (14)
We expect the local geometry in the original space to be preserved
_
in 0 1 0 1
_
their low-dimensional manifold. Therefore, we fix the matrix W = ðwij Þ y1 ð1Þ y2 ð1Þ  yE ð1Þ x i ð1 + ðE  1ÞτÞ
B . .. C B
.. C B .. C
and find the low-dimensional embedding by solving B . C
B . .  . C B . C
_B C B C
 2 ψi B y
B 1 ðhÞ y 2 ðhÞ  yE ðhÞ C B
C=B x i ðLÞ C,
C ð14Þ
_ N 
X X
N  B . C B C
 _  B . .. .. C B .. C
Y = arg min Y i  wij Y j  , ð8Þ @ . .  . A @ . A
Y  
i=1 j=1
y1 ðLÞ y2 ðLÞ  yE ðLÞ x i ðL + ðE  1ÞτÞ
_ _
where ψi ðxÞ = ϕi ψi ðxÞ. Equation (14) suggests a mapping from the fea-
where Y i represents the points in the low-dimensional manifold. Then,
ture manifold to the final coordinate of the reconstructed manifold.
the bottom E nonzero eigenvectors (from Eq. (8)) provide the low-
All the elements on M 0 (see the left matrix in Eq. (14)) are obtained
dimensional embedding M 0
via manifold learning algorithms, whereas partial components on M xi
2 3 are unknown, i.e., x i ðL + τÞ,    ,x i ðL + ðE  1ÞτÞ(see the right matrix in
y1 ð1Þ  yE ð1Þ
6 . 7 Eq. (14)), and they represent the future dynamics of the variable x i ðtÞ.
6
M 0 = 4 .. .. ... 7 ð9Þ
_
. 5, Once ψi is identified, one can find at most ðE  1Þτ-forward dynamics of
y1 ðLÞ ... yE ðLÞ a selected variable x i ðtÞ.
In a dynamical system, each time series variable can be used to
reconstruct a low-dimensional embedding (i.e., M xi , where
where yðtÞ 2 Y i . Consequently, each N-dimensional observation X i is
i = 1,2,    ,N).
_
This approach thus enables the construction of N map-
mapped to an E-dimensional point Y i .
pings (i.e., ψi ,i = 1,2,    ,N) from M 0 to the final coordinate ofM xi , which
(3) An N-dimensional system is embedded into an E-dimensional
yields multistep predictions for all units in high-dimensional dynamical
space from two different approaches, which then suggests a one-to-
systems. In this work, we _
use the Gaussian process regression algo-
one mapping between M xi and M 0 .
rithm to identify every ψi .
(4) We validate the performance by randomly dividing the
ψi : M 0 ! M xi ,ψi ðY ðtÞÞ = X~ i ðtÞ,i = 1,2,    ,N, ð10Þ observed series (x i ð1 + ðE  1ÞτÞ,x i ð2 + ðE  1ÞτÞ,    ,x i ðLÞ) into a training
set and a test set (i.e., cross-validation). The correlation between the
where Y ðtÞ 2 M 0 ,X~ i ðtÞ 2 M xi . From Eq. (10), we infer observed values and predicted values (ρ) and the normalized root
0 1 0 1 mean square error (RMSE) are applied to measure the performance.
y1 ð1Þ y2 ð1Þ  yE ð1Þ x i ð1Þ x i ð1 + τÞ    x i ð1 + ðE  1ÞτÞ
B . .. C B C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ffi
B . .. C B .. .. .. C P
B . .  . C B . .  . C covðx,x ~Þ 1
ðx  x~ Þ
2
B C B C n ð15Þ
ψi B    yE ðhÞ C B C: ρ= , RMSE = ,
B y1 ðhÞ y2 ðhÞ C = B x i ðhÞ x i ðh + τÞ  x i ðLÞ C ηx ηx~ ηx
B . .. .. C B C
B . C B .. .. .. C
@ . .  . A @ . .  . A where x and x~ are the original and predicted data, respectively. ηx
y1 ðLÞ y2 ðLÞ  yE ðLÞ x i ðLÞ x i ðL + τÞ    x i ðL + ðE  1ÞτÞ represents the standard deviation of the series x.
ð11Þ
Benchmark model systems
Since a one-to-one mapping between the feature manifold and the The coupled Lorenz system is defined as Eq. (16), where the ith
reconstructed manifold is held, it is possible to train a mapping from (i = 1,2,    ,N) subsystem is coupled with the (i−1) subsystem via c. To
the feature manifold to each coordinate of the reconstructed mani- make the system closed, we set i−1 as N for i = 1.
fold. In particular, for longer horizon predictions, we aim to find the
mapping from the feature manifold to the final coordinate of the x_ i = σðtÞðyi  x i Þ + cx i1 ,
reconstructed manifold (i.e., xðt + ðE  1ÞτÞ,t = 1,2,    ,L). These pro- y_ i = ax i  yi  x i z i , ð16Þ
cesses can be also explained mathematically, as follows.
z_ i = bz i + x i yi ,
Based on Eq. (11), we deduce the form (12)
0 1 0 1 whereσðtÞ is the time-varying parameter. a and b are set to be generic
x i ð1Þ x i ð1 + τÞ    x i ð1 + ðE  1ÞτÞ x i ð1 + ðE  1ÞτÞ values, i.e., a = 28,b =  8=3.
B . .. .. C B .. C For the time-invariant case (σðtÞ  10), Eq. (16) depicts an ordinary
B . C B C
B . .  . C B . C
_B C B C Lorenz system. Particularly, we obtain a 3-dimensional Lorenz system
ϕi B
B x i ðhÞ x i ðh + τÞ  x i ðLÞ C=B
C B x i ðLÞ C: ð12Þ
C when N = 1 and c = 0. We define a 90-dimensional coupled Lorenz
B . .. .. C B .. C
B . C B C system when N = 30 and c = 0:1.
@ . .  . A @ . A
For the time-varying case, we set σðtÞ to be increased (from an
x i ðLÞ x i ðL + τÞ    x i ðL + ðE  1ÞτÞ x i ðL + ðE  1ÞτÞ
initial value of 10) by 0.2 after every ten time intervals,
i.e., σðtÞ = 10 + 0:2ðtj10Þ.

Nature Communications | (2024)15:2242 10

Article https://doi.org/10.1038/s41467-024-46598-w

To generate the discrete data, we set the initial state as 0.1, and the 17. Fan, J. F. et al. Statistical physics approaches to the complex Earth
output series has a length of 1500. The first 100 data points are ignored system. Phys. Rep. 89, 1–84 (2021).
to avoid transient dynamics. We embed the 3-dimensional ordinary 18. Park, H. J. & Friston, K. J. Structural and functional brain networks:
Lorenz system into a 2-dimensional space, where E = 2,τ = 10. For the From connections to cognition. Science 342, 1238411 (2013).
90-dimensional ordinary and time-varying Lorenz systems, the 19. Fan, J. Q. & Lv, J. C. A selective overview of variable selection in high
embedding dimension and time lag are E = 11 and τ = 1, respectively. dimensional feature space. Stat. Sin. 20, 101–148 (2010).
The diffusion map algorithm is used to find their low-dimensional 20. Liu, X. M. et al. Network resilience. Phys. Rep. 971, 1–108 (2022).
representations. 21. Gao, J., Barze, B. & Barabási, A. L. Universal resilience patterns in
complex networks. Nature 530, 307–312 (2016).
Reporting summary 22. Robertson, D. E. & Wang, D. J. Bayesian approach to predictor
Further information on research design is available in the Nature selection for seasonal strearnflow forecasting. J. Hydrometeorol. 13,
Portfolio Reporting Summary linked to this article. 155–171 (2012).
23. Ben Taieb, S., Bontempi, G., Atiya, A. F. & Sorjamaa, A. A review and
Data availability comparison of strategies for multi-step ahead time series fore-
The details of real-world datasets are listed in SI Appendix Chapter 1.2. casting based on the NN5 forecasting competition. Expert Syst.
All real-world datasets are available at https://github.com/wt1234wt/ Appl. 39, 7067–7083 (2012).
FRMM-framework. The details of datasets from model systems are 24. Zhang, Y. R., Zhang, Y. L. & Haghani, A. A hybrid short-term traffic
given in Methods. flow forecasting method based on spectral analysis and statistical
volatility model. Transp. Res. Part C Emerg. Technol. 43,
Code availability 65–78 (2014).
The related codes are available at https://github.com/wt1234wt/ 25. Whitney, H. Differentiable manifolds. Ann. Math. 37, 645–680
FRMM-framework. (1936).
26. Seung, H. S. & Lee, D. D. The manifold ways of perception. Science
References 290, 2268–2269 (2000).
1. Clauset, A., Larremore, D. B. & Sinatra, R. Data-driven predictions in 27. Busch, E. L. et al. Multi-view manifold learning of human brain-state
the science of science. Science 355, 477–480 (2017). trajectories. Nat. Comput. Sci. 3, 240–253 (2023).
2. Subrahmanian, V. S. & Kumar, S. Predicting human behavior: The 28. Tenenbaum, J. B., de Sivlar, V. & Langford, J. C. A global geometric
next frontiers. Science 355, 489–489 (2017). framework for nonlinear dimensionality reduction. Science 290,
3. Perretti, C. T., Munch, S. B. & Sugihara, G. Model-free forecasting 2319–2323 (2000).
outperforms the correct mechanistic model for simulated and 29. Sugihara, G. et al. Detecting causality in complex ecosystems.
experimental data. PNAS 110, 5253–5257 (2013). Science 338, 496–500 (2012).
4. Wang, W. X., Lai, Y. C. & Grebogi, C. Data based identification and 30. Takens, F. Detecting strange attractors in turbulence. Mathematics
prediction of nonlinear and complex dynamical systems. Phys. Rep. 898, 366–381 (1981).
644, 1–76 (2016). 31. Coifman, R. R. & Lafon, S. Diffusion maps. Appl. Comput. Harmonic
5. Munch, S. B. et al. Constraining nonlinear time series modeling with Anal. 21, 5–30 (2006).
the metabolic theory of ecology. PNAS 120, e2211758120 (2023). 32. Floryan, D. & Graham, M. D. Data-driven discovery of intrinsic
6. Box, G. E. P. & Pierce, D. A. Distribution of residual autocorrelations dynamics. Nat. Mach. Intell. 4, 1113–1120 (2022).
in autoregressive integrated moving average time series models. J. 33. Roweis, S. T. & Saul, L. K. Nonlinear dimensionality reduction by
Am. Stat. Assoc. 65, 1509–1526 (1970). locally linear embedding. Science 29, 2323–2326 (2000).
7. Sugihara, G. Nonlinear forecasting for the classification of natural 34. Belkin, M. & Niyogi, P. Laplacian eigenmaps for dimensionality
time series. Philos. Trans. Phys. Sci. Eng. 348, 477–495 (1994). reduction and data representation. Neural Comput. 15,
8. Ye, H. & Sugihara, G. Information leverage in interconnected eco- 1373–1396 (2003).
systems: Overcoming the curse of dimensionality. Science 353, 35. Zhang, Z. Y. & Zha, H. Y. Principal manifolds and nonlinear dimen-
922–925 (2016). sion reduction via local tangent space alignment. SIAM J. Sci.
9. Rebentrost, P., Mohseni, M. & Lloyd, S. Quantum support vector Comput. 26, 313–338 (2004).
machine for big data classification. Phys. Rev. Lett. 113, 36. Rasmussen, C. & Williams, C. Gaussian processes for machine
130503 (2014). learning (MIT Press, Cambridge, MA) (2006).
10. Hochreiter, S. & Schmidhuber, J. Long short-term memory. Neural 37. Gade, S. V. et al. Impact of the ensemble Kalman filter based cou-
Comput. 9, 1735–1780 (1997). pled data assimilation system on seasonal prediction of Indian
11. Pathak, J., Hunt, B., Girvan, M., Lu, Z. & Ott, E. Model-free prediction summer monsoon rainfall. Geophys. Res. Lett. 49, e2021GL097184
of large spatiotemporally chaotic systems from data: A reservoir (2022).
computing approach. Phys. Rev. Lett. 120, 024102 (2018). 38. Hersbach, H. et al. The ERA5 global reanalysis. Q. J. R. Meteorol.
12. Gauthier, D. J., Bollt, E., Griffith, A. & Barbosa, W. A. S. Next gen- Soc. 1316, 1999–2049 (2020).
eration reservoir computing. Nat. Commun. 12, 5564 (2021). 39. Hassan, M. & Wendling, F. Aiming for high resolution of brain net-
13. Tao, P. et al. Predicting time series by data-driven spatiotemporal works in time and space Electroencephalography Source Con-
information transformation. Inf. Sci. 622, 859–872 (2023). nectivity. IEEE Signal Process. Mag. 35, 81–96 (2018).
14. Ma, H. F., Leng, S. Y., Aihara, K., Lin, W. & Chen, L. N. Randomly 40. Jao, P. K., Chavarriaga, R. & Millan, J. D. EEG-based online regulation
distributed embedding making short-term high-dimensional data of difficulty in simulated flying. IEEE Trans. Affect. Comput. 14,
predictable. PNAS 43, E9994–E10002 (2018). 394–405 (2023).
15. Chen, P., Liu, R., Aihara, K. & Chen, L. N. Autoreservoir computing 41. Battiston, S. et al. Complexity theory and financial regulation. Sci-
for multistep ahead prediction based on the spatiotemporal infor- ence 351, 818–819 (2016).
mation transformation. Nat. Commun. 11, 4568 (2020). 42. Avila, A. M. & Mezic, I. Data-driven analysis and forecasting of
16. Chen, C., Li, R., He, Z. & Chen, L. N. Predicting future dynamics from highway traffic dynamics. Nat. Commun. 11, 2090 (2020).
short-term time series using an Anticipated Learning Machine. Natl 43. Wu, T. et al. Universal window size-dependent transition of corre-
Sci. Rev. 7, 1079–1091 (2020). lations in complex systems. Chaos 33, 023111 (2023).

Nature Communications | (2024)15:2242 11

Article https://doi.org/10.1038/s41467-024-46598-w

44. Tong, Y. Y. et al. Earthquake alerting based on spatial geodetic data Additional information
by spatiotemporal information transformation learning. PNAS 120, Supplementary information The online version contains
e2302275120 (2023). supplementary material available at
45. Wu, T. et al. A novel framework for direct multistep prediction in https://doi.org/10.1038/s41467-024-46598-w.
complex systems. Nonlinear Dyn. 111, 9289–9304 (2023).
46. Scheffer, M. et al. Anticipating critical transitions. Science 338, Correspondence and requests for materials should be addressed to
344–348 (2012). Xiangyun Gao, Feng An, Jianxi Gao or Jürgen Kurths.
47. Grziwotz, F. et al. Anticipating the occurrence and type of critical
transitions. Sci. Adv. 9, eaba4558 (2023). Peer review information Nature Communications thanks Pei Chen,
48. Krakovska, A., Mezeiova, K. & Budacova, H. Use of false nearest Dibakar Ghosh, and the other, anonymous, reviewer(s) for their con-
neighbors for selecting variables and embedding parameters for tribution to the peer review of this work. A peer review file is available.
state space reconstruction. J. Complex Syst. 2015, 932750 (2015).
49. Fraser, A. M. & Swinney, H. L. Independent coordinates for strange Reprints and permissions information is available at
attractors from mutual information. Phys. Rev. A. 33, http://www.nature.com/reprints
134–1140 (1986).
Publisher’s note Springer Nature remains neutral with regard to jur-
Acknowledgements isdictional claims in published maps and institutional affiliations.
X.Y.G. is supported by the National Natural Science Foundation of China
(No.72371229, 71991481, 71991485), and by the Fundamental Research Open Access This article is licensed under a Creative Commons
Funds for the Central Universities (3-7-8-2023-01). H.Z.A. is supported by Attribution 4.0 International License, which permits use, sharing,
the National Natural Science Foundation of China (No.71991481, adaptation, distribution and reproduction in any medium or format, as
71991480], and by the Basic Science Center Project of the National long as you give appropriate credit to the original author(s) and the
Natural Science Foundation of China (No. 72088101). J.X.G. is supported source, provide a link to the Creative Commons licence, and indicate if
by the National Science Foundation (No. 2047488), and by the Rensse- changes were made. The images or other third party material in this
laer- IBM AI Research Collaboration. article are included in the article’s Creative Commons licence, unless
indicated otherwise in a credit line to the material. If material is not
Author contributions included in the article’s Creative Commons licence and your intended
T.W., X.Y.G., F.A., J.X.G. and J.K. designed the research. T.W., X.Y.G. and use is not permitted by statutory regulation or exceeds the permitted
F.A. performed the experiments and wrote the first draft of the paper. use, you will need to obtain permission directly from the copyright
T.W., X.T.S., S.G. and Z.S. analyzed data. T.W., X.Y.G., F.A., H.Z.A., J.X.G. holder. To view a copy of this licence, visit http://creativecommons.org/
and J.K. reviewed and edited the manuscript. licenses/by/4.0/.

Competing interests © The Author(s) 2024

The authors declare no competing interests.

Nature Communications | (2024)15:2242 12