Brain imaging and neuropsychological assessment of individuals recovered from a mild to moderate SARS-CoV-2 infection

Significance In this case–control study, we demonstrate that non-vaccinated individuals recovered from a mild to moderate severe acute respiratory syndrome coronavirus type 2 (SARS-CoV-2) infection show significant alterations of the cerebral white matter identified by diffusion-weighted imaging, such as global increases in extracellular free water and mean diffusivity. Despite the observed brain white matter alterations in this sample, a mild to moderate SARS-CoV-2 infection was not associated with worse cognitive functions within the first year after recovery. Collectively, our findings indicate the presence of a prolonged neuroinflammatory response to the initial viral infection. Further longitudinal research is necessary to elucidate the link between brain alterations and clinical features of post-SARS-CoV-2 individuals.


Methods
Image processing

Estimation of cortical thickness
Surface-based morphometry was conducted in the Computational Anatomy Toolbox for SPM (CAT12) 5 for cortical surface reconstruction and estimation of mean cortical thickness employing the projection-based thickness method 6 , as well as topology correction 7 and spherical map-

Preprocessing
QSIPrep 0.14.2 1 was also used for preprocessing of diffusion-weighted MRI (dMRI). MP-PCA denoising as implemented in MRtrix3's dwidenoise 9 was applied with a 5-voxel window. After MP-PCA, Gibbs unringing was performed using MRtrix3's mrdegibbs. 10 Following unringing, B1 field inhomogeneity was corrected using dwibiascorrect from MRtrix3 with the N4 algorithm. 3 FSL (version 6.0.3:b862cdd5)'s eddy was used for head motion correction and eddy current correction. 11 Eddy was configured with a -space smoothing factor of 10, a total of 5 iterations, and 1000 voxels used to estimate hyperparameters. A linear first level model and a linear second level model were used to characterize eddy current-related spatial distortion.space coordinates were forcefully assigned to shells. Field offset was attempted to be separated from subject movement. Shells were aligned post-eddy. Eddy's outlier replacement was run. Data were grouped by slice, only including values from slices determined to contain at least 250 intracerebral voxels. Groups deviating by more than 4 standard deviations from the prediction had their data replaced with imputed values. Final interpolation was performed using the jac method.
A deformation field to correct for susceptibility distortions was estimated based on fMRIprep's fieldmap-less approach. 12 The deformation field is results from co-registering the b0 reference to the same-subject T1w-reference with its intensity inverted 13 . Registration was performed with antsRegistration (ANTs 2.3.1), and the process regularized by constraining deformation to be nonzero only along the phase-encoding direction and modulated with an average fieldmap template. Based on the estimated susceptibility distortion, an unwarped b=0 reference was calculated for a more accurate co-registration with the anatomical reference. Several confounding time-series were calculated based on the preprocessed DWI: framewise displacement using the implementation in Nipype (following the definitions by 14 ). The head-motion estimates calculated in the correction step were also placed within the corresponding confounds file. Slicewise cross correlation was also calculated. The DWI time-series were resampled to ACPC, generating a preprocessed DWI run in ACPC space with 2mm isotropic voxels.
Many internal operations of QSIPrep use Nilearn 15 and Dipy 16 . For more details of the pipeline, see the section corresponding to workflows in QSIPrep's documentation.

Diffusion tensor imaging and free-water imaging
Fractional anisotropy (FA) and mean diffusivity (MD) were derived from diffusion tensors which were modelled based on preprocessed dMRI using a least-squares fit. 17,18 Further, we employed free-water imaging, a two tensor model, modelling an extracellular compartment of isotropic diffusion, as well as a cellular compartment characterized by hindered/restricted diffusion. 19 Thus, by means of a regularized non-linear fit, free-water, and free-water corrected diffusion tensors were estimated for each study participant from which FA of the tissue compartment was calculated (FAT). 19 Fixel-based analysis pipeline MRtrix3 (v.3.0.2) 20 was utilized to estimate fiber density (FD), fiber cross section (FC), fiber density and cross section (FDC), and complexity (CX) at the voxel-level.
First, the preprocessed DWI was upsampled to a voxel size of 1.25 x 1.25 x 1.25 mm 3 , after which multi-tissue fiber response functions were estimated using the dhollander algorithm. 21 Fiber orientation distributions (FODs) were subsequently estimated via constrained spherical deconvolution 22 (CSD) using an unsupervised single-shell-optimized multi-tissue method (MRtrix3Tissue (https://3Tissue.github.io). 23,24 FODs were intensity-normalized using mtnormalize 25 after which a study-specific unbiased FOD template based on 20 healthy controls and 20 post-SARS-CoV-2 individuals was created with MRtrix3's population_template function. Next, individual FOD images and brain masks were non-linearly registered to the white matter FOD template. Transformed brain masks were used to compute a template mask, i.e. the intersection of all subject masks in template space. In the next step, fixels (= fiber populations within a voxel) were segmented from the FOD template within the template mask, resulting in a template fixel mask which was further refined to respect crossing fibers while excluding false positives (empirically derived crossing fiber and false positive thresholds: 0.06 and 0.18, respectively).
Next, following probabilistic whole-brain tractography based on the FOD template 26 (angle 22.5, maxlen 250, minlen 10, power 1, 20 x 10 6 streamlines, cutoff 0.06) and sphericaldeconvolution informed filtering of whole-brain tractograms 27 , deep learning based tract segmentation was performed with TractSeg 28 . The resulting tract segmentations were utilized to extract averaged diffusion indices for 72 major white matter tracts which served as features for logistic regression models predicting group membership.
Moreover, the transformed individual FOD images were segmented to derive fixels and their apparent FD. 29 Then, fixels of all subjects in template space were reoriented based on the local transformation at each voxel in the warps used previously. Subsequently, each subject's fixels were assigned to template fixels enabling statistical analysis of common, i.e., corresponding, fiber populations. FC was derived from non-linear warps generated during registration of individual FODs to template space after which the logarithm of FC (Log. FC) was calculated to ensure a zero centered normal distribution. 29 FDC was calculated as the product of FC and FD. Based on the whole-brain streamlines tractogram, a fixel-fixel connectivity matrix was computed which was then used for smoothing the fixel metrics FD, Log. FC and FDC.
In order to derive fixel metrics on the voxel-level, they were averaged across all fixels within a voxel using MRtrix3's fixel2voxel function. Moreover, CX, a metric of crossing-fiber organization, was calculated. 30 In an effort to allow for comparisons with conventional diffusion tensor imaging and free-water imaging markers, a study specific FA template was created. Therefore, previously derived individual non-linear warps from native FOD to FOD template space were used to register FA maps to FOD template space. These FA maps in FOD template space were then averaged and the resulting study-specific FA template served as the registration target for nonlinear transformations of FA images from native space to template space utilizing ANTs' SyN registration. 31 The resulting transformations were subsequently applied to the remaining maps of diffusion tensor imaging and free-water imaging metrics.

Tract-based spatial statistics
In order to derive skeletonized maps of each of the estimated diffusion parameters, we conducted tract-based spatial statistics (TBSS) 32,33 utilizing the above described study-specific FA template as the registration target. Briefly, individual FA images in template space got eroded to exclude non-brain voxels on the outer edge of the image. Next, a valid mask containing only the intersection of all subjects' brains was derived and used to mask the average of all previously eroded FA images. This mean FA image was subsequently used to derive a white matter skeleton which was thresholded at FA > 0.25. Next, all individual FA images were projected onto the mean FA skeleton. The resultant projection vectors were used to skeletonize all of the remaining diffusion metrics including those from free-water imaging and fixel-based analysis pipelines. Finally, diffusion markers were averaged across the entire white matter skeleton for further statistical analysis.

Peak width of skeletonized mean diffusivity
Peak width of skeletonized mean diffusivity (PSMD) was calculated based on standard procedures 34 adapted in terms of the non-linear registration step for which we used ANTs' SyN registration. 31 PSMD is calculated as the difference between the 95 th and 5 th percentile of MD values on the white matter skeleton in standard (MNI) space. A mask supplied by the developers was used to exclude white matter areas susceptible to partial volume effects of cerebrospinal fluid (https://github.com/miac-research/psmd/blob/main/skeleton_mask_2019.nii.gz).

Fluid-attenuated inversion recovery MRI
White matter hyperintensity segmentation FSL's Brain Intensity AbNormality Classification Algorithm (BIANCA) 35 with LOCally Adaptive Threshold Estimation (LOCATE) 36 were applied on FLAIR images and T1w images for white matter hyperintensity (WMH) segmentation.
The training dataset for the supervised k-nearest neighbor algorithm (BIANCA) comprised nearly 100 WMH masks manually segmented on FLAIR images by two independent raters. The manual segmentations of both raters were inclusively added into one binary mask for each participant and served as the training dataset for BIANCA 35 and LOCATE 36 . For the training of BIANCA 35 , the following images were used as input: 1) a brain-extracted FLAIR image; 2) a transformation matrix based on a linear registration from FLAIR space to standard MNI space (MNI152NLin2009cAsym template) utilizing FSL's FLIRT tool; 37,38 3) a T1w image rigidly registered to FLAIR space with AntsRegistration; 39 4) the manually segmented WMH masks. We used 3D patches, selected the non-lesion points from "no border", and chose 2.000 training points and 10.000 non-lesion points as segmentation parameters. Based on the initial training of BIANCA 35 , the testing dataset was automatically segmented using the same input images and parameters but without manual segmentations. Participants included in the training dataset were segmented in a leave-one-out validation by defining the query subject parameter in BIANCA. 35 The raw output masks of BIANCA 35 were used as input for LOCATE 36 . In addition to the BIANCA 35 masks, LOCATE 36 further received as input brain-extracted FLAIR and rigidly registered T1-weighted images in FLAIR space, as also used in BIANCA 35 . Moreover, further inputs were a ventricle distancemap created based on Freesurfer v.7.1 output, 40 the manual segmentations for training and a brain mask in FLAIR space. After training with LOCATE 36 , participants included in the training dataset were again segmented with a leave-one-out validation. The remaining participants in the testing dataset were automatically segmented.
After application of the BIANCA 35 and LOCATE 36 algorithms, the segmentations were further refined using Freesurfer v.7.1 parcellations 40 to exclude non-white matter regions. Specifically, a dilated cortical ribbon mask, an eroded ventricle mask, and parcellations from the corpus callosum and basal ganglia were used as exclusion masks. Lesion clusters were filtered for a minimum cluster size of 5 voxels as defined by the 6-connectivity. Finally, the lesion load was retrieved after normalizing for intracranial volume, as calculated by Freesurfer v.7.1. 40 Quality assurance Quality assurance (QA) of MRI data was conducted both quantitatively and qualitatively. First, a neuroradiologist reviewed all imaging data for pathologies. Further, for raw data, quantitative QA measures were derived for T1w and dMRI data utilizing MRIQC 41 , Freesurfer 42 and QSI-Prep. 1 Qualitative QA of raw imaging data was subsequently performed for outliers in framewise displacement and number of slices with signal dropouts (dMRI), as well as cortical thickness, brain volumes and coefficient of joint variation (T1w), in each case defined by ± 2 standard deviations from the mean. The quality of FLAIR images was assessed visually. Last, all derivatives of neuroimaging pipelines were visually assessed in order to ensure appropriate processing.

Machine learning prediction
To further evaluate the predictive capacity of derived imaging markers, they separately served as input to a comparative supervised machine learning pipeline. Therefore, average cortical thickness within Desikan-Killiany atlas parcels 43 was computed and diffusion markers were averaged within predefined anatomical fiber tracts from TractSeg outputs. 28 Per marker, multivariate logistic regression models were trained to predict whether a participant has COVID-19. Elastic net penalties 44 were applied for model regularization and SAGA served as the underlying optimization algorithm 45 . As a binary categorization task was performed, models were scored using accuracy which is also the metric we report. Accuracy scores provide an intuitive measure facilitating between-marker comparison and evaluation of a marker's diagnostic merit beyond abstract effect sizes derived from inferential statistics. Model training, corresponding parameter optimization and evaluation were conducted in a 10-fold nested cross-validation setup (L1-ratios=0.1,0.5,0.7,0.9,0.95,0.99,1, nCs=10) to prevent data leakage and consequent overfitting. Optimization procedures were repeated 100 times with different random split regimens in the cross-validation to make sure that prediction results were not biased by a single arbitrary split. 46 To assess whether prediction performance was statistically significant median accuracy scores obtained from aforementioned analysis were compared to the accuracy distribution of null prediction models where group membership was randomly permuted (npermutations = 1000). The prediction analysis was performed using scikit-learn (v1.0.2). 15

Results
Main analyses Figure S1. Matching results visualized as a balance plot Figure S1 shows the standardized mean differences between the healthy control and post-SARS-CoV-2 groups for each matching variable before (unadjusted, red) and after matching (adjusted, turquoise). The closer the standardized mean difference is to zero, the more similar the groups are. Each matching variable is depicted separately, i.e., from top to bottom, age, sex, years of education, hypertension, hyperlipidemia, smoking and diabetes.    Table  S3.   Table S4.

Correlations of clinical measures with mean diffusivity
Regression analyses with free-water, mean diffusivity and age Table S5. Results of exploratory regression analyses with age.  Table S5.