An updated systematic review and meta-analysis of brain network organization in focal epilepsy: Looking back and forth

Abnormalities of the brain network organization in focal epilepsy have been extensively quantified. However, the extent and directionality of abnormalities are highly variable and subtype insensitive. We conducted meta-analyses to obtain a more accurate and epilepsy type-specific quantification of the interictal global brain network organization in focal epilepsy. By using random-effects models, we estimated differences in average clustering coefficient, average path length, and modularity between patients with focal epilepsy and controls, based on 45 studies with a total sample size of 1,468 patients and 1,021 controls. Structural networks had a significant lower level of integration in patients with epilepsy as compared to controls, with a standardized mean difference of -0.334 (95% confidence interval -0.631 to -0.038; p-value 0.027). Functional networks did not differ between patients and controls, except for the beta band clustering coefficient. Our meta-analyses show that differences in the brain network organization are not as well defined as individual studies often propose. We discuss potential pitfalls and suggestions to enhance the yield and clinical value of network studies.


Introduction
A question that had received considerable attention is how and to what extent we can explain the detrimental consequences of epilepsy for normal brain functioning. These are illustrated by the high incidence of cognitive, neurodevelopmental, and psychiatric comorbidities in patients with epilepsy, even in those with seizures originating from a single epileptic focus (Kanner, 2016;Nickels et al., 2016). To this end, an increasing number of studies have been investigating global functional and structural brain integrity in patients with epilepsy. Particularly network analysis, a branch of graph theory, is used here and allows the schematic representation of the brain (Rubinov and Sporns, 2010;Stam, 2014). Both functional and structural networks can be computed. Functional networks can be described as graphs composed of nodes, representing brain areas, that are linked by edges, denoting functional synchrony (i.e. temporal coherence in activity) between these areas. Structural networks are also composed of nodes and edges, with the edges here representing the physical connections between brain areas (Rubinov and Sporns, 2010). Functional data is usually derived from time series observations using functional MRI (fMRI) and neurophysiological techniques, including electroencephalography (EEG) and magnetoencephalography (MEG) whereas structural data arises from structural MR imaging. Network quantification can be done using a variety of metrics, of which the path length, clustering coefficient, and modularityindicators of a network's level of integration, segregation, and hierarchy (Rubinov and Sporns, 2010) are most commonly used (Table 1). Converging evidence suggests that the healthy brain is characterized by a cost-efficient network configuration with a distinct balance of integration, segregation, and centrality ( Fig. 1) (Rubinov and Sporns, 2010;Stam, 2014).
Based on previous network studies we know that the interictal network configuration in epilepsy is altered (Engel et al., 2013;Richardson, 2012;van Mierlo et al., 2019). These alterations may contribute to the repertoire of clinical manifestations and (cognitive) impairments seen in both generalized and focal epilepsies (Yang et al., 2018). Results are, however, contradictive regarding the extent and directionality of the network abnormalities, which hampers a clear understanding of the network nature of epilepsy. How this holds true for focal epilepsies specifically remains an open question, despite previous aggregation efforts ( van Diessen et al., 2014;Yang et al., 2018). We therefore updated and extended a previous review on the brain network organization in focal epilepsy (van Diessen et al., 2014). Since its publication, a considerable number of new studies have been published, also providing epilepsy type-specific brain network data. By combining all available quantitative whole-brain network data, we aim to provide a more accurate and focal epilepsy type-specific characterization of the functional and structural brain network organization. This could clarify the current understanding of the network nature of focal epilepsy and its comorbidities.

Search strategy
Studies were identified by searching the online PubMed (NCBI) database. We used terms related to epilepsy, neurophysiological and imaging techniques, and network analysis to build the search query. The last search was conducted on March 23rd, 2020. No limiting filters were applied. Details on the search strategy are presented in Supplementary  Table 1.

Selection criteria and definitions
Studies were included if they compared interictal, resting-state global brain networks of focal epilepsy patients with control subjects using the average clustering coefficient, average path length, or modularity (Table 1), irrespective of acquisition technique and network construction method. Focal epilepsy was defined as temporal, frontal, parietal, and occipital lobe epilepsy, benign epilepsy with centrotemporal spikes, and any other unspecified focal epilepsy. All measures should be given as mean or median values with uncertainty estimates. Studies reporting data on epilepsy patients who underwent brain surgery before network construction were excluded, as were studies only reporting combined data on focal and generalized epilepsy types (Allen et al., 2019;Babajani-Feremi et al., 2018;Chiosa et al., 2019;Garcia-Ramos et al., 2016a;Garcia-Ramos et al., 2017).

Data extraction
One of the authors [GS] performed the literature search and screened all titles and abstracts. From potentially eligible articles, full-text versions were reviewed for inclusion and cross-referenced. Data extraction was carried out by one author [GS] using a data collection form and checked at random by another author [WMO]. Disagreements were resolved by discussion with a third author [EvD].
When multiple studies used an identical patient population, we included the study reporting the most comprehensive data. If studies provided outcomes only as graphs, data were manually extracted from these using GraphClick. This approach has proven to be reliable and reproducible (Boyle et al., 2013;Flower et al., 2016). Corresponding authors were contacted in case of missing network data. The following data items were extracted from each study:

General study information
First author, journal, year of publication.

Study population
Group sizes, sex distribution, specific epilepsy diagnosis, age at data acquisition, age at epilepsy onset, duration of epilepsy, anti-epileptic drug usage.

Acquisition technique
Data acquisition modality, sampling frequency (for neurophysiological studies), field strength (for imaging studies).

Network properties
Network size, binary or weighted network, presence or absence of metric normalization.

Network metrics
The mean and standard deviation (SD) of the average path length, average clustering coefficient, and modularity were extracted. If means and SDs were not provided, medians were interpreted as means, and standard errors (SE) and interquartile ranges (IQR) were converted into SDs using the following formulas: Standard error: SD = SE • √ group size Interquartile range: SD = IQR / 1.35 Since studies were allowed to differ in acquisition technique and network construction method, we predefined what metric data to extract. More exactly, if data was provided for a range of sparsity levels, thresholds, or densities, values corresponding to the highest sparsity or Table 1 Definitions of network metrics.

Metric Description
Path length Measure of integration. The average shortest path length is defined as the average minimal number of edges between all possible pairs of nodes in a network. The average shortest path length is equal to the inverted global efficiency (Stam and Reijneveld, 2007). Clustering coefficient Measure of segregation. The average clustering coefficient is defined as the ratio of the number of existing edges between a node's neighbors and the maximum possible number of edges between these neighbors, averaged for all nodes in a network. The average clustering coefficient reflects the average prevalence of clustered connections around individual nodes in a network (Stam and Reijneveld, 2007).

Modularity
Measure of segregation and hierarchy. The modularity is defined as the degree to which the brain is subdivided into modules: groups of nodes with a maximum possible number of within-group edges and a minimum possible number of between-group edges (Rubinov and Sporns, 2010).
Nodes reflect the different brain regions and edges reflect the functional or structural connections between nodes. The healthy brain network is characterized by a balance of (1) path length (red) and clustering coefficient (orange) resulting in a normal order network (horizontal axis) and (2) varying levels of connectedness (vertical axis). This balance gives rise to a hierarchical modular network, in which each module (blue) is composed of smaller components and simultaneously is part of a larger component (diagonal axis). This figure is partly based on Fig. 1 of Stam (2014).
threshold and the lowest density were used. From studies reporting data on both binary and weighted networks, we only extracted weighted data. Only normalized data were used if normalized and non-normalized data were provided (van Wijk et al., 2010;Zalesky et al., 2010). If studies reported data on different network sizes, we extracted data for the largest network. Network metrics reported for subgroups only were combined into a single value by weighted mean calculation. The average path length was defined as the reciprocal of the global efficiency if only the latter was reported. We did not take into account minor differences in the definitions of the network metrics of interest, such as the difference between the average clustering coefficient calculated as the global clustering coefficient or the local averaged clustering coefficient, as it will not influence the eventual results from the meta-analyses. For the neurophysiological studies we were primarily interested in theta frequency band as it is the most consistently reported frequency band altered in focal epilepsy (Bartolomei et al., 2006;Horstmann et al., 2010;van Dellen et al., 2012;van Diessen et al., 2013b). To explore possible differences in the other frequency bands, we extended the main meta-analyses (i.e. including all patients) to the delta, alpha, and beta frequency band. We did not include the gamma frequency band as converging evidence suggest that neurophysiological surface recordings above 20 Hz are contaminated by myogenic artifacts (Whitham et al., 2007). To support future updates of our study, all included studies' characteristics and network data are available via the Open Science Framework: https://osf.io/q7ad4/ (DOI: 10.17605/OSF.IO/Q7AD4).

Quality assessment
Study quality was assessed using an adapted Newcastle-Ottawa Scale (Stang, 2010;Wells et al., 2013). This scale is developed for the evaluation of non-randomized studies and includes the following domains: patient and control selection, comparability of the study groups, and the ascertainment for the exposure of interest, interpreted as the assessment of the brain network organization for our review. The scale uses a star system, wherein higher-quality studies get more stars, with a maximum of eight (Supplementary Table 2).

Methods of analysis
Meta-analyses were performed in R statistical software using the metaphor package (Viechtbauer, 2010). Effect sizes, expressed as standardized mean differences (SMD), for the average path length, average clustering coefficient, and modularity were estimated by applying random-effects models using restricted maximum likelihood estimation to the data (DerSimonian and Kacker, 2007). We used random-effects models since these tend to yield more conservative results by taking into account both intra-study and inter-study variability (DerSimonian and Kacker, 2007). The SMD was used as a summary statistic as we expected scale variability of the metrics between studies due to methodological differences. The SMD has the advantage of estimating effect sizes in a scale-free way (Cummings, 2011). The contribution of each study to the pooled SMD was given by its weight, based on the inverse variance of the study's effect estimate.
Based on previously reported inconsistencies between different recording modalities, we analyzed structural and functional data separately (Garcia-Ramos et al., 2016b;Jiang et al., 2017). As functional MRI data and neurophysiological data (EEG/MEG) have distinct temporal and spatial properties from which connectivity measures are inferred (Bullmore and Sporns, 2009), we ran combined and separate analyses on these two data types. We estimated metric effect sizes for all focal epilepsies as one group, and for temporal lobe epilepsy (TLE) and benign epilepsy with centrotemporal spikes (BECTS) separately as these groups were large enough to perform sub-analyses on. For an overview of all different analyses performed, see Supplementary Fig. 1. All primary analyses included all eligible studies, irrespective of the year of publication. We ran separate meta-analyses for studies published after 2013 that were not included in our previous meta-analysis (van Diessen et al., 2014). For each analysis, we detected, discussed, and removed outliers using funnel plots and the influence diagnostics of the R metaphor package (Viechtbauer and Cheung, 2010).
Heterogeneity between studies was assessed using the I 2 statistic. In contrast to Cochran's Q, the I 2 statistic does not directly depend on the number of studies included and thus might overcome excessive test power in large meta-analyses (Deeks et al., 2020;Higgins et al., 2003). I 2 describes the percentage of variation across studies due to heterogeneity rather than chance, with heterogeneity being quantified as 'low', 'moderate', or 'high', with I 2 around 25 %, 50 %, and 75 %, respectively (Higgins et al., 2003;Melsen et al., 2014). Effect size estimates were considered significant at a p-value < 0.05.

Meta-regression
Meta-regression was performed to explore sources of heterogeneity in the SMD estimates of the functional and structural average clustering coefficient, average path length, and modularity. We regressed the individual study SMDs against the (1) mean duration of epilepsy, (2) percentage of patients (in a study) with a structural brain lesion causing epilepsy, and (3) percentage of patients using anti-epileptic drugs, irrespective of the number of drugs. These factors were previously identified as potential modifiers of brain networks; see for an overview (van Diessen et al., 2013a). Regressions were performed if data from at least five studies were available.

Study selection and quality
We identified 2,992 articles with our search. From 153 articles full text was screened and reference lists were inspected. Fifty-seven studies met the inclusion criteria for our review. Twelve of seventeen contacted corresponding authors did not provide us the requested additional data; their studies were therefore not eligible for inclusion. Eventually, 45 studies were included in our meta-analysis (Fig. 2). These studies together included 966 patients and 962 controls for the functional analyses and 502 patients and 400 controls for the structural analyses. Information on the study population and technical characteristics of the included studies are shown in Tables 2 and 3. Data from 26 studies were extracted using GraphClick.
The quality of the included studies was variable with a range in Newcastle-Ottawa scores between three and eight stars (Supplementary Table 2). Lower scores were mainly the result of incomplete reporting. Based on the funnel plots  and influence diagnostics (not shown) we excluded outliers from the final analyses. Excluded articles are noted in the forest plots' legends.

Modularity
Since only two studies reported structural data on network modularity (Garcia-Ramos et al., 2019;Jiang et al., 2017), we only aggregated functional data from seven studies (Doucet et al., 2015;Garcia-Ramos et al., 2016b;Jiang et al., 2017;Pedersen et al., 2015;Ridley et al., 2015;Vaessen et al., 2013;van Dellen et al., 2012). There was no difference in functional network modularity between patients with focal epilepsy and controls, with an estimated SMD -0.096 (95 % CI − 0.588 to 0.396; p-value 0.702), based on 188 patient and 162 control networks (Fig. 5). Studies were highly heterogeneous with I 2 being 78.7 %. Separate analyses for fMRI and neurophysiological data were not performed as only one study provided neurophysiological data. For this reason, we were also not able to analyze data from other frequency bands than the theta band.

Benign epilepsy with centrotemporal spikes (BECTS)
A total of six functional (Adebimpe et al., 2016(Adebimpe et al., , 2015Besseling et al., 2014;Choi et al., 2019;Ji et al., 2017;Xiao et al., 2015) and two structural studies (Besseling et al., 2014;Garcia-Ramos et al., 2019) reported data on BECTS patients. We only aggregated functional data. The functional average clustering coefficient and average path length did not significantly differ between patients and controls. Modularity data was not available. Heterogeneity between studies was low to moderate, with I 2 being 0.0 % and 43.3 % for the path length and clustering coefficient, respectively. Results are summarized in Supplementary Table 5 and Supplementary Fig. 7.

Recent studies
Results of meta-analyses ran with studies published after 2013 were largely comparable with the results of the analyses including all studies, but additionally yielded a significant decrease in the functional average clustering coefficient in the epilepsy group and TLE subgroup. Heterogeneity between study data remained moderate to high, except for the functional and structural clustering coefficient data in the TLE subgroup (Supplementary Table 6). Analyses on BECTS studies were not repeated, since all BECTS studies included were published after 2013.

Study group sizes
Study groups sizes, in terms of the number of patients and controls included, modestly increased over years ( Supplementary Fig. 8).

Meta-regression analyses
A significant relation was found between the standardized mean difference of the structural path length and the percentage of patients using anti-epileptic drugs in the total epilepsy group (p-value 0.003, based on five studies). No other significant relations were found between the network metrics of interest and the mean epilepsy duration, the percentage of patients using anti-epileptic drugs, and the percentage of patients having a structural brain lesion in the epilepsy group, nor in the TLE subgroup (Table 5). Brain lesions included: tumors, hippocampal damage, gliosis, vascular malformations, cortical dysplasia, atrophy, and aspecific lesions. None of the included papers provided information on the size of the brain lesion. Therefore, we were not able to take this into account in our regression analyses. Due to a lack of data, no separate meta-regression analyses were performed for the BECTS group.

Discussion
In this systematic review with meta-analyses, we pooled data of all available studies on the functional and structural interictal brain network organization in patients with focal epilepsy. In contrast to most individual studies reporting significant network alterations, aggregation of these results yielded only modest effect sizes. We found a significantly increased structural path length in patients with focal epilepsy as one A Focal epilepsy is noted when studies did not further specify the epilepsy type or included a group of patients with various focal epilepsies. B In mean (SD) or median (min-max). C Percentage of patients in a study having a structural brain lesion causing epilepsy. D Percentage of patients in a study using anti-epileptic drugs, irrespective of the number of drugs. E Study also reported structural network data (see Table 3).
group, and in the subgroup of those with TLE, pointing at a lower level of structural integration. When looking at functional networks, we did not find any network differences between patients and controls, except for the clustering coefficient in the beta frequency band, which was significantly lower in epilepsy. No network metric differences at all were found for patients with BECTS. The individual studies' data were heterogeneous. An accurate quantification of the interictal brain network organization in epilepsy can provide valuable insights into the development of epilepsy and its comorbidities in addition to conventional neuroimaging and neurophysiological techniques (Garcia-Ramos et al., 2016a;Kemmotsu et al., 2014;Paldino et al., 2017;Vlooswijk et al., 2011). In a previous systematic review and meta-analysis, we concluded that the interictal brain network organization in patients with focal epilepsy is characterized by a lower level of integration and a higher level of segregation as compared to controls, with similar results for functional and structural networks. Modularity was not assessed. In the present effort, we were able to include a more than tripled number of studies (n = 45) and participants (n >1400) and found comparable results for the network's integration, albeit only in the structural analyses. We were not able to reproduce the significantly increased average clustering coefficient in focal epilepsy patients using fMRI and theta band data. Instead, we found a significantly decreased clustering coefficient in the beta band, that has not been evaluated previously (van Table 3 Characteristics of included studies reporting structural network data. AED: anti-epileptic drugs; BECTS: benign epilepsy with centrotemporal spikes; CT: cortical thickness; DSI: diffusion spectrum imaging; DTI: diffusion tensor imaging; DWI: diffusion weighted imaging; F: female; M: male; NR: not reported; SIFT: spherical-deconvolution informed filtering of tractograms; T: tesla; TLE: temporal lobe epilepsy; -: not applicable. A Focal epilepsy is noted when studies did not further specify the epilepsy type or included a group of patients with various focal epilepsies. B In mean (SD) or median (min-max). C Percentage of patients in a study having a structural brain lesion causing epilepsy. D Percentage of patients in a study using anti-epileptic drugs, irrespective of the number of drugs. E Study also reported functional network data (see Table 2). Diessen et al., 2014). Some studies, however, selectively reported frequency band data (e.g. only beta or theta band data) which might have given rise to reporting bias. Whether the beta band differences truly represent a pathophysiological mechanism thus remains unclear.

Challenges of network studies
Despite the inclusion of more sophisticated network metrics in recent network studies, the majority of studies are still primarily focused on the path length and clustering coefficienttwo highly correlated metrics (Li  The following studies were outliers and excluded from analysis: Bartolomei et al. (2006), Quraan et al. (2013), Wang and Meng (2016) (functional) and Lemkaddem et al. (2014), Garcia-Ramos et al. (2019) (structural). et al., 2011) -for network characterization. Since these metrics do not take into account the essential modular and centrality characteristics of the brain network, informative network data could be neglected. Some of the included studies reported centrality measures, but these were too few to be pooled in a meta-analysis. A related issue comprises the subjective methodological decisions made throughout the network analysis pipeline that can directly influence network metrics and thereby the studies' conclusions (Bastiani et al., 2012;van Diessen et al., 2015;van Wijk et al., 2010).
Most studies use traditional statistical approaches to investigate differences in brain network topology between patients and controls. Previous work by our group revealed that in certain cases no single brain network metric is significantly altered, but instead a combination of metrics (van Diessen et al., 2013c). This not only advocates the inclusion of more network metrics in the analyses but also the use of more sophisticated statistical approaches. Here, machine and deep learning algorithms might be of use (Paulus et al., 2019;Roy et al., 2019). Another limitation is the group size of most studies. Although group sizes modestly increased over time, the chance of studies being unpowered remained high, which may hamper the interpretation of results.

Opportunities and limitations of meta-analyses
Meta-analyses provide a unique opportunity to overcome the limitations of smaller and explorative studies but face limitations. Group analyses often are restricted to the mean values reported per study and are of limited use when the variability between studies is high. We also observed this in our analyses, despite correcting for heterogeneity using random-effects models and the SMD as a summary estimate. High heterogeneity could be an explanation for not having found robust changes in network metrics in patients with epilepsy as compared to controls. A solution to the highly variable between-study-summary estimates is the pooling and re-analysis of individual participant data (IPD). This IPD meta-analysis approach is known to provide more reliable estimates than the standard aggregated data meta-analysis but requires the collection of the original datasets across studies, which is very time consuming and depends on the active participation of corresponding authors in the sharing of data (Simmonds et al., 2005;Stewart and Parmar, 1993;Tudur Smith et al., 2016).

Where do we go from here?
While focal epilepsy was traditionally thought of as a local brain disorder, growing evidence exists that brain abnormalities extend beyond the epileptic zone. This is clearly illustrated by two recent studies from the ENIGMA consortium that reported widespread white matter abnormalities and brain atrophy to be present in focal epilepsy Larivière et al., 2020). Our meta-analyses add to this evidence from a network perspective by showing global alterations in the organizational structure of epileptic brain networks. Since we only assessed the global network organization, it remains unclear to what extent the global alterations are the result of regional network disturbances in and near the epileptogenic zone. Previous studies suggested connectivity patterns and network disturbances to be partly lateralized, with both hemispheres being affected but the most significant disturbances being present in the hemisphere of the epileptogenic focus (Englot et al., 2016;Pourmotabbed et al., 2020). The mechanism(s) underlying global interictal network disturbances in focal epilepsy are poorly understood and might be related to underlying white matter damage (Kim et al., 2008;Otte et al., 2012) or recurrent network inhibition resulting from seizures, as conceptualized in the 'network inhibition hypothesis' (Englot et al., 2016). This latter hypothesis is supported by studies reporting a significant association between brain network disorganization and duration of epilepsy Park et al., 2018;van Dellen et al., 2009). We were not able to confirm this hypothesis with our meta-regression analyses. Due to missing data, we could not take into account seizure type and frequency, which might  be as important as net epilepsy duration. Sub-analysis for specific epilepsy types revealed distinct and more homogenous network alterations in TLE and BECTS. Further in-depth characterization of the brain network organization in different epilepsy types might enable network analysis to be clinically useful for predicting disease severity and patient counseling.
Other clinical applications of network analysis in epilepsy have also been proposed. It might support early-stage epilepsy diagnosis in both children and adults (Douw et al., 2010;van Diessen et al., 2013c) and is proposed to be useful in identifying patients at risk for developing cognitive impairments Vlooswijk et al., 2011). Furthermore, network analysis might be helpful in epilepsy management, as previous studies found that network dynamics and connectivity patterns could predict patients' responsiveness to certain anti-epileptic drug treatment (Anderson et al., 2020), epilepsy surgery planning (Barron et al., 2015;Chiang et al., 2015), and the outcomes of epilepsy surgery (Di et al., 2019;Gleichgerrcht et al., 2018;Wilke et al., 2011). Although these studies together frame the clinical potential of network analysis at the group level, none of the proposed applications has made it towards individual clinical use yet.

Conclusion
Increased sample sizes and methodological standardization are needed to better grasp the specific network alterations in epilepsy. This is in line with the conclusions of a recent systematic review on network alterations in particularly generalized epilepsy (Pegg et al., 2020). The increasing level of experience with network analysiscoming with time might already have started to pay off with regard to methodological harmonization, since we observed less between-study heterogeneity when aggregating only the more recent studies' results. The slightly increasing sample sizes over time might also have positive impact here. Secondly, the inclusion of a greater variety of metrics to characterize the brain network organization might better underpin the complex brain network nature and its alterations in epilepsy. The choice for a combination of network metrics requires careful consideration and should, ideally, not be researcher-dependent. As for the application of advanced statistical methods to imaging or neurophysiological data, deep learning algorithms might be of use here. Machine and deep learning approaches could, for instance, help to select from a predefined set of metrics which of these are most distinguishing between people with and without epilepsy (van Diessen et al., 2013c). One step further and transcending network analysis, strategies can also be built to classify patients with epilepsy and controls with raw imaging and neurophysiological datasets where no data is predefined or no labels are available (Roy et al., 2019). Thirdly, since the brain is in a constant state of change, characterization of the brain network organization does not only require the calculation of static network metrics but may benefit from taking the dynamic properties into account (Pedersen et al., 2017;Rosch et al., 2018). Combining the dynamics of brain regions and their interconnectedness could lead to a more comprehensive understanding of epilepsy-associated network abnormalities. Taken together, we call for the deepening of network analysis and a continuous critical evaluation of its yield as a first step to enhance the field's clinical value for the individual epilepsy patient.

Author contributions
EvD, GS, and WO designed and initiated the study. Data collection and analyses were performed by EvD, GS, and WO. The manuscript was written by EvD, GS, KB, and WO. All authors read and approved the final version of the manuscript.

Role of the funding source
The study was supported in part by the Ming Fund (GS) and a clinical research fellowship from the University Medical Center Utrecht (EvD). The funders had no role in study design, data collection and methodology, decision to publish, or preparation of the manuscript.

Declaration of Competing Interest
None.

Acknowledgements
The study was supported in part by the Ming Fund (GS) and a clinical research fellowship from the University Medical Center Utrecht (EvD).  The funders had no role in study design, data collection and methodology, decision to publish, or preparation of the manuscript.