Integrating patients in time series clinical transcriptomics data

Abstract Motivation Analysis of time series transcriptomics data from clinical trials is challenging. Such studies usually profile very few time points from several individuals with varying response patterns and dynamics. Current methods for these datasets are mainly based on linear, global orderings using visit times which do not account for the varying response rates and subgroups within a patient cohort. Results We developed a new method that utilizes multi-commodity flow algorithms for trajectory inference in large scale clinical studies. Recovered trajectories satisfy individual-based timing restrictions while integrating data from multiple patients. Testing the method on multiple drug datasets demonstrated an improved performance compared to prior approaches suggested for this task, while identifying novel disease subtypes that correspond to heterogeneous patient response patterns. Availability and implementation The source code and instructions to download the data have been deposited on GitHub at https://github.com/euxhenh/Truffle.


Introduction
Transcriptomics data has been collected and profiled in clinical and drug response studies for over a decade (Meyer et al. 2013).In most cases, researchers profile bulk expression, though more recently single-cell data was also profiled in such studies (Wang et al. 2020).The main goal of these studies is to reconstruct networks and systems that are activated in response to the disease, drug, or vaccine, over time (Almon et al. 2003, Huang et al. 2009).
A major challenge in the analysis of data from clinical trials is the fact that different individuals may display different response dynamics (Bar-Joseph et al. 2012, Ding et al. 2022).Even if the same biological process is activated, based on baseline differences (related to age, gender, prior disease history, etc.), these individuals may respond faster or slower to the same treatment.Furthermore, same-day visits do not correspond to the same disease state which makes it challenging to rely on the measured time points for integrating data across these patients.Another challenge is the heterogeneous responses from different individuals.While a single response trajectory is possible, often we observe a (small) number of endotypes."Endotypes" are subtypes of a disease characterized by different pathogenic mechanisms (L€ otvall et al. 2011, Czarnowicki et al. 2019, Battaglia et al. 2020) which can have an impact on the specific optimal treatment.Each of the endotype groups may respond differently to the same treatment and so the overall set of patients cannot be directly integrated when studying treatment or vaccine response.
Several methods have been developed to address the first challenge (aligning patients) (Listgarten et al. 2004, Lin et al. 2008).These often use expectation-maximization (EM) like methods.In these approaches, genes are represented as continuous curves and individuals are assigned to different time points along these (Bar-Joseph et al. 2003).Such methods have been widely applied (Behnke et al. 2010, Czarnewski et al. 2019) but they still suffer from several drawbacks.First, the continuous expression assumption may be problematic when sampling rates are sparse (genes can change a lot between two consecutive measurements) and second, they cannot reconstruct trajectories for multiple subsets of patients but rather assume a homogeneous response among all patients.
Another direction that was explored, especially in the single-cell space, is that of trajectory inference.Unlike the EM methods, these approaches assume the presence of multiple states in the data and allow for multiple subsets or branching.These methods range from linear or tree-based, to more recent adaptations of RNA velocity (Saelens et al. 2019, Lange et al. 2022).However, most of these methods assume no relationships between cells or samples.Only a few methods have focused on the case when samples come from different time points as is often the case with clinical trials data.For example, Tempora (Tran and Bader 2020) assigns temporal scores to each cluster of cells which are used to determine the direction of the edges.Psupertime fits a series of ordinal logistic regression models that separate time points while trying to find a small number of genes that influence the resulting order (Macnair et al. 2022).However, these single-cell methods assume a very large number of samples (in the thousands or tens of thousands) which is not available for most clinical studies including the ones analyzed in this paper.In addition, they usually do not explicitly map the different subgroups within the data, leaving it for subsequent, post-processing, analysis.
In this work, we present Trajectory Inference via Multicommodity Flow with Node Constraints (Truffle), a method that performs pseudotime ordering of samples in short time series data.Truffle is based on the multi-commodity flow algorithm (Leighton et al. 1995) which generalizes minimum cost flow problems to include multiple source and sink nodes.Each sample in our data can be seen as either a source or a sink node and we are interested in recovering directed paths between these that minimize a cost function (typically some distance in gene space).The advantage of Truffle is that these trajectories can be constrained to satisfy timing restrictions and to pass through other nodes which correspond to intermediate disease states not present in the patient specific time series.Endotypes are then determined by constructing a state diagram for different subsets of patients.Truffle allows for the possibility of recovering contrasting endotypes since trajectories are inferred for each patient rather than for the entire dataset.
We tested Truffle on several microarray and bulk RNA-seq datasets.As we show, Truffle can accurately identify relevant disease trajectories and pathways, improving upon prior methods for clinical time series data and methods for singlecell data.A number of novel trajectories identified by Truffle suggest new subsets of patients that can benefit from precision medicine.
Raw gene counts were downloaded from NCBI GEO for the two RNA-seq datasets (psoriasis and COVID-19).Only protein-coding genes that had >0:25 counts per million (CPM) in at least 1% of the samples were kept.In the case of duplicated gene identifiers, the gene with the highest mean expression was considered.Datasets were then normalized for their guanine-cytosine (GC) content and trimmed mean of M-values (TMM) was performed (Robinson and Oshlack 2010).If batch information was present, ComBat was used to extract batch-corrected expression values (Johnson et al. 2006).Only samples with disease/treatment were used for pseudo-ordering.For microarray data, in the case of multiple probesets belonging to a protein-coding gene, only the one with the highest expression was kept.The Crohn's dataset was pre-normalized by Robust Multichip Analysis (RMA).
We removed symptomatic COVID À 19 − from the COVID-19 data and kept only the patients who tested positive for the disease.

Assignment of disease states through clustering
To obtain disease states, we clustered the samples.We followed a standard practice that is also adopted by other computational tools such as Seurat (Hao et al. 2024).We first ran principal component analysis (PCA) to obtain low dimensional embedding vectors which were then used to construct a fuzzy simplicial set as done by Uniform Manifold Approximation and Projection (UMAP) (McInnes et al. 2018).We adjusted the number of neighbors based on the total number of samples-using 15 for Crohn's, 20 for COVID-19, and 5 for psoriasis.Larger numbers resulted in highly connected graphs.This connectivity graph is the input for both Leiden clustering (Traag et al. 2019), and multicommodity flow (below).
To assign states to biological processes, we performed gene set enrichment analysis (GSEA) (Subramanian et al. 2005) using the prerank function of GSEApy (Fang et al. 2023).Genes were ranked based on the following score: gene score i ¼ − log 10 ðadj: pÀ valueÞ � log 2 ðFCÞ where in the first term, adjusted p values were obtained from a two-sided Kolmogorov-Smirnov test (Massey 1951) comparing the diseased and healthy sets of patients, and the second term is the log fold-change in gene expression between the two sets.We rely on the gene ontology (GO) biological processes marker set for the enrichment analysis in this work (Ashburner et al. 2000).

Multi-commodity flow with node capacity constraints
The multi-commodity flow problem with node capacity constraints is defined as follows.Consider a directed graph G ¼ ðV; EÞ, where an edge ðu; vÞ 2 E has an associated cost c u;v .We are given a set of K commodities K :¼ ½K�.The i th commodity is defined by a source and sink node ðs i ; t i Þ.
Multi-commodity flow can be used to model patient trajectories.Assume for simplicity patients with only two visits each.In this setup, each patient corresponds to one commodity, and the two visits represent its source s and sink t.The objective is to recover a smooth disease trajectory between these two endpoints.If the data contains patients with diverse disease states, we can assume that some of the samples will lie "in between" s and t.The shortest path between these two nodes in the neighbors graph captures this smooth transition.By setting edge and node capacities we force the algorithm to look for robust paths (defined here as paths with similar state transitions even though they share no edges).Finally, if a patient has more than two time points, we consider each transition separately.For example, a time series a !b !c is split into two commodities a !b and b !c.Specifically to use multi-commodity for trajectory inference, we use the following constrains.For every commodity i, we wish to learn separate functions f i : E ! f0; 1g that satisfy the following constraints: 1) Max edge capacity: the total amount of commodity that passes over an edge does not exceed its capacity 8ðu; vÞ 2 E : 2) Flow conservation: flow must fully exit source nodes and enter sink nodes.For all i 2 K: Given a node capacity N>0, we also consider the following constraint: 3) Max node capacity: the total amount of commodity that passes through a node does not exceed its capacity Along with flow conservation, constraint three guarantees limits on both incoming and outgoing flow.This variant of multi-commodity flow with node capacity constraints has also been explored before (Charikar et al. 2019).The integer problem has been shown to be NP-complete (Even et al. 1975), however, its fractional form (setting the codomain of f to be [0, 1]) can be solved in polynomial time through linear programming.We use the open source Python optimization library pyomo (Bynum et al. 2021) and the glpk solver (Oki 2012).It is worth noting that faster commercial solvers exist (Meindl and Templ 2013) (Fig. 1).
In the general formulation of the problem, each commodity can have a demand D, and each edge can have a capacity C (Leighton et al. 1995).Since a priori we do not have any preference for individuals, we set D ¼ 1 for all commodities.We set C ¼ 1 for psoriasis and Crohn's datasets.For the COVID-19 data, the problem was infeasible for C ¼ 1, so we used C ¼ 2. Enforcing edge and node capacities prevents outliers and errors in the data from having a large impact.An example has been provided in Supplementary Fig. S1.

Obtaining flow satisfying solutions
We learn f by optimizing the following target function As we are concerned with smooth trajectories, this is initialized as the Euclidean distance between the PCA embeddings for nodes u and v.
Note that for any given commodity defined by source s i and target t i , most of the edges "far away" from s i and t i will not be picked by the solver.We can incorporate this observation into our problem by considering only edges that belong to any path s i !t i of length ≤' for some '.This reduces the runtime for large datasets without compromising the optimality of the solution.For the smaller datasets, we found that the solution to this modified problem was similar to the original one.For the COVID-19 data, we set ' ¼ 4. Unreachable commodities were removed (17%).

Trajectory inference from optimal flow paths
After obtaining a path for each patient, we aggregate this information in the form of a state-transition matrix.In this work, we estimate initial and final state probabilities from the data, although domain expertise or priors determined from larger knowledge bases can be also used.Finally, we can then compute the most likely trajectories by performing random walks of a desired length.This is preferred over simply counting the occurrence of each path since in that case we could miss trajectories which are not identical, but show the same trend.For example, the paths 0−5−2−7 and 0−5−3−2−7 are different, but likely correspond to a similar disease trajectory.Our setup would assign a high probability to transitions 0−5 and 2−7.

STEM analysis of learned trajectories
To determine groups of genes that follow similar transcriptional programs, we perform Short Time-series Expression Miner (STEM) analysis (Ernst and Bar-Joseph 2006).We performed STEM normalization on gene expression values and used the default number of profiles (50), except for paths of length 2 where the maximum possible number is 16.Larger values for the number of profiles resulted in many redundant profiles that were nearly identical.For psupertime only, we reduced the "Minimum Absolute Expression Change" to 0, Integrating patients in time series clinical transcriptomics data since psupertime normalized expression values were in a much smaller range than for the other two methods.

Results
We developed a method to perform pseudotime ordering of multiple short times series clinical data based on optimal flow algorithms (Fig. 1).Our method takes as input gene expression data from multiple subjects along with their specific time point, and tries to reconstruct trajectories that describe distinct disease endotypes.As a proof of concept, we first performed a simulation study with randomly generated data.Truffle accurately recovered the simulated trajectories in this study (Supplementary Fig. S2).To further validate our method, we used clinical data for psoriasis, COVID-19, and Crohn's disease (Table 1).We compare our method against prior work developed for similar tasks including Tempora, psupertime, as well as a baseline that assigns endotypes based solely on clustering analysis.The set GO Biological Processes was used for Tempora.

Truffle recovers trajectories that indicate regeneration and reduction of inflammation in patients with psoriasis
We tested Truffle on bulk RNA data from psoriasis patients treated with secukinumab.The data spans 12 weeks and most patients have data for all four time points (Fig. 2a and  b).Leiden clustering identified six states (Fig. 2c).Cluster 0 predominantly consists of pre-treatment samples (50%) and contains no samples from week 12. Judging by the PASI scores (Fig. 2d), this cluster represents severe chronic plaque psoriasis.GO analysis shows significant upregulation of genes involved in the regulation of immune response (FDR ≤0:001) and defense response to virus & bacterium (FDR � 0, Fig. 2f) when compared to healthy samples.We also see significant upregulation for keratinocyte differentiation (FDR ≤0:001) which is a hallmark of moderate-severe disease states (Ma et al. 2023).Other immune-related processes such as Neutrophil Chemotaxis, Antimicrobial Humoral Response, and Regulation Of Interferon-Beta Production were also up-regulated in this cluster (Supplementary Fig. S3d).In contrast, for cluster 1, approximately 70% of the samples are from week 12 and there are no samples assigned to this cluster from the pre-treatment week.The PASI scores for cluster 1 were also the lowest among all clusters (an average of 2.3).This cluster is enriched for intermediate filament and supramolecular structure organization, and keratinocyte differentiation is no longer significant.Downregulation of processes related to regulation of gene expression is also seen as a result of drug action, along with a reduced immune response.
We first looked at the most common cluster transitions using patients' samples timeline without cost constraints.We found that three patients transitioned from state 0 ! 1, and two remained at state 4. All the remaining transitions were exclusive to only one patient.Next, we ran Truffle to uncover smoother response trajectories.Figure 3 shows the state diagram identified by Truffle as well as the top three paths.The transition 0 ! 1 was supplemented with two intermediate states, 5 and 3. GO analysis (Fig. 2f) shows that state 5 is characterized by a downregulation of defense response mechanisms when compared to state 0, while serving as an intermediary for a number of downregulated terms in state 1.On the other hand, state 3 is characterized by an upregulation of extracellular matrix organization which plays a role in tissue regeneration.Among the baselines, Tempora was able to recover paths of length 1 only (Fig. 4a).However, it correctly identified state 1 as a terminal state, but also 3 and 5. Psupertime identified 294 genes which vary coherently with time.GO analysis shows that these genes are enriched for intermediate filament and supramolecular fiber organization, as well as epidermis development.However, no significant terms involving defense response were found for the psupertime results.
Finally, we performed STEM analysis on the top three trajectories identified by Truffle.Profiles involving upregulation of epidermis development and downregulation of defense response overlapped across all three trajectories.Trajectories 0−5−1 and 4−5−1 contained decreasing profiles which were and Supplementary Fig. S3c].These two trajectories differ in their initial state only.While states 0 and 4 are both enriched for defense response, state 4 shows a downregulation of terms such as cytoplasmic translation and other biosynthetic processes.

Truffle identifies different immune responses to COVID-19
We repeated the analysis with samples from a larger dataset of COVID-19 patients collected at days 0, 3, and 7. Clustering analysis identified 10 states (Fig. 6c).State 8 consisted of day 0 samples, and showed the highest acuity scores (Fig. 6d and e).State 0 showed significant upregulation of inflammatory response and other defense mechanisms when compared to healthy samples (FDR � 0, Supplementary Fig. S5).State 1 was similarly enriched for "Defense Response to Virus," but not for inflammation.About 20% of all patients ended in state 2, which differed from healthy samples only in it being significantly enriched for Antimicrobial Humoral Response and Defense Response To Bacterium (FDR � 0).This suggests that this is a milder state than the previous two, also confirmed by acuity scores where cluster 2 is the only one containing no samples with acuity 4 or 5 (Fig. 6e).Across all three time points, most patients (10) moved from state 0 to state 2. This was also the top trajectory captured by Truffle (factoring in initial and terminal probabilities for each state, Fig. 6f).In contrast, this trajectory was not recovered by Tempora (Fig. 6g).
On the other end, for T 3 , STEM assigned >9000 genes to a strictly increasing profile (profile 41, Supplementary Fig. S4a).This profile was also enriched for processes related  to sensory perception of smell, but this time we see an upregulation of related genes across all four temporal steps.Profile 2 (T 3 ) and profile 9 (T 4 ) indicate downregulation of immune response.Profile 9 is gradual.Looking at GO enrichment of the final state of T 4 (cluster 3), we observe a return to baseline (healthy) for various defense response processes and downregulation of gene regulation activities (Supplementary Fig. S5).Tempora, on the other hand, identified only two paths of length ≥2.These were Q 1 :¼ 6−2−3 and  Integrating patients in time series clinical transcriptomics data Q 2 :¼ 8−7−9−2−3.Three significant STEM profiles were determined for Q 1 , none of which was significantly enriched for any GO process (FDR ¼ 0:05).For Q 2 , STEM returned 11 significant profiles.Among these, only three were enriched for GO processes (Supplementary Fig. S4b).Profile 10 was enriched for sensory perception of smell, and profile 7 was enriched for the only term "Positive Regulation of NF-kappaB Transcription Factor Activity."Meanwhile, profile 37 showed an initial increment, followed by a monotone decrement of processes related to signaling.Finally, psupertime identified 462 relevant genes.GO analysis using these genes returned only one process: "Hydrogen Peroxide Catabolic Process" (FDR ¼ 0:007).

Truffle identifies two contrasting response mechanisms to ustekinumab in patients with
Crohn's disease Finally, we tested Truffle on microarray data from patients with Crohn's disease treated with ustekinumab (VanDussen et al. 2018).The data was collected at weeks 0, 8, and 44.Clustering analysis revealed eight distinct states.States 1 and 4 were not statistically different from healthy samples.States 0, 3, and 6 expressed genes enriched for inflammatory response, while cluster 2 showed a downregulation of the process (Supplementary Fig. S6b and c).
The top Truffle trajectories of length 2 were C 1 :¼ 3−4−1 and C 2 :¼ 2−5−0.C 1 transitions from a state with inflammation into two healthy states, suggesting that patients along this path saw improvement from the drug.In contrast, for C 2 we see an activation of immune response in its final state (cluster 0).Indeed, about 14 patients were clustered under state 0 at week 44, suggesting that they showed partial response to the drug.STEM analysis of C 1 returned several decreasing profiles which were enriched for inflammatory response.In contrast, C 2 was assigned increasing profiles enriched for immune response and activation of T cells (Supplementary Fig. S6d).Thus, Truffle was able to recover two contrasting endotypes for patients in this study.

Discussion
Several trajectory inference methods have been developed to date and these differ in representation power and assumptions made (Saelens et al. 2019).Most of the work has focused on single cell with much less focus on data collected in clinical studies.Here we focus on studies that profile a small number of time-points in multiple patients.To analyze such data, we developed Truffle which respects the time ordering of samples for a given patient, and obtains patient journeys through the disease/treatment process.Truffle is based on multi-commodity flow by splitting short time series into source and target nodes.These are then connected through a path that travels through other intermediate nodes in order to generate a smooth path.We tested Truffle on several time series datasets and compared it to two other methods developed for similar tasks.
For the psoriasis dataset, all patients display a significant health improvement after treatment with secukinumab as indicated by their PASI scores and GO analysis of the terminal state.Since patients respond differently to the treatment, we sought to understand different endotypes within the patient population.Clustering analysis does not lead to accurate grouping of disease subtypes.Some of the other methods were able to capture the improvement either by identifying a healthy final state (Truffle, Tempora) or by showing enrichment for healing biological processes (psupertime).However, Tempora identified only paths of length 1, thus providing lower resolution into the drug response progression, while psupertime does not provide details into different response mechanisms or endotypes due to its linearity assumption.Only Truffle was able to capture temporal dynamics of the treatment process among different patients and obtain different endotypes.For example, Truffle recovered two paths which end in a healthy state but travel through different Both show the downregulation of IL-27 and its pathway genes.Reduction in expression of type I & II interferons (IFNs) and/or tumor necrosis family (TNF) receptors, which are regulators of IL-27, has been previously observed as part of the recovery (Povroznik and Robinson 2020).Furthermore, IL-27 was previously reported to promote the onset of psoriasis (Shibata et al. 2010).However, they also differ in other pathways.One of these trajectories was characterized by an upregulation of extracellular matrix organization (ECM) and downregulation of intermediate filament organization (IFO), while for the other trajectory we observed the opposite.Prior work has shown that activation of ECM is related to the severity of psoriasis (Wagner et al. 2021).We hypothesize that the upregulation of ECM may be an intermediary stage of slow responders.Results show that a subset of patients quickly attained normalization of keratinocyte differentiation (Figs 2 and 3 clusters 1, 3, 5).Such patients can be deemed as super/fast responders to therapy.These patients can be further investigated to better tailor personalized therapy.
For the COVID-19 dataset, prior methods failed to recover smooth trajectories with any significant GO terms.Tempora recovered trajectories that oscillate between time points, which makes them hard to interpret, and psupertime returned only one significant GO process, likely because this linear method was forced to combine heterogeneous subtypes in its trajectories.Truffle identified several trajectories, including ones which showed a downregulation of defense response over time and others where this response was reinstated at day 7.This was confirmed by a reduction of sensory perception of smell during this time step.
While the applications we presented are mainly focused on immunology, we believe that Truffle can also be applied to oncology time series data and that it can also be integrated with time series data from other sources including electronic health records or claims databases.
While successful, Truffle has a few limitations.The datasets we used in this study contained at most 650 samples.The open-source linear solver we used to optimize a graph of this size may not scale to graphs with several thousands of samples.In this case, several simplifications to the problem may need to be introduced, such as limiting the set of edges a commodity can be transported over.For the specific datasets we evaluated, Truffle took 0.12 s to run for the small psoriasis dataset and 22 s for the larger COVID-19 dataset (' ¼ 4) (tests performed on a MacBook Pro with an M3 Pro Max chip).In addition, faster commercial solvers can also be used.
To conclude, Truffle is a method for integrating patient data in time series transcriptomics studies.It is able to both, identify patient trajectories and subgroups within a population.Truffle is available as an open-source software from the link in the abstract.

Figure 1 .
Figure 1.Schematic illustration of Truffle.For each patient, our flow algorithm returns a trajectory that passes through intermediate nodes for a smoother response.These trajectories are then aligned with the clustering results to obtain a state diagram.Finally, by estimating state initial and final probabilities from the data, we can compute and study the top directed trajectories.

Figure 2 .
Figure 2. Clustering analysis of the psoriasis dataset.(a) and (b) Distribution of visits across patients.(c) UMAP plot of cluster assignments.(d) Boxplots of PASI scores for each cluster.(e) Relative frequency of visits by cluster.(f) Top GO terms for each cluster against healthy samples.We used a KS test to rank the genes.A (�) symbol means the category was statistically significant [ð��Þ � q � 0 and ð�Þ � q ≤ 0:05].

Figure 3 .
Fig.S5).State 1 was similarly enriched for "Defense Response to Virus," but not for inflammation.About 20% of all patients ended in state 2, which differed from healthy samples only in it being significantly enriched for Antimicrobial Humoral Response and Defense Response To Bacterium (FDR � 0).This suggests that this is a milder state than the previous two, also confirmed by acuity scores where cluster 2 is the only one containing no samples with acuity 4 or 5 (Fig.6e).Across all three time points, most patients (10) moved from state 0 to state 2. This was also the top trajectory captured by Truffle (factoring in initial and terminal probabilities for each state, Fig.6f).In contrast, this trajectory was not recovered by Tempora (Fig.6g).Next, we studied the top trajectories identified by Truffle at varying levels of resolution.The top trajectories of lengths 3 and 4 were T 1 :¼ 0−1−5−2, T 2 :¼ 0−1−5−4, and T 3 :¼ 0−1−2−5−4, T 4 :¼ 0−2−5−4−3, respectively.For

Figure 4 .
Figure 4. Trajectories uncovered by Tempora and psupertime for the psoriasis dataset.(a) Transition graph identified by Tempora.Five trajectories of length 1 were identified.(b) Separation of time points by psupertime.The y axis is the density of each time point and the x axis is the temporal ordering.(c)The top five genes identified as relevant by psupertime.These correspond to the genes with the largest absolute coefficients.(d) The top GO terms for all the relevant genes (294).Subfigures (b) and (c) were generated using psupertime.

Figure
Figure Selected STEM profiles for the top three Truffle trajectories in the psoriasis dataset.(a-c) Two selected profiles for each of the three trajectories.In (c, right) "IL-27 Mediated Signaling Pathway" obtained a very high combined score (1e6), hence, was removed from the plot for clarity.The full list of profiles can be found in the supplement.

Figure 6 .
Figure 6.Clustering and trajectory analysis for the COVID-19 dataset.(a) and (b) Distribution of visits and distribution of visit counts per patient.(c) Clustering 650 samples from 304 patients.(d) Relative frequency of visits and (e) acuity scores per cluster.(f) A pruned diagram of top state transitions identified by Truffle.Pruning was performed by taking the fewest top edges that amount to ≥50% of a node's outgoing weight.(g) The tree learned by Tempora.Final states are 3, 1, and 5. (h) The top genes that vary with time according to psupertime (plot obtained from psupertime).

Table 1 .
Clinical data used in this study.a a All three datasets contain missing values.We show both the number of patients who tested positive (þ) and the number of healthy control patients (−).b Pretreatment week.i152Hasanajet al.