Given a set of high-dimensional nonlinear multi-omics time series states F = ( F 1 , F 2 , … , F T ) , where T represents the total step, the time series state at t can be described as F t = ( f 1 t , f 2 t , … f n t ) ′ , where n represents the dimension of the time series state, “ ' ”, as the transpose of a vector. Our goal is to predict the future state of the multi-omics time series. Next, we detail how to develop an EKATP as follows.

Since an EKATP is based on the Autoencoder framework, we employ Eq. 1 to define the objective function for Autoencoder ( L i d ). L i d = ‖ F ^ t − F t ‖ MSE (1)

Here, F ^ t is the reconstructed high-dimensional time series state according to encoder ( χ e ) and decoder ( χ d ) of Autoencoder (Supplementary Presentation S1). | | ⋅ | | MSE denotes the mean squared error (MSE), which presents the expected value of the square of the difference between the predicted value and the true value. This loss function term enables us to construct an Autoencoder model that satisfies χ d ∘ χ e ≈ i d , the identity.

According to the description in the delay embedding theory (Supplementary Presentation S1), we employ χ e of the Autoencoder to approximate the delay embedding Φ , mapping the high-dimensional nonlinear input time series state F t back to the low-dimensional time series state Y t by Eq. 2, Y t = ( y t , y t + 1 , ... , y t + L − 1 ) ′ = χ e ( F t ) . (2) where L represents the dimension of the low-dimensional time series state. Similarly, the inverse mapping χ d of mapping χ e is used to approximate the conjugate form of delay embedding Φ , mapping the low-dimensional time series state back to the high-dimensional time series state by Eq. 3. F ^ t = ( f ^ 1 t , f ^ 2 t , … f ^ n t ) ′ = χ d ( Y t ) . (3)

Based on the Koopman theory discussed by Supplementary Presentation S1, we construct the finite dimensional linear matrix C (and matrix D ) to compute the forward (and backward) low-dimensional time series state. Equation 4 shows how to realize the forward prediction for low-dimensional time series state Y t to obtain Y t + 1 . Y t + 1 = C Y t . (4)

Equation 4 can be expanded by Eq. 5. [ y t + 1 y t + 2 y t + 3 … y t + L − 1 y t + L ] = [ 0 1 … 0 0 0 0 … 0 0 0 0 … 0 0 … … … … … 0 0 … 0 1 a 1 a 2 … a L − 1 a L ] [ y t y t + 1 y t + 2 … y t + L − 2 y t + L − 1 ] . (5)

Here, a i is the estimated parameter that needs training, and a 1 ≠ 0 . Equation 6 shows how to realize the backward prediction for a low-dimensional time series state Y t to obtain Y t − 1 . Y t − 1 = D Y t . (6)

Equation 6 can be expanded by Eq. 7. [ y t − 1 y t y t + 1 … y t + L − 3 y t + L − 2 ] = [ b 1 b 2 … b L − 1 b L 1 0 … 0 0 0 1 … 0 0 … … … … … 0 0 … 0 0 0 0 … 1 0 ] [ y t y t + 1 y t + 2 … y t + L − 2 y t + L − 1 ]   . (7)

Here, b i is the estimated parameter that needs training, and b L ≠ 0 . Our goal is to optimise the parameters of the linear matrix C (and D ) of Eqs 5, 7 by model training.

We make the k -steps forward prediction by Eq. 8 and backward prediction by Eq. 9 for the state of the low-dimensional time series Y t . After that, χ d is used to map the low-dimensional predictive time series state back to the high-dimensional predictive time series state by Eq. 10, Y t + k = C k Y t . (8) Y t − k = D k Y t . (9) F ^ t ± k = χ d ( Y t ± k ) . (10) where Y t + k and Y t − k represent the low-dimensional state after k steps of forward and backward prediction, respectively. F ^ t ± k represents the predictive high-dimensional nonlinear state.

Equations 11, 12 define the loss function of forward prediction ( L f w d ) and backward prediction ( L b w d ) to minimize the difference between the high-dimensional predictive value and true states at each step, respectively. L fwd = 1 k ∑ s = 1 k | | F ^ t + s − F t + s | | MSE . (11) L bwd = 1 k ∑ s = 1 k | | F ^ t − s − F t − s | | MSE . (12)

Equation 13 defines the loss function ( L i d y ) to minimize the difference between the predictive low-dimensional state obtained by the C and D matrices and defines such a low-dimensional state that is mapped from the true high-dimensional state by mapping χ e . L idy = 1 k ∑ s = 1 k [ ‖ C s χ e ( F t ) − χ e ( F t + s ) ‖ MSE + ‖ D s χ e ( F t ) − χ e ( F t − s ) ‖ MSE ] . (13)

Additionally, we employ loss function ( L c o n ) by Eq. 14 to train the parameters a i and b i in the matrices C and D , respectively. L con = 1 k ∑ s = 1 k [ ‖ χ d ( D s C s Y t ) − F t ‖ MSE + ‖ χ d ( C s D s Y t ) − F t ‖ MSE ] (14)

Equation 15 optimizes the key parameters for the EKATP by minimizing L . L = λ i d L i d + λ f w d L f w d + λ b w d L b w d + λ i d y L i d y + λ c o n L c o n . (15)

Here, λ i d , λ fwd , λ bwd , λ idy and λ con are user-defined hyperparameters.

We usually employ the chaotic Lorentz system (Lorenz, 1963) to describe a gene expression time series with a low-dimensional manifold (Sauer et al., 1991) by Eq. 16, { x t + 1 = x t + h ( η ( y t − z t ) ) y t + 1 = y t + h ( x t ( ρ − z t ) − y t ) z t + 1 = z t + h ( x t y t − β z t ) , (16) where η and ρ represent the Prandtl and Rayleigh numbers, respectively. β is related to geometry, and t represents time. h represents the level of the complicated nonlinear system. When h is greater, the nonlinear relationship between genes becomes more complicated.

Since gene expression time series contains considerable noise, we employ white Gaussian noise (Li et al., 2017) to simulate the noise by Eq. 17, { x ˜ = x + ε x y ˜ = y + ε y z ˜ = z + ε z . (17) where x ˜ , y ˜ and z ˜ represent data with noise. ε x , ε y and ε z represent the white Gaussian noise for x , y and z by normal distributions N ( 0 , σ 2 ) with a zero mean and a standard deviation σ . The standard deviation σ is referred to as the noise intensity.

Here, we describe how to obtain a high-dimensional gene expression time series with a low-dimensional manifold as follows. First, we generate the three-dimensional time series V = ( V 1 , V 2 , … , V T ) ∈ ℝ 3 ( T is the total step), which is listed in Supplementary Tables S1.1, S1.2, S1.3. Next, we develop a random orthogonal transformation (Anderson et al., 1987) matrix P ∈ ℝ 96 × 3 . Finally, we map the state of a 3-dimensional time series onto the state of a 96-dimensional time series by Eq. 18 to simulate a high-dimensional gene expression time series F = ( F 1 , F 2 , … , F T ) ∈ ℝ 96 with a 3-dimensional manifold, which is listed in Supplementary Tables S1.4, S1.5, S1.6. F = P V . (18)

To prove the accuracy and robustness of the EKATP, we generate a small-scale system containing T = 1,050 steps and choose the last 50 steps to visualize the predictive power of the EKATP.

Figure 2 shows the predictive error in the range of 50 steps under different initial conditions and environments. Detailed information is listed in Supplementary Tables S1.7, S1.8; Supplementary Presentation S4.

Figures 2A,C demonstrates that the EKATP not only has less of a predictive error than the existing methods under a clean environment ( σ =0.00) but also has a stable predictive error when the complexity h increases from 0.003 to 0.006. In particular, with the increase in predictive steps, the predictive error of the EKATP increases slower than that of the existing methods.

Figures 2B,D shows that the EKATP not only has less of a predictive error than previous methods under a noisy environment ( σ =0.01) but also has a predictive error that slightly fluctuates when h increases from 0.003 to 0.006. Moreover, after 25 steps, the predictive error of the EKATP increases much slower than that of the existing methods.

Figure 2 indicates that the EKATP has greater predictive accuracy and robustness than excitation methods in clean and noisy environments.

To further prove the generalizability of the EKATP, we generate a large-scale system containing T = 15,000 steps under the condition of h = 0.003 and σ = 0.00. After that, we randomly choose three different time periods to train and test the model as follows, the procedure of which is detailed in Supplementary Table S1.9.

First, since the 3-dimensional time series state and 96-dimensional time series state are diffeomorphic (Sauer et al., 1991), which is indicated by the data preprocessing procedure, it implies that the mapping between these two time series is reversible. Here, we map the 96-dimensional gene expression predictive results onto a 3-dimensional space by orthogonal inverse transformation (Anderson et al., 1987) to visualize the predictive result of the EKATP.

Figures 3A,B,C demonstrates that the predictive results of the EKATP are close to the true value for different periods of a time series. Figure 3 shows that the EKATP can accurately predict the gene expression time series at different periods, implying that it has a strong generalizability, even in a very complicated nonlinear environment.

We always use a nonlinear pendulum model (Hirsch, 1974) with oscillatory behaviour to describe a proteomics time series with a low-dimensional manifold (Sauer et al., 1991) by Eq. 19, { d 2 θ d t 2 + g l sin ⁡ θ = 0 θ ( t 0 ) = h . (19) where   l , g and t denote the length, gravity and time, respectively. h denotes the initial value of θ , which represents the level of the complicated nonlinear system. When h is greater, the nonlinear relationship between proteins becomes more complicated.

Since a considerable amount of noise exists in a protein time series, we employ white Gaussian noise (Li et al., 2017) to describe it by Eq. 20, { θ ˜ = θ + ε θ θ ˙ ˜ = θ ˙ + ε θ ˙ , (20) where θ ˜ and θ ˙ ˜ represent data with noise. ε θ and ε θ ˙ represent the noise Gaussian terms for θ and θ ˙ by normal distributions N ( 0 , σ 2 ) with a zero mean and a standard deviation σ .

Here, we describe how to obtain a high-dimensional proteomics time series with a low-dimensional manifold. First, we generate the 2-dimensional time series V = ( V 1 , V 2 , … , V T ) ∈ ℝ 2 , which is listed in Supplementary Tables S2.1, S2.2. Next, we develop a random orthogonal transformation (Anderson et al., 1987) matrix P ∈ ℝ 64 × 2 . Finally, we map the state of a 2-dimensional time series onto the state of a 64-dimensional time series by Eq. 18 to simulate a high-dimensional proteomics time series F = ( F 1 , F 2 , … , F T ) ∈ ℝ 64 with a 2-dimensional manifold, which is listed in Supplementary Tables S2.3, S2.4.

To prove the accuracy and robustness of the EKATP, we generate a system containing T = 1,600 steps and choose the last 1,000 steps to visualize the predictability for the EKATP.

Figure 4 shows that the EKATP can effectively predict a proteomic time series under clean and noisy environments within 1,000 steps, the details of which are listed in Supplementary Tables S2.5, S2.6; Supplementary Presentation S4.

Figures 4A,B shows that the EKATP not only has less of a predictive error under a clean environment ( σ =0.00) than the existing methods but also maintains a smaller predictive error when h increases from 0.8 to 2.4. Moreover, the predictive error of the EKATP increases much slower than that of the existing methods when the predictive step increases.

Figures 4C,D demonstrates that the EKATP has less of a predictive error under a noise environment ( σ =0.03) than the existing methods. When h increases from 0.8 to 2.4, the predictive error of the EKATP remains stable. In particular, with the increase in predictive steps, the predictive error of the EKATP increases much slower than that of the existing methods.

Figures 4E,F indicates that the EKATP not only has less of a predictive error under a noise environment ( σ =0.08) than the existing methods but also has a predictive error of the EKATP that remains stable when h increases from 0.8 to 2.4. In particular, when the predictive steps are long enough (after 500 steps), the predictive error of previous methods increases much faster than that of the EKATP.

Figures 4A,C,E shows that the predictive error of the EKATP remains stable when the noise intensity σ increases from 0 to 0.08 under complexity h = 0.8. Figures 4B,D,F shows that the predictive error of the EKATP remains stable when the noise intensity σ increases from 0 to 0.08 under complexity h = 2.4.

Figure 4 demonstrates that the predictive accuracy and robustness of the EKATP outperforms the existing methods under clean and noisy environments.

Since Figure 4 shows that KAE has a better predictive effect than the other existing methods, we use it to compare the predictive performance with the EKATP by visualizing the predictive trajectory.

Indicated by our data preprocessing procedure, since the 2-dimensional time series state and 64-dimensional time series state are diffeomorphic (Sauer et al., 1991), the mapping between these two time series is reversible. Here, we map the 64-dimensional protein time series predictive results onto a 2-dimensional space by orthogonal inverse transformation (Anderson et al., 1987) to visualize the predictive time series trajectory. Figure 5 shows the predictive trajectories of the KAE and EKATP within 1,000 steps under the initial conditions of h = 2.4 and σ = 0.03, which show that the predictive protein time series trajectory of the EKATP (Figure 5B) is much closer to the true trajectory than that of the KAE (Figure 5A). Figure 5 further indicates that the predictive accuracy and robustness of the EKATP is better than that of the KAE.

To further prove that the EKATP has strong generalizability, we randomly selected 20 pieces of different protein time series data for model training and analysis. The details are listed in Table 1; Supplementary Table S2.7.

After we employ 20 different proteomics time series datasets to test the KAE and EKATP, Table 1 shows the predictive error of the KAE and EKATP at 1,000 steps under different initial noise and complexity ( h ) conditions, which demonstrates that the EKATP has less of a statistically significant minimum, maximum, average and variance of the predictive error than the KAE under each noise and complexity ( h ) condition (p-value <0.05) (Gao et al., 2017; Li et al., 2017; Gao et al., 2021). Table 1 implies that the EKATP has statistically significant predictive power for different time series datasets.

We usually employ a nonlinear biological fluid system (Noack et al., 2003) to describe the high-dimensional metabolic time series with a low-dimensional manifold (Sauer et al., 1991) for the flow behavior of biological fluids simulation by Eq. 21, { x ˙ = γ x − ω y + A x z y ˙ = ω x + γ y + A y z z ˙ = − λ ( z − x 2 − y 2 ) , (21) where γ , ω and A determine the size of the fluid. λ determines the speed of the dynamics of z . The different initial values of x , y and z determine the different nonlinear complexities of the metabolomics time series. We use the initial conditions ζ 1 ( x =0, y = -0.01,   z = 0) and ζ 2 ( x =0.01, y = -0.1,   z = 0.5) to generate a high-dimensional metabolomics time series with low complexity h 1 and high complexity h 2 , respectively.

Since the metabolomics time series contains considerable noise, we employ white Gaussian noise (Li et al., 2017) to describe it by Eq. 22, { x ˜ = x + ε x y ˜ = y + ε y z ˜ = z + ε z . (22) where x ˜ , y ˜ and z ˜ represent data with noise. ε x , ε y and ε z represent the white Gaussian noise for x , y and z by normal distributions N ( 0 , σ 2 ) with a zero mean and a standard deviation σ .

Fig. 6 Comparison of the RNN, LSTM, DAE, KAE and EKATP. The abscissa represents the time step, and the ordinate represents the predictive error. (A) The initial conditions are h 1 and σ = 0.000. (B) The initial conditions are h 1 and σ = 0.001. (C) The initial conditions are h 2 and σ = 0.000. (D) The initial conditions are h 2 and σ = 0.001.

Here, we show how to generate a high-dimensional metabolomics time series with a low-dimensional manifold. First, we build up the 3-dimensional time series V = ( V 1 , V 2 , … , V T ) ∈ ℝ 3 , which is listed in Supplementary Tables S3.1; S3.2. Next, we develop a random orthogonal transformation (Anderson et al., 1987) matrix P ∈ ℝ 96 × 3 . Finally, we map the state of the 3-dimensional time series onto the state of the 96-dimensional time series by Eq. 18 to simulate a high-dimensional metabolic time series F = ( F 1 , F 2 , … , F T ) ∈ ℝ 96 with the 3-dimensional manifold, which is listed in Supplementary Tables S3.3, S3.4.

To demonstrate the accuracy and robustness of the EKATP, we generate a system containing T = 900 steps and choose the last 100 steps to visualize the predictive result of the EKATP. Figure 6 shows the predictive results of the metabolic time series under different initial conditions and environments for the last 100 steps. Detailed information is listed in Supplementary Tables S3.5, S3.6; Supplementary Presentation S4.

Figures 6A,C demonstrates that the EKATP has less of a predictive error under a clean environment ( σ =0.000) than the existing methods. When the complexity of h increases, the predictive error of the EKATP remains stable. With the increase in the predictive step, the predictive error of the existing methods increases rapidly, while the predictive error of the EKATP remains small.

Figures 6B,D suggests that the EKATP not only has less of a predictive error under a low noise intensity environment ( σ =0.001) than the existing methods but also has a predictive error of the EKATP that remains stable when h increases. In particular, when the predictive steps are long enough, the predictive error of the EKATP increases much slower than that of the existing methods.

Figure 6 implies that the EKATP has better predictive accuracy and robustness than the existing methods under clean and weakly noisy environments.

Since a metabolomics time series usually has strong noise intensity (Mak et al., 2015), we use the EKATP to predict a high-dimensional metabolomics time series under strong noise intensities to prove its robustness. Because the 3-dimensional time series state and the 96-dimensional time series state are diffeomorphic (Sauer et al., 1991), the mapping between these two time series is reversible. Thus, after we map the 96-dimensional metabolic time series predictive results onto a 3-dimensional space by orthogonal inverse transformation (Anderson et al., 1987), Figure 7 shows the predictive time series trajectories by the EKATP under different intensities of noise. We select the last 100 steps to validate the predictive power of the EKATP as in the previous setup (Supplementary Table S3.6). The results demonstrate that although the true data become gradually messy when we increase the noise intensity σ , the predictive time series trajectory of the EKATP is still very close to the true data to a certain extent (Figures 7A,B,C,D,E,F), which implies that the EKATP still has a satisfactory predictive performance when we increase the noise intensity.

Moreover, we use Eqs 23, 24 to calculate the Pearson correlation coefficient (PCC) (Abar et al., 2017) and the root mean squared error (RMSE) (Abar et al., 2017) between predictive and true data under different noise intensities.

Here, V t and V ^ t represent the true and predictive data at time t . μ and μ ^ represent the average value for true and predictive data, respectively. p represents the predictive step size. PCC = ∑ t = 1 t = p ( V ^ t − μ ^ ) ( V t − μ ) ∑ t = 1 t = p ( V ^ t − μ ^ ) 2 ∑ t = 1 t = p ( V t − μ ) 2 . (23) RMSE = 1 p ∑ t = 1 t = p | | V ^ t − V t | | 2 . (24)

Figure 8A shows that the PCC value of the EKATP gradually decreases when we increase the noise intensity σ , but the overall value is relatively high. Figure 8B indicates that with the increase in noise intensity σ , although the RMSE value of the EKATP gradually increases, it is still relatively small. Thus, we conclude that the EKATP can effectively avoid noise interference and is robust enough under a very strong noise intensity condition.