Longitudinal patterns of cytokine expression at the individual level in humans after laparoscopic sleeve gastrectomy

Abstract The study of the human response to injury has been hampered by the inherent heterogeneity in the models and methods used. By studying a standard injury longitudinally, using individual patient‐level analysis, we endeavoured to better describe its dynamics. We analysed clinical variables, clinical laboratory and plasma cytokines from 20 patients at five time points. Clustering analysis showed two prototype patterns of cytokine behaviour: a concordant type, where cytokines behave the same way for all patients (notably IL‐0 and TNFα), and a variable type, where different patterns of expression are seen for different patients (notably IL‐8, IL‐6 and IL‐1RA). Analysis of the cytokines at the individual patient‐level showed a strong four‐way correlation between IL‐1RA, GCSF, MIP‐1β and MCP‐1. As it holds for most patients and not just on average, this suggests that they form a network which may play a central role in the response to gastro‐intestinal injuries in humans. In conclusion, the longitudinal analysis of cytokines in a standard model allowed the identification of their underlying patterns of expression. We propose that the two prototype patterns shown may reflect the mechanism that separates the common and individual aspects of the injury response.


| 6623
TRAHTEMBERG ET Al. our incomplete understanding of the connection between local injury and the systemic inflammatory response. 2 While several previous studies have analysed cytokines and their correlation with the systemic inflammatory response, sepsis and trauma, several shortcomings have raised questions regarding the reliability of some of these studies, such as the study of heterogeneous populations and injuries that then underwent varying iatrogenic interventions. Some studies used a few cytokines or only single mediator classes; cytokines act in a network; thus, analyzing only a small number of cytokines is of limited value. Another important reason for the lack of consensus in the human cytokine literature relates to the significant variability in the sampling times relative to the injury, given the fact that the acute cytokine responses are rapid and by 48 hours most of the changes are over. 5,6 Also, the changes from baseline within cytokine networks appear to be very important, 1,7 making the sampling strategy crucial in order to observe these changes. All these problems raise concerns that many of these studies would be limited in their ability to reflect the pathological process at hand (for a representative sample see Ref. [8][9][10][11][12]. While these past efforts have been pioneering, the analysis of such complex scenarios with the methods used has led to limited success when trying to collectively generalize and synthesize them. In order to try and deal with the heterogeneity of the immune response in humans, previous efforts have tried to tackle this by looking for response subpopulations with the use of large scale investigations followed by multivariate analysis in stable conditions, 13 with the use of in vitro approaches, 14 or using pattern recognition for defined disease 'snapshots'. 15 Still, acute, dynamic diseases in vivo in humans are much more difficult to study. Important advances have been made trying to account for the individual variability inherent in such cases by relying on data from clinical trials (for example see Ref. 16). While randomization allows these trials to assume heterogeneity is equally partitioned, the heterogeneity itself is very high (patients, treatments, disease progression, sampling time, etc).
On another front, several animal studies have demonstrated that analysis of time-dependent cytokines patterns can provide important insights into the mechanisms of inflammation. 17 However, being single genotype and phenotype studies, these studies lack the variability encountered in humans, and their utility for clinical application is not clear. [18][19][20] In order to overcome some of the problems described, we set out to study a standard injury in humans. Laparoscopic sleeve gastrectomy (LSG) represents such a standard tissue injury model of the inflammatory response to abdominal gastro-intestinal injury. It reproduces the same inflammatory stimulus in a relatively homogenous group of patients who begin at a steady state. In essence, it allows to study a clinical event with conditions that approach a controlled laboratory experiment. Previous studies have examined the long term changes in inflammatory markers following bariatric surgery, [21][22][23] yet none of these studies examined the changes in inflammatory markers and cytokines in the perioperative period, nor was the analysis done in a personalized manner. Also, the number of cytokines previously examined was limited, and they were not clinically correlated. In this study, we aimed to analyse the cytokine patterns and their development over time, at the individual level, as well as their correlation to the patients' clinical characteristics and outcome.

| Experimental design
We performed a prospective, descriptive and analytic study in a group of 20 patients with morbid obesity who underwent laparoscopic sleeve gastrectomy (LSG) as a model for the inflammatory response to abdominal injury. It was performed at the Hadassah Medical Center in Jerusalem, Israel, over a nine-month period.
Informed consent was obtained from all participants prior to inclusion in the study. The study and protocol were approved by specifically, multiple unsuccessful attempts of weight loss and a body mass index (BMI) greater than 40 kg/m 2 , or greater than 35 kg/m 2 if there are associated comorbidities. Exclusion criteria included: refusal or withdrawal of consent, failure to obtain blood samples and baseline haemoglobin level less than 9 gr/dL. Previous studies have shown that even among patients after major trauma or severe illness, most of the dynamics in cytokine levels occurred in the first 48 hours 5,24-27 ; therefore, sampling must be performed early. Accordingly, the times were chosen to try and sample the expected changes in mediator blood levels, together with the logistic constraints of drawing blood from patients. A 6-mL blood sample was obtained from each patient at each of the following five time points: (a) pre-operatively, immediately after insertion of an intravenous line (defined as 0 hours); (b) immediately postoperatively, upon arrival at the recovery room; (c) three hours after the 2nd sample; (d) the morning of post-operative day one (POD 1) (24 hours); and (e) the morning of POD 2 (48 hours). Every blood sample included 2 tubes: an EDTA tube for Luminex assays and a heparin tube used for clinical biochemistry. Blood samples used for Luminex assay were immediately stored at 4°C. The samples were centrifuged within 20 minutes (5 minutes centrifugation time at 300 g relative centrifugal force at a temperature of 4°C.), and the plasma was frozen. Samples 1 to 3 were centrifuged and frozen at −20°C in the recovery room with on-site equipment, and then transferred on the same day to the laboratory for permanent storage at −80°C. Samples 4 and 5, which were taken in the ward, were directly centrifuged and stored at −80°C in the laboratory. In order to minimize degradation and standardize results, all samples underwent up to 2 freeze and thaw cycles, with each specific cytokine measured after the same number of cycles for all samples. Cytokine data are missing for one patient due to technical problems. Clinical data were recorded every time a blood sample was obtained, including systolic blood pressure (SBP), diastolic blood pressure (DBP), heart rate (HR), oxygen saturation (SaO 2 ) and fraction of inspired oxygen (FiO 2 ), body temperature, respiratory rate (RR) and visual analog scale (VAS, 0-10) pain scale.
Patients were observed for surgical complications throughout the admission, such as surgical site and IV cannula site infections.

| Surgical technique
All surgeries were performed by the same surgeon (RE), using an identical technique. The average surgical time (defined from start to end of anaesthesia care) was 71 minutes, with a SD of 21 minutes (range 47 to 135 minutes). This SD was driven by two cases where there were unexpected delays after starting anaesthesia care before the actual surgery began, when excluded, the SD drops to 14 minutes and the average to 66 minutes. The surgical approach was laparoscopic with the use of four trocars. Gastrectomy was performed after releasing the major gastric curvature, from a distance of 5 cm from the pylorus to the gastro-oesophageal junction, adjusted to a 42-Fr intra-gastric bougie catheter. The staple line was tested for homeostasis and after that for leakage using intra-gastric methylene blue, after pylorus obstruction. The specimen of the major gastric curvature was removed through the trocar incision. Haemostasis was secured, and all trocars removed under vision. An intraperitoneal suction drain was placed in all patients and was routinely removed 3 days after surgery.

| Data collection and storage
A paper case report form was kept for each patient, containing clinical information including demographic data, identification information, age, gender, past medical history and current active medical problems, as well as medications and vital signs throughout the perioperative study period. This and all other study data were then recorded digitally using de-identified serial numbers into an Excel

| Cytokine measurements
The Human High Sensitivity Cytokine Premixed kit (R&D systems, Inc) was used to determine the level of IL-1β, IL-2, IL-6, IL-8, IL-10 and Histone H3 was measured using ELISA (Mybiosource.com). All assays were performed according to manufacturer's instructions. The assay plates were read using a Luminex analyzer (MagPix, Luminex Corporation).

| Clinical biochemistry
Serum levels of sodium, potassium, creatinine, urea, lactate dehydrogenase (LDH), ferritin, lactate and creatine phosphokinase (CPK), as well as a complete blood count (CBC) were obtained using the hospital's clinical laboratory services.

| Statistics
Power analysis: Biancotto et al 28 investigated the inter-individual variation of cytokine levels. Based on these data and the expectation of a fourfold increase of cytokine levels (based on our own preliminary data), the minimal sample size needed to detect a statistically significant difference (P < .001) with a power of 0.95 is 10 individuals. The low threshold for the alpha error is necessary due to correction for multiple tests, and the sample size was doubled to account for biological and methodological heterogeneity.
The raw data were prepared as follows: (a) values that were below or above the limit of detection of the assay were set to the corresponding limit of detection. (b) Outliers were defined as 3 interquartile ranges (IQR) above the 75th or below the 25th quartile sampling times as fixed effects. Statistical significance was tested using the Wald test for the random effects' covariance estimates and the F test for the fixed effects. Comparison between sampling times in the model was performed using Tukey's HSD test. Clustering was done using the self-organizing map which is an unsupervised clustering algorithm that bears resemblance to k-means clustering and takes in account the shape of the data being clustered and not just the proximity between the points (ie tries to preserve the topological properties of the original data). 29 This is useful in order to cluster together individuals with similar trajectories over time. Clustering algorithms require pre-specification of the number of clusters to be fitted. In order to decide how many clusters are underlying the data, we chose the largest possible number of clusters that yielded more than one patient per cluster, as a cluster of one patient is trivial statistically. In most cases, this provided three clusters; therefore, for simplicity and clarity, throughout the article we used 3 clusters F I G U R E 1 Vital signs and clinical laboratory trends. The distribution of the different variables is shown over the five sampling times for all patients (n = 20). In the left column, they are shown as box plots. The box is bound by the 3rd quartile at the top and the 1st at the bottom, encompassing the interquartile range; the line inside the box represents the median; the whiskers are drawn extending 1.5 times the interquartile range from the top and bottom of the box; outliers beyond this range are shown as dots. In the right column, all the measurements are shown as dots, with a fitted curve passing through the averages, connecting them; the whiskers encompass the 95% confidence interval of the mean. The units are as follows: blood pressure in mmHg, heart rate (HR) in beats per minute, temperature (Temp) in degrees Celsius, glucose in mg/dL and white blood cells count (WBC) in 10 9 /L F I G U R E 2 Distribution of IL-8, IL-10 and MMP-3 over time. The distribution of the different cytokines is shown over the five sampling times for all patients (n = 19). In the left column, they are shown as box plots. The box is bound by the 3rd quartile at the top and the 1st at the bottom, encompassing the interquartile range; the line inside the box represents the median; the whiskers are drawn extending 1.5 times the interquartile range from the top and bottom of the box; outliers beyond this range are shown as dots. In the right column all of the measurements are shown as dots, with a fitted curve passing through the averages, connecting them; the whiskers encompass the 95% confidence interval of the mean. The values of the cytokines have been row normalized using the median of each patient's five samples, in order to compare between patients and across different experiment batches for all parameters. For clustering, the data were prepared using the Johnson transformation. For the correlation analysis at the individual level, we tested all cytokine pairs, separately for every patient, using Spearman's ρ. Then, the P-values of all patients for a given pair were combined using Fisher's method. In order to control for multiple testing, we used Benjamini and Hochberg's false discovery rate correction with a Q level of 0.05. This was repeated for the clinical laboratory and vital signs together versus the cytokines.

| RE SULTS
The demographic characteristics of the study population were as fol-  Figure S1. We performed mixed model analyses of the variables to discern which component of the responses was due to variabilities in the patients (a random effect) and which component was due to a common response to the procedure, expressed by the progression over sampling times (a fixed effect). Table 1 Figure S2.
As the whole cohort analysis could not explain the cytokine results, we set out to analyse the results differently, at the individual level. As can be seen in Figure 3   F I G U R E 3 SBP, WBC and creatinine clusters. The parameters were plotted so that the 'Y' axes represent the absolute values measured and the 'X' axes represent the five sampling times. Every line represents a single patient, and all patients were included (n = 20). The red, green and blue plots represent the three different clusters found, numbered 1 to 3, for every parameter. Each parameter was clustered separately; therefore, the patient allocation in the clusters is not the same for the different parameters (ie the patients in cluster 1 for the SBP are not necessarily the same as in cluster 1 of the WBC). The cluster means represent the averaged response of each cluster. The SBP is measured in mmHg, the WBC in 10 9 /L and the creatinine in µmol/L |   (Figure 4 and Figure S4).
Using this classification, we found no consistent trends among any given patient, including when analyzing separately pro-inflammatory and anti-inflammatory cytokines. The small sample size precludes the analysis of subgroups of response types.
We also set out to determine the associations between the results. Table 2 Table 3 shows the statistically significant associations between the cytokines themselves, using the same technique as for Table 2. This is a description of the immune reaction elicited by the surgery. Of special note is the quadrilateral correlation between IL-1RA, GCSF, MIP-1β and MCP-1; each one is related to the other pairwise, creating a network.

| D ISCUSS I ON
One of the challenges of personalized medicine lies in identifying molecular targets that show variations between different patients, which are amenable for intervention, and that are likely to affect the clinical course of the patient. Almost by definition, this requires F I G U R E 4 IL-10, TNFα and IL-8 clusters. The parameters were plotted so that the 'Y' axes represent the median-normalized cytokine values, per patient (as in Figure 1), and the 'X' axes represent the five sampling times. Every line represents a single patient, and all patients were included (n = 19). The red, green and blue plots represent the three different clusters found, numbered 1 to 3. Each cytokine was clustered separately; therefore, the patient allocation in the clusters is not the same for the two cytokines (ie the patients in cluster 1 for IL-10 are not necessarily the same as in cluster 1 of IL-8). The cluster means represent the averaged response of each cluster. Normalization per patient allows comparing the fold change in expression over time between the different patients but does not show absolute levels | 6631 analysis of the data at the individual level. In this investigation, we tried to advance that goal by studying the cytokine, vital signs and laboratory responses to a standard injury (LSG) and analyzing the results from a personalized perspective.
When studying the vital signs and clinical laboratory results, the regression model was able to explain a substantial part of the variability, yet provided little insight into the processes at hand.
The cluster analysis, in contrast, uncovered a new layer of information that would have been lost otherwise, reflecting the phys- We would like to discuss the cytokine clustering results in three respects:

| IL-10 and TNFα
Their similar kinetics, but with opposite changes, play out in a way that clearly outlines our understanding of the expected immune response to injury. While these results might seem unsurprising given the current body of knowledge, 32,33 it should be noted that this knowledge comes largely from animal and basic science models. 34,35 In this study, we have been able to show IL-10 and TNFα's kinetics in detail, in response to an actual injury in humans. This demonstrates the actual occurrence in vivo of the behaviour that is expected, based on pre-clinical or partial models. Although others have shown similar kinetics in other human models, these have been averaged results; our analysis was performed at the individual level. This means that the changes we show in IL-10 and TNFα are not averaged courses, but rather that all patients respond in the same way; this carries more weight when trying to establish a universal description of the response to injury. In that sense, these results support the hypothesis that the anti-inflammatory arm is an integral part of the immune response to injury, beginning together with the inflammatory arm and not in reaction to it. 36

| Individual-level analysis
The past three decades have seen many studies on cytokines. Many surgery for all patients together (Figure 2), the conclusion would be that there is considerable variability, and, in average, there is no change between those two points. This is very different from the conclusion drawn from the results of our individual-level analysis, which shows that there are 3 different subgroups of responses ( Figure 4).

| Concordant-type vs personal-type patterns
When studying immune responses to in vivo, one of the major challenges is that, while the macroscopic responses mounted are largely similar, there is significant inter-personal variability at the molecular and cellular level. How can these be reconciled? Recent results, although in vitro, would seem to suggest that genetics can only explain a minority of the heterogeneity. 40 One conclusion from our results is that this heterogeneity seems to have boundaries, as the cytokines are still constrained in their expression dynamics. They do not behave randomly, as their patterns fall within certain forms, as shown in Figure 4 and Figure S4. the environment (such as season of the year). 42 Our results are another way to understand the mechanism of the coexistence of heterogeneity on top of a common injury response.
The correlation analysis presented was also performed in a way that preserves each patient's dynamics, as well as assessing these correlations over time, not just at single time points. We show that there is a core of four cytokines (IL-1RA, GCSF, MIP-1β and MCP-1) that stand out in their reciprocal correlation. It has been shown that cytokines do not act on their own but rather in synergism, forming networks among themselves (as well as with cells and other mediators). 43,44 Other works that analysed cytokine networks found subsets of the same correlations between the same cytokines as we have found. 24 All together, these results suggest that these four cytokines may play an important role in the response to gastro-intestinal injury and bear specific investigation in future studies.
Our work has some limitations, including a small sample size that precluded the performance of more correlations; the use of clinical measurements and plasma biomarkers only; and the use of obese patients only. Nonetheless, it should be noted that the range of BMI in our cohort was very large, from 35 to 64 kg/m 2 . The fact that we still saw consensual cytokine responses from the slightly obese to the severely morbid obese lends reassurance that this response is not affected by obesity.

ACK N OWLED G EM ENTS
We would like to thank Prof. Micha Mandel of the Hebrew University of Jerusalem Department of Statistics for his support in the data analysis plan, data analysis and interpretation.

CO N FLI C T O F I NTE R E S T
The authors declare that no conflict of interest exists related to this work.

AUTH O R CO NTR I B UTI O N
UT and FD contributed equally to this work. All authors participated in the design of the work and the preparation and review of the manuscript. FD collected most of the biological samples and clinical data, with some help from UT. UT performed the laboratory experiments, primary data analysis and first draft of the manuscript. UT, FD, PVvH and SS performed the data analysis and interpretation.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data sets generated and analysed during the current study are available from the corresponding author upon reasonable request.