Introduction

Severe traumatic brain injury (TBI) is a leading cause of death and disability worldwide and across all ages [1], and the mortality rate seems to be unchanged over the past 25 years [2]. TBI affects 50 to 60 millions of new cases per year, with 2.5 million occurring in Europe, and this is the reason why it has been declared a priority for public health policy [3]. Efforts should be put into investments in research and other disciplines, such as clinical management and prevention policies. The clinical challenge is represented by the fact that the patient group is extremely heterogeneous and the pathology is highly dynamic [4]. Therefore, treatment protocols and prediction models are difficult to assess.

The critical care management in the acute phase is focused on reducing the probability and impact of secondary insults, which develop over time as a consequence of raised intracranial pressure (ICP) and/or reduced cerebral perfusion pressure, among other mechanisms [4]. Therefore, continuous ICP monitoring is recommended by international guidelines as a standard of care in all surviving patients with severe TBI to provide information for ICP-directed therapy [5, 6].

Along with ICP, a great deal of neuromonitoring techniques and imaging modalities can be used to improve the understanding of intracranial pathophysiology, which, as mentioned before, is highly heterogeneous and dynamic [4]. Ultimately, more precise targets for therapies could be suggested with this integrated approach. In this perspective, the interaction between the brain and other organs has been suggested as one of the mechanisms that could potentially explain the complexity of this pathology.

In particular, much attention has been given in recent years to the study of interactions existing between the brain and heart [7,8,9]. The dynamical interplay between the two organs is thought to ensure physiological functions and to be involved in pathological conditions, too [7]. For example, in the work by Valenza et al. [7], the authors describe episodes of paroxysmal sympathetic hyperactivity, which often happen when there is a severe axonal injury. In the work by Valenza et al. [7], the authors describe paroxysmal sympathetic hyperactivity in the postresuscitation syndromes after serious anoxic-ischemic brain insults. An interesting pattern of frequency of these paroxysms has been noticed, as described in detail in the work by Valenza et al. [7]. These observations show the importance that the brain–heart coupling has in pathological events. As a consequence, efforts have been made toward exploring analytical methodologies that could tackle this phenomenon, with the final goal of developing metrics describing the brain–heart interaction [8, 10].

Nevertheless, this area is still wide open for investigation, particularly in TBI. In the work by Gao et al. [9], the authors presented an interesting analysis of interaction between brain and heart measures, showing the presence of interaction and Granger causality between ICP, mean arterial pressure and heart rate (HR). In our previous studies, we derived an HR–ICP measure that we denominated “brain–heart crosstalks,” defined as transient elevations of ICP and HR that occur simultaneously [11]. In that work, we presented multiple novelties. For the first time, to the best of our knowledge, the measure of brain–heart crosstalks was defined. Moreover, we presented a novel sliding window method to detect the presence of these events. We studied the crosstalks defined as so in a pediatric population and subsequently in a single-center study with an adult cohort [12], where we also conducted a pilot analysis of the relationship of our novel interaction metric with mortality.

In a further work, we modeled the coupled HR–ICP system as a multilayer network [13] that can be imagined as a framework in which different channels of the same overall modeled structure are included [14]. In this framework, each channel is represented by a layer and each node can maintain different neighbors and characteristics across different domains. In recent years, multilayer networks have been fruitfully applied to a variety of fields such as, for example politics, medicine, economics, social interaction, and time series [14]. Its strength is the capability of modeling relationships across variables (i.e., layers), and therefore their interaction and integration. In our past work [13], to the best of our knowledge, we used for the first time a multilayer network approach to model the relationship between ICP and HR during brain–heart crosstalks events.

In this work, we aimed to integrate several brain–heart crosstalks measures and examine their role in the context of mortality prediction models, also taking into account clinical, demographics, and monitored features that collectively reflect the severity of primary injury. For this purpose, we used the Collaborative European Neurotrauma Effectiveness Research in TBI (CENTER-TBI) high-resolution intensive care unit (ICU) cohort of patients with TBI [15, 16]. Our main hypothesis was that the integrated brain–heart crosstalk measures are statistically related to outcome in patients with TBI.

Methods

Material

The material used in the present work included data collected in 226 patients with TBI as part of the high-resolution CENTER-TBI cohort [15, 16] . This subset consists of patients treated with no external ventricular drain. In the cohort analyzed, 50 patients died. The patient outcome was recorded at 6 months after injury. In Table 1 we present the variables we used as input features. Patients were recruited prospectively between the beginning of 2015 and the end of 2017 from 21 centers across Europe. All patients were admitted to the ICU for their TBI during the course of the study, with high-frequency digital signals recorded from their ICU monitors during the course of their ICU stay.

Table 1 Low resolution variables used as input features

All the patients had invasive ICP monitoring according to the Brain Trauma Foundation guidelines [6].

ICP was acquired using an intraparenchymal strain gauge probe (Codman ICP MicroSensor; Codman & Shurtleff Inc., Raynham, MA) or parenchymal fiber optic pressure sensor (Camino ICP Monitor, Integra Life Sciences, Plainsboro, NJ, US). Arterial blood pressure was obtained through either radial or femoral arterial lines connected to pressure transducers (Baxter Healthcare Corp. CardioVascular Group, Irvine, CA).

All signals were recorded using digital data transfer or digitized via an A/D converter (DT9801; Data Translation, Marlboro, MA), where appropriate, sampled at frequency of 100 Hz or higher, using the ICM + software (Cambridge Enterprise Ltd., Cambridge, UK) or Moberg CNS Monitor (Moberg Research Inc., Ambler, PA, USA) or a combination of both [17]. Signal artifacts were removed using both manual and with automated methods. Post-acquisition processing was conducted using ICM + (https://icmplus.neurosurg.cam.ac.uk). HR was determined by calculating the Fourier transform and finding the fundamental frequency of the arterial blood pressure waveforms over a 10-s window, updated every 10 s. 10-s moving average filter was applied. For each patient, we considered the whole monitored period available that corresponds to the whole duration of ICP monitoring, with initiation of recording within 24 h of the injury.

Low resolution data, such as demographic characteristics, admission, and injury related variables, were retrieved as per version 2.0 of the CENTER-TBI database. In Table 2, we present the summary of the variables describing the dataset used in the present work comparing distributions between the survivors and nonsurvivors.

Table 2 Overview of the variables used in the model

Methods

Brain–Heart Crosstalks Detection and Analysis

In our own previous work [11], we proposed an algorithm to detect the presence of brain–heart crosstalks given two time series of ICP and HR monitored in patients with TBI.

Briefly, the algorithm is based upon a sliding window approach. A window of 10 min of observations is considered jointly in HR and ICP. If a simultaneous increase in HR and ICP, of at least 20% with respect to the minimum value of the window frame, followed by a subsequent decrease in both, is detected, then a brain–heart crosstalk is counted. Such measure was defined to facilitate the study of patterns of interaction between HR and ICP and used as a proxy of brain and heart interaction. Full details of the algorithm can be found in [11]. In the present study we applied our own algorithm just described, to detect the total number of brain–heart crosstalks for adult patients with TBI.

Since the length of observations varied for each patient (depending on various factors, for instance the duration of hospitalization), the raw number of detected brain–heart crosstalks was normalized by the total length of the analyzed time series (ctnp). Figure 1 shows an example of brain–heart crosstalk, with the event marked by the blue rectangle.

Fig. 1
figure 1

Example of one brain–heart crosstalk. In the figure we present time trends for more than 10 min of observations for HR and ICP. The blue rectangle denotes the presence of a simultaneous increase of HR and ICP. The event was detected using the sliding window approach. (HR heart rate, ICP intracranial pressure)

In our later work [13], the behavior of ICP and HR time series, during brain–heart crosstalks events, was further investigated using multilayer network modeling. In particular, the system ICP-HR was modeled as a multiplex network, where each layer represents one of the time series. The conversion between time series and graphs works as follows: nodes represent time points, and an edge between two nodes is present when the criterion of the “horizontal visibility” is met [18]. Briefly, two nodes are connected if they “can see” each other horizontally, i.e., if all the values separating the two in time are of lower, or equal, magnitude. Because in our case we have the same sets of nodes for each layer, we used a special type of multilayer network named multiplex, where each node in a layer is connected to the corresponding one in the other layers.

An example of a multiplex network model, in which one layer represents ICP and the other represents HR, is shown in Fig. 2.

Fig. 2
figure 2

Figure showing the multilayer network model. Each layer is a time series. The top layer here represents HR, whereas the lowest layer represents ICP. Each node is a time stamp, and connections between two nodes are defined according to the horizontal visibility criterion [13]. (HR heart rate, ICP intracranial pressure)

Using the complex network approach, we obtained two network measures: the interlayer mutual information (\({mi}_{ct})\) and the average edge overlap \({(\omega }_{ct}\)). The average edge overlap \({\omega }_{ct}\) quantifies the coherence of the overall graph, and the higher it is, the higher the coherence of the graph layers. The interlayer mutual information \({mi}_{ct}\) models the degree distribution of the two graphs. Because one measure reflects the edges, and the other focuses on the nodes, the two measures complement each other, and they report a complete picture of the joint behavior of the two analyzed time series. They were therefore chosen as representatives of the graphs corresponding to the detected brain–heart crosstalk event windows. Complete details on the whole multiplex modeling can be found in our previous work [13].

The three defined measures of crosstalks, i.e., the count \({ct}_{np}\), the interlayer mutual information \({mi}_{ct}\) and the edge overlap \({\omega }_{ct}\), were added to the set of predictive features in the mortality model.

Mortality Model Predictions

In order to build a predictive mortality model we applied an elastic-net logistic regression machine learning approach [19, 20] obtaining a model for each of the four age groups in which the population was segmented, as described in the Results section. The machine learning approach used, is a slight modification of the standard logistic regression, enriched with a lasso and a ridge regression to overcome multicollinearity issues [19, 20]. It also performs an automatic features selection, embedded in the optimization resolution problem. Two parameters need to be set for the elastic-net model: \(\lambda\) and \(\alpha\), which are respectively used to tune the lasso and ridge contribution in the optimization function. In our experiment \(\alpha\)=0.5 and \(\lambda\) was defined through a leave one out cross validation procedure, in which turn n-1 number of observations are used as training set and the remaining one is used as test set. The \(\lambda\) values selected for each model as well as the performances of each model (in terms of percentage of null deviance explained) are specified in Table 3, at the top of each specific section discussing the model for each age range. Each model was fed with the set of input features described in Table 1 together with the newly introduced brain–heart crosstalk indicators.

Table 3 \(\lambda\), percentages null deviance and coefficients selected by the elastic-net model are shown in the figure

Statistical Analysis

The brain–heart crosstalks algorithm as well as the statistical analysis was conducted using R software [21] version 3.6.0. The point biserial correlation test was performed to verify the correlation between brain–heart crosstalks measures and mortality. Pearson correlation coefficient was used to test correlation between input features. Kruskal–Wallis test was used to verify the presence of statistically significant difference between distributions of brain–heart crosstalks variables in surviving and nonsurviving patients with TBI.

Results

Crosstalks and Preliminary Statistical Analysis

A mean number of 57 brain–heart crosstalks, per patient, with a high standard deviation of 58, were detected. The value of brain–heart crosstalks was normalized by the length of the time series observed per patient, and its correlations with the other predictive variables are reported in Fig. 3. The apparent “squared patterns” of correlations, starting from the left upper corner of the correlation matrix, are likely to correspond to similar features, for instance lab blood samples analysis or CT features. These subgroups appear to be naturally correlated, although other correlations, even if milder, appear in the rest of the plot.

Fig. 3
figure 3

Pearson correlation of the predictive variables for the adult cohort. The matrix is symmetric, and we included here all the clinical, lab, imaging, as well as brain–heart crosstalks variables. (CRP C reactive protein; CT computed tomography; \(ct_{np}\) brain–heart crosstalks; DAI diffuse axonal injury; \(mi_{ct}\) mutual information; SOL shift of the midline structures; ωct average edge overlap)

The relationship between mortality and brain–heart crosstalks was first examined using the distribution analysis of brain–heart crosstalks and related network measures, comparing survivors and nonsurvivors. Figure 4 shows the boxplots of the distributions of ctnp, ωct and \(mi_{ct}\) for the entire cohort (n = 226).

Fig. 4
figure 4

Distribution of the brain–heart crosstalks measures for survivors and nonsurvivors. The p value of the Kruskal–Wallis test performed between the two distributions is also shown in the figure. ( \(ct_{np}\) brain–heart crosstalks)

The point biserial correlation between mortality and ctnp in adults is − 0.13 with a significant p value of 0.04. This seems to suggest an inverse relationship existing between mortality and brain–heart crosstalks in the cohort analyzed. Same can be said for ωct which present − 0.19 as a point biserial correlation (with significant p value of 0.002) and for \(mi_{ct}\) with − 0.09 (with a p value slightly more than the significant threshold, and equal to 0.13).

Moreover, as it can be appreciated from the ctnp boxplot in Fig. 4, the visual difference between the survivors and nonsurvivors distribution of the brain–heart crosstalks measures was confirmed statistically. The Kruskal–Wallis test to compare the distributions of the three measures returned significant p values for all of them: 0.02 for \({mi}_{ct}\), 0.005 for ctnp and 0.006 for ωct.

Mortality Model

The age distribution of the patients included in this study ranged between 16 and 85 years old and is presented in Fig. 5 panel a. Age has been shown to have a significant impact toward TBI prognostic and outcome prediction [22], also because of its influence on the likelihood of comorbidities. To decrease the heterogeneity of patients analyzed in the mortality models, the whole dataset was divided into four age category subgroups. In panel b of Fig. 5 we show the number of patients with the associated age range after the splitting process, whose composition is homogeneous within groups and heterogeneous across groups. In other words, the division into four age groups was performed in a way in which we could obtain four groups of homogeneously aged patients, with comparable sample size, in order to be statistically possible the comparison between groups. In this way we could still include age as a predictor variable in each of the four models, while avoiding using it as a confounding factor, when considering the entire population with no age distinction. In Table 3, we show the results for the models fitted for each age group, including their deviance ratio explained, in the header of each table. Coefficients for all the features in the model are presented. The best fitting models appears to be for the age populations range aged 16–29, 50–64, and 65–85. Indeed, the deviance ratio is high which means that the predictors included in the model explain with very high reliability the mortality outcome. Interestingly, the negative relationship between mortality and ctnp in all the age populations is confirmed as well as the one between mortality and the average edge overlap during crosstalks. In so far as the other populations are concerned, the best fit of the model appears to be for the age range 65–85. We also performed experiments with and without the newly introduced brain–heart crosstalks measures. The deviance ratio actually does improve for the case of age range 16–29 (from 49 to 79%) and 65–85 (from 93 to 99%). In Table 3, we only reported the value of the significant coefficients. If a dot is in the corresponding cell, it means the coefficient is not significant. We also performed further experiments using only the IMPACT (International Mission for Prognosis and Analysis of Clinical Trials [23]) variable scores as predictors and adding subsequently also brain–heart crosstalks measures to the model to evaluate if the IMPACT variables would have benefited from the introduction of our brain–heart crosstalks variables. In particular, we used the standard IMPACT variables of age, motor score, pupils, hypoxia, hypotension, CT variables, glucose, and hemoglobin for predicting mortality using the logistic regression approach and subsequently we added the brain–heart crosstalks measures to the IMPACT scores and evaluated performances of the models so obtained in predicting mortality of the patients. In all the cases the deviance ratio of the model improved, except for the group 50–64 where it remained the same.

Fig. 5
figure 5

Distribution of age in the adult cohort. A Histogram of age distribution in the high-resolution CENTER-TBI cohort on top. B Age ranges groups and their relative numerosity, in the lower part of the figure

To investigate further the correlation between brain–heart crosstalks and mortality, the point biserial correlation coefficients between mortality and the brain–heart crosstalks measures were also obtained. The results, summarized in Fig. 6, shows the point biserial correlation coefficients between the network measures and mortality, exhibiting negative sign in almost all of them. Moreover, in Fig. 6 we present the boxplots of the three newly introduced brain–heart crosstalks measures for survivors and nonsurvivors in the four age groups. In Table 4 we report the p values of the Kruskal–Wallis test performed for the brain–heart crosstalks measures of the four age groups, comparing survivors and nonsurvivors.

Fig. 6
figure 6

Table with the point biserial correlations tests between brain–heart crosstalks measures and mortality for the four age groups (top). In the bottom of the figure, we present the boxplots for the \(ct_{np}\), \(\omega_{ct}\) and \(mi_{{ct{ }}}\) for the survivors and nonsurvivors in the four age groups defined. (\(ct_{np}\), brain–heart crosstalks; \(mi_{ct}\), mutual information; ωct, average edge overlap)

Table 4 Table showing the Kruskal–Wallis tests between survivors and nonsurvivors of each age group analyzed

Discussion

In this article, we investigated the relationship between mortality and events of interaction between brain and heart in patients with TBI, named brain–heart crosstalks.

Brain–heart crosstalks were introduced in our previous work [11]. This novel metric detects events of simultaneous increases of ICP and HR.

Brain–heart crosstalks measures were further extended, through a complex network modeling [13], which led to the computation of two measures for the ICP-HR system: average edge overlap (\(\omega_{ct} )\) and mutual information (\(mi_{ct}\)), which give indication of the behavior of the system during brain–heart crosstalks, as described more in details in the methods section.

From the current study, we could see that the results suggested an inverse relationship between brain–heart crosstalks measures and mortality. The point biserial correlation coefficient was always negative, for the three network measures and mortality. This was true for the case in which we computed it using the 226 patients and when segmenting the population in four age groups. In fact, for our mortality model we divided the population in age groups in the attempt of decreasing the heterogeneity brought by this clinical variable. Important differences in terms of comorbidities, type of lesions and outcome prospective are related to age, particularly in the face of the TBI population aging in western countries.

The first indication of a negative relationship between mortality and brain–heart crosstalk measure, was further confirmed by the Kruskal–Wallis tests between survivors and nonsurvivors distribution of those measures in the whole dataset. In addition, when segmenting the original population, at least one of the brain–heart crosstalks measures was statistically significantly different between survivors and nonsurvivors for each age group, except for the age group 50–64. However, for this age group in the logistic regression model, the coefficient associated to average edge overlap, and the normalized numbers of crosstalks \({ct}_{np}\), was statistically significant and with a negative sign. Logistic elastic-net prediction models appeared to have a good fitting to the data, as the percentage of deviance ratio explained is sufficiently high. This was true except for the age group 30–49. These age-related findings are difficult to explain based on the dataset studied and warrant further investigation in a larger group.

Another important aspect to be considered is that the sign of the \({ct}_{np}\) coefficient was consistently negative for all the four age groups. Given these results, we can confirm that the higher the number of brain–heart crosstalks, the lower is the probability of mortality. This result is confirmed independently of the age group to which the individuals belong to.

Similar behavior was exhibited by the average edge overlap, while for the mutual information two age groups showed positive relationship between mortality and the indicator. This evidence will need further investigation in the future.

Overall, our results suggest that the newly introduced variables of brain–heart crosstalks, might be considered as biomarkers of a healthy brain–heart interaction and therefore the assessment of these biomarkers could play a role in understanding the pathological mechanisms in TBI. The complexity of such a disease is indeed difficult to explain with simply measuring mean ICP values and multimodality physiological monitoring techniques have been developed to improve our understanding and are implemented in the clinical practice. In this perspective, considering brain–heart crosstalks along with other physiology measures in integrated protocols may be a promising approach.

For what concerns the other variables of the prediction models fitted, age came as a statistically significant element only for the two older age groups, but not for the two younger ones. This is in line with findings suggested by the literature, where there is evidence that older age might be a factor worsening the outcome of TBI, because of the comorbidities associated along with other factors [22]. Moreover, from the fitted models, as we can see in Table 3 for the older populations, almost all the variables have statistically significant coefficients, unlike younger population groups. This is consistent with the fact that the model for the two older population groups showed a better fitting than the younger population groups. Possibly, the higher prevalence of cardiologic preexisting disease may explain the relevance of the brain–heart interplay in the elderly.

As a final consideration our analysis shed light on the possibility of integrating brain–heart crosstalks measures into the existing TBI prognostic models. A very well-known prognostic model is the IMPACT score [23], where admission features as for instance age, motor score, pupils, hypoxia, hypotension, CT variables, glucose, and hemoglobin are used to evaluate patient’s condition at the time of hospitalization. Some work has been done for what concerns more complex prognostic models development, but to the best of our knowledge, none has ever taken into account variables expressing the interaction between brain and heart [24]. Therefore, our analysis opens the possibility for implementing the TBI prognostic models with brain–heart interaction measures, and hopefully increasing their precision.

Conclusions

We found an inverse relationship between the number of brain–heart crosstalks events and mortality in patients with TBI. We also applied two further indicators of brain–heart crosstalks from a complex network modeling of the ICP-HR system, which showed similar inverse behavior with mortality. The inverse relationship was also confirmed when fitting the input features into a mortality prediction model.

We are currently investigating how brain–heart crosstalks relate to other physiological variables (e.g., autonomic nervous activities) or clinical variables (treatments a patient receives, airways suctioning, presence of intracranial hypertension or impaired autoregulation). In so far as limitations are concerned, we believe that the division into four age groups was necessary to avoid confounding effects due to patients’ wide age differences. However, we are aware that this reduces the number of observations in each group when fitting the related models thus potentially resulting in fitting bias. Therefore, we plan to extend the study to a larger cohort of patients to explore further the model proposed. Currently we have started a preliminary analysis concerning the evidence of a causality existing between HR and ICP during brain–heart crosstalks, which could lead to another important contribution toward the use of brain–heart crosstalks in clinical practice. For future work, it would also be interesting understanding the possible relationship between different events (such as only raised ICP) with mortality.

Furthermore, we would like to extend the use of the crosstalks as a real-time digital biomarker in patients with TBI, through the use of an ICM+-Python plugin interface [17] running at the bedside. This would ultimately facilitate the use our new metric in the individualized management of this patients. We are in the process of implementing in Python the code developed initially in R.