Research Article

Topological network analysis of patient similarity for precision management of acute blood pressure in spinal cord injury

Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, United States
Division of Physical Medicine and Rehabilitation, Department of Orthopaedics and Rehabilitation, University of New Mexico School of Medicine, United States
San Francisco Veterans Affairs Healthcare System, United States
Computational Research Division, Lawrence Berkeley National Laboratory, United States
Lawrence Berkeley National Laboratory, United States
Rehabilitation Research Center, Santa Clara Valley Medical Center, United States
Department of Psychiatry and Behavioral Science, and University of Minnesota, United States
Institute for Health Informatics, University of Minnesota, United States
Department of Radiology and Biomedical Imaging, University of California, San Francisco, United States
Department of Emergency Medicine, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, United States
Department of Physical Medicine and Rehabilitation, Santa Clara Valley Medical Center, United States
Department of Neurosurgery, Stanford University, United States
Department of Anesthesia and Perioperative Care, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, United States

Nov 16, 2021

https://doi.org/10.7554/eLife.68015

Open access
Copyright information

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record updated: December 3, 2021 (Go to version)
Version of Record published: December 2, 2021 (This version)
Accepted Manuscript updated: November 17, 2021 (Go to version)
Accepted Manuscript published: November 16, 2021 (Go to version)
Accepted: October 23, 2021
Received: March 2, 2021

1. Of interest
A logic-incorporated gene regulatory network deciphers principles in cell fate decisions

Gang Xue, Xiaoyi Zhang ... Zhiyuan Li

Research Article Apr 23, 2024
Further reading

Abstract
Editor's evaluation
eLife digest
Introduction
Methods
Results
Discussion
Data availability
References
Article and author information
Metrics

Abstract

Background:

Predicting neurological recovery after spinal cord injury (SCI) is challenging. Using topological data analysis, we have previously shown that mean arterial pressure (MAP) during SCI surgery predicts long-term functional recovery in rodent models, motivating the present multicenter study in patients.

Methods:

Intra-operative monitoring records and neurological outcome data were extracted (n = 118 patients). We built a similarity network of patients from a low-dimensional space embedded using a non-linear algorithm, Isomap, and ensured topological extraction using persistent homology metrics. Confirmatory analysis was conducted through regression methods.

Results:

Network analysis suggested that time outside of an optimum MAP range (hypotension or hypertension) during surgery was associated with lower likelihood of neurological recovery at hospital discharge. Logistic and LASSO (least absolute shrinkage and selection operator) regression confirmed these findings, revealing an optimal MAP range of 76–[104-117] mmHg associated with neurological recovery.

Conclusions:

We show that deviation from this optimal MAP range during SCI surgery predicts lower probability of neurological recovery and suggest new targets for therapeutic intervention.

Funding:

NIH/NINDS: R01NS088475 (ARF); R01NS122888 (ARF); UH3NS106899 (ARF); Department of Veterans Affairs: 1I01RX002245 (ARF), I01RX002787 (ARF); Wings for Life Foundation (ATE, ARF); Craig H. Neilsen Foundation (ARF); and DOD: SC150198 (MSB); SC190233 (MSB).

Editor's evaluation

The major strengths of this paper are the use of a combination of relatively novel approaches/applications to the identification of important predictors for recovery after spinal cord surgery. These include various data reduction techniques such as dissimilarity matrices and a subject-centered topological network analysis to identify phenotypes.

https://doi.org/10.7554/eLife.68015.sa0

eLife digest

Spinal cord injury is a devastating condition that involves damage to the nerve fibers connecting the brain with the spinal cord, often leading to permanent changes in strength, sensation and body functions, and in severe cases paralysis. Scientists around the world work hard to find ways to treat or even repair spinal cord injuries but few patients with complete immediate paralysis recover fully.

Immediate paralysis is caused by direct damage to neurons and their extension in the spinal cord. Previous research has shown that blood pressure regulation may be key in saving these damaged neurons, as spinal cord injuries can break the communication between nerves that is involved in controlling blood pressure. This can lead to a vicious cycle of dysregulation of blood pressure and limit the supply of blood and oxygen to the damaged spinal cord tissue, exacerbating the death of spinal neurons. Management of blood pressure is therefore a key target for spinal cord injury care, but so far, the precise thresholds to enable neurons to recover are poorly understood.

To find out more, Torres-Espin, Haefeli et al. used machine learning software to analyze previously recorded blood pressure and heart rate data obtained from 118 patients that underwent spinal cord surgery after acute spinal cord injury.

The analyses revealed that patients who suffered from either low or high blood pressure during surgery had poorer prospects of recovery. Statistical models confirming these findings showed that the optimal blood pressure range to ensure recovery lies between 76 to 104-117 mmHg. Any deviation from this narrow window would dramatically worsen the ability to recover.

These findings suggests that dysregulated blood pressure during surgery affects to odds of recovery in patients with a spinal cord injury. Torres-Espin, Haefeli et al. provide specific information that could improve current clinical practice in trauma centers. In the future, such machine learning tools and models could help develop real-time models that could predict the likelihood of a patient’s recovery following spinal cord injury and related neurological conditions.

Introduction

Spinal cord injury (SCI) may result in motor, sensory, and autonomic dysfunction in various degrees, depending on the injury severity and location. In the USA, the incidence of SCI is around 18,000 new cases per year, with a total prevalence ranging from 250,000 to 368,000 cases (National Spinal Cor Injury Statistical Center, 2021). The dramatic life event of SCI imposes a high socioeconomic burden, with an estimated lifetime cost between $1.2 and $5.1 million per patient (National Spinal Cor Injury Statistical Center, 2021).

Prior retrospective observational single-center studies in humans suggest that lower post-surgery mean arterial pressure (MAP) predicts poor outcome (Cohn et al., 2010; Hawryluk et al., 2015; Saadeh et al., 2017; Ehsanian, 2020), which have resulted in clinical guidelines focused on avoidance of hypotension by maintaining MAP >85 mmHg during the first 7 days of injury (Aarabi et al., 2013). The rational for MAP augmentation to avoid hypotension is based on the hypothesis that decreased spinal cord prefusion leads to ischemia and additional tissue lost (Mautes et al., 2000; Soubeyrand et al., 2014). Importantly, experimentally raising MAP during acute SCI in animals by using vasopressors increases the risk for hemorrhage and consequent tissue damage (Soubeyrand et al., 2014; Streijger et al., 2018; Guha et al., 1987). In acute cervical patients with SCI, spinal cord hemorrhage correlates with poor prognosis for neurological recovery (Miyanji et al., 2007). These findings collectively suggest that hypo- and hypertension must be accounted for in MAP management, but there remains a gap in evidence from clinical studies that definitively informs MAP guidelines (Saadeh et al., 2017).

One of the challenges resulting in the lack of strong evidence for MAP management in patients with acute SCI is the heterogeneity of injury. The heterogeneity of SCIs results in data complexity that benefit from modern analytic tools. Using the machine intelligence approach of topological data analysis, we have previously shown that hypertension during SCI surgery (ultra-acute phase) predicts long-term functional recovery in rodent models (Nielson et al., 2015). These prior cross-species findings motivated the present multicenter study where we apply a data-driven workflow in patients with ultra-acute SCI surgical records from two different Level 1 trauma centers. By employing machine intelligence tools, we show that deviation from an optimal MAP range during surgery predicts lower likelihood of neurological recovery and suggest new targets for therapeutic intervention.

Methods

Retrospective data extraction and cohort selection

Operating room (OR) records from n = 94 patients (98 surgical records, 3 patients with multiple surgeries) from the Zuckerberg San Francisco General Hospital (ZSFG, from 2005 to 2013) and n = 33 patients (33 surgical records) from the Santa Clara Valley Medical Center (SCVMC, from 2013 to 2015) that underwent spinal surgery were collected retrospectively. For ZSFG, monitoring data was extracted from the values manually recorded by the anesthesiologist at 5 min intervals (Q5). For SCVMC, monitoring data was extracted from the Surgical Information Systems (SIS) (Alpharetta, GA) at 1 min intervals (Q1). Demographics and outcome variables were extracted from an existing retrospective registry. AIS (American Spinal Injury Association [ASIA] Impairment Scale) grade at admission (first complete AIS upon admission to the hospital before surgery) and discharge (latest complete AIS grade after surgery before discharge from hospital) were estimated using the available ISNCSCI exams (International Standards for Neurological Classification of SCI) and the neurosurgery, trauma surgery, emergency department, and physical medicine and rehabilitation physical exam notes. To ensure compatibility between centers on the estimated AIS grades, one independent physician conducted the estimates for all the patients in each center (SM for SCVMC and JT for ZSFG) and one independent ISNCSCI certified physician (WW) extracted the AIS grades for all the patients (across centers). In case of conflict between grades, both physicians established a consensus. From the total 131 surgical records, three records were excluded for not having monitoring recorded for both MAP and HR, three were excluded because surgeries were not related to SCI, and seven surgeries were excluded from three patients that were submitted to more than one surgery. The final cohort for exploratory topological data analysis included 118 patients with complete MAP and heart rate (HR) monitoring. For confirmatory regression analysis, 15 patients were excluded because AIS grade could not be extracted either at admission and/or discharge (Figure 2—figure supplement 1). AIS improvement was defined as an increase of at least one AIS grade from admission to discharge. The final list of extracted variables included: MAP and HR continuous monitoring, age, length of surgery in minutes, days from surgery to hospital discharge, estimated AIS grade at admission, estimated AIS grade at discharge and AIS improvement (‘yes’, ‘no’). All data was de-identified before pre-processing and analysis. Protocols for retrospective data extraction were approved by Institutional Research Board (IRB).

Cohort statistics

The differences in the AIS improvers/non-improvers population characteristics were tested at the univariate level using R (see software below). For continuous numerical variables (age, length of surgery, and days from surgery to discharge), the group mean differences were tested using unpaired Student’s t-test (two-sided test). For categorical variables (AIS admission, AIS at discharge, and dichotomized neurological level of injury [NLI]), their levels were compared using Fisher’s exact test (two-sided test). p-Values are presented in Table 1.

Table 1

Cohort demographics split by AIS (American Spinal Injury Association [ASIA] Impairment Scale) improvement.

	AIS improve. N/A(n = 15)	AIS improve. NO(n = 61)	AIS improve. YES(n = 42)	Univariate p-value
Age (years)				0.12
Mean (SD)	46.0 (17.6)	45.3 (19.1)	51.4 (19.7)
Median [min, max]	45.5 [19.0, 87.0]	47.0 [18.0, 82.0]	55.0 [18.0, 86.0]
Missing	1 (6.7%)	2 (3.3%)	1 (2.4%)
AIS admission				0.013
A	1 (6.7%)	33 (54.1%)	18 (42.9%)
B	0 (0%)	5 (8.2%)	8 (19.0%)
C	0 (0%)	5 (8.2%)	11 (26.2%)
D	0 (0%)	14 (23.0%)	5 (11.9%)
E	0 (0%)	4 (6.6%)	0 (0%)
Missing	14 (93.3%)	0 (0%)	0 (0%)
AIS discharge				<0.0001
A	0 (0%)	35 (57.4%)	0 (0%)
B	0 (0%)	5 (8.2%)	5 (11.9%)
C	1 (6.7%)	4 (6.6%)	15 (35.7%)
D	0 (0%)	14 (23.0%)	17 (40.5%)
E	1 (6.7%)	2 (3.3%)	5 (11.9%)
Missing	13 (86.7%)	1 (1.6%)	0 (0%)
Surgery duration (min)				0.66
Mean (SD)	433 (167)	392 (146)	407 (181)
Median [min, max]	432 [121, 725]	389 [120, 728]	343 [151, 950]
Missing	1 (6.7%)	2 (3.3%)	1 (2.4%)
Surgery to discharge (days)				0.33
Mean (SD)	9.50 (2.12)	18.8 (20.6)	23.4 (23.8)
Median [min, max]	9.50 [8.00, 11.0]	11.0 [1.00, 128]	14.5 [4.00, 120]
Missing	13 (86.7%)	4 (6.6%)	2 (4.8%)
Dichotomized neurological level of injury at admission				0.054
Cervical	2.00 (13.3%)	36 (59%)	33 (78.6%)
Non-cervical	13.00 (86.7%)	25 (41%)	9 (21.4%)

Topological network extraction and exploration of monitoring data workflow

Monitoring data pre-processing

Two datasets were generated, one for each center. To ensure compatibility, both datasets were harmonized. Given the difference in the sampling frequency of the monitoring data (Q5 for ZSFG and Q1 SCVMC) between centers and protocols for data collection, SCVMC monitoring data was first pre-processed. Briefly, electronic data was exported from the SIS SQL database, de-identified and imported into MATLAB version 2016b (Mathworks Inc, Natick, MA) for filtering. A custom MATLAB script generated by the SCVMC team implemented filtering criteria established by clinicians and researchers to correct for invalid data (e.g., motion artifacts and injections). Thus, MAP values under 10 and above 200 mmHg as well as point-to-point changes greater than 40 mmHg were filtered, as these instances were found to represent data artifacts. This process was validated by comparing clinical curated data and the extracted data from the script with an accuracy of 99.1%. After filtering, SCVMC Q1 monitoring data was downsampled to Q5 by taking the average of five consecutive Q1 intervals for compatibility with ZSFG data. Given that the duration of the monitoring for each patient was different and the continuous time-series data is not aligned between patients (without time dependency on monitoring values), the empirical cumulative distribution function (CDF) for each time-series and each patient was computed. To account for the different scales between MAP and HR, a bin width for CDF was set as a 1% of the range of each measure, producing 100 CDF bins for MAP and 100 bins for HR (Figure 2—figure supplement 2). Additionally, the average MAP (aMAP) and HR (aHR) across time for each patient was calculated for posterior analysis.

Similarity between patients

The CDF bins for both MAP and HR were serialized in one vector per patient and the Euclidean distance ( $d (p, q) = \sqrt{\sum_{i = 1}^{n} {(p_{i} - q_{i})}^{2}}$ , where p and q are two patients’ CDF vectors and n is the number of CDF bins) calculated for each pair of patients’ vectors. The Euclidean distance matrix was then processed using dimensionality reduction.

Dimensionality reduction

Our goal for dimensionality reduction was to increase the signal-to-noise ratio by mapping to a lower number of dimensions (d) that contained the major information in the dataset. Dimensionality reduction was achieved by using the Isomap algorithm (Tenenbaum et al., 2000). Isomap is a non-linear dimensionality reduction method that uses multidimensional scaling (MDS) with geodesic distances instead of the Euclidean distances as the classic MDS does, and it has been suggested before for health data (Weng et al., 2005). Isomap performs three steps: (1) construct an NNG (near neighbor graph), (2) estimate the geodesic distances from the graph (shortest paths), (3) compute MDS embedding with the geodesic distances. The algorithm takes one parameter (k or e) to set the threshold for the NNG (we used k, the number of near neighbors for NNG). For selecting k, we considered the minimal k the smallest one that produced a connected NNG, in our case k = 3. Then two criteria were established for both k optimization and d selection, distance preservation (RV, residual variance) and topological persistent homology preservation as described (Rieck and Leitte, 2015; Paul and Chalup, 2017). We considered Isomap solutions for k = 3–7 (Figure 2—figure supplement 1). The RV was computed as (Tenenbaum et al., 2000): $1 - R^{2} (\hat{D_{m}}, D_{y})$ where $R$ is the standard linear correlation coefficient taken over all entries of ${\hat{D}}_{m}$ and $D_{y}$ . ${\hat{D}}_{m}$ is the input distance matrix that the algorithm is trying to estimate the real dimensions of (k-geodesic distance matrix for k-Isomap). $D_{y}$ is the Euclidean distance matrix of the low-dimensional solution. No major differences were observed in RV between the solutions for different k, except for the first dimension where RV increases as k increases. Isomap persistence diagrams were obtained using Vietoris-Rips filtration (Paul and Chalup, 2017) for ${\hat{D}}_{m}$ and $D_{y}$ for different d solutions (Figure 2—figure supplement 1). Then the topological zero-dimensional and one-dimensional Wasserstein (power = 2) distances (WD0 and WD1, respectively) were computed between ${\hat{D}}_{m}$ and $D_{y}$ . In persistent homology, dimension 0 measures zero-dimensional holes in the data (the connectivity of the datapoints, i.e., the number of connected sets) and topological dimension 1 measures one-dimensional holes, namely loops. We sought to select the solution (determine d and k) that minimized the WD0 and WD1, indicative of the optimal solution preserving the major topological information (Rieck and Leitte, 2015; Paul and Chalup, 2017). Given that k = 6 and k = 7 showed the lowest WD0 and WD1, we considered k = 6 as the final solution (Figure 2—figure supplement 1). A d = 4 (four dimensions kept in the k = 6 Isomap) was chosen for being the one at the ‘elbow’ of the RV, the one that minimized WD0 in k = 6 Isomap and presented a good compromised WD1.

Network analysis

A network from the k = 6 d = 4 Isomap solution was created for visual representation of the connectivity of patients (similarity) in the low-dimensional space. In this network, nodes represent each patient and edges the connection of two patients that are similar in the Isomap solution. An adjacency (whether two nodes are connected) matrix was obtained by computing a k-NNG for the low-dimensional space. The cutoff threshold for adjacency was set to the minimal k that produced a full-connected network. For network clustering, the walktrap algorithm was used (Pons and Latapy, 2005) as implemented in the igraph R package. Walktrap takes a single parameter, the number of random steps the algorithm uses to determine if nodes are in the same cluster or not. To select the optimal number of steps we computed walktrap solutions of a set of random steps (1–100) and chosen the first solution which maximized modularity (Q) as implemented in igraph R package (Clauset et al., 2004). In network analysis, modularity can be interpreted as the proportion of within cluster compared to the between clusters connectivity (edges). This solution was 47 random steps, producing a dendrogram tree of connectivity which maximal modularity cut the tree in 11 clusters (Figure 2—figure supplement 1). Then the network was contracted for visual representation of a cluster network of patients, where nodes represented the clusters and edges connected clusters that had at least one edge in the similarity network. Both the similarity and the cluster networks were used for exploratory network analysis and hypothesis generation by mapping patient features and visual inspection (Figure 2—figure supplement 1). We used the assortativity coefficient (Ar) to explore the possibility that the network was capturing the association between patients and the time of MAP out of a range (Figure 4; see time MAP out of range). The Ar was calculated using the igraph implementation in R (Newman, 2003) and it can be interpreted as the Pearson coefficient (−1–1) between nodes connectivity and value of a variable.

Regression analysis

Logistic regression

We first used a logistic regression to model the probability of predicting improvement by aMAP. Visual inspection of the plot (Figure 2) suggested a non-linear relationship between aMAP and the probability of improvement. Consequently, the following logistic models were considered ( $l$ being the log-odds or logit of the probability of improving): the null model with no predictors ( $l = β_{0}$ , the simple model ( $l = β_{0} + β_{1} x$ ), the two-grade polynomial model ( $l = β_{0} + β_{1} x + β_{2} x^{2}$ ), the three-grade polynomial model ( $l = β_{0} + β_{1} x + β_{2} x^{2} + β_{3} x^{3}$ ) and a natural spline model ( $l = β_{0} + f_{1} (x)$ )where $f_{1} (x)$ is the natural spline function with 2 or 3 degrees of freedom (df)). The natural cubic spline was chosen to relax the symmetric constraints of polynomial models given that the visual inspection of the data suggested an asymmetric aMAP range (Figure 2, distribution of aMAP of improvers is skewed to the left). The results of fitting these models and the likelihood comparison between them (by likelihood ratio test) are shown in Table 2. The best fitting model was the two-grade polynomial (Figure 3) and the natural spline (2df) with significant coefficients, confirming our hypothesis. These results were confirmed by leave-one-out cross-validation (LOOCV) (Table 2). To account for the potential confounding effect of AIS grade at admission (given differences between groups, Table 1), aHR, length of surgery (minutes), days from surgery to discharge and age, we fitted the quadratic model with those covariates and LOOCV (Table 3). Considering the independence of the predictors (small correlation coefficients between variables; Table 4), the results of the quadratic term being significant still holds for the covariate model.

Table 2

Logistic regression likelihood ratio test and leave-one-out cross-validation (LOOCV) error (n = 103 patients).

Model	AIC	Residual df	Residualdeviance	Deviance	p-Value	LOOCV error
Null model $(l = β_{0})$	141.26	102	139.26			0.246
Linear model $(l = β_{0} + β_{1} x)$	134.8	101	130.80	8.46(vs. null model)	0.0036**(vs. null model)	0.231
Quadratic model $(l = β_{0} + β_{1} x + β_{2} x^{2})$	128.48	100	122.48	8.32(vs. linear model)	0.0039**(vs. linear model)	0.210
Cubic model $(l = β_{0} + β_{1} x + β_{2} x^{2} + β_{3} x^{3})$	126.97	99	118.97	3.50(vs. quadratic model)	0.061(vs. quadratic model)	0.213
Natural Spline model (df = 2) $(l = β_{0} + f_{1} (x))$	128.29	100	122.29	8.50(vs. linear model)	0.0035**(vs. linear model)	0.210
Natural Spline model (df = 3) $(l = β_{0} + f_{1} (x))$	127.13	99	119.13	3.34(vs. quadratic model)	0.067(vs. quadratic model)	0.213

** p < 0.01.

Table 3

Evaluation of logistic regression (Wald test) and leave-one-out cross-validation (LOOCV) error.

Model: $l = β_{0} + β_{1} x_{1} + β_{2} x_{1}^{2}$ where $x_{1} :$ average MAP (n = 103 patients)
LOOCV: average observed accuracy = 0.66; average kappa statistic = 0.334
Predictor	Coef. estimate (logit)	Std. error	z-Value	p-Value
Intercept	$β_{0}$ = –0.55	0.242	–2.293	0.02183*
Average MAP ( $x_{1}$ )	$β_{1}$ = 8.62	2.944	2.931	0.00338**
Average MAP ( $x_{1}^{2}$ )	$β_{2}$ = –7.601	3.039	–2.501	0.0123*

*p < 0.05; **p < 0.01.

Table 4

Evaluation of logistic regression with covariates (Wald test) and leave-one-out cross-validation (LOOCV).

Model: $l = β_{0} + β_{11} x_{1} + β_{12} x_{1}^{2} + β_{2} x_{2} + β_{3} x_{3} + β_{4} x_{4} + β_{5} x_{5} + β_{6} x_{6} + β_{7} x_{7} + β_{8} x_{8} + β_{9} x_{9}$ , where $x_{1} :$ average MAP; $x_{2}$ : average HR; $x_{3} :$ length of surgery (min); $x_{4} :$ days to AIS discharge (days); $x_{5} :$ age; $x_{6} :$ AIS admission D (‘yes’,’no’); $x_{7} :$ AIS admission C (‘yes’,’no’); $x_{8} :$ AIS admission B (‘yes’,’no’); $x_{9} :$ AIS admission A (‘yes’,’no’); (AIS admission E was set as the reference level for AIS admission variable and is part of the intercept) (final n = 93)
LOOCV: average observed accuracy = 0.688; average kappa statistic = 0.362
Predictor		Coef. estimate (logit)		Std. error		z-Value		p-Value
Intercept		$β_{0}$ = –1.530		121.8		–0.013		0.99
Average MAP ( $x_{1}$ )		$β_{11}$ = 7.398		3.112		2.377		0.017*
Average MAP ( $x_{1}^{2}$ )		$β_{12}$ = –8.053		3.530		–2.281		0.022*
Average HR ( $x_{2}$ )		$β_{2}$ = –2.087		0.0245		–0.851		0.394
Length of surgery ( $x_{3}$ )		$β_{3}$ = 0.0011		0.0015		0.728		0.466
Days to AIS discharge ( $x_{4}$ )		$β_{4}$ = 0.0037		0.0109		0.344		0.730
Age ( $x_{5}$ )		$β_{5}$ = 0.0082		0.013		0.634		0.526
AIS admission D ( $x_{6}$ )		$β_{6}$ = 1.454		1.218		0.012		0.990
AIS admission C ( $x_{7}$ )		$β_{7}$ = 1.645		1.218		0.014		0.989
AIS admission B ( $x_{8}$ )		$β_{8}$ = 1.585		1.218		0.013		0.989
AIS admission A ( $x_{9}$ )		$β_{9}$ = 1.527		1.218		0.013		0.990
Correlation matrix (Spearman)
	Average MAP	AverageHR	Length of surgery		Days to AIS discharge		Age		AIS admission
Average MAP	1
Average HR	–0.126	1
Length of surgery	–0.152	0.101	1
Days to AIS discharge	0.088	–0.059	0.165		1
Age	0.006	–0.245	0.011		0.022		1
AIS admission	0.024	0.003	–0.01		0.258		–0.13		1

*p < 0.05.

Time out of MAP range

We sought to determine a range of MAP in which time outside that range might predict improvement. To consider the time at which MAP was outside a range, we performed an increasing window of MAP for either a symmetric range or an asymmetric one. For the symmetric range, a 1 mmHg range increment at each site of the center (90 mmHg, the mean MAP for improvers) was created. For the asymmetric range, the lower limit was fixed at 76 mmHg and the upper limit was incremented 1 mmHg at the time. The time of MAP (in min) being outside each range was estimated for each patient.

LASSO regression

LASSO (least absolute shrinkage and selection operator) regression (Tibshirani, 1996) was used for selecting a single range of MAP (see time MAP out of range) predictor of the logistic model: $l = β_{0} + \sum_{j = 1}^{p} β_{j} x_{j}$ , where $x_{j}$ is the jth MAP range. LASSO takes as parameter lambda that sets the amount of shrinkage or regularization (using L1-norm penalty). LOOCV was used to determine the lambda that shrunk the models to one predictor or MAP range. The one-predictor solutions (MAP range of 76–104 for symmetric range model, and 76–117 for the asymmetric range model) were used as the solo predictor of AIS improvers in a logistic regression with LOOCV (see above). It is important to note that given the high multicollinearity in the range data, the Q5 time estimation and the low sample size, the LASSO solution should be taken with caution and as an indicator of the MAP range hypothesis rather than a hard rule for medical decision making.

Prediction modeling

Logistic regression (see above) was used to build prediction models for three binary outcome metrics: AIS improvement of at least one grade from admission to discharge, whether patient was AIS grade A at discharge, or whether the patient was AIS grade D at discharge. For each one of the classification tasks, the following predictors were considered: quadratic aMAP (both linear and quadratic terms), aHR, length of surgery (min), days from surgery to discharge, age, AIS grade at admission, and dichotomized NLI. We performed model selection (a.k.a. feature selection) through an exhaustive search of all potential combinations of at least one of the predictors using the glmulti R package (Calcagno, 2020). The most parsimonious models were selected to be the one minimizing the small-sample corrected Akaike information criteria (AIC) for each task. We then investigated the performance of each one of the most parsimonious models using LOOCV and adjusting the classification threshold to balance prediction sensitivity and specificity. Briefly, each model was trained n (patient) times with an n−1 training sample and tested the performance in the remaining sample. A vector of n probabilities of predictions was then used to measure the LOOCV model performance. Model fitting and prediction performance were conducted using the caret R package (Kuhn et al., 2019). Receiving operating curves (ROC) and area under the curve (AUC) for the LOOCV prediction were obtained using the ROCR R package (Sing et al., 2005).

Software

All data wrangling, processing, visualization, and analysis was performed using the R programming language (R version 3.5.1) (R core team, 2019) and RStudio (RStudio version 1.2.1335) (Team, 2018) in Windows 10 operating system, with the exception of the Q1 OR measures form SCVMC that were preprocessed in MATLAB before downsampling to Q5 in R. The most relevant R functions and packages (beyond the installed with R) used and the references for each function/package and methods are reported in the following table. For more details, see the source code available.

R packages used

Package	Version	Usage	Reference
igraph	1.2.4.1	Building, manage and analyze networks	Csárdi and Nepusz, 2006
dplyr	0.8.3	Data cleaning and wrangling	Wickham et al., 2018
ggplot2	3.2.1	Data visualization and plotting	Wickham, 2016
vegan	2.5–5	For Isomap implementation	Jari et al., 2019
RColorBrewer	1.1–2	To control and create colormaps	Neuwirth, 2014
TDAstats	0.4.0	Utilities for topology data analysis for persistent homology	Wadhwa et al., 2018
cccd	1.5	For generating NNGs	Marchette, 2015
table1	1.1	Generates table of demographics	Rich, 2018
glmnet	2.0–18	For fitting LASSO	Friedman et al., 2010
glmnetUtils	1.1.2	For fitting LASSO	Ooi, 2019
caret	6.0–84	To fit logistic regression with leave-one-out cross-validation	Kuhn et al., 2019
splines	3.5.1	To fit the spline models	R core team, 2019
VisNetwork	2.0.7	Visualization suit for network graphs using the vis.js JavaScript library	Almende, 2019
stats	3.5.1	Fit generalized linear models	R core team, 2019
glmulti	1.0.8	For model search	Calcagno, 2020
ROCR	1.0–11	For ROC visualization and performance	Sing et al., 2005
reshape2	1.4.3	From wide to long view dataframe formatting	Wickham, 2007

Data and code availability

The final de-identified datasets for analysis are deposited and accessible at the Open Data Commons for SCI (odc-sci.org, RRID:SCR_016673) under DOIs 10.34945/F5R59J and 10.34945/F5MG68. The R code to run all the analysis present in this publication, including visualizations, is available as supplementary material. The code would reproduce the entire analysis and plots when run using the same versions of R, RStudio, and packages specified in this publication. Otherwise results might change.

Results

Exploratory network analysis suggests an upper and lower limit of intra-operative MAP for recovery

Intra-operative monitoring records (MAP, HR) and neurological outcome data were extracted and curated from two Level 1 trauma centers. A final cohort of 118 patients was included (Figure 1a and Table 1). The cohort represents a varied dataset of intra-operative MAP and HR patterns and respective aMAP across time in surgery and aHR across time in surgery values (Figure 1b–c). Using a machine intelligence analytical pipeline (Figure 2a), we extracted a similarity network of patients (Pai and Bader, 2018; Parimbelli et al., 2018) from a low-dimensional space embedded using a non-linear algorithm, Isomap (Tenenbaum et al., 2000), on a distance matrix derived from the MAP and HR records and then performed topological network extraction using persistent homology metrics (Rieck and Leitte, 2015; Figure 2a and Figure 2—figure supplement 1). The results of this dimensionality reduction suggested that four dimensions are enough to capture most of the variance and the topological structures of the original data (Figure 2—figure supplement 1c-e). Clustering the network of patients through a random-walk algorithm, Walktrap (Pons and Latapy, 2005), revealed 11 different clusters where patients were regarded to share intra-operative hemodynamic phenotypes (Figure 2 and Figure 2—figure supplement 1f-h). Importantly, this workflow was unsupervised: only the OR hemodynamic time-series was used to derive patient clustering, and therefore any association captured by the network must be dependent on hemodynamic patterns. Exploratory network analysis showed a gradient distribution of patients by their aMAP (Figure 2b–d) and aHR (Figure 2e–g), confirming that the network captured a valid representation of the raw high-dimensional dataset. We then investigated the association of the clusters to patient recovery as defined by whether the patient improved at least one AIS grade A–D (Roberts et al., 2017) between time before surgery and time of discharge from the hospital. Mapping the proportion of patients with AIS improvement onto the similarity network (Figure 2h–j) revealed that patients with recovery localized to clusters associated with a middle range of MAP (Figure 2k). Those clusters also showed a higher proportion of less severe AIS grades at discharge (AIS C, D, and E) than other clusters (Figure 2—figure supplement 2). In contrast, clusters of patients showing an extreme variation of MAP were highly enriched with patients with no AIS recovery and patients with more severe AIS grades at discharge (AIS A and B, Figure 2—figure supplement 2). This analysis suggested that there is a limited range of MAP during surgery associated with neurological recovery.

Figure 1

Download asset Open asset

High-frequency monitoring operating room (OR) data.

Flowchart of retrospective study and cohort selection criteria (a). A final cohort of 118 patients were identified and values of mean arterial pressure (MAP) (b) and heart rate (HR), (c) by time (bins of 5 min; Q5) retrospectively extracted from patients’ records. Colormaps represent the MAP (mmHg; green marks normotensive MAP, while blue and red marks hypotension and hypertension, respectively) and HR (beats per min, bpm; dark yellow lowest to purple highest) at each Q5 time, depicting the temporal fluctuation of each measure for each patient (row). The average MAP (aMAP, right plot in b) and average heart rate (aHR, right plot in c) were computed.

Figure 2 with 2 supplements see all

Download asset Open asset

Topological network analysis of intra-operative monitoring.

Intra-operative mean arterial pressure (MAP) and heart rate (HR) sampled every 5 min (Q5) were curated, processed, and formatted in a unique data matrix (a) (Figure 2—figure supplement 1). The similarity matrix between patients was computed and a four-dimensional subspace extracted using Isomap (Figure 2—figure supplement 1). A network was constructed where nodes represent patients and edges the connection of pairs of patients under a specified threshold of similarity (see Methods). The network was clustered and collapsed (Figure 2—figure supplement 1) by using the walktrap algorithm conveying in 11 clusters. These networks captured both the average MAP (aMAP) (**b–d**) and the average HR (aHR) (**e–g**) in a gradient fashion. Similarly, at least one AIS grade gain at discharge (‘yes’, ‘no’) was mapped over the network (**h–j**, gray: 15 AIS grades could not be extracted). Clusters of higher proportion of patients with recovery had an aMAP in a middle range, while clusters with higher proportion of patients without recovery presented extreme aMAPs (k). The mean cluster aHR showed a less apparent relationship with the proportion of AIS improvers (l).

MAP has a non-linear association with probability of recovery

The exploratory network analysis revealed that clusters with higher proportion of patients that increased AIS of at least one grade were associated with having a middle range aMAP (Figure 2 and Figure 2—figure supplements 1 and 2) and that clusters of patients with aMAP on the extremes contained fewer improvers. We hypothesized that there might be a non-linear relationship between intra-operative MAP and the probability of AIS grade improvement. To confirm this hypothesis, logistic regression models with LOOCV were used (Figure 3, Table 2). We fitted a null model (no predictors) as well as linear, polynomial, and cubic models of aMAP (Figure 3, Table 2) to test the non-linearity of the hypothesis. The linear model showed a significant improvement over the null model with a positive association, suggesting that the higher the aMAP, the higher the probability of AIS grade improvement. However, polynomial logistic regression demonstrated a significant quadratic fit (Table 2) with lower LOOCV error than the linear model, indicating that a quadratic form of aMAP better predicts the probability of improvement. Notably, the cubic model did not show significant improvement over the quadratic one. Exploratory network analysis suggested an asymmetrical function of AIS improvement with respect to aMAP (Figure 2k); we therefore also tested spline models to relax the symmetry of polynomial models. A spline model of degree 2 resulted in a significant fit over the linear model (Table 2) while a spline model of degree 3 resulted in a similar fit as compared to the cubic model. There was no evidence from which to choose between the spline model of degree 2 and the quadratic model. Accordingly, we did not pursue the asymmetric model further, although we explore an asymmetric MAP range below.

Figure 3

Download asset Open asset

Non-linear relationship of average mean arterial pressure (aMAP) with the probability of improving at least one AIS grade.

Logistic regression models were fitted to study the potential non-linearity of the aMAP predictor as suggested by the exploratory analysis. Six different models were studied: naïve, linear, quadratic, cubic, spline of degree 2, and spline of degree 3. The estimated leave-one-out cross-validation (LOOCV) error for each model showed that both the quadratic and the spline of degree 2 have the minimal cross-validation error (a). This suggests that the linear model did not capture all the potential explainable variance of the response variable by aMAP, while the cubic and spline of degree 3 were probably overfitting the model (Table 2). (b) shows the logit function (blue line) and standard error (gray ribbon) of the quadratic model over the fitted values (points).

Factors influencing MAP association with recovery

We sought to explore additional patient characteristics that might explain or affect MAP association with recovery. To test whether other factors could be responsible for the observed non-linear association, we first compared the quadratic model with aMAP as a predictor alone, a model that also includes several covariates (aHR, length of surgery, days from surgery to discharge, age, and AIS grade at admission), and a model with only the covariates. The significance of the quadratic fit holds even after accounting for the covariates (Table 3, 4), and none of the terms in the covariates-only model had a significant coefficient (Table 5). These results indicate that even in the presence of the other factors, aMAP is still non-linearly associated with AIS grade improvement at discharge.

Table 5

Evaluation of logistic regression covariates only (Wald test) and leave-one-out cross-validation (LOOCV).

Model: $l = β_{0} + β_{2} x_{2} + β_{3} x_{3} + β_{4} x_{4} + β_{5} x_{5} + β_{6} x_{6} + β_{7} x_{7} + β_{8} x_{8} + β_{9} x_{9}$ , where $x_{2}$ : average HR; $x_{3} :$ length of surgery (min); $x_{4} :$ days to AIS discharge (days); $x_{5} :$ age; $x_{6} :$ AIS admission D (‘yes’,’no’); $x_{7} :$ AIS admission C (‘yes’,’no’); $x_{8} :$ AIS admission B (‘yes’,’no’); $x_{9} :$ AIS admission A (‘yes’,’no’); (AIS admission E was set as the reference level for AIS admission variable and is part of the intercept) (final n = 93)
LOOCV: average observed accuracy = 0.612; average kappa statistic = 0.17
Predictor	Coef. estimate (logit)	Std. error	z-Value	p-Value
Intercept	$β_{0}$ = –1.585	138.2	–0.011	0.991
Average HR ( $x_{2}$ )	$β_{2}$ = –0.0209	0.029	–0.911	0.362
Length of surgery ( $x_{3}$ )	$β_{3}$ = 0.0016	0.00141	1.151	0.250
Days to AIS discharge ( $x_{4}$ )	$β_{4}$ = 0.0105	0.0106	0.993	0.320
Age ( $x_{5}$ )	$β_{5}$ = 0.0052	0.012	0.424	0.672
AIS admission D ( $x_{6}$ )	$β_{6}$ = 1.511	1.382	0.011	0.991
AIS admission C ( $x_{7}$ )	$β_{7}$ = 1.715	1.382	0.012	0.991
AIS admission B ( $x_{8}$ )	$β_{8}$ = 1.643	1.382	0.012	0.990
AIS admission A ( $x_{9}$ )	$β_{9}$ = 1.574	1.382	0.011	0.991

Patients with more severe injuries are more likely to suffer hemodynamic dysregulations (Lehmann et al., 1987). Hence, we studied whether the relationship of MAP and AIS improvement was maintained in the subcohort of patients with an AIS grade of A at admission. We first filtered the data for the subcohort and then fitted a full model as above but without the AIS grade at admission covariate. The resulting model showed the linear aMAP coefficient to be significant and the quadratic term close to significant (p = 0.14) with the second biggest coefficient (Table 6). A likelihood ratio test between a linear model with covariates and a quadratic model with covariates resulted in p-values = 0.07. On the other hand, in the full model with covariates fitted to the entire cohort, none of the AIS grades at admission had significant coefficients, which suggested that the non-linear relationship of MAP with neurological recovery was sustained across injury severity in that model. This apparent divergence in results might be explained by the reduction in power for the AIS A cohort model.

Table 6

Evaluation of logistic regression in American Spinal Injury Association (ASIA) Impairment Scale (AIS) A at admission cohort (Wald test) and leave-one-out cross-validation (LOOCV).

Model: $l = β_{0} + β_{11} x_{1} + β_{12} x_{1}^{2} + β_{2} x_{2} + β_{3} x_{3} + β_{4} x_{4} + β_{5} x_{5}$ , where $x_{1} :$ average MAP; $x_{2}$ : average HR; $x_{3} :$ length of surgery (min); $x_{4} :$ days to AIS discharge (days); $x_{5} :$ age (final n = 51)
LOOCV: average observed accuracy = 0.63; average kappa statistic = 0.197
Predictor	Coef. estimate (logit)	Std. error	z-Value	p-Value
Intercept	$β_{0}$ = –0.931	3.433	–0.271	0.786
Average MAP ( $x_{1}$ )	$β_{11}$ = 10.79	5.014	2.153	0.031*
Average MAP ( $x_{1}^{2}$ )	$β_{12}$ = –6.73	4.591	–1.468	0.142
Average HR ( $x_{2}$ )	$β_{2}$ = –0.016	0.035	–0.468	0.639
Length of surgery ( $x_{3}$ )	$β_{3}$ = 0.0039	0.0026	1.504	0.132
Days to AIS discharge ( $x_{4}$ )	$β_{4}$ = 0.0067	0.014	0.477	0.633
Age ( $x_{5}$ )	$β_{5}$ = –0.012	0.020	–0.599	0.549

*p < 0.05.

Next, given that the level of the cord injury can be related to different degrees of hemodynamic dysregulation (Lehmann et al., 1987), we studied the effect of the NLI at admission on the association of MAP and patient recovery. Our cohort was very heterogeneous on the NLI, with most patients having cervical injuries and the rest distributed along the mid and lower segments of the cord (Table 7). Thus, we divided the population into two categories: cervical and non-cervical patients. Running the same full model with just the cervical patients resulted in similar results as compared to the full model on the entire cohort, maintaining the quadratic aMAP significance (Table 8). In the non-cervical cohort, we did not find a significant association of the quadratic aMAP to recovery (Table 9). We then performed additional analyses to determine whether this difference in aMAP relationship to recovery between cervical and non-cervical patients was due to differences in the likelihood of recovery between the two NLI populations. A univariate analysis suggested that the proportion of improvers and not improvers in the cervical and non-cervical population were marginally different (Table 1). Moreover, a logistic regression predicting AIS grade improvement by NLI categorization indicated that non-cervical patients were significantly less likely to recover ( $β$ = –0.93, p = 0.041). While these results suggest that a quadratic aMAP is important for predicting AIS grade recovery in cervical patients, the lack of significant results in the non-cervical patients must be interpreted with caution due to the reduced number of cases, the heterogeneous distribution, and the low number of improvers in the group.

Table 7

Neurological level of injury cases.

	Cervical (n = 71)							Non-cervical (n = 32)
NLI	C2	C3	C4	C5	C6	C7	C8	T2	T3	T4	T5	T6	T7	T8	T9	T10	T11	T12	S1	S5
Cases	3	3	24	28	4	8	1	1	3	3	1	2	1	3	1	3	2	4	6	2

Table 8

Evaluation of logistic regression in Cervical cohort (Wald test) and leave-one-out cross-validation (LOOCV).

Model: $l = β_{0} + β_{11} x_{1} + β_{12} x_{1}^{2} + β_{2} x_{2} + β_{3} x_{3} + β_{4} x_{4} + β_{5} x_{5} + β_{6} x_{6} + β_{7} x_{7} + β_{8} x_{8}$ , where $x_{1} :$ average MAP; $x_{2}$ : average HR; $x_{3} :$ length of surgery (min); $x_{4} :$ days to AIS discharge (days); $x_{5} :$ age; $x_{6} :$ AIS admission D (‘yes’,’no’); $x_{7} :$ AIS admission C (‘yes’,’no’); $x_{8} :$ AIS admission B (‘yes’,’no’); (AIS admission A was set as the reference level for AIS admission variable and is part of the intercept, no AIS admission E was present in this cohort) (final n = 93)
LOOCV: average observed accuracy = 0.688; average kappa statistic = 0.362
Predictor	Coef. estimate (logit)	Std. error	z-Value	p-Value
Intercept	$β_{0}$ = 2.747	3.018	0.91	0.363
Average MAP ( $x_{1}$ )	$β_{11}$ = 7.594	3.056	2.485	0.013*
Average MAP ( $x_{1}^{2}$ )	$β_{12}$ = –7.528	3.358	–2.242	0.025*
Average HR ( $x_{2}$ )	$β_{2}$ = –0.055	0.034	–1.608	0.108
Length of surgery ( $x_{3}$ )	$β_{3}$ = 0.0014	0.0019	0.720	0.472
Days to AIS discharge ( $x_{4}$ )	$β_{4}$ = 0.0022	0.012	0.182	0.855
Age ( $x_{5}$ )	$β_{5}$ = 0.0079	0.016	0.482	0.630
AIS admission D ( $x_{6}$ )	$β_{6}$ = –0.747	0.87	–0.840	0.730
AIS admission C ( $x_{7}$ )	$β_{7}$ = 0.745	0.80	0.925	0.355
AIS admission B ( $x_{8}$ )	$β_{8}$ = 0.301	0.88	0.346	0.401

*p < 0.05.

Table 9

Evaluation of logistic regression in non-cervical cohort only (Wald test) and leave-one-out cross-validation (LOOCV).

Model: $l = β_{0} + β_{11} x_{1} + β_{12} x_{1}^{2} + β_{2} x_{2} + β_{3} x_{3} + β_{4} x_{4} + β_{5} x_{5} + β_{6} x_{6} + β_{7} x_{7} + β_{8} x_{8} + β_{9} x_{9}$ , where $x_{1} :$ average MAP; $x_{2}$ : average HR; $x_{3} :$ length of surgery (min); $x_{4} :$ days to AIS discharge (days); $x_{5} :$ age; $x_{6} :$ AIS admission D (‘yes’,’no’); $x_{7} :$ AIS admission C (‘yes’,’no’); $x_{8} :$ AIS admission B (‘yes’,’no’); $x_{9} :$ AIS admission A (‘yes’,’no’); (AIS admission E was set as the reference level for AIS admission variable and is part of the intercept) (final n = 93)
LOOCV: average observed accuracy = 0.688; average kappa statistic = 0.362
Predictor	Coef. estimate (logit)	Std. error	z-Value	p-Value
Intercept	$β_{0}$ = –1.883	352.4	–0.005	0.996
Average MAP ( $x_{1}$ )	$β_{11}$ = –0.206	4.713	–0.044	0.965
Average MAP ( $x_{1}^{2}$ )	$β_{12}$ = –8.064	7.643	–1.055	0.291
Average HR ( $x_{2}$ )	$β_{2}$ = –0.0002	0.0649	0.004	0.997
Length of surgery ( $x_{3}$ )	$β_{3}$ = 0.0018	0.0054	0.336	0.737
Days to AIS discharge ( $x_{4}$ )	$β_{4}$ = 0.076	0.0613	1.240	0.215
Age ( $x_{5}$ )	$β_{5}$ = –0.0047	0.051	–0.921	0.357
AIS admission D ( $x_{6}$ )	$β_{6}$ = 1.727	3.524	0.005	0.996
AIS admission C ( $x_{7}$ )	$β_{7}$ = 3.557	5.782	0.005	0.996
AIS admission B ( $x_{8}$ )	$β_{8}$ = 1.738	3.524	0.005	0.995
AIS admission A ( $x_{9}$ )	$β_{9}$ = 1.686	3.524	0.005	0.996

Finally, we sought to determine whether the probability of recovery associated to MAP could be influenced by the time the patient is in the hospital. For that, we break down the potential causal pathway between MAP dysregulation, AIS improvement, and days from surgery to discharge. We first fitted a logistic regression model with AIS improvement as response and days to discharge as the only predictor. This resulted in a non-significant p-value of p = 0.32, suggesting that days to discharge does not associate with probability of improvement. Second, we fitted a linear model with days to discharge as a response and quadratic aMAP (both linear and quadratic terms) as predictors. This resulted as a significant coefficient of the quadratic term (p = 0.047), although the model was not significant (p = 0.13 for the F statistic) and the adjusted R² was small (0.0217). We also investigated whether days to discharge interacts with MAP and quadratic MAP to predict AIS improvement, with no significant results on the interaction (interaction days to discharge with aMAP: linear term p = 0.61; quadratic term p = 0.91). These suggest that these two factors do not moderate each other. Finally, eliminating days to discharge from the full covariate model predicting AIS improvement does not have a major effect on the model fit. A likelihood ratio test between both models shows a non-significant change in variance explained (p = 0.729) with a deviance difference of ~0.1%. All together indicates that the non-linear relationship between aMAP and AIS improvement is independent of the days from surgery to discharge.

Intra-operative MAP range from 76-[104-117] mmHg predicts recovery

Since aMAP can obscure episodes of high deviation from average (Hawryluk et al., 2015) and has a non-linear relationship with recovery, we hypothesized that there might be a range of intra-operative MAP that better predicts AIS grade improvers. To test this hypothesis, we analyzed the amount of time MAP was out of a specific range (Figure 4). Since our modeling suggested both a symmetric and asymmetric range, we performed two different analyses. First, starting at a MAP of 90 mmHg, we implemented an algorithm to iteratively expand the MAP range symmetrically 1 mmHg higher and lower and calculate the time MAP was outside the range (Figure 4a). Exploratory analysis of the similarity network indicated a high association between the time out of a MAP range of 73–107 mmHg with the topological distribution of patients (Figure 4b and Figure 4—figure supplement 1). To validate this range and the associated lower and upper MAP thresholds, we used a logistic model with LASSO regularization with the predictors being the time outside of each MAP range as in Figure 4b. This allowed us to systematically reduce the number of relevant predictors until only one remained (non-zero coefficient). Interestingly, the logistic LASSO regression with LOOCV revealed that a MAP range from 76 to 104 mmHg was optimal in our dataset since it produced the most reproducible prediction of recovery (average LOOCV prediction accuracy of 61.16%; Figure 4c and Figure 2—figure supplement 2, Table 9). Next, we studied the possibility of an asymmetric range by fixing the lower limit to 76 mmHg and increasing the upper limit by 1 mmHg at the time (Figure 4d). The association of the patient distribution in the network plateau at a range of 76–116 mmHg (Figure 4e and Figure 4—figure supplement 1) and the logistic LASSO found the range 76–117 mmHg be the most predictive of recovery (average cross-validation prediction accuracy of 57.28%; Figure 4f and Figure 4—figure supplement 2a, Table 10). While both the exploratory analysis and the logistic LASSO produced similar ranges, the later analysis is performed through statistical modeling rather than descriptive associations, and therefore we further discuss the results of the LASSO.

Figure 4 with 3 supplements see all

Download asset Open asset

Range of mean arterial pressure (MAP).

To find the optimal MAP range, a moving MAP range was computed and the time of MAP outside range calculated (a and d, example of the same patient for symmetric and asymmetric map range, respectively). Calculating the assortativity coefficient (Ar) of the network for each range revealed that the distribution of patients in the network was most highly associated with the range 73–107 mmHg for the symmetric range (b, Figure 4—figure supplement 1), and 76–116 mmHg for the asymmetric range study of upper limit threshold (e, Figure 4—figure supplement 1). A logistic LASSO (least absolute shrinkage and selection operator) regression (Figure 4—figure supplement 2a, Figure 4—figure supplement 3) was used as a confirmatory analysis and to obtain the MAP range that most highly predicts AIS grade recovery. For the symmetric range, the time of MAP outside the 76–104 mmHg (c, Table 10) was found to be the ‘last-standing’ predictor during LASSO regularization (Figure 4—figure supplement 2a), suggesting that greater duration outside this range is associated with lower probability of neurological recovery. In the case of the asymmetric range (Figure 4—figure supplement 3), the last non-zero coefficient was for the range 76–117 mmHg (f, Table 11).

Table 10

Least absolute shrinkage and selection operator (LASSO) solution and logistic regression of most predictive symmetric range with leave-one-out cross-validation (LOOCV).

Model: $l = β_{0} + β_{1} x_{1}$ , where $x_{1} :$ time of MAP outside range 76–104 mmHg (n = 103 patients)
LOOCV: average observed accuracy = 0.61; average kappa statistic = 0.158
Predictor	Coef. estimate (logit)	Std. error	z-Value	p-Value
Intercept	$β_{0}$ = 0.368	0.333	1.103	0.269
Time MAP out 76–104 ( $x_{1}$ )	$β_{1}$ = –0.006	0.0026	–2.566	0.0103*

*p < 0.05.

Table 11

Least absolute shrinkage and selection operator (LASSO) solution and logistic regression of most predictive asymmetric range with leave-one-out cross-validation (LOOCV).

Model: $l = β_{0} + β_{1} x_{1}$ , where $x_{1} :$ time of MAP outside range 76–117 mmHg (n = 103 patients)
LOOCV: average observed accuracy = 0.5728; average kappa statistic = 0.102
Predictor	Coef. estimate (logit)	Std. error	z-Value	p-Value
Intercept	$β_{0}$ = 0.2881	0.287	1.002	0.316
Time MAP out 76–117 ( $x_{1}$ )	$β_{1}$ = –0.00788	0.0027	–2.828	0.00468**

**p < 0.01.

Altogether, the findings indicate that the time of MAP outside a measurable normotensive range during surgery is associated with lower odds of recovering at least one AIS grade. Our analysis suggests the optimal range for recovery is between 76–104 and 76–117 mmHg. Notice that while range 76–104 mmHg has higher predictive utility than 76–117 mmHg (mean LOOCV accuracy of 61.16 % vs. 57.28%), the difference in variance of the probability of AIS improvement explained by these two predictors is minimal (<4% difference in RV). Therefore, from a modeling perspective, we broadly conclude that the upper limit of the MAP range is probably anywhere between 104 and 117 mmHg.

Building a predictive model of outcome

Finally, we wanted to study the prediction utility of a model including the analyzed features together with other patient characteristics. We focused on three classification tasks: a model predicting AIS improvement of at least one grade at discharge, a model predicting AIS A at discharge, and a model predicting AIS D at discharge. We chose to predict AIS A and D instead of a multiclass prediction of the AIS at discharge in concordance to our previous studies (Kyritsis et al., 2021) and because of the low representation of other grades in our dataset (Table 1). For each of the three classification tasks, we performed an exhaustive search of all possible additive models with at least one of the predictors of interest: quadratic aMAP, aHR, length of surgery, days from surgery to discharge, age, AIS grade at admission, dichotomized NLI (cervical, non-cervical), time of MAP out of range 76–104, and time of MAP out of range 76–117. We selected the parsimonious model as the model that minimized the small-sample corrected AIC (Table 12). Next, for the selected best model for each task, we performed LOOCV performance evaluation and prediction threshold calibration balancing prediction sensitivity and specificity (Figure 5). The model predicting AIS improvement had a cross-validated AUC of 0.74, the model predicting AIS A at discharge had a cross-validated AUC of 0.88, and the model predicting AIS D at discharge had a cross-validation AUC of 0.84. Other metrics of classification performance can be seen in Table 12. Both the parsimonious model predicting AIS improvement and the one predicting AIS A at discharge included quadratic aMAP as an important predictor. The model predicting AIS A also included the time of MAP out of range 76–117 mmHg. The model predicting AIS D did not include any of the MAP associated terms, suggesting that patients discharged with AIS D can be predicted without considering their MAP during OR. In fact, training the same model but with the inclusion of the quadratic aMAP term resulted in slightly worse prediction performance (AUC 0.84 vs. 0.83). Training the models predicting AIS improvement and AIS A at discharge but without a MAP component (quadratic MAP term or time of MAP out of range) reduced the model performance considerably (AUC, AIS improvement: 0.74 vs. 0.52; AIS A discharge: 0.88 vs. 0.78).

Figure 5

Download asset Open asset

Leave-one-out cross-validation (LOOCV) performance of prediction models.

We built three prediction models, one to predict American Spinal Injury Association (ASIA) Impairment Scale (AIS) improvement of at least one grade at discharge (AIS impro., a), one to predict AIS A at discharge (A at disch., b) and one to predict AIS D at discharge (D at disch., c). The sensitivity and specificity for each model was computed out of the prediction probability of LOOCV, where each leave-one-out patient is predicted with the model that was trained without that patient. For each model, the classification threshold was set at the probability that balances sensitivity and specificity (dashed red line). The receiving operation curve (ROC) and area under the curve (AUC) for the three models are presented in d.

Table 12

Best prediction models of outcome.

Model predicting AIS improvement: $l = β_{0} + β_{11} x_{1} + β_{12} x_{1}^{2} + β_{2} x_{2} + β_{3} x_{3} + β_{4} x_{4} + β_{5} x_{5}$ Model predicting AIS A: $l = β_{0} + β_{11} x_{1} + β_{12} x_{1}^{2} + β_{2} x_{2} + β_{3} x_{3} + β_{4} x_{4} + β_{5} x_{5} + β_{6} x_{6} + β_{7} x_{7}$ Model predicting AIS D: $l = β_{0} + β_{2} x_{2} + β_{3} x_{3} + β_{4} x_{4} + β_{5} x_{5} + β_{8} x_{8} + β_{9} x_{9}$ where $x_{1} :$ average MAP; $x_{2}$ : AIS admission A (‘yes’, ‘no’); $x_{3}$ : AIS admission B (‘yes’, ‘no’); $x_{4}$ : AIS admission C (‘yes’, ‘no’); $x_{5}$ : AIS admission D (‘yes’, ‘no’); $x_{6}$ : NLI non-cervical; $x_{7}$ : Time MAP out 76–117; $x_{8}$ : Length of surgery; $x_{9}$ : Age; (AIS admission E and NLI cervical were set as the reference levels for the corresponding variable and are part of the intercept). All metrics are on LOOCV prediction (n = 93)

	Model AIS improv.	Model AIS A	Model AIS D
Predictor	Coef. estimate (logit)	Coef. estimate (logit)	Coef. estimate (logit)
Intercept	$β_{0}$ = –16.24	$β_{0}$ = 20.466	$β_{0}$ = 1.558
Average MAP ( $x_{1}$ )	$β_{11}$ = 7.374	$β_{11}$ = 27.031
Average MAP (Cohn et al., 2010) ( $x_{1}$ )	$β_{12}$ = –8.215	$β_{12}$ = –17.138
AIS admission A ( $x_{2}$ )	$β_{2}$ = 15.54	$β_{2}$ = –22.814	$β_{2}$ = 2.324
AIS admission B ( $x_{3}$ )	$β_{3}$ = 16.1818	$β_{3}$ = –20.38	$β_{3}$ = 0.41
AIS admission C ( $x_{4}$ )	$β_{4}$ = 16.752	$β_{4}$ = –19.01	$β_{4}$ = –2.591
AIS admission D ( $x_{5}$ )	$β_{5}$ = 14.828	$β_{5}$ = 0.217	$β_{5}$ = –2.624
NLI non-Cervical ( $x_{6}$ )		$β_{6}$ = –1.228
Time MAP out 76–117 ( $x_{7}$ )		$β_{7}$ = 0.017
Length of Surgery ( $x_{8}$ )			$β_{8}$ = –0.0044
Age ( $x_{9}$ )			$β_{9}$ = 0.03
Model performance metric	Metric value	Metric value	Metric value
Accuracy (95% CI)	0.73 (0.629, 0.818)	0.82 (0.735, 0.898)	0.806 (0.71, 0.881)
AUC	0.743	0.88	0.87
Kappa	0.45	0.629	0.573
Sensitivity	0.71	0.812	0.793
Specificity	0.74	0.836	0.812
Positive predicted value	0.658	0.72	0.657
Negative predicted value	0.788	0.89	0.896

Altogether, this suggests that models can be built for predicting AIS improvement or AIS A at discharge and that such the model performance critically depends on MAP during OR. Conversely, we did not find evidence that predicting AIS D at hospital discharge is dependent on intra-operative MAP.

Discussion

Acute hypotension is common in patients with SCI due to neurogenic shock (Lehmann et al., 1987; Krassioukov et al., 2007) and autonomic dysregulation (Lehmann et al., 1987), probably contributing to post-traumatic spinal ischemia (Streijger et al., 2018; Hall and Wolf, 1987), which is known to cause impaired neurological recovery in animal models (Fehlings et al., 1989). Level 4 evidence from a small single-center case series study in the 1990s suggested that MAP augmentation to 85–90 mmHg during the first 5–7 days after injury was linked to neurological recovery in acute SCI (Levi et al., 1993; Vale et al., 1997). These results are the basis of clinical guidelines for avoidance of hypotension in acute SCI management (Aarabi et al., 2013). However, while numerous clinical studies support MAP augmentation, the arbitrary, recommended MAP goal has been controversial (Cohn et al., 2010; Hawryluk et al., 2015; Saadeh et al., 2017). Recent analysis of high-frequency ICU monitoring data (Hawryluk et al., 2015) and systematic meta-analysis of post-surgery management (Saadeh et al., 2017) suggest that the MAP threshold to avoid is actually lower (~75 mmHg) than the current recommendation of 85 mmHg, and that MAP management might be effective at shorter duration (< 5 days post-injury) than the 7-day goal (Saadeh et al., 2017). The present study represents a multicenter, data-driven, and cross-validated re-evaluation in a different setting (during surgery as compared with prior ICU studies).

Our analysis support that there must be a MAP range during surgery at which neurological recovery is maximized, providing further evidence that MAP management for maintaining normotension might be more beneficial for patient outcome than MAP augmentation for hypotension avoidance alone (Ehsanian, 2020; Nielson et al., 2015). The low boundary of 76 mmHg found in our ultra-early analysis further supports previous suggestions for lowering the intervention threshold (Cohn et al., 2010; Hawryluk et al., 2015; Saadeh et al., 2017). On the other side, we find an upper boundary to MAP management between 104 and 117 mmHg, above which the probability of improvement is reduced. Thus, the proposal for MAP augmentation with vasopressors to increase spinal cord perfusion (Saadoun and Papadopoulos, 2016) has a limit since spinal hyper-perfusion pressure can be detrimental (Saadoun and Papadopoulos, 2016). The physiological rational is that high blood pressure induced by vasopressors can translate to increased risk of hemorrhage in the injured spinal cord, exacerbating tissue damage (Soubeyrand et al., 2014; Streijger et al., 2018; Guha et al., 1987). Moreover, the use of some vasopressors might cause more complications in patients (Inoue et al., 2014) while also potentially contributing to intra-spinal hemorrhage. In fact, recent results in acute experimental SCI suggest controlling for hemodynamic dysregulation through a cardiac-focused treatment instead of using standard vasopressors such as norepinephrine (Williams et al., 2020). Specifically, the authors demonstrated that dobutamine can correct for hemodynamic anomalies and increase blood flow to the spinal cord while reducing the risk of hemorrhage compared to norepinephrine. Furthermore, hypertension during surgery in rodent SCI has been associated with lower probability of recovery (Nielson et al., 2015), probably related to higher tissue damage. Our findings together with previous work (Ehsanian, 2020) also translate these animal study results to humans, indicating that prolonged periods of hypertension early after injury can be a predictor of poor neurological recovery in patients with SCI.

An important finding of our study is the indication that level of injury and injury severity modify the association of MAP with neurological recovery. We observed that normotensive MAP during surgery predicts AIS improvement in patients with cervical SCI but not in patients with lower injuries (thoracic, lumbar, and sacral). While the heterogeneity of our population and low sample size for patients with non-cervical SCI sets limitations on interpreting the results, our finding raises a relevant question regarding precision management of patients with SCI. Patients with cervical SCI present more frequently with hemodynamic and cardiac abnormalities than patients with thoracolumbar SCI, increasing the need for treatment (Lehmann et al., 1987). This is due to sympathetic dysregulation in upper cord injuries, which reduces sympathetic tone likely causing reduced heart contractility, bradycardia, and hypotension (Lehmann et al., 1987; Myers et al., 2007; Teasell et al., 2000). This is particularly true for individuals with severe cervical injury (Lehmann et al., 1987). In that context, our results may indicate that those patients with cervical SCI that are more difficult to maintain within a normotensive MAP are probably less likely to improve in neurological function. Alternatively, it could also be the case that more aggressive MAP management treatment is performed in these patients during their course in the hospital, which could increase the chances of aggravating secondary cord injury. Hemodynamic instability early after injury could serve as a prognostic physiology-based biomarker in a subset of the population, providing a potential tool for precision medicine in SCI. Hence, we have established basic prediction models around non-linear features of MAP that could serve as a benchmark for future machine learning development.

Another relevant contribution of this work is the analytical workflow. First, we demonstrate that topology-based analytics can undercover associations for hypothesis generation during exploratory analysis in a cross-species validation. Our group has previously used a similar workflow in data from animal models (Nielson et al., 2015) suggesting that hypertension is a predictor of neurological recovery and providing rational for the present study. Hence, our work constitutes a successful story of translating machine intelligence analytical tools from animals to humans. Second, we provide further illustration that patient similarity networks are useful and interpretable representations of multidimensional datasets that capture associations during exploratory analysis that can then be validated by network-independent confirmatory analysis. Third, we successfully combine Isomap, a non-linear dimensionality reduction method, with topology-based metrics to evaluate embedding solutions. Fourth, our method for finding the MAP range could be expanded and deployed in other settings. Lastly, our workflow captures representations of multidimensional time-series of different lengths into a network that is actionable.

Limitations of this study include the retrospective nature of the analysis, the relatively small sample size (although large for SCI), and the use of an estimated ordinal scale (AIS grade) as an indicator of neurological recovery. An important consideration is the difficulty of determining AIS grade early after injury. Moreover, other factors not considered in this analysis such as MAP levels before or after surgery or the use of vasopressors might influence the results. Future research with more granular data should address these and other important questions.

Data availability

Source data has been deposited to the Open Data Commons for Spinal Cord Injury (odc-sci.org; RRID:SCR_016673) under the accession number ODC-SCI:245 (https://doi.org/10.34945/F5R59J) and ODC-SCI:246 (https://doi.org/10.34945/F5MG68).

The following data sets were generated

1. Torres-Espin A
2. Haefeli J
3. Ehsanian R
4. Torres D
5. Almeida C
6. Huie J
7. Chou A
8. Dirlikov B
9. Suen C G
10. Nielson JL
11. Kyritsis N
12. Duong-Fernandez X
13. Thomas LH
14. Hemmerle T
15. Morozov DD
16. Sanderson D
17. Talbott N
18. Manley J
19. Dhall GT
20. Whetstone SS
21. Bresnahan WD
22. Beattie JC
23. McKenna MS
24. Pan SL
25. Ferguson JZ
(2021) Open Data Commons for Spinal Cord Injury
ASIA Impairment Scale, level of injury, intraoperative time series mean arterial pressure and heart rate after spinal cord injury in patients in a multi-site retrospective TRACK-SCI cohort: site 1 of 2.
https://doi.org/10.34945/F5R59J
1. Torres-Espín A
2. Haefeli J
3. Ehsanian R
4. Torres D
5. Almeida C A
6. Huie J
7. Chou A
8. Dirlikov B
9. Suen CG
10. Nielson JL
11. Kyritsis N
12. Duong-Fernandez X
13. Thomas LH
14. Hemmerle DD
15. Morozov D
16. Sanderson N
17. Talbott JF
18. Manley GT
19. Dhall SS
20. Whetstone WD
21. Bresnahan J C
22. Beattie MS
23. McKenna SL
24. Pan JZ
25. Ferguson AR
(2021) Open Data Commons for Spinal Cord Injury
Intraoperative time series mean arterial pressure and heart rate after spinal cord injury in patients in a multi-site retrospective TRACK-SCI cohort: site 2 of 2.
https://doi.org/10.34945/F5MG68

References

1. Aarabi B
2. Hadley MN
3. Dhall SS
4. Gelb DE
5. Hurlbert RJ
6. Rozzelle CJ
7. Ryken TC
8. Theodore N
9. Walters BC
(2013) Management of acute traumatic central cord syndrome (ATCCS)
Neurosurgery 72 Suppl 2:195–204.
https://doi.org/10.1227/NEU.0b013e318276f64b
- PubMed
- Google Scholar
Software
1. Almende BV
(2019) Benoit Thieurmel and Titouan Robert. visNetwork: Network Visualization using ‘vis.js, version 2.0.7
R Package.

http://datastorm-open.github.io/visNetwork/
Software
1. Calcagno VG
(2020)
Model Selection and Multimodel Inference Made Easy, version 1.0.7

Gmulti.
(2004) Finding community structure in very large networks
Physical Review E 70:066111.
https://doi.org/10.1103/PhysRevE.70.066111
- Google Scholar
1. Cohn J
2. Wright J
3. McKenna S
4. Bushnik T
(2010) Impact of Mean Arterial Blood Pressure During the First Seven Days Post Spinal Cord Injury
Topics in Spinal Cord Injury Rehabilitation 15:96–106.
https://doi.org/10.1310/sci1503-96
- Google Scholar
Software
1. Csárdi G
2. Nepusz T
(2006) The igraph software package for complex network research, version 1
BibSonomy.

http://igraph.sf.net
1. Ehsanian R
(2020) Exploration of surgical blood pressure management and expected motor recovery in individuals with traumatic spinal cord injury: This article has been corrected since Advance Online Publication and a correction is also printed in this issue
Spinal Cord 58:377–386.
https://doi.org/10.1038/s41393-019-0370-5
- Google Scholar
(1989) The effect of nimodipine and dextran on axonal function and blood flow following experimental spinal cord injury
Journal of Neurosurgery 71:403–416.
https://doi.org/10.3171/jns.1989.71.3.0403
- PubMed
- Google Scholar
(2010) Regularization Paths for Generalized Linear Models via Coordinate Descent
Journal of Statistical Software 33:1–22.
https://doi.org/10.18637/jss.v033.i01
- PubMed
- Google Scholar
1. Guha A
2. Tator CH
3. Piper I
(1987) Effect of a calcium channel blocker on posttraumatic spinal cord blood flow
Journal of Neurosurgery 66:423–430.
https://doi.org/10.3171/jns.1987.66.3.0423
- PubMed
- Google Scholar
1. Hall ED
2. Wolf DL
(1987) Post-traumatic spinal cord ischemia: relationship to injury severity and physiological parameters
Central Nervous System Trauma 4:15–25.
https://doi.org/10.1089/cns.1987.4.15
- Google Scholar
1. Hawryluk G
2. Whetstone W
3. Saigal R
4. Ferguson A
5. Talbott J
6. Bresnahan J
7. Dhall S
8. Pan J
9. Beattie M
10. Manley G
(2015) Mean Arterial Blood Pressure Correlates with Neurological Recovery after Human Spinal Cord Injury: Analysis of High Frequency Physiologic Data
Journal of Neurotrauma 32:1958–1967.
https://doi.org/10.1089/neu.2014.3778
- PubMed
- Google Scholar
(2014) Medical and Surgical Management after Spinal Cord Injury: Vasopressor Usage, Early Surgerys, and Complications
Journal of Neurotrauma 31:284–291.
https://doi.org/10.1089/neu.2013.3061
- PubMed
- Google Scholar
Software
1. Jari H
2. Peter R
3. Minchin R
4. O’Hara b
5. Gavin LS
6. Stevens H
7. Peter Solymos M
(2019) Eduard Szoecs and, & Helene Wagner. Vegan: Community Ecology Package, version 2.5-7
R Package.

https://cran.r-project.org
(2007) Assessment of autonomic dysfunction following spinal cord injury: Rationale for additions to International Standards for Neurological Assessment
Journal of Rehabilitation Research and Development 44:103–112.
https://doi.org/10.1682/JRRD.2005.10.0159
- PubMed
- Google Scholar
Software
1. Kuhn M
2. Contributions from Jed Wing
3. Weston S
4. Williams A
5. Keefer C
6. Engelhardt A
7. Cooper T
8. Mayer Z
9. Kenkel B
10. the R Core Team
11. Benesty M
12. Lescarbeau R
13. Ziem A
14. Scrucca L
15. Tang Y
16. Candan C
17. Hunt T
(2019) Classification and Regression Training, version 6.0-90
Caret.

https://github.com/topepo/caret/
1. Kyritsis N
2. Torres-Espín A
3. Schupp PG
4. Huie JR
5. Chou A
6. Duong-Fernandez X
7. Thomas LH
8. Tsolinas RE
9. Hemmerle DD
10. Pascual LU
11. Singh V
12. Pan JZ
13. Talbott JF
14. Whetstone WD
15. Burke JF
16. DiGiorgio AM
17. Weinstein PR
18. Manley GT
19. Dhall SS
20. Ferguson AR
21. Oldham MC
22. Bresnahan JC
23. Beattie MS
(2021) Diagnostic blood RNA profiles for human acute spinal cord injury
The Journal of Experimental Medicine 218:e20201795.
https://doi.org/10.1084/jem.20201795
- PubMed
- Google Scholar
(1987) Cardiovascular abnormalities accompanying acute spinal cord injury in humans: incidence, time course and severity
Journal of the American College of Cardiology 10:46–52.
https://doi.org/10.1016/S0735-1097(87)80158-4
- PubMed
- Google Scholar
1. Levi L
2. Wolf A
3. Belzberg H
(1993)
Hemodynamic parameters in patients with acute cervical cord trauma: description, intervention, and prediction of outcome

Neurosurgery 33:1007–1016.
- PubMed
- Google Scholar
Software
1. Marchette DJ
(2015) cccd: Class Cover Catch Digraphs, version 15
R Package.

https://CRAN.R-project.org/package=cccd
(2000) Vascular events after spinal cord injury: contribution to secondary pathogenesis
Physical Therapy 80:673–687.
https://doi.org/10.1093/ptj/80.7.673
- PubMed
- Google Scholar
(2007) Acute cervical traumatic spinal cord injury: MR imaging findings correlated with neurologic outcome--prospective study with 100 consecutive patients
Radiology 243:820–827.
https://doi.org/10.1148/radiol.2433060583
- PubMed
- Google Scholar
1. Myers J
2. Lee M
3. Kiratli J
(2007) Cardiovascular disease in spinal cord injury: an overview of prevalence, risk, evaluation, and management
American Journal of Physical Medicine & Rehabilitation 86:142–152.
https://doi.org/10.1097/PHM.0b013e31802f0247
- PubMed
- Google Scholar
Book
1. National Spinal Cor Injury Statistical Center
(2021)
Spinal Cord Injury Facts and Figures at a Glance

National Spinal Cor Injury Statistical Center.
- Google Scholar
Software
1. Neuwirth E
(2014) RColorBrewer: ColorBrewer Palettes, version Version 11-2
R Package.

https://CRAN.R-project.org/package=RColorBrewer
1. Newman MEJ
(2003) Mixing patterns in networks
Physical Review E 67:ArXivcond–Mat0209450.
https://doi.org/10.1103/PhysRevE.67.026126
- Google Scholar
1. Nielson JL
2. Paquette J
3. Liu AW
4. Guandique CF
5. Tovar CA
6. Inoue T
7. Irvine K-A
8. Gensel JC
9. Kloke J
10. Petrossian TC
11. Lum PY
12. Carlsson GE
13. Manley GT
14. Young W
15. Beattie MS
16. Bresnahan JC
17. Ferguson AR
(2015) Topological data analysis for discovery in preclinical spinal cord injury and traumatic brain injury
Nature Communications 6:8581.
https://doi.org/10.1038/ncomms9581
- PubMed
- Google Scholar
Software
1. Ooi H
(2019) glmnetUtils: Utilities for ‘Glmnet, version 112
R Package Version.

https://github.com/hongooi73/glmnetUtils
1. Pai S
2. Bader GD
(2018) Patient Similarity Networks for Precision Medicine
Journal of Molecular Biology 430:2924–2938.
https://doi.org/10.1016/j.jmb.2018.05.037
- PubMed
- Google Scholar
(2018) Patient similarity for precision medicine: A systematic review
Journal of Biomedical Informatics 83:87–96.
https://doi.org/10.1016/j.jbi.2018.06.001
- PubMed
- Google Scholar
1. Paul R
2. Chalup SK
(2017) A study on validating non-linear dimensionality reduction using persistent homology
Pattern Recognition Letters 100:160–166.
https://doi.org/10.1016/j.patrec.2017.09.032
- Google Scholar
1. Pons P
2. Latapy M
(2005) Computing communities in large networks using random walks
Computer and Information Sciences 3733:284–293.
https://doi.org/10.1007/11569596_31
- Google Scholar
Software
1. R core team
(2019) A language and environment for statistical computing R Foundation for Statistical Computing
R Core Team.

https://www.R-project.org/
Software
1. Rich B
(2018) table1: Tables of Descriptive Statistics in HTML, version 1.4.2
R Package Version.

https://github.com/benjaminrich/table1
1. Rieck B
2. Leitte H
(2015) Persistent Homology for the Evaluation of Dimensionality Reduction Schemes
Computer Graphics Forum 34:431–440.
https://doi.org/10.1111/cgf.12655
- Google Scholar
(2017) Classifications In Brief: American Spinal Injury Association (ASIA) Impairment Scale
Clinical Orthopaedics and Related Research 475:1499–1504.
https://doi.org/10.1007/s11999-016-5133-4
- PubMed
- Google Scholar
1. Saadeh YS
2. Smith BW
3. Joseph JR
4. Jaffer SY
5. Buckingham MJ
6. Oppenlander ME
7. Szerlip NJ
8. Park P
(2017) The impact of blood pressure management after spinal cord injury: a systematic review of the literature
Neurosurgical Focus 43:E20.
https://doi.org/10.3171/2017.8.FOCUS17428
- PubMed
- Google Scholar
1. Saadoun S
2. Papadopoulos MC
(2016) Spinal cord injury: is monitoring from the injury site the future?
Critical Care 20:308.
https://doi.org/10.1186/s13054-016-1490-3
- PubMed
- Google Scholar
(2005) ROCR: visualizing classifier performance in R
Bioinformatics 21:3940–3941.
https://doi.org/10.1093/bioinformatics/bti623
- Google Scholar
1. Soubeyrand M
2. Dubory A
3. Laemmel E
4. Court C
5. Vicaut E
6. Duranteau J
(2014) Effect of norepinephrine on spinal cord blood flow and parenchymal hemorrhage size in acute-phase experimental spinal cord injury
European Spine Journal 23:658–665.
https://doi.org/10.1007/s00586-013-3086-9
- PubMed
- Google Scholar
1. Streijger F
2. So K
3. Manouchehri N
4. Gheorghe A
5. Okon EB
6. Chan RM
7. Ng B
8. Shortt K
9. Sekhon MS
10. Griesdale DE
11. Kwon BK
(2018) A Direct Comparison between Norepinephrine and Phenylephrine for Augmenting Spinal Cord Perfusion in a Porcine Model of Spinal Cord Injury
Journal of Neurotrauma 35:1345–1357.
https://doi.org/10.1089/neu.2017.5285
- PubMed
- Google Scholar
Software
1. Team R
(2018) Integrated Development for R, version 1.1
RStudio, Inc.

http://www.rstudio.com/
(2000) Cardiovascular consequences of loss of supraspinal control of the sympathetic nervous system after spinal cord injury
Archives of Physical Medicine and Rehabilitation 81:506–516.
https://doi.org/10.1053/mr.2000.3848
- PubMed
- Google Scholar
(2000) A Global Geometric Framework for Nonlinear Dimensionality Reduction
Science 290:2319–2323.
https://doi.org/10.1126/science.290.5500.2319
- PubMed
- Google Scholar
1. Tibshirani R
(1996) Regression Shrinkage and Selection via the Lasso
Journal of the Royal Statistical Society 58:267–288.
https://doi.org/10.1111/j.2517-6161.1996.tb02080.x
- Google Scholar
1. Vale FL
2. Burns J
3. Jackson AB
4. Hadley MN
(1997) Combined medical and surgical treatment after acute spinal cord injury: results of a prospective pilot study to assess the merits of aggressive medical resuscitation and blood pressure management
Journal of Neurosurgery 87:239–246.
https://doi.org/10.3171/jns.1997.87.2.0239
- PubMed
- Google Scholar
(2018) TDAstats: R pipeline for computing persistent homology in topological data analysis
Journal of Open Source Software 3:860.
https://doi.org/10.21105/joss.00860
- PubMed
- Google Scholar
1. Weng S
2. Zhang C
3. Lin Z
4. Zhang X
(2005) Mining the structural knowledge of high-dimensional medical data using isomap
Medical & Biological Engineering & Computing 43:410–412.
https://doi.org/10.1007/BF02345820
- PubMed
- Google Scholar
1. Wickham H
(2007) Reshaping Data with the reshape Package
Journal of Statistical Software 21:1–20.
https://doi.org/10.18637/jss.v021.i12
- Google Scholar
Software
1. Wickham H
(2016) ggplot2: Elegant Graphics for Data Analysis
Springer-Verlag.

http://ggplot2.org
Software
(2018) A Grammar of Data Manipulation, version 1.0.7
Dplyr.

https://dplyr.tidyverse.org
1. Williams AM
2. Manouchehri N
3. Erskine E
4. Tauh K
5. So K
6. Shortt K
7. Webster M
8. Fisk S
9. Billingsley A
10. Munro A
11. Tigchelaar S
12. Streijger F
13. Kim K-T
14. Kwon BK
15. West CR
(2020) Cardio-centric hemodynamic management improves spinal cord oxygenation and mitigates hemorrhage in acute spinal cord injury
Nature Communications 11:5209.
https://doi.org/10.1038/s41467-020-18905-8
- PubMed
- Google Scholar

Article and author information

Author details

Abel Torres-Espín

Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States

Contribution
Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review and editing

Contributed equally with
Jenny Haefeli

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-9787-8738
Jenny Haefeli

Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States

Contribution
Conceptualization, Formal analysis, Investigation, Methodology, Writing – review and editing

Contributed equally with
Abel Torres-Espín

Competing interests
No competing interests declared
Reza Ehsanian

Division of Physical Medicine and Rehabilitation, Department of Orthopaedics and Rehabilitation, University of New Mexico School of Medicine, Albuquerque, United States

Contribution
Writing – review and editing, Data interpretation

Competing interests
No competing interests declared
Dolores Torres

Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States

Contribution
Data curation, Writing – review and editing, Data curation

Competing interests
No competing interests declared
Carlos A Almeida

Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States

Contribution
Data curation, Writing – review and editing

Competing interests
No competing interests declared
J Russell Huie
1. Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States
2. San Francisco Veterans Affairs Healthcare System, San Francisco, United States
Contribution
Writing – review and editing, Data interpretation

Competing interests
No competing interests declared
Austin Chou

Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States

Contribution
Writing – review and editing, Data interpretation

Competing interests
No competing interests declared
Dmitriy Morozov

Computational Research Division, Lawrence Berkeley National Laboratory, Berkeley, United States

Contribution
Formal analysis, Methodology, Visualization, Writing – review and editing, Data interpretation

Competing interests
No competing interests declared
Nicole Sanderson

Lawrence Berkeley National Laboratory, Berkeley, United States

Contribution
Formal analysis, Methodology, Visualization, Writing – review and editing, Data interpretation

Competing interests
No competing interests declared
Benjamin Dirlikov

Rehabilitation Research Center, Santa Clara Valley Medical Center, San Jose, United States

Contribution
Data curation, Writing – review and editing

Competing interests
No competing interests declared
Catherine G Suen

Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States

Contribution
Data curation, Writing – review and editing, Data interpretation

Competing interests
No competing interests declared
Jessica L Nielson
1. Department of Psychiatry and Behavioral Science, and University of Minnesota, Minneapolis, United States
2. Institute for Health Informatics, University of Minnesota, Minneapolis, United States
Contribution
Writing – review and editing, Data interpretation

Competing interests
No competing interests declared
Nikos Kyritsis

Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States

Contribution
Writing – review and editing, Data interpretation

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-7801-5796
Debra D Hemmerle

Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States

Contribution
Data curation, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-2796-6107
Jason F Talbott

Department of Radiology and Biomedical Imaging, University of California, San Francisco, San Francisco, United States

Contribution
Conceptualization, Investigation, Methodology, Writing – review and editing

Competing interests
No competing interests declared
Geoffrey T Manley

Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States

Contribution
Conceptualization, Resources, Writing – review and editing

Competing interests
No competing interests declared
Sanjay S Dhall

Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States

Contribution
Conceptualization, Resources, Writing – review and editing

Competing interests
No competing interests declared
William D Whetstone

Department of Emergency Medicine, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States

Contribution
Data curation, Methodology, Validation, Writing – review and editing

Competing interests
No competing interests declared
Jacqueline C Bresnahan
1. Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States
2. San Francisco Veterans Affairs Healthcare System, San Francisco, United States
Contribution
Conceptualization, Resources, Writing – review and editing

Competing interests
No competing interests declared
Michael S Beattie
1. Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States
2. San Francisco Veterans Affairs Healthcare System, San Francisco, United States
Contribution
Conceptualization, Resources, Writing – review and editing, Data interpretation

Competing interests
No competing interests declared
Stephen L McKenna
1. Department of Physical Medicine and Rehabilitation, Santa Clara Valley Medical Center, San Jose, United States
2. Department of Neurosurgery, Stanford University, Stanford, United States
Contribution
Conceptualization, Data curation, Methodology, Resources, Writing – review and editing, Data interpretation

Competing interests
No competing interests declared
Jonathan Z Pan

Department of Anesthesia and Perioperative Care, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States

Contribution
Conceptualization, Data curation, Methodology, Writing – review and editing, Data interpretation

For correspondence
jonathan.pan@ucsf.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-5814-3707
Adam R Ferguson
1. Weill Institute for Neurosciences; Brain and Spinal Injury Center (BASIC), Department of Neurological Surgery, University of California, San Francisco; Zuckerberg San Francisco General Hospital and Trauma Center, San Francisco, United States
2. San Francisco Veterans Affairs Healthcare System, San Francisco, United States
Contribution
Conceptualization, Formal analysis, Funding acquisition, Methodology, Supervision, Writing – original draft, Writing – review and editing, Data interpretation

For correspondence
adam.ferguson@ucsf.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-7102-1608
The TRACK-SCI Investigators

Competing interests
No competing interests declared
1. MS Beattie
2. JC Bresnahan
3. JF Burke
4. A Chou
5. CA de Almeida
6. SS Dhall
7. AM DiGiorgio
8. X Doung-Fernandez
9. AR Ferguson
10. J Haefeli
11. DD Hemmerle
12. JR Huie
13. N Kyritsis
14. GT Manley
15. S Moncivais
16. C Omondi
17. JZ Pan
18. LU Pascual
19. V Singh
20. JF Talbott
21. LH Thomas
22. A Torres-Espin
23. P Weinstein
24. WD Whetstone

Funding

National Institute of Neurological Disorders and Stroke (R01NS088475)

Adam R Ferguson

National Institute of Neurological Disorders and Stroke (R01NS122888)

Adam R Ferguson

National Institute of Neurological Disorders and Stroke (UH3NS106899)

Adam R Ferguson

U.S. Department of Veterans Affairs (1I01RX002245)

Adam R Ferguson

U.S. Department of Veterans Affairs (I01RX002787)

Adam R Ferguson

Wings for Life

Abel Torres-Espín
Adam R Ferguson

Craig H. Neilsen Foundation

Adam R Ferguson

Department of Defense (SC150198)

Michael S Beattie

Department of Defense (SC190233)

Michael S Beattie

Foundation for Anesthesia Education and Research (A123320)

Jonathan Z Pan

Department of Energy (DE-AC02-05CH11231)

Dmitriy Morozov

National Institute of Neurological Disorders and Stroke (U24NS122732)

Adam R Ferguson

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

Supported in part by NIH/NINDS: R01NS088475 (ARF); R01NS122888 (ARF); UH3NS106899 (ARF); U24NS122732 (ARF); Department of Veterans Affairs: 1I01RX002245 (ARF), I01RX002787 (ARF); Wings for Life Foundation (ATE, ARF); Craig H Neilsen Foundation (ARF); DOD: SC150198 (MSB); SC190233 (MSB); DOE: DE-AC02-05CH11231 (DM).

Ethics

Human subjects: This study constitutes a retrospective data analysis. All data was de-identified before pre-processing and analysis. Protocols for retrospective data extraction were approved by Institutional Research Board (IRB) under protocol numbers 11-07639 and 11-06997.

Version history

Received: March 2, 2021
Accepted: October 23, 2021
Accepted Manuscript published: November 16, 2021 (version 1)
Accepted Manuscript updated: November 17, 2021 (version 2)
Version of Record published: December 2, 2021 (version 3)
Version of Record updated: December 3, 2021 (version 4)

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.