The dried roots and rhizomes of Rheum tanguticum Maxim. ex Balf. (No. SUCM-20190106) were collected from Gannan Autonomous Prefecture, Gansu Province, China. They were identified as the four-year-old authentic rhubarb by Professor YongGang Yan (Shaanxi University of Chinese Medicine, Xi’an, China), and the voucher specimens were deposited in the Herbarium of Shaanxi University of Chinese Medicine. Quality control of rhubarb can be found in Additional file 1.

Male Sprague Dawley rats were randomly divided into eight groups: normal group, model group, and six doses of rhubarb groups (n = 8). The constipation model was established to observe the effect of rhubarb. The defecation characteristics and pathological sections of colons were employed for purgative effect and toxicity evaluation. Serum samples were collected for metabolomics analysis. The animal experiments were approved by the Animals Ethics Committee of Shaanxi University of Chinese Medicine (No. 2020078). Detailed experimental design and procedures can be found in Additional file 1.

The metabolic profiling of serum samples was performed by an Acquity™ UPLC system (Waters, Milford, USA) coupled to a Synapt™ Q-TOF mass spectrometer. Raw UPLC-Q-TOF/MS data were processed by MassLynxTMv4.1 software (Waters, Milford, MA, USA) for peak detection, noise removal, filtering and alignment to generate a data matrix that was composed of retention time, m/z value and normalized ion intensity for each peak area.

Dose–response modeling mainly consists of four steps: dose design, multiple comparison between groups, degree value calculation of the changed features, and fitting of dose–response curves. The overall workflow for dose–response metabolomics with therapeutic effects and adverse reactions was shown in Fig. 1.

Overall workflow for dose–response metabolomics with therapeutic effects and adverse reactions

For multiple comparisons, one way analysis of variance is performed to compare the relative abundance (RA) of each detected feature. If there are no significant differences (P < 0.05 or P < 0.01, user defined P value) between the administration group and the model group, the features are considered meaningless. If significant differences exist between the model group and the normal group, as well as the administration group and the model group, it can be divided into three situations according to the regulation degree of TCM: “moderate positive reaction”, “excessive positive reaction” and “negative reaction”. If there are no difference between the model group and the normal group, while significant differences exist between the administration group and the normal group/the model group, it means that the modeling did not cause change in metabolites but the administration caused metabolites disturbance, which are assigned as “unexpected reaction” induced by TCM. Multiple comparisons between the groups and the formula are shown in Fig. 2.

Multiple comparison of metabolic features between groups

Then, degree value of the changed features is calculated by applying Eqs. (1)–(6). Degreei in Eq. (1) has a common premise that TCM can adjust the metabolites disorder caused by modeling, and there are three situations: when 0 < Degreei ≤ 1, it indicates that TCM can adjust the level of metabolites to near normal or normal, which is defined as “moderate positive reaction” and is the so-called “efficacy”, in addition, the “moderate positive reaction” value is expressed as the sum of the deviation degree and can be calculated by Eq. (2); when Degreei > 1, it means that TCM overregulates the metabolites to higher than normal levels, which is assigned as “excessive positive reaction”, whose value is expressed as the sum of the deviation degree and can be calculated by Eq. (3); when Degreei < 0, it shows that TCM not only fails to adjust the level of metabolites to the normal level, but also can regulate them to the opposite direction, which is defined as “negative reaction” and is also a type of adverse reactions, whose value is expressed as the sum of the deviation degree and can be calculated by Eq. (4). When Degreei in Eq. (5) is less than 0, it indicates that TCM administration causes the down-regulation of metabolites, at the same time, when Degreei in Eq. (5) is more than 0, it shows the up-regulation of metabolites, both of which mean that TCM can lead to the metabolites to deviate from the normal level. It is defined as “unexpected reaction” and is expressed as the sum of the absolute value of the deviation degree, which can be calculated by Eq. (6).

$$}_i = \frac - RA_ }} - RA_ }}, \, \left( \right)$$

(1)

$$}\,\,\,0 < \frac - RA_ }} - RA_ }} \le 1,\,}\,\left( } \right) = \mathop \sum \limits_^q \frac - RA_ }} - RA_ }},\,\left( \right)$$

(2)

$$}\,\,\frac - RA_ }} - RA_ }} > 1,\,}\,\left( }\,}\,}} \right) = \mathop \sum \limits_^p \frac - RA_ }} - RA_ }}, \, \left( \right)$$

(3)

$$}\,\frac - RA_ }} - RA_ }} < 0,\,}\,\left( }\,}} \right) = \mathop \sum \limits_^0 \frac - RA_ }} - RA_ }},\left( \right)$$

(4)

$$}_i = \frac - RA_ }} }},\left( \right)$$

(5)

$$}\,\left( } \right)\, = \mathop \sum \limits_^m \left| - RA_ }} }}} \right|,\left( \right)$$

(6)

RAAi means the relative abundance of the feature in administration group, and RANi means the relative abundance of the feature in normal group, then RAMi means the relative abundance of the feature in model group.

The final step is to fit the dose–response curves, and dose–response models are regression models where the independent variable (x) refers to the dose and the dependent variable (y) refers to the response (efficacy or adverse reactions), in which y is the value derived from the metabolomics features according to the Eqs. (2), (3), (4) and (6). The positive value of y indicates the efficacy, and the negative value of y represents the adverse reactions. The basic Trendline package in R. Plot was used to construct the dose–response curves, the regression lines and confidence intervals were drawn after, and the models built in the ‘trendline’ function were used to display the regression equations, R-squares and P-values. The following R code was run:

$$> }\left( },},} = }}},}.} = },} = },} = } \right)$$

“exp3P” indicates the formula of y = a×\(^}\) + c; level of confidence interval was 0.95 by default; R2 indicates the R-Squared value of regression model; p indicates the p-value of regression model; AIC or BIC indicates the Akaike’s information criterion or Bayesian information criterion for fitted model, the smaller the AIC or BIC, the better the model; RSS indicates the residual sum of squares of regression model.