fNIRS data collection and analyses
REGRESSION

Data acquisitionData acquisition was performed using previously described methodsIn addition to the DTOFs, the Flow1 system is capable of capturing parallel streams of data including: data from 5 inertial measurement units (IMU) distributed throughout the headset as well as temperature sensor data for both sources and detectors that are acquired at all times. Data from auxiliary sources, if connected, will be collected as well (Fig. For the fifteen participants who completed all four study visits, technical issues occurred during some sessions leading to partial (3 of 30 dosing sessions) or full (1 of 30 dosing sessions) data loss (Supplementary Fig.  Data preprocessing and feature extraction: relative hemoglobin concentrationsWe implemented a channel selection procedure using both a histogram shape criteria as well as the presence of physiological signalsNext, histograms from the selected channels were used to compute the moments of the DTOFs–specifically the sum, mean, and variance momentsThe relative changes in preprocessed moments were then converted to relative changes in absorption coefficients for each wavelength using the sensitivities of the different moments to changes in absorption coefficients (moments Jacobians, i.e. derivative of each moment with respect to a change in absorption coefficient), following the procedure recently described inFinally, we converted the changes in absorption coefficients at each wavelength to changes in oxyhemoglobin and deoxyhemoglobin concentrations (HbO and HbR, respectively) using the extinction coefficients for the two wavelengths and the modified Beer–Lambert law (mBLL)Extraction of physiological signalsPulse rate, pulse rate variability, and breathing rate extracted from TD-fNIRS data: Within each session, we used the sum moment (total counts) of within-module channels (SDS = 10 mm) from the prefrontal plates (prior to performing short channel regression). Data from channels with high scalp-coupling indexPhysiological measurements from auxiliary physiological sensors: We applied a fast Fourier transform (FFT) on the time series of COAbsolute chromophore concentrationsThe DTOF is the convolution between the time-resolved temporal point spread function (TPSF) and the Instrument Response Function (IRF). As such, the mean and variance moments of the TPSF can be obtained by subtracting the IRF moments from the DTOF moments. We leveraged calibration data for Flow1 that characterized the IRF for each laser in the system as a function of both temperature and optical power, which allowed us to estimate the IRF throughout the recording session by mapping the instantaneous temperature and power setting to the calibrated IRF value. Absolute moments were then computed by subtracting estimated IRF moments (i.e. mean and variance) from DTOF moments recorded during the experimental sessions.We then used the analytical expression that converts the mean and variance moments to coefficients of absorption for a semi-infinite homogeneous mediumFractional Amplitude of Low Frequency Fluctuations (fALFF)fALFF quantificationThe time series of data from each channel (for each chromophore) was converted to the frequency domain using an FFT. The ratio of the power in the 0.01–0.08 Hz frequency range was calculated relative to the full 0–0.25 Hz frequency range. For time-varying fALFF analysis, we used a window size of 5 min with a stride of 10 s and fALFF was calculated for each window (Figs.  RMEQ total score predictionTo predict participants’ total RMEQ scores, we used three features of PR, HbO fALFF, and HbR fALFF (all time-varying features normalized by their initial value in a given session). We used a Leave One subject Out (LOO) cross-validation combined with linear regression. On each fold, both saline and ketamine sessions of a given participant were left out as the validation set and the RMEQ scores were regressed against the data from all other subjects; predictions on the validation set are shown in Fig. To test whether the performance of the model was above chance level, we used a shuffling procedure. Here, the RMEQ scores (“targets”) were randomly shuffled across sessions so that the samples for each session were assigned a new RMEQ score from a randomly selected session. The entire LOO CV training procedure was then repeated using these shuffled targets, yielding a coefficient of determination and mean MSE. This shuffling procedure was repeated 1,000 times. We then calculated a p-value for the better-than-chance significance of our models’ results by computing the percentage of shuffles which yielded a higher coefficient of determination or lower MSE.To tease apart the contributions of the neural and physiological features, partial correlations were also computed. Here, both fALFF features were used as covariates to PR, while only PR was used as a covariate for each of the fALFF features. As we are interested in exploring whether the neural signal contains additional useful information on top of physiological signals, and since the two fALFF features are highly correlated with each other, including one as a covariate for the other would artificially lower its partial correlation.RMEQ prediction time-seriesTo extract time-varying RMEQ total score predictions, the model (full model, using all features) was evaluated at every time point during the ketamine session. Prediction on the validation set for each participant is shown in Fig. Functional connectivityConnectivity computationThe time series of each channel was further processed by applying a high-pass filter using a bank of discrete cosine basis functions, followed by a low-pass filter using an acausal finite impulse response (FIR) filter. The pass band was 0.01–0.1 Hz, as is customary in the functional connectivity literature. Functional connectivity was computed for all pairs of channels as the Pearson correlation coefficient between their time courses, for each chromophore. For the baseline sessions, we selected a 400 s segment of data ending before the breathing challenge, from 100 to 500 s, to compute connectivity. This segment of data can be considered to be a resting session.For the dosing sessions, we selected a 400 s segment of data, from 1100 to 1500 s, to compute connectivity. We chose this segment to be near the end of the session (maximal expected neural effect), and of the same duration as the baseline segment.Global brain connectivity (GBC)For each channel in the prefrontal plates, a connectivity strength score was obtained by averaging all positive connections to other channels measured over the whole head. A connectivity strength map (limited to the prefrontal plates) was thus computed for each session.Statistical ComparisonsWe used the baseline session as a baseline for all connectivity analyses. We reasoned that the metric of interest was a change in the state of participants, and using a within-visit reference allowed us to focus on this change. For functional connectivity comparisons, we computed the median of the within-network connections and between-network connections separately for each network and each participant, then used a 