Biomarker models

During ad libitum energy intake in the 9H and 5H conditions, participants consumed more energy than needed to maintain body weight [17], consistent with findings from other studies of insufficient sleep with ad libitum energy intake [13, 48–50]. To increase the specificity of our biomarker models for insufficient sleep, the compounds most impacted by this increased energy intake were excluded from all biomarker models. To identify the compounds most impacted by increased energy intake we exclusively analyzed samples in 9H versus BL, within condition order A, to remove the potential confounding effects of increased energy intake combined with insufficient sleep that occurred in 5H. Within the 9H versus BL comparison, compounds with an area under the receiver operator characteristic plot (AUROC) ≥ 0.70, detected using the biomarker module in Metaboanalyst [51], were excluded from all biomarker models—95 compounds met this criterion. Of these 95 compounds, three (propionylcarnitine, phosphatidylcholine(PC)-[36:3], and sphingomyelin(SM)-[d34:1]) were previously found to be impacted by insufficient sleep or total sleep deprivation [30, 33, 39].

Because findings from untargeted metabolomics studies show ~15%–20% of the human plasma metabolome is circadian regulated [32, 34, 35], and nearly half of the metabolites in a targeted analysis showed 24-hour time-of-day patterns [29], we trained biomarker models using two planned datasets, one that included and one that a priori excluded compounds with 24-hour time-of-day patterns. We also performed an exploratory analysis using only the compounds with 24-hour time-of-day patterns (Supplementary Table S1). MetaCycle [52] (v1.1.0) for R [53] (v3.5.3) was used to identify compounds with 24-hour time-of-day patterns in at least one condition. LOG2 normalized metabolomics data from the BL, 9H, and 5H conditions, using the T1, T5, T9, T13, T17, and T21 samples, were loaded separately into MetaCycle and the “meta2d” function was used with “cycMethod” set to “JTK” and “LS.” The resulting p-values for each compound for each condition were combined using the minP method [54]. Statistical significance for detecting a 24 hours time-of-day pattern was set at the 10% False Discovery Rate corrected level, similar to our previous proteomics findings [28]. Based on these analyses, we identified 985 compounds with 24-hour time-of-day patterns in at least one condition.

For all biomarker models, plasma samples were divided into training and validation datasets. After filtering and QC analyses, 153 plasma samples were included in the final biomarker model analyses. Within the BL, 9H, and 5H conditions, the plasma samples were randomly divided into training (2/3 samples) and validation (1/3 samples) datasets (Supplementary Figures S2 and S3), balanced by time-points, using the caret [55] package (v6.0) for R (v3.5.3). We assigned 2/3 samples to the training dataset, as opposed to dividing the samples evenly between the training and validation datasets, to maintain higher statistical power in the training dataset during the discovery phase of our biomarker analyses. The training dataset consisted of 104 samples with 57 from men and 47 from women. The validation dataset consisted of 49 plasma samples with 20 from men and 29 from women. Because participants had 9 hours per night sleep opportunities in the BL and 9H conditions, samples from BL and 9H were defined as adequate sleep and 5H samples defined as insufficient sleep.

Compound selection for biomarker models was conducted using least absolute shrinkage and selection operator (LASSO) in the GLMNET [56] package (v2.0–16) for R (v3.5.3), and included annotated and unannotated compounds. Biomarker models were trained individually on datasets that included and excluded the 985 compounds with 24-hour time-of-day patterns using the “cv.glmnet” function with “nfolds” set to the number of samples for leave one sample out cross-validation (LOOCV), “type.measure” set to “mae,” “family” set to “binomial” for logistic regression, “alpha” set to “1.0” for LASSO, and “nlambda” set to “500.” Lambda requires tuning; thus, the mean value of lambda across 500 searches that provided the minimum cross-validated error within the training set was defined as the optimal lambda value for each final model. The predicted class (e.g. adequate sleep) in the validation datasets was determined using the “predict” function with “type” set to “class” and “s” set to the optimal lambda, additionally, the logistic regression linear predictor was determined using the “predict” function with “type” set to “link” and “s” set to the optimal lambda.

We also developed linear support vector machine (SVM) biomarker models using the e1071 package (v1.7-0.1) [27] for R (v3.5.3). SVM models were developed using subsets of the compounds identified by LASSO. Specifically, compound subsets for training the SVM models were created using the “LASSO frequency” function within the Biomarker module in MetaboAnalyst [51] that provides the frequency of LASSO compound selection using Monte Carlo cross-validation. SVM models were trained on the resulting subsets consisting of compounds with ≥100%, 90%, and 80% LASSO selection frequencies. The optimal compound subset for each SVM model was identified based on the compound subset that produced the highest Matthew’s Correlation Coefficient [57, 58] (MCC) in the validation dataset using the insufficient and adequate sleep classifications. MCC is invariant to the balancing of samples between classes [57] and ranges from −1 to 1, with −1 representing perfect negative correlation (all samples misclassified), and 1 representing perfect positive correlation (all samples correctly classified). Thus, MCC was a priori selected as the primary performance metric of each biomarker model. Biomarker models were trained individually on datasets that included and excluded compounds with 24-hour time-of-day patterns using the “svm” function with “cross” set to the number of samples in the training set for LOOCV, “kernel” set to “linear,” “scale” set to “FALSE,” and “type” set to “C-classification.” Additionally, the argument “class.weights” was set to the inverse of the number of samples in the insufficient and adequate sleep classifications in the training dataset to account for the larger number of adequate sleep samples. The “cost” function for SVM requires tuning. The “tune.svm” function was used to identify the optimal “cost” value for each SVM model with “scale” set to “FALSE,” “cost” set to the range 2^–10 to 2², “cross” set to “10” for 10-fold cross-validation, “class.weights” set to the inverse of the number of samples for insufficient and adequate sleep classifications in the training dataset as described above, and “kernel” set to “linear.” Ten-fold cross-validation was used to tune the “cost” parameter, and the tuning procedure was repeated 25 times for each SVM model separately. The mode of the optimal “cost” value from the 25 tuning repetitions was selected as the final “cost” value for each SVM model. The predicted class and SVM decision values for the validation datasets were determined using the “predict” function with “decision.values” set to “TRUE” and “scale” set to “FALSE.”