The 64-channel EEG data were recorded from subjects with the two-tone auditory oddball paradigm (see Methods). Mathematically, we suppose that the dynamics of brain signals under the standard stimulus evolves as follows: x ( t + 1 ) =   f ( x ( t ) ,   λ ) , (1) where x = ( x 1 , x 2 , … , x n ) T ∈ R n is a vector containing signals from all n electrodes, f is assumed to satisfy the existence and uniqueness condition, and λ is a bifurcation parameter. Because we assume the invariance of the brain dynamics, the signals y for the deviant stimulus (the deviant tone in two-tone auditory oddball paradigm, see Methods), which is a small perturbation from x , follow the same equation y ( t + 1 ) = f ( y ( t ) , λ ) . Thus, the evolution of the difference process z ( t )   =   y ( t ) − x ( t ) can be approximated as follows: z ( t + 1 ) ≜ y ( t + 1 ) − x ( t + 1 ) ≈ A ( λ ) ( y ( t ) − x ( t ) ) +     ξ ( t ) = A ( λ ) z ( t ) + ξ ( t ) , (2) where the matrix A ( λ ) ∈ R n × n is the Jacobian matrix of f ( ⋅ ,   λ ) , and ξ ( ⋅ ) is inevitable noise in the brain and assumed as a small additive noise term which is independent of A ( λ ) . When the parameter λ → λ 0 , where λ 0 is a codimension-one catastrophic bifurcation point of the system (Kéfi et al., 2013), the largest-modulus eigenvalue of A ( λ ) will influence system stability. According to the DNM theory (Liu et al., 2012), the network between nodes can be partitioned into a DNM group and a non-DNM group. The DNM group leads the criticality of the system, and deviations of signals and correlations between any nodes in the DNM group will increase sharply around a critical point. Utilizing these characteristics, we constructed the DNMnet and compared the dynamical differences between brains of HCs and those of patients with PD. The DNMnets under duration deviant (dD) experiments and frequency deviant (fD) experiments are shown in Figures 1A–D. Details can be found in the Methods.

Mismatch negativity (MMN) (Erickson et al., 2016) is one of the most well-known patterns in PD, whose amplitude of neural activities is reduced compared to the healthy brain. MMN is associated with cognitive impairments (Chung et al., 2017) and can serve as a biomarker for early interventions (Bodatsch et al., 2015) and development of novel treatments (Perez et al., 2017; Kantrowitz et al., 2018). However, mechanisms underlying altered connectivity remain unknown in the clinical field. Using the DNM theory, we found that the MMN amplitude is explained by the difference process z ( t ) of the brain signal. Furthermore, the nodes in the DNMnet with a large deviation in the signal are all located in the fronto-central area, which corresponds to the region with significant MMN patterns (Figures 1A–D). The DNM theory can provide a dynamical explanation of the MMN pattern, and in return, the MMN pattern provides an evidence for the criticality of the normal brain under the DNM framework.

A DNMIndex was designed to measure the criticality based on the characteristics extracted from the DNMnet as follows: DNMIndex = Π ( σ ( x ) ) ⋅ Π ( | ρ | ) , (3) where “Π ( ⋅ ) ” stands for the geometric mean of a vector, “ σ ( ⋅ ) ” is the component-wise standard deviation, “ x ” denotes the brain signals after preprocessing (corrected epochs) of nodes in the DNMnet, “ | ⋅ | ” is the component-wise absolute value function, and “ ρ ” denotes the correlations between nodes with edges in the DNMnet. The DNMIndex will be higher at the critical state than at non-critical states. Using the DNMIndex, we found the brain of HC is in a significantly more critical state than those of patients with PD and UHR individuals in both dD and fD cases (Figures 1E–F).

Mutual information between each pairs of nodes over time and its entropy are important metrics for isotropy and the criticality of a network. Based on the nodes in the DNMnet, we define the mutual information over time from node x to node y as follows: I x → y = ∬ p ( x t ,   y t + Δ t ) log p ( x t , y t + Δ t ) p ( x t ) p ( y t + Δ t ) d x d y (4) where { x t | t   =   1 ,   2 ,   … } and { y t | t =   1 ,   2 ,   … } are the observed time series from node x and y ; Δ t is a time interval; p ( x t ,   y t + Δ t ) is the joint probability density function; p ( x t ) and p ( y t + Δ t ) are the respective marginal density functions; and the integral is over the entire ( x t ,   y t + Δ t ) space. According to the DNM theory, if x and y belong to the DNM group, their correlation will become strong near the critical point; thus, I x → y grows rapidly. However, if either x or y belongs to the non-DNM group, mutual information will remain stable and bounded. Thus, near the critical point, some rapidly growing mutual information will change their distribution over the entire DNMnet. We will use two indices to detect the anisotropy of mutual information in the network—the distribution entropy (DE) and network entropy (NE)—whose definitions are shown in the Methods. With the time interval Δ t = 20 ms, we calculated the DE and NE of the HC, UHR, and PD groups (Figures 2A–B for the dD case, and Supplementary Figures S7A–B for the fD case). The HC group always exhibited significantly higher DE and NE than the PD and UHR groups, while the difference between the PD and UHR groups was not significant. These results support the critical brain hypothesis that the HC group should be at a critical state, whose entropy is high possibly to ensure adaptivity. In addition, correlations between positive symptoms (clinical score) of PD and DE/NE were calculated (Supplementary Figure S11), which also implied that the worse the positive symptoms are, the more low-entropy and stable the brain is. In Supplementary Figure S12, we also showed that the orders of DE/NE for three groups were not sensitively influenced by the parameter Δ t .

In the biological sense, the high entropy at the critical state reflects a possibility that the brain can send information to different brain regions in time to deal with various stimuli.

The causality pattern is another evidence for criticality, which can be detected by the conditional transfer entropy (CTE) between different nodes, which is given as CT E x → y = CMI ( x t ,   y t + Δ t | x ¯ t ) = ∫ ∫ ∫ p ( x t ,   y t + Δ t , x ¯ t ) log p ( x t , y t + Δ t | x ¯ t ) p ( x t | x ¯ t ) p ( y t + Δ t | x ¯ t ) d x d y d x ¯ t , (5) where CMI denotes the condition mutual information, x and y are two variables to detect causality, x ¯ includes all other variables except x, and Δ t is the time interval. CTE is different from mutual information because it excludes indirect influences and is usually used to detect the direct causality. If the brain is at a stable state, causalities between different electrodes should be apparently present, while for critical states, random-like activities decrease the possibility of the appearance of regular causality patterns. With the time interval Δ t = 20 ms, we can obtain heatmaps of the mean CTE between every pair of 64 electrodes for the PD, UHR, and HC groups (Figures 2C–E for the dD case, and Supplementary Figures S7C–E for the fD case). The larger the value of CTE is, the stronger is the direct causality between variables. We found that the PD had a stronger causality network compared with the UHR/HC, which implies that psychosis makes the brain fall into a more stable state. In contrast, it is difficult for the normal brain to form general causality patterns at a critical state, as shown by the heatmap of the HC group. The probability density functions of CTE for each group are plotted in Figure 2F for a better intuition. In Supplementary Figure S12, we also showed that the orders of CTE for three groups were not sensitively influenced by the parameter Δ t .

For EEG data, we performed analysis of covariance (ANCOVA) with age, premorbid IQ, and antipsychotic dose as covariates because these variables were significantly different among groups (Supplementary Table S1). For the duration deviant experiments, all these covariates showed no significant effects (F_1,91 < 2.58, p > 0.11). Therefore, these factors did not affect our findings in the duration deviant experiments. For the frequency deviant experiments, premorbid IQ showed significant effects, whereas age and antipsychotic dose showed no significant effects (F_1,91 < 0.75, p > 0.39). Therefore, we performed ANCOVA with premorbid IQ as a covariate in the frequency deviant experiments.

ANCOVA of the DNMIndex in the frequency deviant condition revealed significant effects of groups (F_2,95 = 6.61, p = 0.01) and premorbid IQ (F_1,95 = 11.08, p < 0.001). Post-hoc analyses revealed that the HC group had a significantly higher mean value than the PD (p < 0.001) and the UHR (p = 0.004) groups, while the difference between the PD and UHR groups was not significant (p = 1.00).

ANCOVA of DE in the frequency deviant condition revealed significant effects of groups (F_2,95 = 5.72, p = 0.005) and premorbid IQ (F_1,95 = 6.80, p = 0.01). Post-hoc analyses revealed that the HC group had a significantly higher mean value than the PD group (p = 0.003), while the difference between the HC and UHR groups (p = 0.33) and the difference between the UHR and PD groups (p = 0.48) were not significant. ANCOVA of NE in the frequency deviant condition revealed significant effects of groups (F_2,95 = 5.44, p = 0.006) and premorbid IQ (F_1,95 = 8.43, p = 0.005). Post-hoc analyses revealed that the HC group had a significantly higher mean value than the PD group (p = 0.004), while the difference between the HC and UHR groups (p = 0.58) and the difference between the UHR and PD groups (p = 0.31) were not significant.

These findings revealed the lower DNMIndex, lower distribution entropy, and lower network entropy of the frequency deviant experiments in the PD group even after controlling premorbid IQ. In UHR, the DNMIndex was lower, whereas DE and NE showed insignificant difference for the frequency deviant experiments after controlling premorbid IQ.

The level of the criticality of the brain contained in the DNMnet can provide a quantitative measure to weigh the risk of mental disease. From the EEG samples, we applied an ensemble classifier based on the DNMnet (see Methods). AUCs and cross validations were used to test the training accuracy and to ensure the reliability of the classifier (Supplementary Figure S8). With the final risk output, we found that the PD group exhibited the highest risk, and the HC group showed the lowest risk, whereas the UHR group evinced a medial risk (Figure 3A for the dD case and Figure 3B for the fD case). This risk analysis indicates that DNM results could help distinguish between the three groups and has a potential to guide prepsychotic diagnosis.

A total of 49 HC subjects, 24 UHR individuals, and 29 patients with PD participated in this study, which contained participants of the Integrative Neuroimaging Studies for Schizophrenia Targeting Early Intervention and Prevention (IN-STEP) (Koike et al., 2013). Detailed participant information is provided in the Supplementary Notes and Supplementary Table S1. The Research Ethics Committee of the Faculty of Medicine, The University of Tokyo, approved this study (approval No. 629, 2226). We conducted this study in accordance with the Declaration of Helsinki. Written informed consent was obtained from all the participants.

A 64-channel Geodesic EEG System (Electrical Geodesics Inc, Eugene, OR) was used to acquire EEG data. Electrodes were referenced to the vertex, and impedances were kept below 50 kΩ. The sampling rate was 500 Hz, and the analog filter bandpass was set at 0.1–100 Hz. The locations of the EEG electrodes are shown in Supplementary Figure S1A.

The two-tone auditory oddball paradigm with 2000 stimuli was performed for each subject when obtaining EEG data. For the duration deviant (dD) experiments, 90% stimuli were standard tones (1,000 Hz, 50 ms) and 10% were deviant tones (1,000 Hz, 100 ms). For the frequency deviant (fD) experiments, 90% stimuli were standard tones (1,000 Hz, 50 ms) and 10% were deviant tones (1,200 Hz, 50 ms). All stimuli were 80 dB SPL with a 1 ms rise/fall time. The stimulus onset asynchrony was 500 ms. The oddball paradigms were counter-balanced, and tones were presented binaurally through earphones while participants watched a silent cartoon.

Original 64-channel signal files (Supplementary Figure S1B) were input into MATLAB (9.3.0) and were further preprocessed using the EEGLAB (v14_1_1b) package (Delorme and Makeig, 2004) (Supplementary Figures S1C–E). Detailed preprocessing procedures are provided in the Supplementary Material. Finally, corrected epochs (deviant epochs with the mean signal of standard epochs being subtracted for each sample) were utilized for further analyses (Supplementary Figure S1F).

Corrected epochs were obtained by subtracting the event-related potential (ERP) waveforms in response to the standard stimuli from those in response to the deviant stimuli. We defined the peak latency as the most negative peak between 100 and 250 ms relative to the onset, and the MMN amplitude was calculated as the mean signal around the peak latency (we used the window of 135–205 ms for the duration MMN and the window of 100–200 ms for the frequency MMN) (Nagai et al., 2013). The MMN patterns for the 64 electrodes in the three groups (PD, UHR, and HC) under dD and fD are shown in Supplementary Figures S3 and S4, respectively. Using the MMN amplitude for each channel as the feature, we can compute p-values of the group difference under the Wilcoxon rank-sum test. After FDR correction, we identified 12 electrodes (Nos. 3, 4, 5, 9, 17, 18, 22, 30, 43, 54, 55, and 58) in the dD and 12 electrodes (Nos. 3, 4, 5, 8, 9, 17, 18, 30, 43, 54, 55, and 58) in the fD that exhibited significant difference (adjusted p < 0.05) between the HC and UHR/PD groups, albeit insignificant difference (adjusted p > 0.05) between the UHR and the PD groups. It should be noticed that all electrodes with significant MMN difference were located in the fronto-central area.

The neuronal signal under standard stimuli is supposed to obey the dynamics in Eq. 1, whereas the corrected epochs after preprocessing follow Eq. 2. Criticality is defined as when the dynamics is around a codimension-1 local bifurcation (Chen et al., 2012). The network of electrodes can be partitioned into a DNM group and a non-DNM group (Chen et al., 2012) by measuring the fluctuations of signals and their correlations between electrodes for each sample. The DNMnet is the leading network for complex systems that drives the system toward or away from critical states. Two conditions for constructing the DNMnet were used in this study:

1. Any member of the nodes in DNMnet is highly fluctuating around a critical point.

The one-sided paired-sample t-test and the one-sided Wilcoxon rank-sum test were used for selected nodes and edges, respectively. We selected nodes and edges with significant statistical difference between HC and PD groups. Significance was set at p < 10 − 3 for nodes and p < 0.005 for edges after FDR correction. Supplementary Tables S2, S3 and Supplementary Figures S5, S6 show the statistical properties of the nodes and edges in DNMnets. Specifically, 20 nodes (the one-sided paired-sample t-test with p < 10 − 3 after FDR correction) were extracted from the dD experiment dataset (Supplementary Figure S5). 12 edges with significantly increasing correlations (the one-sided Wilcoxon rank-sum test with   p < 0.005 after FDR correction) were also detected (Supplementary Figure S6). No edge with significantly decreasing correlation (the one-sided Wilcoxon rank-sum test with significance level α = 0.05 ) in HC was found. We denoted the 20 nodes and 12 edges as the DNMnet for the dD case (Figure 1B). The same analysis was also applied to the fD experiment dataset, from which 17 nodes (the one-sided paired-sample t-test with p < 10 − 3 after FDR correction, Supplementary Table S2) and 37 edges (the one-sided Wilcoxon rank-sum test with p < 0.005 after FDR correction, Supplementary Table S3) were selected to form the DNMnet for the fD case (Figure 1D). No significant decreasing correlation (the one-sided Wilcoxon rank-sum test with significance level α = 0.05 ) in HC was found in the fD experiment dataset, either.

We remark that the Lyapunov Exponent is also a possible measure for detecting criticality. We note that the dynamical network we used here refers to the leading subnetwork (Chen et al., 2012). In the celebrated paper on synchronization (Arenas et al., 2008), the dynamical network is considered as the whole system of the interacting dynamical units. In this manuscript, most of the time we use the notation the DNMnet to avoid the confusion.

When calculating the DNMIndex using Eq. 3, the corrected epochs were used as the signal. Nodes and edges in the DNMnet were considered.

We denote a network as G   =   { V ,   E } , where V   =   { x 1 ,   x 2 ,   … ,   x n } are n nodes and E   =   { e i , j | i ,   j ∈   { 1 ,   2 ,   … ,   n }   and   i ≠   j } are N e   =   n ( n − 1 ) different directed edges. Using the mutual information in Eq. 4, the distribution entropy (DE) is defined as the normalized entropy of all mutual information on edges, as follows: Ent d   =   − 1 log ⁡ N e ∑ i = 1 n ∑ j = 1 j ≠ i n p i j ⁡ log ⁡ p i j , (6) where p i j = I x i → x j ∑ i = 1 n ∑ j = 1 ,   j ≠ i n I x i → x j (7) is the normalized mutual information over the entire network. Mutual information I x i → x j is calculated under Gaussian approximation as I x → y =   − 1 / 2 ⋅ log ( 1 − ρ x t y t + Δ t 2 )   , where ρ xy is the Pearson correlation coefficient between x and y. In contrast, the network entropy (NE) is defined as follows Ent n =   − 1 n ∑ i = 1 n ∑ j = 1 j ≠ i n q i j ⁡ log ⁡ q i j , (8) where q i j = I x i → x j ∑ j = 1 ,   j ≠ i n I x i → x j (9) is the normalized mutual information exiting each node. We use the full network between nodes in DNMnet to calculate the DE and NE. If the brain works around a critical state, some increasing mutual information will expand their distribution density on the network. Thus, DE and NE are larger at a critical state than at non-critical states. High entropy also reflects a possibility that the brain can send information to different brain regions in time to deal with various stimuli.

For the CTE defined in Eq. 5, every pair between 64 electrodes is considered to construct the causality patterns in (Figures 2C–E). CTE is also calculated under Gaussian approximation as CTE x → y =   − 1 / 2 ⋅ log ( 1 − ρ x t y t + Δ t | x ¯ t   2 ) , where ρ xy | z is the partial correlation between x and y conditional on z.

Using the knowledge of the DNMnet after dynamical network analysis, we applied an ensemble classification on the EEG data to estimate the risk index for psychosis. A flowchart of the risk estimation is shown in Supplementary Figure S2.

We assume that there are n k corrected deviant epochs after preprocessing for sample k, which may belong to one of the PD, HC, or UHR groups. PD and HC epochs are used as the training set. Features for training comprise the elements in the covariance matrix for each sample. We assigned each epoch a label: 1 for PD and 0 for HC. Ensemble classification was then applied, and cross validation was performed. Using the classification tree, we can test all PD, HC, and UHR epochs. Finally, the risk for person k is defined as the mean output (the label value) of the epochs belonging to this sample, and this is defined as R k   = 1 n k ∑ i = 1 n k O i , (10) where R k is the risk for sample k, n k is the number of corrected epochs in sample k, and O i is the classification output of the ith epoch. R k is always between 0 and 1. To verify the accuracy and compare different classifiers, we used a risk accuracy index. If there are n HC samples in HC, n UHR samples in UHR, and n PD samples in PD, we sort the risks for all samples in an ascending order from smallest to largest. The risk accuracy for the HC group is set as the percent of HC samples in the first n HC smallest R k s. The risk accuracy for the UHR group is set as the percent of UHR samples in the medial n UHR medium R k s. The risk accuracy for the PD group is set as the percent of PD samples in the last n PD largest R k s. Furthermore, the risk accuracy for the disease group is set as the percent of UHR and PD samples in the last n UHR + n PD as the largest R k . The training accuracy of the classifier is defined as 1 − loss . Thus, the DNM classifier has more than 80% training accuracy (Supplementary Table S4). The risk accuracy indices for the test set were also calculated and are listed in Supplementary Table S4.

We used SPSS (IBM Corp., New York, United States) and MATLAB 9.3.0 (MathWorks Inc.) for statistical analyses. For demographic and clinical data, we performed a chi-squared test, independent t-tests, and analysis of variance (ANOVA) for comparison among groups. Bonferroni correction was performed in post-hoc analyses of ANOVA. For EEG data, we performed paired-sample t-tests and Wilcoxon rank-sum tests to compare difference in the MMN, DNM measures, and risks between different groups. When we quantified an alternative hypothesis of a one-side inequality, such as some index with HC > PD, we used the one-side Wilcoxon rank-sum tests. While we quantified an alternative hypothesis of a two-side inequality, such as some index with UHR≠PD (i.e., UHR < PD or UHR > PD), we used the two-side Wilcoxon rank-sum tests. False discovery rate (FDR) was also controlled for multiple comparisons (Benjamini and Hochberg, 1995). Analysis of covariance (ANCOVA) with age, premorbid IQ, and antipsychotic dose as covariates were also conducted, which showed that they did not affect the difference of measures for criticality between groups.