## 2 Materials and Methods
### Machine Learning Pipeline  
We analyzed the entire set of features using a machine-learning pipeline. We split the available samples of PLR curves into a series of train/test sets via stratified 5-fold cross-validation: in this way, we ensured that in each fold an equal number of samples of each class was considered by the model for both the training and validation.  
We rescaled the set of features using their median values, normalizing according to the 1 st and 3 rd quartiles, i.e., a robust scaling algorithm: in this way we minimized the dependency from possible outliers. The parameters for the features standardization were estimated from the training set and propagated to the related test set, avoiding cross contaminations and biases to the model.  
The processed set of features was used as input for a penalized classification model, trained for the automated prediction of the two classes (ON vs. Healthy). We used a Ridge classifier for this task, performing a further 10-fold cross-validation on the training set for the estimation and tuning of the best hyperparameters.  
The pipeline performances were estimated by monitoring the classification accuracy score (\(ACC\)) and the Matthews Correlation Coefficient (MCC), given by:  
\[ACC=\frac{TP+TN}{TP+TN+FP+FN} \tag{15}\]  
\[MCC=\frac{TP\times TN-FP\times FN}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}} \tag{16}\]  
where _TP_, _TN_, _FP_, and _FN_ represent the number of True Positives, True Negatives, False Positives, and False Negatives, respectively. The use of the \(MCC\) score is recommended for binary classification tasks for the quantification of the correct predictions (accuracy) in relation to the numerosity of the two classes [36].  
We repeated the entire proposed pipeline for 200 different cross-validation procedures, monitoring the score metrics along them. In this way we could test the robustness of the proposed pipeline in relation to different train/test splits, quantifying its efficiency according to the distribution of the obtained scores.  
Feature ranking was assessed by fitting the pipeline on the whole dataset and computing the accuracy drop after a random permutation of each variable. For each variable, 500 different random permutations were simulated to get the average and the corresponding error bars. Since the importance of each variable is contained in the interval [0, 1], to avoid unrealistic confidence intervals outside that range and be more precise, we determined the 95% confidence interval by bootstrapping. For each variable importance, we resampled with replacement 5000 times the original data. Variables with higher importance are then the ones having the most impact on the model prediction and the error bars give information about the statistical significance of each variable, with a _p_-value of 5%.  
To assess the performance of the classification at the patient level, we classified a subject as ON if it had at least 2 curves classified as pathologic. We checked the robustness of this condition by using a different number of curves for the classification of the patients. For instance, we considered as well a subject as a case with at least 4 curves classified as pathologic, which is equivalent to a majority vote. The rationale behind our first choice is that one subject was classified as a pathological case if they have at least two PLR curves classified as pathological because pupil response issues may not arise simultaneously on both eyes and may not be detected at the same time by both protocols.  
The entire pipeline was written in Python v3.10 using the scikit-learn library [37] for the feature computation as well as for the implementation of the pre-processing steps and the machine-learning model.