End-to-end Stroke imaging analysis, using reservoir computing-based effective connectivity, and interpretable Artificial intelligence

2.1 Data and preprocessing

The dataset was previously collected by the School of Medicine of the Washington University in St. Louis and complete procedures can be found in [10]. They collected MRI data and behavioral examinations of stroke patients and healthy controls. The imaging data comprise structural and functional MRI from controls and patients suffering from hemorrhagic and ischemic stroke. Acquisitions were done within the first two weeks of the stroke onset (i.e., acute). Structural scans include T1-weighted, T2-weighted, and diffusion tensor images. Functional images include a resting state paradigm. Scanning was performed with a Siemens 3T Tim-Trio scanner. Briefly, we closely followed the pre-processing steps outlined in [27]. Following a quality control of fMRIPrep outputs, 104 stroke subjects and 26 control subjects were qualified for further analysis. For our purposes, it suffices to say that structural scans were used in combination with functional acquisitions to co-register all participants into a common template. Gray matter signal was finally obtained after artifact removal and parcellated into 100 regions of interest (ROIS) [15, 38]. For every subject and patient, these 100 time series (i.e., one for each ROI) were fed into the pipeline outlined below to obtain the subject-specific effective connectivity maps.

The dataset is not public but it is available upon request to the original authors [10]. The used code is instead available at the URL https://github.com/Wotaker/Effective-Connectivity-Reservoir-Computing.

2.2 Reservoir computing

Reservoir computing networks (RCN), despite being known for more than two decades, have been largely eclipsed by other frameworks. A reservoir network is a set of artificial neurons that are randomly connected between themselves thus forming a recurrent architecture [26, 31]. Sometimes this is also called echo-state network since the internal dynamics of the reservoir (or “echo state”) maintain information about the system’s input history. In this framework, an input series u_t is fed into this high dimensional dynamical system of N units through a non-linear activation function, where Wⁱⁿ is an N x (N_in +1) matrix of random weights including biases, N_in is the dimensionality of the multivariate input, and f_in is the non-linearity. At each time step t the former projection is used to drive the reservoir units r_t. The current state of each unit is a combination of the past states as well as the current input, where W is an N x N matrix of random weights, and λ is the leakage that controls the importance of the reservoir’s history to the current time stamp t. The final component of the reservoir is a set of readout weights W^out that extract information from the hidden state and map onto specific predictions. That is, The predictions y_t might be of arbitrary dimension N^out and, importantly, are linear w.r.t. to the reservoir states. Within this paradigm, only that readout weights W^out are trained via incremental linear regression optimization [45, 28], with α being a regularization parameter, R is the matrix obtained after concatenating all the reservoir states, and Y contains the outputs. Once again, the readout weights contain a set of N_out biases.

Noteworthy, as opposed to other architectures suited for time series forecasting, only a reduced set of output weights needs to be trained, thus increasing its computational efficiency. The random weights Wⁱⁿ are drawn from a uniform distribution bounded between −1 and 1. The recurrent connections are drawn from a standard normal distribution and are later scaled by the spectral radius. The latter largely ensures that the network possesses the echo-state property, although there is recent evidence disagreeing with this aspect [52, 13].

Briefly, the main idea behind reservoir-like computing is that a given input pushes the reservoir to specific locations in a high-dimensional manifold [26, 31]; the output weights are then optimized to retrieve information from the nearby regions. Were the input to move the reservoir away to other points, the output weights would not be able to recover meaningful information hence completely missing the prediction. Further evidence suggests that RCNs supersede deep learning-based models for temporal series prediction even on the verge of chaos [40]. Richer approaches aim to train the reservoir connections themselves and have been proven to be useful in understanding the dynamical properties of cortical networks [35], offering an interesting framework for similar use cases. The parameter values used in our experiments can be found in Table 1.

2.3 Reservoir computing networks to map causal interactions in lesioned brains

Traditionally, effective connectivity in neuroimaging can be estimated in different ways, as dynamic causal modeling [20], GC [23], continuous-time implementations [21], or information theory [50]. Granger-like interpretations are often preferred due to their relative computational costs and implementation, though they are not exempt from controversy [24] thus justifying alternative approaches.

An unrelated proposal relies on the properties of the state-space of the dynamical system to reconstruct asymmetric mappings between delayed embeddings of each component of the system [46]. That is, it leverages Taken’s theorem to find the optimal neighborhood as well as the exact delay at which the reconstruction is optimal. Recent extensions [48, 51, 7] have incorporated non-linear methods as well as reducing the number of ad-hoc parameters. Most prominently, reservoir computing has proven to be an efficient and accurate alternative to automatize the process almost in its entirety [25].

Let’s consider the relationship between two one-dimensional variables, x and y, where it hypothesizes that the delay at which interactions take place is not smaller than the sampling rate (e.g., Time of Repetition in functional MRI). The prediction skill, denoted by ρ_x→y (τ), is defined as the Pearson correlation between the true time series, y(t + τ), and the predicted series ŷ(t + τ) from the input x(t).

Noteworthy, the Pearson correlation between the true and reconstructed series (ρ) is used to estimate directedness, though other metrics like mean squared error could also be used. Directionality can still be assessed using the same hypothesis testing mechanisms [48]. Moreover, the time series are fed into the reservoir all-at-once, letting the network project all of them. The neighboring points in the variable’s embedding are then remapped to the target embedding via the training of the output weights. It should noted that this represents a deviation from more canonical usages [25, 13]. To investigate the causal relationship between variables, we first calculate both ρ_x→y (τ) and ρ_y→x(τ) in a given temporal domain. We then examine the values of τ at which either ρ_x→y (τ) or ρ_y→x(τ) reaches its peak value [51, 25]. To streamline the subsequent description, we introduce the following notation:

Empirically, directionality is then defined as follows [46]:

• if τ_x→y is positive, and τ_y→x is negative, we say that x causes y;

• if τ_x→y is negative, and τ_y→x is positive, we say that y causes x;

• if both τ_x→y and τ_y→x are negative, we say that x and y causes each other.

Despite seeming counterintuitive, information of y is present in earlier observations of x and, consequently, that current information of the cause x is useful to predict future observations of the consequence y (see [46] for a comprehensive explanation). In certain systems, predictability scores peak at negative lags τ < 0 for both directions, being the height of the peaks informative of the coupling strength [25]. However, the existence of this bidirectionality does not necessarily invalidate the former statements [17].

It was quickly noted that in large and noisy networks, such as the brain, it is unlikely that the predictability scores in Eq. 5 reach clear and distinct peaks. Functional signals are notoriously noisy [43], and indeed prediction with this approach is challenging [4]. A solution to this issue relies on assessing the minimal requirements that are needed to suggest causal interactions [17]. For that, the difference between prediction scores should be evaluated and contrasted with proper surrogate predictions [34, 39, 30]. That is, which can be interpreted as an indication of the potential causality direction (Table 2). The scores in Eqs. 5 and 7 can be contrasted against the 95% confidence interval obtained from a surrogate testing procedure [17]. It has been shown that the requirements to define causality can be compressed into a reduced set of δ-scores [17], for directed interactions, and for bidirectional interactions. is the p-value after testing the alternative hypothesis H₁ against the surrogate data (Fig. 2). For instance, is a p-value for the hypothesis that x influences y. The values of the δ-scores range from 0 to 1, with higher values indicating greater confidence in the existence of a causal relationship with a coupling delay of τ between the examined variables.

• Download figure
• Open in new tab

Figure 2.

Predictability scores from an the same chaotic system defined in [46, 17]. Solid lines show the predictability in Eq. 5 between embeddings. Shaded regions show 1 standard error of the mean. Transparent lines show the predictability of the surrogate system, which is used to define the expected level of chance against which the hypotheses are tested. In this academic example, it can be said that strong asymmetric interactions between two time series exist at different temporal lags.

Then, for a given lag τ, a matrix A_τ collects the δ-scores, where each element [x, y] represents the causal relationship from node signal x to ROIs signal y,

The effective connectivity (RC) matrix A_τ is a final representation of the effective connectivity network of every subject; it is directed, non-symmetric, and can incorporate bidirectional causality connections. For our experiments, for every possible interaction x → y, we trained 20 different reservoirs and tested against 100 shuffled targets, strictly following what was outlined in [17]. Furthermore, only uni-directional connections were kept from the adjacency matrix in Eq. 10.

In our experiments, we investigated the classification of pathological groups with the effective connectivity matrices used as features (Fig. 3 TOP), and we also compared those to the effective connectivity matrices obtained by Granger causality, representing one of the state-of-art approaches. As a last step, for each entry A_τ [x, y], we standardized all samples by subtracting the mean connectivity of the control group and dividing by the standard deviation. Finally, these standardized causal relationships (i.e., directed graphs) were fed into two simple graph classifiers to explore and explain the most informative nodes and links to detect stroke occurrence.

• Download figure
• Open in new tab

Figure 3.

Working diagram of causality given by the reservoir computing (TOP) and graph convolutional architecture (BOTTOM).

2.4 Graph convolutional neural networks

Graph convolutional neural networks (GNNs) are a variation of traditional convolutional neural networks which capitalize on graph data representations and can learn non-trivial representations by leveraging the complex topological organization of the data [49]. Intuitively, a graph constitutes a non-Euclidean geometric space where complex relationships between data points can be embedded and forwarded as inputs into a GNN [8]. More formally, a graph 𝒢 = (𝒱, ℰ) is defined as a set of nodes 𝒱 = {1, …, n} and a set of edges ℰ = {(i, j) | i, j ∈ 𝒱} where (i, j) represents a link or interaction between the i-th and j-th nodes. Initially, each node i ∈ 𝒱 is associated with a column feature vector .

Every layer l of a GNN updates the hidden representation of each node by aggregating information from the neighborhoods: where are the new node representations, 𝒩_i is the neighborhood of the i-th node, f_θ denotes a nonlinearity, and F is a permutation-invariant aggregator. Several proposals exist for the aggregation operator, determining the expressive power, interpretability, learning stability, and scalability of the network [49].

The non-symmetric effective connectivity maps derived are also non-attributed, that is, there are no node features to be aggregated in Eq. 11. Although non-attributed graphs are classifiable, they dramatically increase the problem’s difficulty. Fortunately, the Local Degree Profile (LDP) method effectively decreases the challenge by setting the attributes of each node to local neighboring properties [9]. Thus, we computed the in and out degree of each node as well as the minimum, maximum, mean, and standard deviation of the in/out degree of its neighbors. This created a feature vector of dimension 10 that was propagated through the directed adjacency matrix for every subject. The neural network consisted of l = 2 hidden layers and was trained for 150 epochs with a learning rate of 0.005 to minimize the binary cross entropy between the predicted and true classes (Fig. 3 BOTTOM). The metrics were computed with a balanced class weight to account for the different number of samples in each class. The model was tested in a 10-fold cross-validation scheme and used a validation set to test for overfitting.

2.5 Local Topology Profile

A recent extension of the LDP attribution outlined before incorporates other local properties to the already-mentioned descriptors. This Local Topology Profile [1] has been shown to improve the accuracy over its parent version, namely LDP. Following the original proposal, we extended the feature vector with the edge betweenness centrality [22], the overlap between node neighborhoods (i.e., Jaccard index), and the local degree score [29].

However, as an attempt to further reduce the complexity of the workflow, we used the 13 LTP features with a random forest classifier of max depth 2 and a maximum number of features equal to 5. As in the GCN classifier, we used class weights to balance the dataset and used a 10-fold cross-validation scheme. The architecture used in practice is summarized in Figure 3.

2.6 Local Interpretable Machine-Agnostic Explanations

To explain the features allowing the classification we used the LIME (Local Interpretable Model-agnostic Explanations) approach. This technique explains the prediction of any classifier by learning the model locally around the prediction [37]. In our case, this was used to highlight the edges that contributed to the classification performance the most. LIME assigns a coefficient to each edge on the EC matrix based on the contribution to the final classification score.

Positive values were useful in identifying the stroke group, whereas negative values were consistent in identifying the control group. The total explainability values of each ROI were calculated for both groups separately. These values were thresholded with the arbitrary threshold of 0.02 for the stroke group and −0.02 for the control group (because these directions helped the correct decisions). Edges associated with wrong decisions were not studied due to their lack of meaning in neurological terms.

3.1 Effective connectivity maps derived from Reservoir Computing

EC maps were not readily interpretable given the complex interactions expected to occur at different spatial and temporal scales. Consensus stipulates that information transfer is obscured by the hemodynamic response function, which effectively masks the corresponding temporal delay between cause-consequence associations. We computed effective connectivity maps between 100 ROIs at two different delays (Time of Repetition = 1 and 2; see Fig. 4). The average maps showed clear patterns of hemispheric segregation while at the same time exhibiting strong connectivity between homotopic regions. In canonical functional connectivity studies, this a priori segregated structure can be considered as an initial quality assessment of the resulting maps, forming the basis for an accurate description of the functional relationships expected to occur in brain disease.

• Download figure
• Open in new tab

Figure 4.

Group averaged effective connectivity matrices for two different Times of Repetition. Top: −1 TR. Bottom: −2 TRs. The left column is the average of subjects suffering from a stroke located in the left hemisphere. The middle column is the average of subjects suffering from a stroke located in the right hemisphere. The right column is the average of the control group.

Even though stroke occurrence is not entirely random [10, 47], their exact morphologies and functional disconnection patterns are highly variable. We further examined the properties of the directed networks by computing the average directed connectivity for controls, subjects suffering from right-hemispheric stroke, and subjects suffering from a stroke located on the left hemisphere (Fig. 5).

• Download figure
• Open in new tab

Figure 5.

Global effective connectivity alterations between regions located in the same hemisphere (top) and between regions located in different hemispheres (bottom). Error bars depict 1 standard error of the mean. Insets show the average difference between left-left and right-right effective connectivity (top) and between left-right and right-left effective connectivity (bottom). Statistical significance was assessed via a two-sample t-test (‘*’ p<0.05). Global connectivities were obtained by averaging the weight value over the connections belonging to the corresponding hemispheres H.

Global hemispheric connectivity was computed by averaging the EC maps within and between hemispheres. That is, averaging the values in each on of the 4 visible squares in the average EC maps (Fig. 4). Briefly, intra- and inter-hemispheric connectivity was severely altered in all patients, showing a clear break of symmetric communication w.r.t. the control group, especially for right-impaired subjects [27].

3.2 Classification results

The results of the classification are reported in Tables 4 and 3 respectively for the GCN and LTP classifiers. Results are reported for both the proposed method and Granger Causality: Average AUC, accuracy, precision, recall, and F1 are reported. As expected, the LTP (augmented with a random forest classifier) generally increased the classification metrics, although both models are comparable. It should be noted that classifying effective connectivity graphs is a complicated task due to sample heterogeneity [12, 2], and that very similar scores compared to the chance levels (e.g., an increase of 0.2-0.3) are found in the literature [1].

3.3 Node and edge importance in stroke detection

We used the LIME explainability framework on the LTP classifier due to its slightly better performance and higher computational efficiency to highlight the most descriptive ROIs and edges related to stroke onset. Importantly, the explanations were done on top of the EC matrices obtained with the reservoir method and not the granger one. For each node in the EC networks, we summed all the explainability coefficients to assess the contribution of each connection arising in each node to the correct classification (i.e., sum over all columns). Lastly, binarized and thresholded explainability values were projected back to the surface mesh (Fig. 6; see also Methods and [27]). The resulting maps show that regions in visual, dorsal, and ventral attention have the most contribution to the classification performance for stroke subjects, while ventral attention and frontoparietal networks contributed the most to the detection of control subjects.

• Download figure
• Open in new tab

Figure 6.

Interpretation of the LIME explainability outputs for each group. Cortical projection of the total contribution of each ROI (left) and its association with one of the 7 resting-state networks. The Dorsal attention network is distinctively necessary to discriminate the presence of a lesion.