Coagulation phenotypes in sepsis and effects of recombinant human thrombomodulin: an analysis of three multicentre observational studies

Background A recent randomised trial showed that recombinant thrombomodulin did not benefit patients who had sepsis with coagulopathy and organ dysfunction. Several recent studies suggested presence of clinical phenotypes in patients with sepsis and heterogenous treatment effects across different sepsis phenotypes. We examined the latent phenotypes of sepsis with coagulopathy and the associations between thrombomodulin treatment and the 28-day and in-hospital mortality for each phenotype. Methods This was a secondary analysis of multicentre registries containing data on patients (aged ≥ 16 years) who were admitted to intensive care units for severe sepsis or septic shock in Japan. Three multicentre registries were divided into derivation (two registries) and validation (one registry) cohorts. Phenotypes were derived using k-means with coagulation markers, platelet counts, prothrombin time/international normalised ratios, fibrinogen, fibrinogen/fibrin-degradation-products (FDP), D-dimer, and antithrombin activities. Associations between thrombomodulin treatment and survival outcomes (28-day and in-hospital mortality) were assessed in the derived clusters using a generalised estimating equation. Results Four sepsis phenotypes were derived from 3694 patients in the derivation cohort. Cluster dA (n = 323) had severe coagulopathy with high FDP and D-dimer levels, severe organ dysfunction, and high mortality. Cluster dB had severe disease with moderate coagulopathy. Clusters dC and dD had moderate and mild disease with and without coagulopathy, respectively. Thrombomodulin was associated with a lower 28-day (adjusted risk difference [RD]: − 17.8% [95% CI − 28.7 to − 6.9%]) and in-hospital (adjusted RD: − 17.7% [95% CI − 27.6 to − 7.8%]) mortality only in cluster dA. Sepsis phenotypes were similar in the validation cohort, and thrombomodulin treatment was also associated with lower 28-day (RD: − 24.9% [95% CI − 49.1 to − 0.7%]) and in-hospital mortality (RD: − 30.9% [95% CI − 55.3 to − 6.6%]). Conclusions We identified four coagulation marker-based sepsis phenotypes. The treatment effects of thrombomodulin varied across sepsis phenotypes. This finding will facilitate future trials of thrombomodulin, in which a sepsis phenotype with high FDP and D-dimer can be targeted. Supplementary Information The online version contains supplementary material available at 10.1186/s13054-021-03541-5.


Background
The global rate of sepsis-related mortality remains high and the annual age-standardised mortality owing to sepsis is 148.1 deaths per 100,000 of the global population [1]. Recombinant human thrombomodulin (rhTM) has anti-inflammatory and anticoagulation activities [2], and it has been suggested as an adjunct therapy for patients with sepsis, particularly those with sepsis-induced coagulopathy [3]. Nevertheless, a recent phase III randomised controlled trial revealed no beneficial effect of rhTM in patients with sepsis-induced coagulopathy [4]. This result can be explained by the heterogeneity of patients with sepsis and inappropriate criteria of coagulopathy [5] using the prothrombin time/international normalised ratio (PT-INR) and a platelet count based on subgroup analysis of an international phase II trial of rhTM [6].
Sepsis is a highly heterogeneous syndrome with variable aetiology and pathophysiology [7]. Thus, a specific therapy may benefit some, but not all, patients with sepsis. Several recent studies have classified sepsis into several phenotypes with distinct characteristics using cluster analysis [8][9][10], an unsupervised machine learning method that can identify relatively homogeneous groups in a heterogeneous population [11]. Furthermore, these studies indicated that specific therapies conferred benefits only in patients with specific phenotypes of sepsis [8][9][10].
Identifying the target patients receptive to rhTM treatment through grouping based on biomarker cut-offs is challenging. To address this issue, we examined latent sepsis phenotypes in terms of coagulopathy and identified which phenotypes would benefit from rhTM using machine learning approaches.

Study design and settings
Details of the methods and analytical processes in the present study are provided in the Supplemental Digital Content. This was a secondary analysis of the following multicentre registries: the Japan Septic Disseminated Intravascular Coagulation (JSEPTIC-DIC) study (UMIN-CTR ID: UNIN000012543) [12], Tohoku Sepsis Registry (UMIN000010297) [13], and Focused Outcomes Research in Emergency Care for Acute Respiratory Distress Syndrome, Sepsis, and Trauma (FORECAST) sepsis study (UMIN000019742) [14]. All three registries include information on consecutive patients admitted to ICUs for severe sepsis or septic shock [15,16]. Briefly, the JSEPTIC-DIC study retrospectively reviewed data derived from 3195 consecutive patients with severe sepsis or septic shock, aged ≥ 16 years, admitted to 42 ICUs at 40 institutions in Japan, between January 2011 and December 2013 [12]. The Tohoku Sepsis Registry prospectively registered 616 consecutive patients who were admitted to ICUs with severe sepsis, or those who developed severe sepsis after admission to the ICUs or general wards at 10 institutions (three university hospitals and seven community hospitals) in the Tohoku District, Northern Japan, between January 2015 and December 2015 [13]. The multicentre prospective FORE-CAST sepsis study included 1184 consecutive patients aged ≥ 16 years, who were admitted to 59 ICUs in Japan with severe sepsis according to the sepsis-2 criteria [15], between January 2016 and March 2017 [14]. These studies were approved and the need for informed consent was waived by the institutional review boards at the participating hospitals.

Definitions
Ventilator-free days were defined as the number of days on which a patient did not require mechanical ventilation during the initial 28 days following enrolment. The number of ventilator-free days of patients, who died within day 28, was assigned as 0. ICU-free days were calculated similarly.

Statistical analysis Analytical cohorts
We derived sepsis phenotypes from the JSEPTIC-DIC study (n = 3195) and Tohoku Sepsis Registry (n = 499) and validated the phenotypes using the FORECAST sepsis study (n = 1184).

Cluster derivation
We initially assessed the distribution and missingness in phenotyping variables. Non-normal data were log-transformed and scaled. Patients without 28-day mortality information were excluded. Missing data were imputed using the random forest method for each study cohort with the missForest package [17]. Random forest imputation is a nonparametric algorithm that accommodates nonlinearities and interactions and does not require the specification of a specific parametric model [18]. This approach generated single-point estimates by random draws from independent normal distributions centred on conditional means predicted by random forest. Random forest applies bootstrap aggregation of multiple regression trees to reduce the risk of overfitting and combines estimates from many trees [17]. Missingness was imputed using patient characteristics, laboratory data, outcomes, and other covariates, including in-hospital management (Additional file 1: Supplemental documents).
We applied k-means with Euclid distance, which is a basic and widely used machine learning-based clustering approach, to derive sepsis phenotypes [9,11]. We then determined the optimal number of clusters using a consensus clustering approach that assessed both quantitatively and visually to estimate the number of unsupervised classes in a dataset by inducing sampling variability with sub-sampling [19]. In consensus clustering, we evaluated the separation of consensus matrix heatmaps using the elbow method, cumulative distribution function, and cluster-consensus plots. We visually evaluated the clustering using t-Distributed Stochastic Neighbour Embedding (t-SNE) for reducing dimensionality and visualising high-dimensional datasets [20]. We also derived phenotypes using a divisive hierarchical clustering approach as an alternative to k-means, for confirming the cluster consistency. The number of clusters was determined using the dendrogram and the elbow and gap statistic methods [21].

Evaluation of rhTM effects in derived phenotypes
We examined the associations between rhTM and the clinical outcomes in each derived cluster, using a generalised estimating equation to adjust for hospital-level variance. We analysed the associations after adjusting for the potential confounders of age, sex, comorbidities, and sequential organ failure assessment (SOFA) scores. In the derivation cohort, we did not adjust for the management before and after admission to the hospitals because the information on when each management was initiated was not available. We examined the cluster-level effects of rhTM by including the interaction term, rhTM use -xcluster, in the model to examine different effects of rhTM across clusters. In addition, to confirm the robustness of the association of interest, we applied a Bayesian regression model to assess the associations between rhTM and the clinical outcomes for each derived cluster, based on k-means in the derivation cohort [21,22]. Bayesian regression was performed using a Markov chain Monte Carlo procedure with four chains of 2000 iterations per chain. The results are shown as beta coefficients with 95% credible intervals and displayed as odds ratios with 95% credible intervals for simplicity.

Cluster validation and evaluation of rhTM effects
We predicted patient phenotypes in the FORECAST sepsis study as external data based on coagulation markers of clusters in the derivation cohort (derived from JSEP-TIC-DIC and Tohoku Sepsis Registry). Predictions arose from the Euclidean distance from each patient to the centroid of each FORECAST phenotype. In each predicted cluster in the FORECAST sepsis study, we first described the frequency and clinical characteristics of the clusters. Thereafter, we used a generalised estimating equation to account for patient clustering within hospitals to assess associations between rhTM and clinical outcomes in each predicted cluster in the external data. The adjusted variables were age, sex, comorbidities, SOFA scores, and in-hospital management, including renal replacement therapy, and treatment with steroids, intravenous immunoglobulin, antithrombin, and vasopressors. As the FORECAST sepsis data included information on the time of management, we included management before or after admission as a covariate to estimate the effects of rhTM on clinical outcomes. For sensitivity analyses, we used a generalised estimating equation, applying the acute physiology and chronic health evaluation (APACHE II) score and source of infection as adjusted variables, instead of the SOFA score. The source of infection was categorised into respiratory, abdominal, skin and soft tissue, urinary tract, and others. As in the derivation cohort, we analysed the association by Bayesian regression with a Markov chain Monte Carlo procedure with four chains.

Patients in the derivation cohort
We excluded 117 patients without 28-day mortality information in the derivation cohort from the two multicentre registries, leaving 3694 patients who were eligible for analysis (3195 from JSEPTIC-DIC and 499 from the Tohoku Sepsis Registry). Table 1 summarises the patients' characteristics. The median age was 72 years, and 40% of patients were female. Overall, rhTM was administered to 26.2% of patients (the distribution of the proportion of patients in the institutes is shown in Additional file 3: Figure S1). The in-hospital mortality and 28-day mortality rates were 32.1% and 20.4%, respectively.

Derivation of clinical sepsis phenotypes
We assessed the distributions and missingness among the phenotyping variables. Antithrombin activity was the most lacking variable, 51.3% and 45.8% in derivation and validation datasets, respectively ( Table 2). According to clustering using k-means, a four-class model including the phenotype clusters derivation dA, dB, dC, and dD ("d" represents "derivation") may be an optimal fit. The heatmap matrix (Additional file 3: Figure S2), elbow method (Additional file 3: Figure S3), and cumulative distribution function curve (Additional file 3: Figure S4), indicated that the four-class model was optimal, whereas the cluster-consensus plot suggested that two, three, or four clusters were optimal (Additional file 3: Figure S5). The four-class model was supported by the t-SNE plot with clear separation (Fig. 1). Additional file 3: Figure S6 shows a cluster dendrogram obtained using a divisive hierarchical clustering approach. The elbow method showed that a two-or four-cluster model is optimal (Additional file 3: Figure S7), whereas the gap statistic method [21] showed that the four-cluster model was optimal (Additional file 3: Figure S8).
Patients in cluster dA were likely to have a severe physiological status and organ dysfunction (high APACHE II and SOFA scores), coagulopathy (low platelet counts, prolonged PT-INR, low fibrinogen, and extremely high FDP and D-dimer levels), high lactate levels, and high mortality (Table 1). Approximately 90% of patients in this cluster required vasopressors. The characteristics of patients in cluster dB were similar to that in cluster dA in terms of severity, but likely to have abdominal infection with normal white blood cell counts, moderate coagulopathy with moderate FDP and D-dimer levels, and low antithrombin activity. Patients in clusters dC and dD had moderate and mild disease, respectively. Although patients in cluster dC had coagulopathy with high FDP and D-dimer levels, those in cluster dD were likely to have respiratory infection without coagulopathy. The phenotypes were similar according to the four-cluster hierarchical clustering (Additional file 2: Table S1).

Characteristics of phenotypes in the validation cohort
Additional file 2: Table S4 shows the patients' characteristics in each cluster in the validation cohort. The median age was 73 years, 40% of the patients were women, and rhTM was administered to 21.2% of patients. In-hospital and 28-day mortality rates were 23.4% and 19.0%, respectively. These characteristics were similar to those in the derivation cohort but the rate of rhTM treatment and mortality were relatively lower.
We used only coagulation markers to predict clusters in the validation cohort, and the characteristics were similar to those in the derivation cohort ("v" represents "validation"). Similar to the patients in cluster dA, those in cluster vA were likely to have a severe physiological status and organ dysfunction (high APACHE II and SOFA scores), coagulopathy (low platelet counts, prolonged PT-INR, low fibrinogen, and extremely high FDP and D-dimer levels), high lactate levels, and moderately high mortality. Patients in cluster vB had a high mortality rate with moderate coagulopathy and moderate FDP and D-dimer levels. Patients in clusters vC and vD had moderate and mild disease, respectively. Patients in cluster vC had coagulopathy with high FDP and D-dimer levels, whereas those in cluster vD did not have coagulopathy.    Table 3). In contrast, rhTM was not associated with better outcomes in the other clusters. As secondary outcomes, rhTM was associated with increased number of ventilator-free days (6.7 days [95% CI 0.8-12.7 days]) in Cluster vA, but not with the number of ICU-free days or discharge location in any of the clusters (Additional file 2: Table S5). In the sensitivity analyses, rhTM was associated with better outcomes, when APACHE II score and source of infection were applied as adjusted variables, instead of the SOFA score (adjusted RD, − 0.24% [95% CI − 0.45 to − 0.02%] for 28-day mortality; RD − 0.26% [95% CI − 0.47 to − 0.04%] for in-hospital mortality; Additional file 2: Table S6). The associations between rhTM, and 28-day and in-hospital mortalities were consistent with the findings of the Bayesian regression analysis (Additional file 2: Table S3 and Additional file 3: Figure S9).

Discussion
This secondary analysis of the sepsis registries identified four phenotypes with various coagulation features among patients with severe sepsis. Treatment with rhTM was associated with lower in-hospital mortality rates only in the phenotype with severe coagulopathy, characterised by low platelet counts, extremely high levels of FDP and D-dimer (phenotype clusters dA and vA), and severe organ dysfunction. These results were not identified in the other phenotypes. The severity of coagulopathy is defined by the DIC scoring systems, such as the International Society on Thrombosis and Haemostasis (ISTH) scoring system for diagnosing overt DIC [23] and Japanese Association for Acute Medicine DIC scoring system [24], both of which have been applied in many studies. The difference between these systems and machine learning-based clustering is the use of a trivial cut-off. Table 1, Additional file 2: Tables S1, and S4 show that each phenotype cluster included patients with various ISTH DIC scores without a clear cut-off that overlapped with the other clusters. This suggests that clustering based on machine learning can detect novel phenotypes that cannot be identified using the conventional scoring systems.
Recombinant human thrombomodulin has anticoagulation effects and was shown to be beneficial for patients Six coagulation markers (bold font) were used for clustering APACHE acute physiology and chronic health evaluation, DIC disseminated intravascular coagulation, FDP fibrinogen/fibrin degradation product, IQR interquartile range, PT-INR prothrombin time-international normalised ratio, SIRS systemic inflammatory response syndrome, SOFA sequential organ failure assessment, WBC white blood cells Variables with asterisk (*) were potential confounders that were adjusted in a generalised estimating equation. *P between clusters  with sepsis and coagulopathy in observational studies and in a subgroup analysis of a phase II trial [3,6]. The latest phase III trial focused on patients with sepsis with cardiovascular or respiratory dysfunction as well as coagulopathy according to subgroup analysis of the phase II trial [4]. However, the phase III trial did not identify a positive effect of rhTM on survival, suggesting that differentiating a subgroup that may benefit from rhTM is difficult using conventional methods with clear cut-offs. In our study, despite overlapping characteristics, various DIC scores, and differences in severity among clusters, cluster vA (dA) was the only phenotype in which rhTM was associated with better survival outcomes. This suggests that machine learning clustering can identify optimal clinical phenotypes for rhTM treatment. Additionally, the machine learning clustering described herein used only six variables, all of which are general markers that can be measured in most hospitals. Additionally, the results can be available soon after admission before deciding to administer rhTM in an emergency room or ICU.
Other studies using machine learning-based clustering for patients with sepsis also suggest that several specific therapies have beneficial effects only in patients with specific phenotypes. A Toll-like receptor 4 antagonist, protocol-based resuscitation, activated protein C, and fluid input affected each phenotype differently [9,10]. The effectiveness of rhTM also varied across phenotypes in our study. Therefore, selecting an optimal clinical phenotype may be key to the success of a specific therapy for patients with sepsis. Including entire populations with sepsis may explain why previous randomised trials found no beneficial effects of adjunctive therapies [25][26][27][28]. The goal of precision/tailored medicine is to select the optimal therapy for patients, for which machine learning-based clustering can be effective. Although our study does not fully address the definite endotypes of coagulation in sepsis biologically or pathophysiologically, our findings improve the understanding of the true endotypes of sepsis with coagulation.

Limitations
This study had several limitations. We used three registries that included different variables. Therefore, unmeasured confounders, and a lack of information such as the timing of rhTM administration may have biased our findings. Nevertheless, the data included detailed clinical information that is generally used for adjustment, and the results in the validation cohort accounted for management before or after admission. The number of missing variables for clustering was not small; therefore, missingness may have limited our findings. However, missing data imputation using the random forest approach is considered valid and valid imputation reduces bias, even when the proportion of missingness is high [29]. Our data did not include information on the long-term outcomes (i.e. 6-/12-months mortality) and SOFA score at several weeks after admission, although long-term outcomes are also important. Our data did not include the duration of rhTM administration, which was presumably 6 days, according to the generally prescribed dose and duration in Japan. We could not evaluate whether the phenotypes and efficacy of rhTM are applicable for patients with sepsis defined by the Sepsis-3 criteria [30], as three observational studies enrolled patients with sepsis using the Sepsis-2 definition [16], and the datasets did not include SOFA scores before admission. Heparin is commonly used for anticoagulant therapy, worldwide. However, we could not include heparin treatment for sepsis-induced coagulopathy in the model, because two of three registries did not collect the data. Only 5% (167/3195) of the patients were treated with heparin for coagulopathy with sepsis (excluding use for venous thromboembolism) in JSEPTIC-DIC study [12]. We need to develop a model to determine the phenotypes of individual patients to be able to perform a clinical trial in the future. Finally, our data were derived from Japanese patients; therefore, the generalisability of the results may be limited.

Conclusions
The findings derived using machine learning clustering indicated that rhTM can benefit only patients with a severe coagulopathy phenotype. Identifying patients for whom a therapy will have a beneficial effect can lead to precision/tailored medicine in critical care. To achieve this goal, the accuracy of phenotyping should be increased by analysing more patients, and through further validation. A randomised trial focusing on suitable phenotypes determined by effective phenotyping is warranted.