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 longterm functional recovery in rodent models, motivating the present multicenter study in patients.
Methods:
Intraoperative monitoring records and neurological outcome data were extracted (n = 118 patients). We built a similarity network of patients from a lowdimensional space embedded using a nonlinear 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–[104117] 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); DOE: DEAC0205CH11231 (DM).
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 subjectcentered topological network analysis to identify phenotypes.
https://doi.org/10.7554/eLife.68015.sa0eLife 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, TorresEspin, 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 104117 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. TorresEspin, 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 realtime 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 singlecenter studies in humans suggest that lower postsurgery 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 (ultraacute phase) predicts longterm functional recovery in rodent models (Nielson et al., 2015). These prior crossspecies findings motivated the present multicenter study where we apply a datadriven workflow in patients with ultraacute 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 deidentified before preprocessing and analysis. Protocols for retrospective data extraction were approved by Institutional Research Board (IRB).
Cohort statistics
The differences in the AIS improvers/nonimprovers 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 ttest (twosided 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 (twosided test). pValues are presented in Table 1.
Topological network extraction and exploration of monitoring data workflow
Monitoring data preprocessing
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 preprocessed. Briefly, electronic data was exported from the SIS SQL database, deidentified 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 pointtopoint 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 timeseries data is not aligned between patients (without time dependency on monitoring values), the empirical cumulative distribution function (CDF) for each timeseries 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 ($\mathrm{d}\left(\mathrm{p},\mathrm{q}\right)=\sqrt{\sum _{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{p}}_{\mathrm{i}}{\mathrm{q}}_{\mathrm{i}}\right)}^{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 signaltonoise 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 nonlinear 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}\left(\hat{{D}_{m}},{D}_{y}\right)$ 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 (kgeodesic distance matrix for kIsomap). ${D}_{y}$ is the Euclidean distance matrix of the lowdimensional 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 VietorisRips filtration (Paul and Chalup, 2017) for $\hat{D}}_{m$ and ${D}_{y}$ for different d solutions (Figure 2—figure supplement 1). Then the topological zerodimensional and onedimensional Wasserstein (power = 2) distances (WD0 and WD1, respectively) were computed between $\hat{D}}_{m$ and ${D}_{y}$. In persistent homology, dimension 0 measures zerodimensional holes in the data (the connectivity of the datapoints, i.e., the number of connected sets) and topological dimension 1 measures onedimensional 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 lowdimensional 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 kNNG for the lowdimensional space. The cutoff threshold for adjacency was set to the minimal k that produced a fullconnected 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 nonlinear relationship between aMAP and the probability of improvement. Consequently, the following logistic models were considered ($l$ being the logodds or logit of the probability of improving): the null model with no predictors ($l={\beta}_{0}$, the simple model ($l={\beta}_{0}+{\beta}_{1}x$), the twograde polynomial model ($l={\beta}_{0}+{\beta}_{1}x+{\beta}_{2}{x}^{2}$), the threegrade polynomial model ($l={\beta}_{0}+{\beta}_{1}x+{\beta}_{2}{x}^{2}+{\beta}_{3}{x}^{3}$) and a natural spline model ($l={\beta}_{0}+{f}_{1}\left(x\right)$ )where ${f}_{1}\left(x\right)$ 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 twograde polynomial (Figure 3) and the natural spline (2df) with significant coefficients, confirming our hypothesis. These results were confirmed by leaveoneout crossvalidation (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.
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={\beta}_{0}+{\sum}_{j=1}^{p}{\beta}_{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 L1norm penalty). LOOCV was used to determine the lambda that shrunk the models to one predictor or MAP range. The onepredictor 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 smallsample 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 (Supplementary file 1 and Source code 1).
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 leaveoneout crossvalidation  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 deidentified datasets for analysis are deposited and accessible at the Open Data Commons for SCI (odcsci.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 intraoperative MAP for recovery
Intraoperative 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 intraoperative 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 lowdimensional space embedded using a nonlinear 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 1ce). Clustering the network of patients through a randomwalk algorithm, Walktrap (Pons and Latapy, 2005), revealed 11 different clusters where patients were regarded to share intraoperative hemodynamic phenotypes (Figure 2 andFigure 2—figure supplement 1fh). Importantly, this workflow was unsupervised: only the OR hemodynamic timeseries 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 highdimensional 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.
MAP has a nonlinear 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 nonlinear relationship between intraoperative 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 nonlinearity 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.
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 nonlinear 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 covariatesonly model had a significant coefficient (Table 5). These results indicate that even in the presence of the other factors, aMAP is still nonlinearly associated with AIS grade improvement at discharge.
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 pvalues = 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 nonlinear 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.
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 noncervical 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 noncervical 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 noncervical 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 noncervical population were marginally different (Table 1). Moreover, a logistic regression predicting AIS grade improvement by NLI categorization indicated that noncervical patients were significantly less likely to recover ($\beta $ = –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 noncervical 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.
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 nonsignificant pvalue 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^{2} 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 nonsignificant change in variance explained (p = 0.729) with a deviance difference of ~0.1%. All together indicates that the nonlinear relationship between aMAP and AIS improvement is independent of the days from surgery to discharge.
Intraoperative MAP range from 76[104117] mmHg predicts recovery
Since aMAP can obscure episodes of high deviation from average (Hawryluk et al., 2015) and has a nonlinear relationship with recovery, we hypothesized that there might be a range of intraoperative 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 (nonzero 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 crossvalidation 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.
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, noncervical), 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 smallsample 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 crossvalidated AUC of 0.74, the model predicting AIS A at discharge had a crossvalidated AUC of 0.88, and the model predicting AIS D at discharge had a crossvalidation 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).
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 intraoperative 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 posttraumatic 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 singlecenter 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 highfrequency ICU monitoring data (Hawryluk et al., 2015) and systematic metaanalysis of postsurgery 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 postinjury) than the 7day goal (Saadeh et al., 2017). The present study represents a multicenter, datadriven, and crossvalidated reevaluation 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 ultraearly 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 hyperperfusion 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 intraspinal hemorrhage. In fact, recent results in acute experimental SCI suggest controlling for hemodynamic dysregulation through a cardiacfocused 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 noncervical 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 physiologybased biomarker in a subset of the population, providing a potential tool for precision medicine in SCI. Hence, we have established basic prediction models around nonlinear 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 topologybased analytics can undercover associations for hypothesis generation during exploratory analysis in a crossspecies 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 networkindependent confirmatory analysis. Third, we successfully combine Isomap, a nonlinear dimensionality reduction method, with topologybased 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 timeseries 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 (odcsci.org; RRID:SCR_016673) under the accession number ODCSCI:245 (https://doi.org/10.34945/F5R59J) and ODCSCI:246 (https://doi.org/10.34945/F5MG68).

Open Data Commons for Spinal Cord InjuryASIA Impairment Scale, level of injury, intraoperative time series mean arterial pressure and heart rate after spinal cord injury in patients in a multisite retrospective TRACKSCI cohort: site 1 of 2.https://doi.org/10.34945/F5R59J

Open Data Commons for Spinal Cord InjuryIntraoperative time series mean arterial pressure and heart rate after spinal cord injury in patients in a multisite retrospective TRACKSCI cohort: site 2 of 2.https://doi.org/10.34945/F5MG68
References

Management of acute traumatic central cord syndrome (ATCCS)Neurosurgery 72 Suppl 2:195–204.https://doi.org/10.1227/NEU.0b013e318276f64b

SoftwareModel Selection and Multimodel Inference Made Easy, version 1.0.7Gmulti.

Finding community structure in very large networksPhysical Review E 70:066111.https://doi.org/10.1103/PhysRevE.70.066111

Impact of Mean Arterial Blood Pressure During the First Seven Days Post Spinal Cord InjuryTopics in Spinal Cord Injury Rehabilitation 15:96–106.https://doi.org/10.1310/sci150396

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 issueSpinal Cord 58:377–386.https://doi.org/10.1038/s4139301903705

Regularization Paths for Generalized Linear Models via Coordinate DescentJournal of Statistical Software 33:1–22.https://doi.org/10.18637/jss.v033.i01

Effect of a calcium channel blocker on posttraumatic spinal cord blood flowJournal of Neurosurgery 66:423–430.https://doi.org/10.3171/jns.1987.66.3.0423

Posttraumatic spinal cord ischemia: relationship to injury severity and physiological parametersCentral Nervous System Trauma 4:15–25.https://doi.org/10.1089/cns.1987.4.15

Assessment of autonomic dysfunction following spinal cord injury: Rationale for additions to International Standards for Neurological AssessmentJournal of Rehabilitation Research and Development 44:103–112.https://doi.org/10.1682/JRRD.2005.10.0159

Diagnostic blood RNA profiles for human acute spinal cord injuryThe Journal of Experimental Medicine 218:e20201795.https://doi.org/10.1084/jem.20201795

Cardiovascular abnormalities accompanying acute spinal cord injury in humans: incidence, time course and severityJournal of the American College of Cardiology 10:46–52.https://doi.org/10.1016/S07351097(87)801584

Hemodynamic parameters in patients with acute cervical cord trauma: description, intervention, and prediction of outcomeNeurosurgery 33:1007–1016.

Cardiovascular disease in spinal cord injury: an overview of prevalence, risk, evaluation, and managementAmerican Journal of Physical Medicine & Rehabilitation 86:142–152.https://doi.org/10.1097/PHM.0b013e31802f0247

BookSpinal Cord Injury Facts and Figures at a GlanceNational Spinal Cor Injury Statistical Center.

Mixing patterns in networksPhysical Review E 67:ArXivcond–Mat0209450.https://doi.org/10.1103/PhysRevE.67.026126

Patient Similarity Networks for Precision MedicineJournal of Molecular Biology 430:2924–2938.https://doi.org/10.1016/j.jmb.2018.05.037

Patient similarity for precision medicine: A systematic reviewJournal of Biomedical Informatics 83:87–96.https://doi.org/10.1016/j.jbi.2018.06.001

A study on validating nonlinear dimensionality reduction using persistent homologyPattern Recognition Letters 100:160–166.https://doi.org/10.1016/j.patrec.2017.09.032

Computing communities in large networks using random walksComputer and Information Sciences 3733:284–293.https://doi.org/10.1007/11569596_31

Persistent Homology for the Evaluation of Dimensionality Reduction SchemesComputer Graphics Forum 34:431–440.https://doi.org/10.1111/cgf.12655

Classifications In Brief: American Spinal Injury Association (ASIA) Impairment ScaleClinical Orthopaedics and Related Research 475:1499–1504.https://doi.org/10.1007/s1199901651334

ROCR: visualizing classifier performance in RBioinformatics 21:3940–3941.https://doi.org/10.1093/bioinformatics/bti623

Cardiovascular consequences of loss of supraspinal control of the sympathetic nervous system after spinal cord injuryArchives of Physical Medicine and Rehabilitation 81:506–516.https://doi.org/10.1053/mr.2000.3848

Regression Shrinkage and Selection via the LassoJournal of the Royal Statistical Society 58:267–288.https://doi.org/10.1111/j.25176161.1996.tb02080.x

TDAstats: R pipeline for computing persistent homology in topological data analysisJournal of Open Source Software 3:860.https://doi.org/10.21105/joss.00860

Mining the structural knowledge of highdimensional medical data using isomapMedical & Biological Engineering & Computing 43:410–412.https://doi.org/10.1007/BF02345820

Reshaping Data with the reshape PackageJournal of Statistical Software 21:1–20.https://doi.org/10.18637/jss.v021.i12
Decision letter

Arduino A MangoniReviewing Editor; Flinders Medical Centre, Australia

Matthias BartonSenior Editor; University of Zurich, Switzerland

Marcel KoppReviewer

Aaron PhillipsReviewer; University of Calgary, Canada
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Decision letter after peer review:
Thank you for submitting your article "Topological network analysis of patient similarity for precision management of acute blood pressure in spinal cord injury" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Marcel Kopp (Reviewer #2).
The Reviewers and Editors have discussed the reviews with one another, and this decision letter is to help you prepare a revised submission.
Essential revisions:
1) Please provide information regarding what covariate adjustment was used in the confirmatory logistic regression models.
2) Can the authors provide an explanation of why they chose the specific forms of clustering to identify patient phenotypes? Other, perhaps simpler and more common unsupervised machine learning algorithms could have been used.
3) Are the results sensitive to the defined outcome of improvement of at least one AIS grade? What happens if this is increased to e.g. 2 grades?
4) Two different approaches to analysis were used – i.e. essentially clustering of some form and also logistic regression (using e.g. quadratic and spline functions). Can the authors comment on whether these 2 approaches can be used interchangeably or whether one would be preferred over to the other to answer the research questions of interest. What advantage does the clustering approach have in terms of the research question?
5) Why was a simple accuracy metric used and not e.g. AUC? Does the accuracy metric account for an imbalance in the outcomes?
6) The LOOCV accuracy was not that high suggesting a lot of other factors might influence outcomes. Is the accuracy really high enough to support the use of MAP being used by clinicians for decision making and/or interventions to control MAP during surgery?
7) What variables were used for the LASSO. The prediction accuracy again seems very low for highdimensional dataset.
8) Only one dataset was used without splitting the data into a training and validation dataset. Are similar results for the topological network analysis obtained if the data is split for training and validation?
9) What was the modularity of the final network and does it suggest significant clustering?
10) Why was days from surgery to discharge used in the logistic regression models? Might it not be considered a mediator rather than a confounder – and how does its exclusion from the model influence the result?
11) The limitations are mentioned but not discussed or justified. This leads to the following questions: a) Why was the lesion level not included in the analysis? and b) Why did the authors only analyze MAD values during surgery? Because the analysis of MAP data from the ICU period published elsewhere showed similar results regarding the lower limit of MAP, wouldn't it be of interest to know how much overlap there is between the populations with critical MAP values during surgery and during stay in the ICU?
12) Introduction: Neither hypertension nor hypotension following acute SCI has been conclusively demonstrated to impact neurological recovery. Instead, guidelines and more recent work are based purely on observations and posthoc regression analyses. While the purported mechanism for repeated hypotensive episodes is clear, readers may benefit from at least a brief description of why both hypertension and hypotension could plausibly be important (aside from the fact that, again, noncausal observations demonstrated a relationship in the author's prior work).
13) AIS scores: Based on Table 1 most patients were discharged within 2 weeks after injury. The neurological exam is not so reliable at this point. This is a big limitation of the current work. Although a six month followup would be ideal to determine whether neurological recovery occurred, the authors should at least mention this.
14) All the analyses seem to have been conducted on an AIS change. The authors should demonstrate that their analysis holds for a more linear measure of recovery (e.g., total motor score).
15) Based on Figure 2 – Supplement 2, it is difficult to ascertain whether clusters contain a higher proportion of individuals that show an AIS improvement, and those individuals tend to have a MAP > 80 and < 100, is due to these clusters having individuals who were less severe to begin with (i.e., C,D,E) and therefore less likely to be hemodynamically unstable. One way to answer this would be to examine AIS A patients in a separate analysis and determine whether these findings hold. Because, the alternative explanation here is that this analysis is effectively finding a proxy for initial injury severity (i.e., more severe, more hemodynamically unstable) – and not that hemodynamic instability per se is the problem. Another analysis that could help complement this work and avoid this confound would be to use total motor score as the outcome instead of AIS conversion.
16) Logistic regression – Based on Figure 2p the overall trend looks more to be that higher MAP = > Pr(δ AIS grade > 0). The exception is only 2 data points on the top end. It is difficult to determine how robust the notion is that there is a 'too high' component to this data. Indeed it seems that a linear model does quite well for this analysis as well. Please see my comment below but this should be addressed as the concept of also having a 'top cutoff' is an extremely important clinical feature here.
17) Time outside MAP – The authors use an approach that systematically increases their window in both directions to find their optimal range of 76104. However, what happens if you then hold 76 and only increase on the upper end? Does this rapidly degrade the relationship? If not, again this would suggest that the evidence for a top end cutoff is not as strong. While I understand the authors briefly looked at this (methods) it seems worth exploring further as this is a critical point for clinical management. I do not see a good reason that the time outside the threshold can not at least be plotted to determine this relationship.
18) Data availability – The code and analysis should be made available to the reviewers. It is impossible to determine the accuracy of this type of analysis without it.
19) Discussion – It seems that the authors should discuss the confound of injury severity being linked with worse hemodynamic instability, and also worse neurological recovery. It would be helpful to include some of the suggested analyses to convince the readers that this confound does not explain the results since it is the most likely alternative explanation.
Reviewer #1:
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 subjectcentred topological network analysis to identify phenotypes. The weaknesses include its relatively modest prediction accuracy and the lack of internal and external validation in the primary network analysis.
Reviewer #2:
The major strength of the paper is the statistically highly advanced analysis based on highresolution data from acute SCI care, i.e. intraoperative mean arterial blood pressure (MAP) and heart rate. The steps of data exploration and analysis and their results are presented transparently. In conjunction with the results of previous studies suggesting that the lower threshold of MAP levels to be avoided in the ICU is about 75 mmHg, the main results of the this study imply that the minimum target MAP may be lower than the currently recommended 85 mmHg also during surgery. The analysis, which combines machine learning algorithms with logistic regression models, may serve as a template for datadriven studies also on other aspects of critical care in the field of SCI.
A weakness of the study is that some of the baseline neurological criteria were not included in the analysis. In particular, the neurological level of injury could be important for the research question, because the degree of blood pressure dysregulation also depends on the lesion level. The authors mention this limitation but do not explain why they accept it. Another limitation is the relatively small sample size of the study. Therefore, the specific results might have limited generalizability. Nevertheless, the study is an essential contribution to the readjustment of MAP threshold recommendations in the very acute stage of SCI and provides key information for the design of future precision medicine studies.
Reviewer #3:
The authors present a nice analysis of the relationship between intraoperative mean arterial pressure and neurological recovery after spinal cord injury, and I appreciate the opportunity to review this work. In general, the results are interesting and largely in line with a growing body of evidence that supports the use of hemodynamic monitoring in the acute phase after injury. There are three primary points with regard to this work that the authors can and should address prior to publication.
First, the exclusive use of AIS conversion as the outcome of neurological improvement is not ideal. Demonstrating that these findings are robust to other more continuous measures of neurological improvement such as motor/sensory scores would go a long way towards demonstrating that this finding is robust.
Second, the authors state in the discussion that MAP management may only be needed for <5 days postinjury. There does not appear to be strong data in the paper to support this point. Either more data should be added that supports this contention or this point should be removed.
Third, although the authors briefly discuss the confound of injury severity being linked to hemodynamic instability, the results are compelling enough at the moment to discount a role of injury severity. Individuals with more severe injuries will necessarily have greater hemodynamic instability, particularly in the hyperacute and acute phase after injury. These same individuals are also less likely to exhibit a conversion of their AIS score (e.g., individuals with AIS A/B). The authors control for this in some of their analyses (e.g., including a coefficient of initial AIS score in their regression models), yet their results seem to indicate that the clusters of individuals they focus on are indeed those with less severe injuries to start with. Demonstrating that injury severity is not a factor would require an analysis of only individuals with AIS A injuries at admission. This would be very compelling and enhance the impact of this work.
Overall, this is an interesting analysis that will be of use to the field.
https://doi.org/10.7554/eLife.68015.sa1Author response
Essential revisions:
1) Please provide information regarding what covariate adjustment was used in the confirmatory logistic regression models.
We adjusted the models by heart rate, length of surgery, days from surgery to discharge, age, and AIS grade at admission. We added text to the Results section that helps to clarify the models and the covariates.
2) Can the authors provide an explanation of why they chose the specific forms of clustering to identify patient phenotypes? Other, perhaps simpler and more common unsupervised machine learning algorithms could have been used.
We have added further explanation to the manuscript. In short, we and others have previously shown that topological analysis drives discovery of biologically meaningful associations with high resolution (Nielson et al., 2015; Nicolau et al., 2011). Here, we use topological network analysis for several reasons. First, physiological data is nonlinear, which conforms to a manifold in the multidimensional space. We use Isomap, a wellestablished nonlinear dimension reduction method, to find the low dimensional embedding of the data (Balasubramanian, 2002; Tenenbaum, 2000). Isomap works by constructing a neighbor network to approximate the geodesic distance between observations. Second, to ensure that Isomap approximates the shape of the manifold we use topological metrics. Third, the optimal Isomap solution had 4 dimensions, which makes it challenging for representing data for exploration. We construct the network of patients as networks are efficient data representation methods for multidimensional data. Therefore, using the walktrap clustering algorithm is a natural progression as it is based on the network and its topology.
Indeed, other methods could be used, however, it has been shown that networks (graphs) allow for more efficient data exploration than other approaches as they are capable of representing multidimensional relationships on a single 2D visual space (Benson et al., 2016). Ultimately, the question of which unsupervised method to use in each case is a scholarly question that requires its own dedicated research. Here, we focused on topological methods as we have successfully used them before to explore hemodynamics after SCI in animal models, a direct benchtobedside translation of both physiological variables and analytical methods (Nielson et al., 2015).
3) Are the results sensitive to the defined outcome of improvement of at least one AIS grade? What happens if this is increased to e.g. 2 grades?
This is a good question. Unfortunately, we do not have enough data to provide an accurate answer, despite having one of the largest cohorts assembled in SCI (N = 118, greater than 80% of acute SCI studies registered in Clinical trials.gov). Though increasing one AIS grade is clinically meaningful, we acknowledge the limitations of our study and agree with the reviewers that an analysis using a more sensitive outcome such as a continuous metric (e.g. motor/sensory scores) would increase the value of our findings. This is an area for future research that will benefit from the ongoing assembly of ever larger cohorts in the field.
4) Two different approaches to analysis were used – i.e. essentially clustering of some form and also logistic regression (using e.g. quadratic and spline functions). Can the authors comment on whether these 2 approaches can be used interchangeably or whether one would be preferred over to the other to answer the research questions of interest. What advantage does the clustering approach have in terms of the research question?
We have now clarified this in the results and discussion. These two methods are complementary rather than interchangeable. The patient network analysis is unsupervised and was used with the intention of exploratory discovery. It is through this method where we observed the potential double bound (upper and lower limits) for arterial pressure in humans, and the dependence of it to AIS at admission. The logistic regression is used in two ways, one to confirm the hypothesis formulated through the network analysis and two to find a precise MAP threshold through LASSO regularization. The fact that the MAP threshold for prediction of improvement captured through the network is very close to the one found by the LASSO indicates that our network methodology captures valid relationships in the data without previous knowledge of outcome (unsupervised). This suggests that our network method can be used for data discovery in future research including more physiological measures. We take the values obtained through LASSO as more confirmatory because it is a modeling approach, while the network analysis is exploratory. We added text in the manuscript clarifying these points
The use of topological analysis and clustering through networks has some advantages for exploratory analysis and discovery in the question at hand. We show that even though we are technically only including 2 variables in the analysis (MAP and HR), the time series of these two conform to a multidimensional manifold that we approximate by Isomap plus topological metrics. The network allows us to represent that multidimensional space in a 2D visual based on the similarity between patients. This is convenient for fast exploration of the multidimensional space, where relationship between patient similarities and any variable can be easily depicted by mapping the variable values to coloring the nodes. This and other related topological data analysis methods have been described as “Machine Intelligence” approaches as they aid the analyst with expert matter knowledge to perform fast discovery (Carlsson, 2009; Wasserman, 2018).
5) Why was a simple accuracy metric used and not e.g. AUC? Does the accuracy metric account for an imbalance in the outcomes?
We have now performed new analysis on prediction modeling. For these, we have calculated AUC and ROC. We have included new text and a new figure reflecting the new results.
6) The LOOCV accuracy was not that high suggesting a lot of other factors might influence outcomes. Is the accuracy really high enough to support the use of MAP being used by clinicians for decision making and/or interventions to control MAP during surgery?
We have included new analysis for the prediction modeling and added a subsection in the results. In these, we show that the prediction accuracy can be improved through model selection, where MAP is an important feature to maintain in the model for good prediction performance, especially for predicting AIS improvement or AIS A at discharge. It is important to note that all our metrics of performance are for LOOCV, which measures performances over unseen data for n (patients) trained models and thus offers an estimation of the generatability of the data and some protection against model overfitting. In that sense, LOOCV accuracy would be in most cases lower than accuracy of a model that has not been measured with crossvalidation, which is the case of most logistic regression in the clinical literature. With the new prediction modeling we included, our model has a LOOCV AUC of 0.75 for predicting AIS improvement at discharge. While we acknowledge the limitations of the study, we think that our prediction models are informative, and indeed are among the first ever to rigorously test ultraacute predictors of outcome in SCI. These models can provide a benchmark for future prediction modeling in the early level1 trauma clinical SCI.
7) What variables were used for the LASSO. The prediction accuracy again seems very low for highdimensional dataset.
We have added the range threshold analysis in its own section in the results and separated the explanation in the methods to address this point. LASSO was used to determine the most predictive MAP range and respective thresholds. The variables included are the time spent out of each one of the MAP ranges. In that sense, we have not used LASSO as supervised predictive ML method but rather as a feature selection method (through the regularization process) to obtain the last standing coefficient and correspondingly the best predicting MAP range. Therefore, we do not expect our use of LASSO to produce the most accurate predictive model but to provide the most predictive MAP range.
8) Only one dataset was used without splitting the data into a training and validation dataset. Are similar results for the topological network analysis obtained if the data is split for training and validation?
We agree with the reviewer about the need to validate the results and have added explicit discussion of this point. Topological network analysis learns the crossvalidated data manifold space, through a specific type of validation that subsamples the subspace of all networks. Splitting between validation and training in the exploratory topological analysis presents two major challenges: (1) splitting the data with our sample size can be problematic since would leave relatively few cases for training and testing with precision; (2) training and testing split is normally performed in supervised models to determine how well a trained model can predict unseen data. The topological network analysis process that we implemented is an unsupervised method used to explore the relationship between patients and physiological variables. In that sense, we are not using the network as a model to be “trained” but as a descriptive statistic of the shape of the data manifold. Nonetheless, the metrics of topological descriptions that we use (persistent homology) find the best configuration (most stable) given the data's many potential configurations (Carlsson, 2009). In that sense, we are validating the topological structure of the network for our specific case.
We conducted further validation by (1) performing confirmatory analysis through hypothesis testing (logistic models) and (2) using leaveoneout crossvalidation (LOOCV). LOOCV is perhaps the best alternative since splitting the data in this case would leave very few cases for training and testing which would ultimately underpower the analysis.
9) What was the modularity of the final network and does it suggest significant clustering?
The modularity of the final network was Q = 0.837. Modularity has a maximum of 1, representing complete clustering above expected in a random graph with the same structure (same number of nodes and degrees). Our final Q is quite high, suggesting the presence of clusters above what should be expected by chance. We have included this information in the corresponding figure legend.
10) Why was days from surgery to discharge used in the logistic regression models? Might it not be considered a mediator rather than a confounder – and how does its exclusion from the model influence the result?
We now include discussion of this important point.
We use time of surgery to discharge in the model as a covariate, to control for the potential effect in outcome related to the time to discharge. Considering time of surgery to discharge as a mediator is indeed a possibility. That would be true if: (a) time to discharge affects probability of improving (the outcome), and (b) MAP dysregulation during surgery affects days from surgery to discharge. The logical progression in the causal pathway would be as follows: higher MAP OR dysregulation indicates/causes changes in time to discharge, and time to discharge then affects probability of improvement (either increase or decrease). Since we see that MAP dysregulation (or proxies of it) actually reduces probability of improvement, then there are two options: (1) MAP dysregulation during OR directly or indirectly associates with longer time to discharge, and higher time to discharge reduce probability of improvement, or (2) MAP dysregulation associates with shorter time to discharge which in turns, is related to reduced probability of improvement. One would expect that the longer the patient is in the hospital, the more time the patient has to improve before discharge. However, the longer a patient stays in the hospital is probably an indicator of their severity. We have investigated this matter using different models to dissect the potential causal pathway.
First, we fitted a logistic model with AIS improvement as response and days to discharge as the only predictor. This resulted in a nonsignificant p value of p = 0.32, suggesting that days to discharge doesn’t 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^{2} was small (0.0217). Visual examination of the data indicates that the density of patients that tend to stay longer in the hospital are the ones with a MAP in the normotensive range ~ 80100. We also investigated whether days to discharge interacts with MAP and quadratic MAP to predict AIS improvement, with no significant results on the interaction, and the inclusion of the interaction do not change the quadratic association of MAP with the response. These suggest that these two factors do not moderate each other. We also performed the same analysis with time MAP out of 76104 mmHg instead of aMAP as a different proxy for MAP dysregulation, with similar results.
Finally, eliminating days to discharge from the full covariate model predicting AIS improvement doesn’t have a major effect on the model fit. A likelihood ratio test between both models shows a nonsignificant change in variance explained (p=0.79) with a deviance difference of ~ 0.1%.
Therefore, given the data, we think there is little evidence for considering an effect of MAP dysregulation on AIS improvement mediated (or moderated) by length of stay in the hospital. We have maintained days from surgery to discharge as covariate because although not significant, we still think it is important to adjust the model for differences in length of hospital stay between patients. We have included this analysis in the manuscript as results. In addition, we understand the limitation of the above analysis where all other variables and their relationship have not been considered. Nonetheless, we want to acknowledge the importance of dissecting the correlational observations as a first approximation to determine the causal effects. We will definitely explore these ideas in the near future as larger cohorts of data become available.
11) The limitations are mentioned but not discussed or justified. This leads to the following questions: a) Why was the lesion level not included in the analysis? and b) Why did the authors only analyze MAD values during surgery? Because the analysis of MAP data from the ICU period published elsewhere showed similar results regarding the lower limit of MAP, wouldn't it be of interest to know how much overlap there is between the populations with critical MAP values during surgery and during stay in the ICU?
We thank the reviewers for pointing to these important questions. We have expanded our analysis to include the level of injury. We found that the association of MAP range of AIS grade improvement is dependent on whether the patient had a cervical injury or not. At this time, we do not have enough samples to tackle this question with higher granularity (e.g. per segment level) as our data include very few cases with middle and lower level injuries. This is a limitation due to the infrequent nature of SCI, the difference in injury incidence per segment level, and the slow accumulation of patient cohorts, even at some of the busiest level 1 trauma centers in the country. We have included a new section in the results and expanded the discussion accordingly.
Regarding the analysis of MAP only during surgery, there are some considerations: First, our previous results in animal models suggest that OR hemodynamics is an important predictor of recovery. Second, this is a retrospective cohort in which systematic ICU monitoring is not available for answering the proposed question. Third, and most importantly, monitoring in the OR is perhaps the best controlled scenario for hemodynamic monitoring in these patients. We acknowledge the importance of the question of whether MAP during OR is sufficient as a proxy for patient dysregulation and how it confounds with hemodynamic dysregulation/treatment during ICU. As part of the TRACKSCI study, we are collecting highfrequency monitoring data in a prospective cohort and we will explore these questions in the future.
12) Introduction: Neither hypertension nor hypotension following acute SCI has been conclusively demonstrated to impact neurological recovery. Instead, guidelines and more recent work are based purely on observations and posthoc regression analyses. While the purported mechanism for repeated hypotensive episodes is clear, readers may benefit from at least a brief description of why both hypertension and hypotension could plausibly be important (aside from the fact that, again, noncausal observations demonstrated a relationship in the author's prior work).
We have included further details in the introduction describing the importance of both hypertension and hypotension from clinical and animal studies. We agree with the reviewers that a causal link between hemodynamic control and recovery in humans remains elusive, and it should be noted that randomized controlled studies are challenging in this context given the ethical concerns of such studies. Nevertheless, some experimental animal studies 'clamping MAP' at different set points strongly suggest causal effects (Nout et al., 2012), and physicianscientists in the field are conducting ongoing studies that appear strongly supports a potential causal role. Such a causal role is theoretically wellfounded, given that in traumatic brain injury and stroke, tissue perfusion/oxygenation, and brain tissue survival are supported by prevention of hypotension, but it is also possible that hypertension can induce 'hemorrhagic conversion' and bleeding into the tissue with devastating effects on lesion expansion. Spinal cord injury is known to involve similar secondary injury processes in both rodent and primate models (Crowe et al., 1997).
13) AIS scores: Based on Table 1 most patients were discharged within 2 weeks after injury. The neurological exam is not so reliable at this point. This is a big limitation of the current work. Although a six month followup would be ideal to determine whether neurological recovery occurred, the authors should at least mention this.
We agree with the limitations and have added some text in the discussion. These data are unique in that they reflect ultraacute data from the level 1 trauma center. This is a moment in the care pathway when a highly heterogeneous patient population is concentrated in one location prior to being dispersed to a broad variety of discharge scenarios (largely depending on insurance status in the US healthcare system). This makes long term followup data difficult to obtain. Nevertheless, as the TRACKSCI prospective study progresses, we are collecting contact information and have funds to track patients out to 12 months follow up together with high granular physiological data both in the OR and ICU. It will take many years to collect a large N, but in the future this dataset will provide a rich source of information to continue researching these important questions.
14) All the analyses seem to have been conducted on an AIS change. The authors should demonstrate that their analysis holds for a more linear measure of recovery (e.g., total motor score).
We agree with the reviewers that AIS change has limitations. We now acknowledge them in the manuscript. We have also added further analysis as a result of this revision. Specifically, we investigated models predicting AIS A and AIS D at discharge in addition to the AIS change analysis we had. Please see the following responses to details on the new analysis.
15) Based on Figure 2 – Supplement 2, it is difficult to ascertain whether clusters contain a higher proportion of individuals that show an AIS improvement, and those individuals tend to have a MAP > 80 and < 100, is due to these clusters having individuals who were less severe to begin with (i.e., C,D,E) and therefore less likely to be hemodynamically unstable. One way to answer this would be to examine AIS A patients in a separate analysis and determine whether these findings hold. Because, the alternative explanation here is that this analysis is effectively finding a proxy for initial injury severity (i.e., more severe, more hemodynamically unstable) – and not that hemodynamic instability per se is the problem. Another analysis that could help complement this work and avoid this confound would be to use total motor score as the outcome instead of AIS conversion.
This is an important consideration. We study the question of injury severity being a driver for the effect in several analyses. First, as suggested by the reviewers, we conducted the model fitting only in patients with AIS A at admission. This resulted in a significant linear MAP term, and a quadratic MAP with the second biggest coefficient and a p = 0.14. Comparing a linear model with the quadratic one in this cohort through a likelihood ratio test resulted in a p = 0.07. This suggests that while the linear term of MAP as predictor of AIS improvement is maintained significantly for AIS A patients, the evidence for a quadratic form is weak in this cohort. However, with the reduction in sample size (n=46) we may be underpowered for detecting such an effect. Second, we have included AIS at admission as a covariate in the full model to adjust for potential differences between different injury severities. In this scenario, the quadratic form of MAP predicting AIS improvement is maintained. All together, this suggests that if nonlinear relationship of MAP and AIS improvement depend on initial injury severity we lack the statistical power to fully investigate this question. We have included this analysis in the document pointing out the caveats of it. Future work should address this, as we agree with the reviewers will be important for guiding precision in clinical recommendations.
16) Logistic regression – Based on Figure 2p the overall trend looks more to be that higher MAP = > Pr(δ AIS grade > 0). The exception is only 2 data points on the top end. It is difficult to determine how robust the notion is that there is a 'too high' component to this data. Indeed it seems that a linear model does quite well for this analysis as well. Please see my comment below but this should be addressed as the concept of also having a 'top cutoff' is an extremely important clinical feature here.
We agree with the reviewers that having a ‘top cutoff” is an extremely important clinical feature. Although the linear model performs better than the unconditional model (no predictors), our results show that the polynomial quadratic model or the spline of degree 2 model outperforms the linear model, both in the variance explained and the LOOCV error. We interpret this as indicating that a nonlinear association between MAP and the prediction of AIS increase explains the data better than the linear model. We have expanded our explanation in the results to help the reader. Also, please see our answer below regarding the upper limit of the cutoff.
17) Time outside MAP – The authors use an approach that systematically increases their window in both directions to find their optimal range of 76104. However, what happens if you then hold 76 and only increase on the upper end? Does this rapidly degrade the relationship? If not, again this would suggest that the evidence for a top end cutoff is not as strong. While I understand the authors briefly looked at this (methods) it seems worth exploring further as this is a critical point for clinical management. I do not see a good reason that the time outside the threshold can not at least be plotted to determine this relationship.
We have included further details on the result section about the potential asymmetry of the range. While it is true that we see similar performance of the polynomial quadratic model (symmetric) and the spline model of degree 2 (asymmetric), we did not have strong evidence to justify either model as being 'more correct', so we chose simplicity as spline models will tend to overfit more than a polynomial. For that reason, we used the symmetric window increase in our study of the threshold range. As suggested by the reviewer, we performed the same analysis but holding the lower limit to 76 and increasing the upper limit by 1. This results in the most predictive range to be 76117. This suggests that our data and methodology is limited to tune the upper level and that it might be between 104 to 117. We have included this analysis in the results and discuss the implication of it.
18) Data availability – The code and analysis should be made available to the reviewers. It is impossible to determine the accuracy of this type of analysis without it.
We are adding the code that reproduces the entire analysis (including figures) as supplementary material, and the data has been uploaded to odcsci.org, a communitybased repository for SCI research data sharing and publication. The data would be available and active upon acceptance with DOIs (10.34945/F5R59J and 10.34945/F5MG68).
19) Discussion – It seems that the authors should discuss the confound of injury severity being linked with worse hemodynamic instability, and also worse neurological recovery. It would be helpful to include some of the suggested analyses to convince the readers that this confound does not explain the results since it is the most likely alternative explanation.
We have performed the suggested analysis when possible and discussed their results. While we agree that there are limitations in our study, we also contend that the findings with the new analyses are important. Although there are caveats to any observational clinical study, this is the largest study todate on this topic, with a sample size that is 4fold larger than the last study in this area. The work increases the body of evidence to guide both the clinical conversation and future prospective research.
https://doi.org/10.7554/eLife.68015.sa2Article and author information
Author details
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 TorresEspí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 (DEAC0205CH11231)
 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: DEAC0205CH11231 (DM).
Ethics
Human subjects: This study constitutes a retrospective data analysis. All data was deidentified before preprocessing and analysis. Protocols for retrospective data extraction were approved by Institutional Research Board (IRB) under protocol numbers 1107639 and 1106997.
Senior Editor
 Matthias Barton, University of Zurich, Switzerland
Reviewing Editor
 Arduino A Mangoni, Flinders Medical Centre, Australia
Reviewers
 Marcel Kopp
 Aaron Phillips, University of Calgary, Canada
Publication 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
© 2021, TorresEspín et al.
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.
Metrics

 607
 Page views

 103
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.