Search for schizophrenia and bipolar biotypes using functional network properties

Abstract Introduction Recent studies support the identification of valid subtypes within schizophrenia and bipolar disorder using cluster analysis. Our aim was to identify meaningful biotypes of psychosis based on network properties of the electroencephalogram. We hypothesized that these parameters would be more altered in a subgroup of patients also characterized by more severe deficits in other clinical, cognitive, and biological measurements. Methods A clustering analysis was performed using the electroencephalogram‐based network parameters derived from graph‐theory obtained during a P300 task of 137 schizophrenia (of them, 35 first episodes) and 46 bipolar patients. Both prestimulus and modulation of the electroencephalogram were included in the analysis. Demographic, clinical, cognitive, structural cerebral data, and the modulation of the spectral entropy of the electroencephalogram were compared between clusters. Data from 158 healthy controls were included for further comparisons. Results We identified two clusters of patients. One cluster presented higher prestimulus connectivity strength, clustering coefficient, path‐length, and lower small‐world index compared to controls. The modulation of clustering coefficient and path‐length parameters was smaller in the former cluster, which also showed an altered structural connectivity network and a widespread cortical thinning. The other cluster of patients did not show significant differences with controls in the functional network properties. No significant differences were found between patients´ clusters in first episodes and bipolar proportions, symptoms scores, cognitive performance, or spectral entropy modulation. Conclusion These data support the existence of a subgroup within psychosis with altered global properties of functional and structural connectivity.


INTRODUCTION
Schizophrenia may include a heterogeneous population with diverse cerebral alterations that may underlie its clinical variability. Several lines of evidence support this contention. A meta-analysis showed increased variability within schizophrenia in relevant structural measurements (Brugger & Howes, 2017). Moreover, structural magnetic resonance imaging (MRI) data could discriminate subgroups within schizophrenia with different clinical and biological characteristics (Lubeiro et al., 2016;Planchuelo-Gómez et al., 2020). Specifically, the existence of a group within the schizophrenia syndrome characterized by higher cortical curvature and lower cortical thickness was revealed (Lubeiro et al., 2016). Moreover, a cluster could be found in schizophrenia and bipolar disorder patients characterized by global cortical thinning associated with cognitive deficits (Planchuelo-Gómez et al., 2020). Neurophysiological data were also found to be different between treatment-resistant (TR) and non-TR schizophrenia patients (Molina et al., 2008). These differences included greater clinical severity in TR patients and a more severe profile of alterations in cerebral anatomical and electrophysiological parameters. Deficit and nondeficit schizophrenia patients showed different structural network properties (Wheeler et al., 2015). Furthermore, white matter abnormalities discriminated between first-episode (FE) patients with or without severe negative symptoms (Sun et al., 2015). Attending to studies of mixed diagnostic groups of psychosis, another study defined three biotypes cutting across the schizophrenia and schizoaffective and bipolar diagnoses based on cognitive, neuroanatomical, and neurophysiological data (Clementz et al., 2016), highlighting the existence of a subgroup characterized by worse cognition and widespread gray matter (GM) deficits. Furthermore, another study revealed a subgroup of schizophrenia patients with widespread volumetric reductions and worse cognitive deficits (Weinberg et al., 2016). Finally, a recent study aimed at identifying subgroups within schizophrenia and bipolar disorder patients based on their neurocognitive profile revealed the existence of a cognitive severely impaired group which showed higher symptom scores, a hypersynchronic basal connectivity state, and lower fractional anisotropy of frontal tracts (Fernández-Linsenbarth et al., 2021). Identifying biotypes in psychosis such as schizophrenia and bipolar disorder may contribute to consider them not as homogeneous entities but as syndromes, that is, a collection of symptoms and signs that may have different substrates. Thus, studies based on this consideration of the likely existing heterogeneity in psychosis could contribute to a better understanding of their neurobiological underpinnings, biomarkers definition, and the development of personalized treatments.
The use of cerebral parameters related to the mental functions could be useful in searching pathophysiologically meaningful biotypes within psychosis. Although the cerebral substrates of mental activity are incompletely understood, some facts can be reasonably assumed and may be useful for the purpose of exploring the presence of those biotypes. Mental activity is likely based on the fast-evolving synchronization of neural assemblies distributed across the brain (Buzsáki & Draguhn, 2004;Varela et al., 2001). The electroencephalogram (EEG) reflects the bioelectrical signal resulting from such synchronization; hence, EEG can be a useful tool for the analysis of the substrates of mental functions. One advantage of the EEG is its high temporal resolution, due to the swiftly evolving synchrony of neural assemblies underlying cognition (Dehaene & Changeux, 2011;Uhlhaas & Singer, 2010).
Moreover, higher mental functions such as those altered in schizophrenia likely implicate many cerebral regions, supporting the consideration of global cerebral networks rather than single-electrode local measurements in analyses aimed at studying possible subtypes of psychosis.
In the present study, our aim is to explore the possibility of discriminating meaningful biotypes within psychosis based on the global network properties of the EEG and their modulation during cognitive activity. In previous studies, we have shown significant differences in these network properties between schizophrenia patients and controls (Gomez-Pilar, de Luis-García, Lubeiro, de la Red, et al., 2018) and between schizophrenia and bipolar patients (Cea-Cañas et al., 2020). The hypothesis here is that these network properties would be more severely altered in a subgroup of patients also characterized by alterations in other biological, cognitive, and/or clinical measurements.
We included bipolar patients due to the clinical and genetic overlapping between these syndromes, considering that these patients are included within the psychosis spectrum. Moreover, in order to discard a major effect of chronicity in the results, we also included first episode patients thus being able to carry out repeated comparisons of these patients within possible resulting subgroups.

Symptoms assessment
Symptoms were scored with the Positive and Negative Syndrome Scale (PANSS) (Kay et al., 1987)

EEG-based brain graphs and connectivity strength calculation
In the construction of EEG-based brain graphs, network nodes are a mathematical representation of the EEG electrodes, whereas the values of the network edges are calculated from the neural coupling between each pair of electrodes (Stam & van Straaten, 2012). This coupling can be estimated with different methods. Here, we selected the phase-locking value (PLV) across successive trials (Lachaux et al., 1999), which is sensitive to low-amplitude oscillatory EEG components (Spencer et al., 2003) in addition to nonlinearities (van Diessen et al., 2015).
The PLV, in turn, can be computed using different methodologies; we used the continuous wavelet transform (CWT) using the convolution of each trial with a scaled and translated version of the complex Morlet wavelet. Thereby, the phase information from each trial is computed (Bob et al., 2008) considering cones of influence to remove edge effects (Torrence & Compo, 1998).
Applying the CWT approach for the performance of filter and phase extraction in one operation, the PLV between two signals, x(t) and y(t), was obtained evaluating the variability of the phase difference across successive trials (Gomez-Pilar, de Luis-García, Lubeiro, de Uribe, et al., 2018;Lachaux et al., 1999): where N t is the number of trials, Δφ xy is the instantaneous phase difference between the signals x and y, k is the time interval, and s is the scaling factor of the mother wavelet.
We generated functional connectivity matrices using the PLV values. Due to the fact that no threshold was applied, these connectivity matrices ranged between 0 and 1; 0 was obtained when two signals had no synchronization, and 1 was obtained when two signals were perfectly synchronized.
We selected two windows from the EEG signal: (i) the prestimulus window, which corresponded to a period of expectation before the stimulus onset from −300 ms to the stimulus onset and (ii) the response window, which is related to the P3b response (150-450 ms after the stimulus onset). Thus, the prestimulus window is located during task performance and is completely different from resting state. This procedure was applied both for the EEG theta band (

Graph parameters
Once the functional connectivity matrices were obtained, the resulting matrices were analyzed by means of different parameters from graph-theory field in order to characterize global connectivity properties of the brain network. Specifically, the present study focused on four parameters of the brain network: segregation (clustering coefficient [CLC]), integration (path length [PL]), small-world index (SW), and CS (Gomez-Pilar, de Luis-García, Lubeiro, de la Red, et al., 2018). These parameters were computed in two windows: prestimulus (300 ms before stimulus onset) and response (150-450 ms from the stimulus onset, centered around the P300 peak). The corresponding difference between the response and prestimulus window was considered as the modulation value, measuring the degree of change of graph parameters across time, that is, chronnectomics. Prestimulus and modulation network values were used in clustering analysis. Complete details of the graph parameters calculations can be found in the Supporting Information.
It is important to note here the bias introduced by volume conduction effects in EEG studies (Brunner et al., 2016), particularly in approaches involving connectivity metrics. These effects hamper the estimation of the connectivity between the real activity sources. Probably the most used technique to tackle this problem is to perform a source analysis, that is, to solve the inverse problem. These approaches, such as low-resolution tomography (LORETA), are able to identify activity sources from which connectivity metrics can be computed.
However, they are not completely bias-free, and the inverse solution problem remains unsolved, especially for nonhigh-density EEG (e.g., Hassan et al., 2014;Michel et al., 2004). A well-known alternative strategy to minimize field spread is based on the assumption that volume conduction affects the connectivity in a similar way in two different experimental contrasts, such as the prestimulus and response conditions (Bastos & Schoffelen, 2016). Therefore, the comparison between these two conditions, as has been done in this study, reduces volume conduction effects.

Spectral entropy
Furthermore, in previous studies, we identified a deficit of brain activity modulation with cognition in schizophrenia patients during a P300 task using the spectral entropy (SE) parameter (Bachiller et al., 2014; Gomez-Pilar, de Luis-García, Lubeiro

Structural data
The structural data were based on the assessment of GM morphometry from T1-weighted data and structural connectivity from diffusionweighted MRI data. MRI acquisition details can be found in the Supporting Information. The acquisition parameters were the same employed previously to characterize psychosis subgroups from MRI data (Planchuelo-Gómez et al., 2020).

MRI processing
Using the T1-weighted images, the segmentation pipeline from FreeSurfer (http://surfer.nmr.mgh.harvard.edu) version 6.0.0 was employed for the automatic cortical parcellation of GM regions (Dale et al., 1999). The average cortical thickness and subcortical GM volume were extracted from regions included in the Desikan-Killiany atlas (Desikan et al., 2006). We restricted our analysis to 14 bilateral cortical regions (Table S1). Moreover, we calculated the GM volume of the hippocampus, thalamus, caudate, putamen, and pallidum, as in our previous work (Lubeiro et al., 2016).

Diffusion tensor imaging data
From the diffusion tensor imaging (DTI), the FA in connections between pairs of regions was assessed, following the processing pipeline described in Lubeiro et al. (2017). Anatomically constrained tractography was obtained using the diffusion-weighted data (Jenkinson et al., 2012), considering the FA as the structural connectivity metric of interest. The evaluated connections from the tractography were focused on regions from the prefrontal cortex (rostral middle frontal and superior frontal gyri) and the limbic system (entorhinal cortex, parahippocampal gyrus, and hippocampus). Connections in which null values were found in a third (or more) of the subjects were discarded. A total of 46 homolateral connections were analyzed.

Cluster extraction
Data were divided in main and replication datasets, according to the recruitment center. The former included 80 patients with chronic schizophrenia, 35 patients with FE, and 34 patients with BD. The latter was composed of 22 patients with chronic schizophrenia and 12 patients with BD. The EEG parameters employed in the cluster extraction were the PL, the SW coefficient (the CLC divided by the PL), and the CS, all three in the prestimulus and modulation windows, for a total of six variables included in the clustering analysis. These variables were normalized with Z-scores before proceeding to the clustering process to avoid a bias caused by the difference between the values from each parameter. We decided to employ exclusively EEG network parameters to focus on the functional connectivity differences within psychosis and to avoid the employment of an excessive number of parameters from diverse sources that may cause overfitting. Moreover, as secondary analysis, we evaluated structural connectivity and cognition in the clusters extracted from the EEG network parameters.
In the main dataset, 26 indices were employed to extract the optimal number of clusters (Supporting Information). A large number of indices was employed to avoid the bias that may be produced by the use of a single or few indices. The majority rule (the value obtained with a higher number of indices) was used as the selection criterion. In case of a draw, the lowest number was chosen. This process was implemented with the NbClust package included in R, where further details about this clustering process can be found (Charrad et al., 2014).
For a specific number of clusters according to the majority rule, the subgroups were extracted using the k-means clustering methods and 50 initial random centroids. The centroids with the best silhouette profile were chosen (Rousseeuw, 1987).
We used the Clustering Large Application (CLARA) method to assess the consistency of the k-means clustering (Kaufman & Rousseeuw, 1990). Briefly, the CLARA method is based on the kmedoids method, which is similar to the k-means method using the median instead of the mean values. Here, the dataset was divided into 50 random subsets and the Euclidean distance was used as metric. The classification results were compared with the k-means clusters.
To determine the most important features characterizing the subgroups, a linear discriminant function was computed. A jackknife procedure was employed to test the classification accuracy from the discriminant function (Severiano et al., 2011). The discriminant scores were separately computed for controls.
In the replication dataset, the previous procedures were also used.
The discriminant function obtained with the main dataset was applied to the patients from the replication dataset and the classification results were compared with the k-means clusters.
As secondary analysis, we repeated the clustering process excluding the BD patients, considering that some of them did not have psychosis and the possible bias that they may introduce into the clustering results.

Statistical analysis of the clusters
Chi-squared and t-tests were used to compare age, sex distribution, parental education level, positive and negative symptoms, illness duration, and treatment doses between patients' clusters.
Analysis of variance (ANOVA) followed by pairwise comparisons with Bonferroni correction was used to compare cognitive performance, SE modulation, structural connective network, and regional thickness between patients' clusters and controls. The p-values from the global ANOVA tests were corrected for multiple comparisons following the Benjamini-Hochberg false discovery rate procedure, grouping these p-values by sets of comparisons (e.g., EEG graph theory parameters). We repeated these comparisons between FE patients in each cluster using Mann-Whitney U-tests to discard a major effect of chronicity in the results, without further correction for multiple comparisons due to the relatively small sample size.
As in our previous studies (Gomez-Pilar, de Luis-García, Lubeiro, de Uribe, et al., 2018), the SE modulation in individual sensors was introduced in a factor analysis and the resulting factor scores were used in subsequent comparisons.

Cluster solutions
The optimal number of clusters in the main dataset according to the majority rule was determined to be 2. A principal component analysis (PCA) was carried out to summarize the information from the EEG data and illustrate the clusters, but it was not used for clusterization.
The reason of not using PCA in the analysis but for illustration is that Bearing in mind the relatively low number of variables in this analysis, we considered that there was no need to summarize the information and lose part of the original information. Figure 3 shows the cluster plot based on the PCA values, and Figure S1 shows the opti- There were no significant differences in patients' diagnoses distribution between clusters, considering the subjects from both datasets together (χ 2 = 0.136, df = 2, p = 0.93). There were no significant differences between clusters in age, illness duration, or sex distribution.
Compared to the k-means classification of the main dataset, the CLARA method obtained similar results, with 96.0% of the subjects classified in the same clusters. Table S2 shows the comparison of the clusters obtained with k-means and CLARA.
Regarding the jackknife analysis, 96.0% of accuracy was obtained with respect to the k-means estimation. The classification results and

F I G U R E 4
Scatter plot of the distribution of (a) prestimulus path length and modulation and (B) pre-stimulus small-worldness and modulation in the identified clusters. Circles represent healthy controls (HC), triangles represent patients from cluster 1 (C1), and squares represent patients from cluster 2 (C2). The ellipsoids have a radius of 1 SD the coefficients of the discriminant function are shown in Tables 1 and   S3. In Table 1, it can be observed that the modulation of SW and PL are the factors with the highest weight. Figure 4 depicts the residuals plot with the value of those factors in the identified clusters. More details regarding the differences of the EEG network parameters are shown in Section 3.3. With respect to the classification of the subjects from the replication sample using the discriminant function from the main dataset, 82.4% of the subjects were classified in the same cluster as using the k-means method. The classification of the replication dataset is shown in Table S4.
The comparison of the discriminant values showed that patients from C2 presented higher scores compared to patients from C1 and controls (p < .0001 in both cases), while lower scores were identified in C1 patients with respect to controls (p < .0001; Figure 5). In

F I G U R E 5
Violin and box plots illustrating the discriminant scores of the patient subgroups and healthy controls from the main dataset. C1, cluster 1; C2, cluster 2; HC, healthy controls; LDA, linear discriminant analysis Figure 5, it is worth noting that the outliers could give a wrong impression about the differences between the two clusters groups and controls. The distribution of C2 scores is displaced to considerably higher values compared to controls, and the distribution of C1 scores to lower values. Similar trends were observed with the scores of the replication dataset, but with no statistically significant differences ( Figure S4).
Regarding the secondary analysis excluding the BD patients, the These values are shown in Table S5.

Cognitive and clinical comparisons
Current treatment doses were not significantly different between clusters (Table 2).
Positive and negative symptoms were not significantly different between C1 and C2 patients, nor between FE patients in these clusters (Table S6).
Cognitive performance according to BACS and WCST was not significantly different between clusters. Global IQ did not differ between clusters. There were no significant differences between C1 and HC in percentage of perseverative errors (WCST) nor between C2 and HC in problem solving and total IQ. For all other cognitive domains, scores were significantly lower in both patients' clusters as compared to HC (Table 2). When considered alone, FE patients in C2 showed lower performance than FE cases in C1 in working memory (U = 55, z = −2.15, p = .031) and, at trend level, in problem solving (U = 61.0, z = −1.7964, p = .10) (Table S6). There were no significant differences in social cognition scores between clusters.

EEG network parameters
A significant effect of the group was found for all variables of the functional network (4.44 < F < 62.48; .01 > p > .0001). In the post hoc com-parisons, prestimulus CS, CLC, and PL were higher in C2 as compared to C1 and HC. SW was lower in C2 compared to C1 and HC (Table 3).
Modulation of CLC and PL parameters was smaller in C2 compared to C1 and HC, indicating that these parameters decreased in the active window in C2 but increased in C1 and HC. SW modulation was higher in C2 (where SW value was higher in the response window) than in C1 and HC (where SW was lower in the response window). Modulation of CS was lower in C2, implying that CS increased more in C1 than in C2 (Table 3). p < .001) were also smaller than C1 FE patients, and modulation of SW (U = 4, z = −3.83, p < .001) was higher in C2 patients (i.e., its values were higher in the response) (Table S7).

Entropy modulation
As in our previous studies, a single factor summarized most of the variance for entropy modulation, with all sensors contributing positively to that factor. Factor scores for spectral entropy modulation values did not differ between clusters (C1 mean = 0.355, SD = 0.518; C2 TA B L E 3 Prestimulus network values and its modulation with P300 task Note: Data are shown as mean (SD). All the ANOVA results survive the correction for multiple comparisons. * p < .05; ** p < .001 as compared to healthy controls. # p < .05; ## p < .001 between patients' clusters.
The previous p-values are adjusted by corrections for multiple comparisons.   (Table 4).

Cortical thickness
Significant effects of group were found for bilateral caudal and ros- in both groups, although this decrease was more widespread in C2. We did not find significant differences in symptoms severity, cognitive performance, illness duration, or antipsychotic doses between clusters.
In our previous results using graph-theory applied to EEG data in schizophrenia we found several baseline alterations. We reported that schizophrenia patients showed higher CS values at baseline compared to HC (Gomez-Pilar, de Luis-García, Lubeiro, de la Red, et al., 2018) and to bipolar patients (Cea-Cañas et al., 2020). Besides, we also described in schizophrenia patients higher CLC values at baseline (Gomez-Pilar et al., 2017). These alterations were found in the schizophrenia patients when considered as a whole and could be expected to be more severe in a subgroup of cases such as the cluster 2 in the present study.
Since SW is the CLC/PL ratio, the higher modulation SW values in C2 is possibly consequence of the higher increase of PL in C1 and HC, One of such substrates may relate to an alteration of the functional network, but our data suggest that even with a normal functional network other factors may hamper cognition and underlie symptoms, which seems coherent with the data supporting the biological heterogeneity of schizophrenia substrates (Arnedo et al., 2015;Molina & Blanco, 2013;Volk et al., 2012). Our results suggest that an alteration in structural connectivity (only found in C2) may be associated with an altered functional network in only a subset of cases, although we cannot conclude that the former causes the latter. Perhaps the other patients' symptoms and cognition may be underpinned by a biochemical disbalance not reflected in the functional architecture but that could be caught by other kind of functional analyses. Furthermore, we assessed functional connectivity related to a task, but not resting-state functional connectivity. The specific behavior associated with a particular task may present different properties compared to resting-state functional connectivity. It has been reported that task-based and restingstate functional connectivity present different network properties (Di et al., 2013), considering additionally the distinct networks for specific tasks and resting-state conditions. Therefore, structural connectivity, cognition, and clinical symptoms may be more related to resting-state or other tasks rather than the assessed task-based functional connectivity of this study.
In support of such possibility, both clusters showed a decreased modulation of their EEG activity during an oddball task as measured with SE. We have described and replicated such a modulation deficit in schizophrenia and bipolar disorder as a possible biomarker for the altered function in this disorder . That deficit was associated with cognitive deficits and negative symptoms . This could explain why clusters in the present study did not differ in these clinical dimensions, since SE modulation is similar between the clusters extracted in this study.
In a previous study that aimed to explore the existence of clusters based on the neurocognitive profile, we reported that patients with more severe cognitive deficits were also characterized by higher prestimulus CS of the EEG network (Fernández-Linsenbarth et al., 2021).
Thus, we could expect a significant cognitive deficit in C2 with a higher CS, but this was not found. Nevertheless, in that study, both cognitive clusters showed a similar pattern of differences as compared to controls in other neurophysiological data, that is, the same decrease in SE modulation and a smaller CS modulation than HC. No other functional network characteristics were different between the two cognitive clusters identified in that study, where both clusters showed a significant deficit in cognition as in the present one. Probably, depending on the clustering criteria (cognition or EEG), the resulting group correlates may differ slightly, although the global pattern would be that both normal and altered EEG network characteristics at baseline and its modulation may be associated with different degrees of cognitive alteration in schizophrenia and bipolar patients. This could also be consistent with the possibility that altered functional network dynamics are only one of many possible factors leading to symptoms and cognitive impairment in these syndromes.
In the present study, C2 patients also showed altered structural connectivity network properties. In a previous report, we did not find a significant correlation between abnormal values of EEG network parameters and DTI-based network parameters in schizophrenia (Gomez-the way these disorders are conceived could also shed light on the traditionally lack of consistent replication of cerebral findings in them. As a future research line, the integration of MRI, cognition, and EEG data should be performed to understand the interactions between structure and function, and their relationship with cognition. Previously, the heterogeneity concerning GM morphometry and cognition has been assessed (Fernández-Linsenbarth et al., 2021;Planchuelo-Gómez et al., 2020). The separate evaluation of each area is important to understand the diverse individual mechanisms related to psychosis and the associations between the different paths of psychosis.
To conclude, it is worth noting that our data support the existence of different subgroups within psychosis and may contribute to considering schizophrenia and bipolar disorder not as homogeneous nosologically entities but as syndromes. This notion can also help address the lack of consistent results in previous literature. Finally, the characterization of subgroups could contribute to understanding underlying pathophysiological mechanisms. These mechanisms may be more easily identified through the analyses of biological parameters characterizing subtypes than by comparing schizophrenia patients and controls.
This could also raise the possibility of developing personalized treatments based on the most relevant altered underpinnings.

Limitations
Our study has limitations. The first limitation is its sample size.
Although a larger sample size would have been desirable, we were able to obtain solid results. Second, biological tests were not available in all subjects. Third, the lack of untreated patients mean that we cannot discard an effect of treatment on EEG values, although we have previously reported nonsignificant effects of antipsychotic, antidepressants, lithium, and benzodiazepines for these values . Finally, although the volume conduction effects were minimized by comparing two experimental contrast (i.e., EEG modulation), the results should be cautiously interpreted.

CONCLUSIONS
In conclusion, we found that an abnormal EEG-based connectivity network is present in approximately half of the patients with schizophrenia and bipolar patients in which significant anatomical changes related to GM cortical thickness and white matter connectivity were also found. These alterations seem independent of chronicity and antipsychotic treatment. EEG network alterations may characterize a biotype across schizophrenia and bipolar diagnoses.

ACKNOWLEDGMENTS
This work was supported by the following grants: "Instituto de Salud Carlos III" grant number PI18/00178, "Gerencia Regional de Salud de Castilla y León" grant number GRS 1721/A/18, and by predoctoral grants from the "Consejería de Educación -Junta de Castilla y León" and the European Social Fund grant numbers VA-183-18 to Inés Fernández-Linsenbarth, 376062 to Álvaro Planchuelo-Gómez, and VA-223-19 to Rosa M. Beño-Ruiz-de-la-Sierra. We appreciate the collaboration of patients and healthy controls in our research.

DATA AVAILABILITY STATEMENT
The datasets that support the findings of this study is available from the corresponding author upon request.

PEER REVIEW
The peer review history for this article is available at https://publons. com/publon/10.1002/brb3.2415