All procedures of experiments were reviewed and approved by the Ethics Review Committee for Animal Experimentation, China Academy of Chinese Medical Sciences, and it was implemented according to the national laws of China as well as local guidelines. A total of 60 Sprague–Dawley male rats (weighing 180–220 g) were used in the study. All rats were randomized into six groups: sham, vehicle, DHI dose-1 (DH1, 1 ml/kg, iv), DHI dose-2 (DH2, 2.5 ml/kg, iv), DHI dose-3 (DH3, 5 ml/kg, iv), and DHI dose-4 (DH4, 10 ml/kg, iv), with 10 animals for each group.

The rat transient middle cerebral artery occlusion (MCAO) model was induced as described elsewhere (Guo et al., 2015; Zhu et al., 2015; Wenjuan et al., 2018). Briefly, rats were anesthetized with 10% chloral hydrate (4 ml/kg, i.p.), and a monofilament with a silicon-coated tip was inserted into the internal carotid artery and advanced until it obstructed the middle cerebral artery. The blood circulation was blocked for 1 h and then restored by removing the monofilament. Sham-operated rats received the same manner as the MCAO rats without ischemia, the external carotid artery was surgically prepared for the insertion of the filament, but the monofilament was not inserted. All animal experiments were conducted in accordance with the ethical principles of animal use and care.

Then, the drug was injected through the caudal vein with 1 ml/kg euvolemia, which was diluted by physiological saline of varying amounts according to the concentrations in the different groups once daily for three continuous days. The sham and vehicle groups were injected with the same doses of saline. DH injections were supplied by ShanXi BuChang Pharmaceutical Corp., Ltd.

Blood samples were taken from the abdominal aorta on the third day after anesthetization, and the brain was removed rapidly.

To evaluate the protective effect of DHI against cerebral ischemia, the cerebral infarct volume, brain weight index, and neurological deficit score were measured.

Triphenyltetrazolium chloride (TTC) staining was used to determine infarction volume. Two-millimeter-thick slices of the cerebrum were immersed in a 2% TTC solution at 37°C for 30 min. Infarct area was measured using the Scion Image program, and the areas of ischemic lesions were calculated using [contralateral hemisphere volume −(ipsilateral hemisphere volume − measured injury volume)].

After 36 h of reperfusion, the brains were swiftly removed. The brain samples and bodies of the rats were weighed immediately. The brain weight index was expressed as brains mass (g)/total rat mass (g).

Neurological deficit scoring (NDS) was measured at first, second, and third after reperfusion by an investigator who was not involved in the study. Consciousness, breathing, smell, vision and hearing, reflexes, motor function, overall activity, orientation, and presence of seizures were scored, according to the neurological deficit score system [0–4 scale; 0, no motor deficit (normal); 1, forelimb weakness and torso turning to the ipsilateral side when held by tail (mild); 2, circling to the contralateral side but normal posture at rest (moderate); 3, unable to bear weight on the affected side at rest (severe); and 4, no spontaneous locomotor activity or barrel rolling (critical)] adapted from a previous study.

Serum samples were collected for targeted metabolomics analysis based on an Agilent 1290 rapid-resolution liquid chromatography system coupled with an Agilent 6550 QTOF/MS system, and the metabolite extraction process was performed as described in previous research with little modification (Dettmer et al., 2011). The targeted metabolomics analysis was used to determine the key nodes in the cerebral ischemia modules adapted in this work, which were described in a previous study (Wenjuan et al., 2018).

Metabolite expression was quantitively identified by the HPLC system for each group. Based on the metabolite expression data, the metabolic networks were constructed using the WGCNA R package (Langfelder and Horvath, 2008). A matrix of pairwise correlations was constructed between all pairs of metabolites across the measured samples using appropriate soft thresholding for each group. The thresholds were selected when the network achieved the best scale-free topology criterion (Langfelder et al., 2008).

After construction of the metabolic networks for the six groups, we characterized them in terms of network entropy, which can be decomposed into contributions from individual nodes (Tan and Wu, 2004; Li et al., 2009):

where I_i is the importance of node i, N is the number of nodes in the network, and k_i is the degree of node i. The network structure entropy is defined as follows:

A method that integrates all metabolite data was proposed to explore the equilibrium state between Yin (downregulated expression) and Yang (upregulated expression). All of the metabolites among the different groups were compared. First, we compared the vehicle group with the sham group. Then, we compared the DHI groups with the vehicle group. The six metabolites’ expression transformation patterns were defined: The metabolites with expression levels that changed from negative to positive were defined as positive transformation (PT). Those with expression levels that changed from positive to negative, on the contrary, were defined as negative transformation (NT). The metabolites with expression that changed from a low positive value to a higher positive value were defined as positive add (PA). Those metabolites with expression values that remained negative yet demonstrated an increasing trend were defined as negative add (NA). Those with high negative values that decreased to lower negative values were named as negative reduce (NR). Metabolites with expression levels that remained positive despite a decreasing trend were defined as positive reduce (PR). Among these, PT, PA, and NA were regarded as yang because of their upregulation character. By comparison, NT, PR, and NR were regarded as yin because of their downregulation character.

All of the metabolites were classified into the above six patterns. The cumulative summation of each class of metabolites was denoted as PT^a, PA^a, NA^a, NT^a, PR^a, and NR^a, respectively. Subsequently, the Yin/Yang ratio and difference were observed. Because PT and NT reflect qualitative changes, a specific comparison between PT^a and NT^a was performed. The formulas are as follows:

In this section, we used Kendall’s coefficient of concordance, or Kendall’s Tau, to evaluate the association between any two different metabolites. Recall that the definition of Kendall’s Tau is as follows:

In the above formula, P₁ and P₂ are two independent points distributed as (X,Y) and the probability P_r(⋅) is defined on ℱ with ℱ being the distribution function of bivariate random vectors (X,Y). The two points P₁ and P₂ on the plane are said to be concordant if joining them produces a positive slope; they are discordant if the slope is negative. Thus, τ=0, implying that X and Y are independent, and τ±1 means that Y = f(X) for some monotone function f. The sample estimator of τ is

Here,

and sgn(u) is the sign function.

Thus, the Kendall’s Tau estimator t_n reflects the proportion of concordant pairs of points in a sample P₁,,P_n minus the proportion of pairs of points that are discordant. We estimated Kendall’s Tau between each two distinct set of metabolites.

To quantitatively describe the synchronization of pairwise nodes, we determined an absolute value greater than 0.5 as a threshold value. We selected pairwise nodes that were in accordance with the threshold for network construction, and an open-source software, Cytoscape 3.6.1, was used to visualize the networks. To identify the synchronous pairwise nodes with a higher degree of synchronization, we defined synchronous pairwise nodes and asynchronous pairwise nodes as those with correlation coefficients greater than 0.7 or less than −0.7, respectively. Functional modules were identified using the Molecular Complex Detection (MCODE) method (Bader and Hogue, 2003).

Kendall’s Tau was estimated between each two distinct metabolites, which were reconstructed from the different groups. If the module had been full-course synchronous during DH1, DH2, DH3, and DH4, it was called DH1–DH4; otherwise, it was called a partial synchronous module, such as DH1/DH2, DH3/DH4, or DH2/DH3. In this process, there existed both corotating changes and reverse changes, which were defined as synchronization and anti-synchronization, respectively. The synchronous module and anti-synchronous module represent the synergistic and antagonistic relationships between the metabolites, respectively.

The pathway enrichment analysis of synchronous modules was employed on the MetaboAnalyst¹ to evaluate the most relevant biological roles involved in the conditions. This analysis was conducted using a modified Fisher’s exact test. We selected all pathways that were significant with a p-value < 0.05 after correcting for multiple-term testing.

Compared with the vehicle group, DHI significantly reduced the infarction volume, increased brain weight index, and improved the neurological deficit score with a dose-dependent relationship. Overall, efficacy increased with increasing dose (Figures 1A–C; Guo et al., 2015; Zhu et al., 2015; Wenjuan et al., 2018). In terms of the brain weight index, the slope of the linear fit between the DHI multi-groups was 0.89. Compared with vehicle, DH3 and DH4 showed a significant difference, but no significant difference was found between vehicle, DH1, and DH2. Between the multiple DHI groups, DH3 and DH4 were significantly different from DH1 to DH2, but no significant difference was found between DH3 and DH4 or between DH1 and DH2. Therefore, we divided the four DHI dose groups into a low-dose category (DH1 and DH2) and a high-dose category (DH3 and DH4). The high-dose category was significantly more effective than the low-dose group.

Based on the fatty acid/amino acid profiling and fusing maps of comprehensive serum metabolome through our previous studies, 41 metabolites were selected as candidate biomarkers related to efficacy of DHI (Wenjuan et al., 2018). Detailed information on these compounds is shown in Supplementary Table 1.

As illustrated in Figure 1D, compared with the sham group, the network entropy value decreased from 3.59 to 3.51 for the vehicle group. This result was reversed when DHI were given. With the DHI concentration increasing, the entropy values also gradually increased, up until it reached the normal level, where it remained relatively constant. Specifically, DH1 entropy value was 3.52, which was slightly higher than the vehicle group and increased by 12.5%. With an increasing dose, the entropy value of the metabolic network gradually increased. In the DH2 group, this value gradually returned to the normal level, which was increased by 94%. The DH3 and DH4 groups tended to maintain a relatively stable, high level. In general, this result is consistent with the trend of the efficacy measured by animal experiments. It also indicated that the 41 metabolites that we selected can represent efficacy phenotypes at a molecular level.

Based on the expression level, 41 metabolite biomarkers in each group were classified into the six patterns of transformation (Figure 2A), as described in section “The Yin and Yang Metabolite Expressional Transformation Patterns” (Figure 2B). Generally, Yin dominates in the cerebral infarction state, which accounts for approximately 66.7% of the total metabolites. This was entirely reversed in the DHI groups which manifested by changes in these metabolites’ expression from Yin to Yang. The same result was found in the Yang aspect. In particular, seven metabolites constantly maintained the PT or NT reverse trend, including hydroxypyruvate, stearic acid, L-carnitine, L-arginine, xanthine, tetrahydrofolate, and L-leucine. Some metabolites, such as myristic acid, gradually reversed from NA to PT. With the dose increasing, some dose-insensitive metabolites only changed into PT in DH4, such as arachidonate and Cys-Gly. Two metabolites, i.e., cystathionine and anandamide, maintained Yin status and did not reverse. In brief, PT and NT transformation power are dominant whether in the pathologic or DHI treatment process. According to the dose–effect interaction, the metabolites’ transformation manifested as PT gradually increased and NT gradually decreased (Figure 2C), i.e., the Yin-to-Yang conversion.

To describe the specific trends of six transforming force patterns, we calculated the cumulative sum of PT, NT, NA, PA, PR, and NR in the different DHI dose groups, respectively. The details are shown in Figure 3. From a quantitative perspective, transformation power represented by PT and NT was still dominant. Compared with the vehicle group, the number of NT metabolites constantly decreased as the DHI dose increased, and the correlation coefficient was 0.76 (Figure 3B). Moreover, the number of PT metabolites was maintained at a high proportion in multiple dose groups (Figure 3A). Although their inclusive metabolites were different, there was no overlap of PT metabolites among the multiple DHI dose groups and the vehicle group (Figure 2D). From the vehicle group to the DHI treatment groups, the number of NA and PR metabolites showed a dose-dependent growth associated with increases in DHI concentration (Figures 3C,D), and the correlation coefficients were 0.82 and 0.21, respectively. There was no correlation among the vehicle group and multiple DHI dose groups, because of the relatively few numbers of PA metabolites (Figure 3E). NR in the DH1, DH2, and DH3 groups were similar to the vehicle group, but in the DH4 group, it was twice as many as the vehicle group (Figure 3F). Generally, this result shows that inhibiting the downregulation power of Yin was mainly contributing to the efficacy of DHI treatment.

As described in section “The Yin and Yang Metabolite Expressional Transformation Patterns,” we calculated the ratio and difference of the Yang/Yin cumulative sum. First, PT, NT, NA, PA, PR, and NR were divided into two Yin/Yang categories, and their differences were compared. The results show that the difference between Yin and Yang gradually increased, and the equilibrium state was broken when DHI was introduced (Figures 4A–C). Second, since PT and NT play a dominant role, we compared these two powers in isolation, and the same trends were obtained. Compared with the PT/NT conversation, the overall Yin/Yang comparison had higher linear fit scores on the ratio and difference, i.e., 0.85 and 0.91, respectively (Figures 4D–F).

In terms of synchronization node pairs, we found that a metabolite pair—oleic acid and hexadecanoic acid—always exists in dose change, and only the negative control node pair of octadecatrienoic acid and L-threonine appears in the conversion dose group.

In the synchronous network, we identified three modules of the DH1–DH4 full-course synchronous network and two, three, and two modules of DH1/DH2, DH3/DH4, and DH2/DH3 synchronous networks, respectively, using the MCODE algorithm (Figure 5 and Supplementary Figures 1, 3, 6). Based on pathway enrichment analysis, the three full-course synchronized modules were enriched to specific signal pathways of linoleic acid metabolism, fatty acid biosynthesis, and phenylalanine, tyrosine, and tryptophan biosynthesis, and so on (Supplementary Figure 5).

With different DHI treatment doses, similar and individual synchronous modules may reveal molecular mechanisms. By comparing the synchronous modules in the DH1/DH2 and DH3/DH4 networks, a conservative allosteric module was obtained, which contained three constant nodes: taurine, L-tyrosine, and L-leucine (Figure 6A). It may represent the basic mechanisms of DHI in the treatment of cerebral infarction. Meanwhile, there are three unique synchronous modules obtained. One was composed of hexadecanoic acid, oleic acid, and linoleate in the DH1/DH2 synchronous network (Supplementary Figure 3). The other two modules in the DH3/DH4 synchronous network were composed of octadecatrienoic acid, phytosphingosine, myristic acid, and glyoxylate, L-valine, and L-glutamate (Supplementary Figure 6). More importantly, a new emerging module of the DH3/DH4 synchronous network was found, which was composed of glyoxylate, L-glutamate, and L-valine. The KEGG enrichment analysis of this specific module indicated that it was mainly involved in the glutamine and glutamate metabolism pathway (Supplementary Figures 4, 7). As for the DH2/DH3 synchronous network, which represents the low to high DHI dose transformation, two synchronous modules were obtained, and they were unique to DH2/DH3 (Figure 6A). These modules mainly participated in the regulation of glycine, serine, threonine, glyoxylate, and dicarboxylate metabolism (Supplementary Figure 5).