Characterization of Microcirculation in Multiple Sclerosis Lesions by Dynamic Texture Parameter Analysis (DTPA)

Objective Texture analysis is an alternative method to quantitatively assess MR-images. In this study, we introduce dynamic texture parameter analysis (DTPA), a novel technique to investigate the temporal evolution of texture parameters using dynamic susceptibility contrast enhanced (DSCE) imaging. Here, we aim to introduce the method and its application on enhancing lesions (EL), non-enhancing lesions (NEL) and normal appearing white matter (NAWM) in multiple sclerosis (MS). Methods We investigated 18 patients with MS and clinical isolated syndrome (CIS), according to the 2010 McDonald's criteria using DSCE imaging at different field strengths (1.5 and 3 Tesla). Tissues of interest (TOIs) were defined within 27 EL, 29 NEL and 37 NAWM areas after normalization and eight histogram-based texture parameter maps (TPMs) were computed. TPMs quantify the heterogeneity of the TOI. For every TOI, the average, variance, skewness, kurtosis and variance-of-the-variance statistical parameters were calculated. These TOI parameters were further analyzed using one-way ANOVA followed by multiple Wilcoxon sum rank testing corrected for multiple comparisons. Results Tissue- and time-dependent differences were observed in the dynamics of computed texture parameters. Sixteen parameters discriminated between EL, NEL and NAWM (pAVG = 0.0005). Significant differences in the DTPA texture maps were found during inflow (52 parameters), outflow (40 parameters) and reperfusion (62 parameters). The strongest discriminators among the TPMs were observed in the variance-related parameters, while skewness and kurtosis TPMs were in general less sensitive to detect differences between the tissues. Conclusion DTPA of DSCE image time series revealed characteristic time responses for ELs, NELs and NAWM. This may be further used for a refined quantitative grading of MS lesions during their evolution from acute to chronic state. DTPA discriminates lesions beyond features of enhancement or T2-hypersignal, on a numeric scale allowing for a more subtle grading of MS-lesions.


Introduction
The search for novel imaging biomarkers in multiple sclerosis (MS) has modified the concepts of neuroimaging from identifying imaging sequelae of demyelination on conventional MR-images towards strategies capable of examining functionality and pathophysiology of the MS brain. A recent technique that contributed to a better understanding of vascular changes in MS, is dynamic susceptibility contrast enhanced (DSCE-) MRI. DSCE-MRII identifies various patterns of impaired perfusion, either in nonenhancing MS lesions or NAWM, compared to healthy controls. On the other hand, inflammatory activity is accompanied by increased perfusion in lesions during the acute phase of the disease, compared to NAWM [1][2][3].
Local perfusion changes in MS lesions are currently interpreted as a consequence of local inflammation-mediated vasodilatation -a phenomenon secondary to hyperemia and blood congestion within the brain parenchyma [4,5]. However, different observations in perfusion imaging studies in MS have challenged the interpretation of abnormal perfusion as a reactive phenomenon to inflammation. The occurrence of demyelinating lesions is not inevitably coupled to the presence of a local preceding inflammatory reaction [6,7], and diffuse NAWM changes in the absence of structural lesions may be the consequence of a down-regulation in cerebral micro-circulation due to astrocyte dysfunction, or secondary axonal damage in the NAWM [8]. Recently, some authors suggested a formation of new outflow routes along plaque formation bypassing obstructed pathways as the reason for local perfusion changes [9].
The contribution of vascular changes to the generation of MS lesions in MS is still a matter of debate. Currently, it remains to be elucidated whether alteration in CNS perfusion in MS is a cause or rather a consequence of disease pathogenesis. From a methodological point of view, current concepts of perfusion imaging have mainly addressed vascular changes that are described either by the time needed for a contrast agent to pass the vasculature (MTT), the total volume of blood within the cerebral vasculature during the passage of a contrast bolus (CBV), or the amount of blood that perfuses the brain per time unit (CBF) [10]. Calculation of CBV, CBF, and MTT from concentration-time curves is based on the indicator dilution methods for non diffusible tracers [11]. While this technique has been clinically applied for more than a decade to investigate perfusion deficits in acute stroke [12] and altered microcirculation in brain tumors [13,14], its application in MS is still new and MS-related perfusion characteristics has not yet been evaluated. A principal limitation of model-based perfusion imaging in MS is the variable blood-brain barrier disruption in active MS lesions requiring complex pharmacokinetic modeling to correct for extravasation of contrast medium from the plasma to the extracellular space.
A potential alternative to leakage modeling is the analysis of the tissue of interest by means of texture analysis of the MR-images at different time points along with the bolus passage of Gadolinium. Texture parameters [15][16][17] enable quantitative analyses of MRimages. Although no strict definition of image textures exists, they are described as complex visual patterns which are composed of spatially organized, repeated fingerprints characterized by individual brightness, size, shape, etc. [18]. In recent years, software developments facilitated the study of texture parameters of medical images, e.g. the MaZda package [19,20]. Texture analysis has been applied for structural imaging in MS [21,22] in breast cancer, liver cirrhosis, brain tumors, epilepsy or acute ischemic stroke [23][24][25]. Kassner et al. demonstrated that statistical or spectral textural features outperformed visual assessment in discriminating between tumors, as well as in discerning subtle anatomic changes associated with a high risk of seizures in patients with epilepsy [23][24][25].
Recently a novel approach has been proposed, namely Dynamic Histogram Analysis (DHA) where a simple form of texture analysis is applied to a time series of DSCE-MRII-images [14]. DHA investigates the time dependency of mean and standard deviation parameters of the voxel intensity values of the tissue of interest. This DHA approach has been successfully applied to differentiate cerebral gliomas in accordance with their histological classification [14]. In this study, a further extension to DHA is proposed, the DTPA. DTPA investigates MR differenceimages, calculated by subtraction of the first steady state baseline image from images during bolus passage and reperfusion. In this paper, we extend the focus of texture analysis to quantitative studies of changes in micro-structural perfusion and leakage. This explorative study aims for investigating the differences in dynamical texture parameters between NELs, ELs, and NAWM during contrast agent passage. We hypothesized that: (a) texture parameters show characteristic tissue-and time-dependent (dynamic) behavior of MS plaques during bolus passage that differ from NAWM; and (b) statistically significant differences in texture parameters can be detected between NAWM, ELs and NELs, sufficient to discriminate the three tissue types from each other.

Patients
Patients with clinical isolated syndrom (CIS) or MS diagnosed according to the 2010 McDonald's criteria, with or without immunomodulatory therapy, were included if at least one EL was present on MRI [26]. In equivocal cases, the inflammatory/ demyelinating origin of the disease had to be demonstrated via biopsy and/or MR spectroscopy. Besides the clinical evaluation and MR imaging, lumbar puncture and blood tests were performed in all patients to exclude alternative diagnoses. Main exclusion criteria were any other diseases that could better explain the patient's symptoms and signs. Patients gave written informed consent prior to the study entry and the study was approved by the local cantonal ethics committee Bern, Switzerland. (c) at early-reperfusion, and (d) at late-reperfusion. The MS lesion shows a non-uniform structure that varies during bolus passage; (e) Concentration curves were computed within 3 TOIs; the first within the brightest circular part of the MS-lesion, the second TOI in the edema surrounding the enhancing part of the lesion; the third TOI in the NAWM on the contralateral site. The concentration curve of the EL reaches its maximum slightly earlier in time than the NAWM due to BBB disruption. Interestingly, the contrast time curve of edema reaches less high concentrations, indicating hypoperfusion of the edema. doi:10.1371/journal.pone.0067610.g001 2.2. MRI parameters DSCE-imaging was performed on our institutes' Siemens MRscanners using standard echo planar imaging (EPI) sequences at 1.5 T and 3.0 T. For our 1.5T scanner a TE of 47 ms and a TR of 1440 ms was used. For our 3T scanner we had a TE of 32 ms and a TR of 1500 msec. At both field strengths a FOV of 230 mm6230 mm was used with a slice thickness of 5 mm. A total of 19 parallel images with a time resolution of 1500 ms were acquired on our 3.0T system, and a total of 12 parallel images with a time resolution of 1440 ms were acquired on our 1.5T scanner. At 3.0T a total of 80 different time points were acquired, whereas on 1.5T a total of 40 different time points were sampled after bolus application. All patients received Gadobutrol (Gadovist TM ) 0.1 ml kg -1 bodyweight. Patients were positioned comfortably in the head coil and padding on either side of the head was used to help immobilization. Further the intravenous line with a long tube was put before acquisition to avoid unnecessary MRI table moving during examination.

Image Processing
The image analysis method we applied focuses on DSCEdifference images: These difference images are computed from the difference between the first steady state EPI image and images recorded during bolus arrival, passage and reperfusion. In mathematical terms, let I m ½x,y,t n (1ƒmƒN; 1ƒnƒN; 0ƒxƒ256; 0ƒyƒ256)be the n-th image of the image time series at slice location m. The difference images D m ½x,y,t n are computed as follows: D m ½x,y,t n ~I m ½x,y,t 0 {I m ½x,y,t n . In case of a non-leaky blood brain barrier (BBB), the contrast of D m ½x,y,t n is such that the brighter the voxel of D m ½x,y,t n , the more contrast agent there is in the voxel at that time point t n . From D m ½x,y,t n all time dependent texture parameter maps (TPMs) are computed, resulting in TPM m ½x,y,t n . The rationale for the computed texture maps in this study is given in section ''Computed TPM types during BLP, IFP, OFP and RFP'' below. To illustrate the time-dependent perfusion behaviour in a EL, figure 1 a-d display typical DSCE EPI images during the four time-intervals that are investigated in this study.
2.3.1. Analysis of texture parameter maps. Figure 2 shows schematically a 11|9 voxel region detail of an MR-image I½x,y,t n at a time point t~t n during bolus passage. Each voxel ½x,y is assigned with a specific voxel intensity distribution; e.g. a voxel ½x,y in the center of NAWM will have a different voxel intensity distribution compared to a voxel ½x,y centered in the gray matter (see Figure 1b). Gray matter, due to its strong vascularisation, will in general be more hypointense during in-, out-and reperfusionphases compared to voxels within NAWM. Furthermore, NAWM constitutes a smaller range of voxel intensity distribution than voxels closer to the interface between two tissue types. This results in a highly variable voxel intensity distribution. The statistical properties of regions with variable center points ½x,y strongly differ within an image as a function of tissue type, and of the position within a tissue (e.g. representing a voxel at the periphery or in the center of the tissue). Extending our previous work [14], we now incorporated further histogram based texture parameters in the analysis of voxel intensity distributions computed from the DSCE-EPI images. Texture parameters computed as a function of the location ½x,y can be used to characterize these (non-) Gaussian voxel intensity distributions. In order to quantify these textures as a function of position ½x,y and time, we developed an application that is called dynamic texture parameter analyzer (DTPA). This application computes a series of different texture parameters for user-defined square voxel neighbour kernel size; in this paper we used a 5|5 voxel range (Figure 2), for every region around ½x,y as a function of time. In this manner, a time series of texture parameter images or maps (TPMs) is obtained. For illustration, an example of the TPMs computed by the application at one time point during the reperfusion phase is given in Figure 3. In this paper we focused on the first four moments of the mean and their variances during baseline, contrast agent inflow, outflow, and reperfusion phase to characterize MS-lesions of different types (NEL, EL), and NAWM.

Data Normalization.
Our study sample has been recorded on three MR scanners runnung at with 1.5 and 3T of the same vendor using the same DSCE-Perfusion sequence. Differences in field strength, patient head size, and head positioning lead to significant differences in RF-coil-load and thus to changes in MRI-voxel intensities I m ½x,y,t n from patient to patient. Additionally, in DSCE-MRI the voxel intensities of the difference image series D m ½x,y,t n do not only depend on the flow properties of the blood vessels, but also on the cardiac ejection fraction. In order to be able to compare the texture parameters between two different patients, normalization for both sources of variance is mandatory. Since both I m ½x,y,t n and D m ½x,y,t n derived texture parameters are studied, a twofold normalization is required: 2.3.3. Corrections for Differences in RF-coil-loads. Normalization for these differences was performed by setting a user defined region encircling a specific tissue of interest (TOI) to a certain value. The normalization constant g 1 is defined as the quotient of average voxel intensity value of voxels within this user-defined region m TOI divided by a user-defined constant t TOI , that resembles the ideal value for this tissue type, so g 1~mTOI =t TOI , and the normalized image time series I m,norm ½x,y,t n ~I m ½x,y,t n =g 1 for all t n . With the selection of the most homogeneous tissue located at a predefined anatomic structure, the inter-subject variance can be minimized further. In this study, normalization was performed by defining a region in NAWM of the first baseline image and with a t TOI~1 000.
2.3.4. Correction for different cardiac ejection fractions. Tissue concentrations C m ½x,y,t n time curves (see Figure 2b) as a function of time can be computed from the formula C m ½x,y,t n c onst:| ln (Im½x,y,t n =Im½x,y,t 0 ). Depending on the cardiac ejection fraction, the duration of inflow and outflow phase may differ substantially, resulting in different peak contrast concentration values that prevent a time-point-by-time-point comparison between patients for D m ½x,y,t n and all its derived time dependent texture parameters. To overcome this problem, we compared identical time intervals of normalized D m ½x,y,t n as described in [14]. The following time intervals were defined (see also Figure 2b For the normalization of inter-patient differences due to differences in cardiac ejection fraction the area-under-the-timedifference-curve D m ½x,y,t n of a user defined TOI is determined for the IFP and OFP and is denoted by a AUC . Within the scaling factor g 2~aAUC =t AUC the dominator t AUC is set to 200, for all images D m,norm ½x,y,t n ~D m ½x,y,t n =g 2 . In the present implementation, g 1 and g 2 are determined for the same region in NAWM. The time points demarking the four time intervals (BLP, IFP, OFP, RPP) are determined automatically by the developed program by numerical analysis of the time dependent bolus passage function. Due to differences in timing caused e.g. by differences in cardiac function, only separately time averaged values over the four defined time intervals as defined above (BLP, IFP. OFP, RPP) can be adequately compared.

Computed TPM types during BLP, IFP, OFP and
RPP. This section shortly describes the computed histogram based on TPMs. Figure 3 illustrates all derived texture parameter maps computed in this study for a single patient with RRMS.
1. Original EPI MR image time series: TPM-ORIG The original Im½x,y,t n EPI image series can be regarded as a texture map itself, since every type of MR-pulse sequences is designed to obtain a specific image contrast, i.e. having specific textural characteristics. DSCE-raw images show limited but substantial contrast differences during the BLP, so that both EL and NEL can clearly be identified (Figure 1a). In case of an intact BBB, the passage of Gadubutrol causes T 2 Ã -shortening effects resulting in a signal decrease relative to the BLP. Counteracting this T 2 Ãshortening effect, voxel intensities I m ½x,y,t n increase during the RPP and OFP in ELs with a disrupted BBB. The leakage effect during the OFP and RPP is illustrated in Figure 3, (denoted as ''ORIG''). The EL (red) can be clearly depicted.
2. Difference image time series: TPM-DIFF The brightness of a voxel is proportional to its Gadobutrol concentration for the difference image series D m ½x,y,t n if the BBB is preserved. In case of BBB disruption T 1 -shortening effects in the tissue result in a signal decrease during the first passage of the contrast agent and reperfusion (see Figure 1e and Figure 3 DIFF). This TPM is further denoted by TPM-DIFF.
3. Local Average Difference Maps: TPM-AVER This TPM, denoted by Ave½x,y,t n , is the moving average filtered version of TPM-DIFF over a 5|5 voxel intensity range. In Figure 3, the TPM denoted with ''AVE'' gives an example of the appearance of this texture, illustrating better SNR, but worse spatial resolution. This TPM was added to provide the same spatial resolution as the other textures, and is denoted by TPM-AVER. 4. Local Variance: TPM-VARI The TPM-VARI Var½x,y,t n time series computed from D m ½x,y,t n is sensitive to edges within the image. The broader the local voxel intensity range, the brighter it appears in this TPM. To some extent this texture is related to the gradient image of D m ½x,y,t n . This texture thus reflects the local heterogeneity of a tissue [14]. In Figure 3, the texture map indicated with ''VARI'' is an example of this texture type, illustrating how an EL appears in this type of TPM. Note that the surrounding edema behaves differently. This TPM is further denoted by TPM-VARI.
5. Local standard deviation: TPM-STDEV The TPM StDev½x,y,t n is derived from Var m ½x,y,t n by taking its square root. This TPM is further denoted by TPM-STDEV. 6. Local Skewness: TPM-SKEW Skew½x,y,t n is the third moment about the mean of D m ½x,y,t n (see Figure 3, ''SKEW''). This TPM-SKEW is a measure of the asymmetry of the voxel intensity distribution, and in this case indicates the asymmetry of contrast agent distribution within the image voxels. An increase in asymmetry reflects an increase in skewness. For positive skewness, there is a preference in the distribution for brighter voxel intensity values, and for a negative skewness there is a preference for darker voxel intensity values.
7. Local Kurtosis: TPM-KURT Kurt½x,y,t n is the fourth moment about the mean of D m ½x,y,t n (see Figure 3, ''KURT''). This texture parameter measures the steepness or peakedness of voxel intensity distribution. The software computes the excess kurtosis, which sets the kurtosis for Gaussian distributed voxel intensities to zero. Voxel intensity distributions which are less peaked than the Gaussian distribution have an excess kurtosisv0 (platykurtic), and voxel distributions more peaked than the Gaussian have an excess kurtosisw0 (leptokurtic). This kurtosis TPM is denoted by TPM-KURT. 8. Local variance-of-variance: TPM-VAVA VarVar½x,y,t n of Var m ½x,y,t n is the variance-of-the-variance TPM, indicated as TPM-VAVA. Since Var m ½x,y,t n is a stochastic variable itself, it also has a variance (see Figure 3, ''VAVA''). Due to the large dynamic range of the TPM-VAVA and the limitation to 512 color discretization display-levels in the current implementation (see Figure 3), low intensity voxel values are all displayed in black, although they all have distinct floating point values. Interestingly, in contrast to the TPM-VARI this texture is not an ''edge detector'' as can be clearly seen in the image.
2.3.6. Determined Texture Parameters. Within the texture maps as displayed in Figure 3, manually shaped regions can be drawn. The basic statistical properties of these regions (NAWM, EL, and NEL) are determined within all TPMs. To clarify the interdependence of the TPMs, Figure 4 displays the TPMs (eight dark gray boxes), and the parameters determined from them. For each TOI (EL, NEL or NAWN) defined in the TPMs and for each time interval, the six statistical parameters (light gray box) are determined, resulting in a total of 8|3|16~144 parameters. The user defines the TOIs in one of the eight TPMs and the TOIs are copied fully automatically into all other TPMs.
2.3.7. Selection Criteria for the TOI. For this pilot study, all TOIs (EL, NEL and NAWM) were defined by a board certified  For every patient at least one image covering one or more lesions was analyzed, and a maximum of 3 EL and 3 NEL lesions was not exceeded. A total of 37 TOIs were defined in NAWM, 27 TOIs in ELs, and 29 TOIs in NELs. In most cases lesions and NAWM were analyzed on the same slice. Since in a few cases only one lesion and the corresponding NAWM were analyzed, the number of NAWM measurements exceeds the number of EL and NEL measurements. Four out of 27 ELs were surrounded by vasogenic edema (see Figure 1). Additionally to the DTPA we calculated the bolus passage (i.e. perfusion) characteristics of every TOI [27].

Workflow
Data analysis included the following workflow steps: (a) launch DTPA analyzer tool; (b) mount examination and load DICOM DSCE image-series (TPMs are computed and displayed); (c) draw contour in NAWM and press normalize button (data are twofold normalized now); (d) define the TOIs by drawing contours around EL, NEL and NAWM contours and press reporting button; all statistical parameters are written to text files that can be easily imported into a statistical analysis application; and (e) per TPM with three different tissue types defined, six parameters were determined for each tissue type, and for each of the four time intervals separately, making a total of 72 parameters, where 72 independent Wilcoxon-tests were performed in order to test for differences between the tissue types.

Statistical evaluation
The statistical software package R-Studio (0.96.331, available from http://www.rstudio.org/) was used for the statistical analysis of the data. In this explorative study we investigated (a.) which dynamic texture parameters discriminate between NAWM, NEL and EL, and (b.) during which time interval (BLP, IFP, OFP or RPP) these tissue types show most discriminative properties. The distribution of all texture parameters was examined in advance, in order to be able to select the correct test statistic. Since the texture parameters that we determined were not normally distributed, we employed a one way ANOVA followed by Wilcoxon sum rank tests. If the p-value of the ANOVAs F-test satisfied the condition p,0.05, the multiple Wilcoxon sum rank tests between NAWM, NEL and EL were performed. The p-values of these tests were corrected for multiple Wilcoxon sum rank testing by setting the ''BH'' option, which stands for the Benjamini, Hochberg, and Yekutieli method to control for false discovery rate. The false discovery rate of this method is a less stringent condition than the family-wise error rate, so these methods are more powerful than e.g. the Bonferroni test.

Clinical data
The results of every patient regarding gender, age, diagnosis, EDSS, disease duration, disease modifying treatment (DMT), acute disease exacerbation and symptomatic EL are presented in Table 1. 18 patients (83.3% female) with a mean age of 38.9 years (range 20-74 years) and a mean disease duration of 7.7 years were enrolled. 15 patients had relapsing-remitting MS (RRMS), two patients secondary progressive MS (SPMS) and one patient CIS. Three of these patients had tumefactive lesions and biopsy was performed in two of them showing predominately macrophage and T-cell mediated demyelination in one and antibody-mediated demyelination in the other patient (Table 1) corresponding to pattern I (patient 5) and II (patient 2) of the classification of Table 3. Overview of all texture parameters for which there exist statistical significant differences between all three TOIs. demyelination by Lassmann and Lucchinetti [28]. In the third patient, MR-spetroscopy and MR-angiography was sufficient to further underline the inflammatory/demyelinating origin of the lesion(s) (patient 1). 12 patients presented with acute disease exacerbation at the time of MRI, with symptomatic ELs in eight patients. None of the patients received steroid treatment before or at the time of MRI. Table 2 summarizes the number of statistical differences obtained for all TPMs as a function of the time interval and texture parameter type: For the different periods IFP, OFP and RPP we calculated the parameters, that discriminated between two tissue types (NEL vs. EL, NEL vs. NAWM or EL vs. NAWM). Measured during all timed periods, the most discriminative TPM is TPM-ORIG, followed by TPM-STDEV and TPM-VARI. Within the four preselected time periods that have been analysed, during the RPP the highest number of discriminative TPMs were detected. During RPP, for the TPM-STDEV 11 different TPs were detected that discriminated differences in the tissue under investigation. During OFP, TPM-STDEV (8 TPs) and TPM-VARI (8 TPs) were most discriminative, while during IFP, TPM-DIFF (9 TPs) yielded the most significant results.

Discriminators between two tissue types
Overall, the differences in tissue types between between EL, NEL and NAWM are predominantly expressed during RPP (62 out of 144 parameters were discriminative). The TPM types with least significant discrimination between the TOIs are the TPM-SKEW, with only (5 TPs), and TPM-KURT with 10 significant different TPs. The program did not determine the skewness and kurtosis during the base line period and are therefore indicated by ''n.a''. doi:10.1371/journal.pone.0067610.t004 Table 3 summarizes all parameters that discriminate NAWM, NEL and EL. This table lists p-values of the One-way ANOVA (p,0.05). In addition, we performed a repeated Wilcoxon sum rank tests corrected for multiple comparisons using the Benjamini-Hochberg method. TPM-DIFF maps, TPM-STDEV and TPM-  VARI were the strongest discriminators between the TOIs, and yielded the highest significance levels during the OFP, followed by the RPP. The TPM-ORIG was less discriminative between the TOIs.

The Time Evolution of Texture Parameters
As an example, Table 4 displays the derived time evolution of all texture parameters for TPM-STDEV derived statistical parameters. Two perfusion imaging parameters discriminated between the different tissue types, as investigated in this study, were mean TPM (p,0.002) and TPM-VAVA (p,0.0001). The strongest effects could be denoted for the mean TPM during the outflow for the VAVA TPM during the reperfusion period (p,0.002). Notably, only the mean TPM discriminated between NAWM and EL during inflow, while all six parameters discriminated the ELs during the outflow period (p,0.05). The best discrimination between EL and NAWM was achieved by the VAVA TPM (p,3.5610 211 ) during outflow and reperfusion. NEL were discriminated from NAWM by three parameters during outflow (mean TPM, p,0.002; TPM-SKEW, p,0.03 and TPM-VAVA (p,4610 24 ). The best discrimination between EL and NAWM and NEL and NAWM was achieved by the TPM-VAVA during reperfusion. For discrimination between NEL and EL, the mean TPM during reperfusion was the best discriminator (p,9.7610 26 ), followed by mean TPM during outflow (2.8610 25 ). Notably, the TPM-VAVA contributed less to the discrimination of NEL and EL during outflow (p,0.03) and reperfusion (p,2.5610 25 ). During the inflow period, the mean TPM was the best discriminator between NAWM and EL (4.4610 25 ), while TPM-VAVA most reliably discriminated between NAWM and NEL (2.5610 25 ). In summary, the mean TPM could be used as a discriminator between EL, NEL and NAWM during inflow and outflow, while TPM-VAVA discriminated during the outflow and reperfusion period. The TPM-VAR, the TPM-SKEW and TPM-KURT failed to discriminate between the different tissue types.

All TPMs Mean Parameters as a Function of Time
In Figure 5, 6 the mean values of all 8 TPMs for different time intervals are displayed. The indicated ranges represent the 95% confidence intervals. Differences between ELs, NELs and NAWM can be observed in e.g. TPM-DIFF and TPM-VARI especially during the OFP and RPP phase. Differences between NELs, ELs and NAWMs can be observed in the TPM-ORIG. The differences between the TOIs from TPM-KURT and TPM-SKEW are less pronounced.

Discussion
In this paper, we suggest a novel strategy towards in-vivo tissue characterization of ELs, NELs and NAWM in multiple sclerosis. The goal of this pilot study was to investigate the time dependent behavior of statistical texture parameters derived from time series of DSCE difference images. While previous studies aimed at either static texture analysis or description of cerebral blood flow, cerebral blood volume and blood-brain barrier permeability, we aimed to analyze specific texture parameter constellations that differentiate between ELs, NELs and NAWM by analyzing DSCE image time series data. We further analyzed whether these texture parameters evolve as a function of time intervals during bolus passage, and whether these differences are sufficiently large to enable a quantitative grading of the perfusion-related texture of MS lesions. Three points stand out from this study: (a) the voxel intensity distributions of tissue types in the DSCE (difference) image series vary as a function of time during bolus passage; (b) histogram based texture parameters reveal characteristic time responses for ELs, NELs and NAWM; and (c) different tissue types can be discriminated by statistically significant difference in histogram based texture parameter signatures.
The key finding of our study -that DTPA detects subtle differences in time dependent contrast agent distributions otherwise obscured by visual inspection of T1w CE images -implicates that quantitative differences in voxel intensity distributions of the examined DSCE difference images are related to alterations in micro-circulation and BBB properties of the investigated tissue types. The BBB constitutes a dynamic interface forming a neurovascular unit that controls the supply of nutrients while shielding neurons and glial cells off from potentially harmful substances. In MS, GD-DTPA diffuses into the nervous system along the leaky BBB, and conventional MRI provides only a gross estimate of tissue damage. Leukocyte infiltration into the perivascular space is mainly driven by BBB disruption and parenchymal inflammation, although some authors demonstrated passive diffusion of hydrophilic molecules and leukocyte recruitment in post-capillary segments [29]. Local blood flow changes precede the plaque formation process and elude conventional MRI [3]. Inversely, persistence of inflammatory activity along reconstituted BBB in tissue that lacks frank BBB disruption is obscured [30]. The poor clinico-radiological association of the T 2 lesion load in MS may thus be influenced by microscopic inflammation that contributes more strongly to disability [31]. Thus, quantitative knowledge about persistent inflammatory activity would enhance the knowledge about tissue damage beyond a deterministic and analysis of ELs, NELs and NAWM.
Texture analysis applied to MR-images enables such an extraction of quantitative information by post-hoc analysis. Different texture parameters have been suggested in the literature: (a) histogram-based parameters; (b) gradient-based parameters; (c) run length-based parameters; (d) co-occurrence matrix-based parameters; (e) auto regressive model-parameters; and (f) wavelet parameters [18]. In the present study, the analysis to histogram based texture parameters (i.e. average, variance, skewness, kurtosis and variance-of-the-variance) is restricted to explore whether differences between ELs, NELs and NAWM can be extracted directly from DSCE (difference) images. To the best knowledge of the authors, texture analysis has not been applied previously to time series of DSCE-MRI (difference) images aiming to characterize tissue response to a bolus of Gd-containing contrast agent as a function of time.
Previous studies have examined the feasibility of texture analysis to differentiate ELs (88%) and NELs (96%) with high sensitivity [32]. ''Coarse'' texture analysis in acute enhancing lesions (ELs) predicted tissue injury based on the severity of structural disorganization within acute lesions [33,34], where ''fine'' texture refers to a regular pattern and ''coarse'' texture corresponds to irregular tissue. Texture analysis of T 2 w lesions predicted poor recovery and mild ongoing tissue injury [35]. The study further indicated that recovery of acute lesions tends to be associated with the degree of coarse texture during enhancement. While these studies aimed at static characteristics related either to persistence or recovery of acute lesions over time, DTPA analysis first aims at the discrimination of how pathological processes influence the effect of bolus passage of Gadobutrolum on the selected subset of histogram TPs as a function of time and tissue type. Notably, discrimination between ELs, NELs and NAWM was performed on raw EPI images (thus incorporating the different static and dynamic tissue properties) during bolus passage and the difference image time series (TPM-DIFF) and its derived TPMs. In the latter, the different structural pathological components that influence the static signal components of the tissues are canceled out by the subtraction of the raw image at t 0 . Due to the large differences in the histogram-based TPMs over different time intervals and their tissue dependent specific response, DTPA provides numeric information on a continuous scale about the amount of BBB disruption in MS-lesions, instead of a classification in enhancing or non-enhancing lesions. Yet, a formal proof for a relationship between inflammation and the DTPA information is currently lacking.
Potential relevance of TPMs for the clinical evaluation of MS lesions:

TPM-ORIG/TPM-DIFF/TPM-AVE
The TPM-ORIG shows monotonously increasing values for ELs during IFP, OFP, and RPP. Yet, this TPM did not discriminate ELs from NELs due to: (a) hypo-and hyper-perfused ELs; and (b) large variations in leakage among ELs; the latter is demonstrated by the TPM-DIFF in Figure 6. The TPM-DIFF further indicates that even NELs continue to have a subtle BBB disruption since the TPM-DIFF is statistically significantly smaller than zero, indicating leakage. This leakage is not observed by visual inspection of T 1 w images, but is clearly observed with DTPA. Therefore the TPM-DIFF might be used as a novel gradual marker for BBB disruption that is otherwise discarded on routine imaging. The difference in TPM-DIFF (or TPM-AVE) parameter between IFP and RPP (indicated by DTPM-DIFF) obtrudes itself as a surrogate marker to grade the MS-lesion state: the higher the value of DTPM-DIFF, the larger is the grade of leakage. The DTPM-DIFF-value allows grading of the leakage of the MS-lesion on a continuous numeric scale, instead of an ordinal scale known from T 1 w imaging (enhancing/non-enhancing).
Further clinical studies are necessary to validate the added value DTPM-DIFF for the predictive evaluation of disease progression.

TPM-STDEV/TPM-VAR/TPM-VAVA
These TPMs are all measures of heterogeneity of the examined tissues. For IFP, OFP and RPP the values of the ELs were , NELs , NAWM. TPM-VAR/TPM-STDEV values increase slowly from IFP towards RPP for ELs, and strongly decrease for NELs and NAWM. For TPM-VAVA we observe reductions as a function of time for all three tissue types. A second DTPA-derived numeric MR marker is proposed from the difference in the numeric values obtained for the TPM-VAVA between IFP and RPP, and is indicated by DTPM-VAVA. The smaller DTPM-VAVA, the more severely inflammatory is the MS-lesion. This also applies to DTPM-STDEV and DTPM-VAR. Again, further research must be done to find out whether these parameters correlate with disease activity.

TPM-SKEW/TPM-KURT
Due to their low SNR (see discussion below) the comparison of the mean values between the different tissues did not reach statistical significance. TPM-SKEW decreased over time in ELs, NELs and NAWM. ELs showed the highest variation in its values including a sign reversal during the late stages. TPM-KURT decreases for NAWM and ELs towards RPP, whereas the ELs increased. For the used 5|5 voxel kernel it seems that based on these two texture types no meaningful surrogate markers can be defined. For higher kernel sizes it is likely to find significant differences between these tissue types as well.

Limitations of the study
This study aimed to introduce DTPA and to demonstrate its feasibility in a clinical setting. The female preponderance in our Figure 7. Direct computation and analysis of time dependent texture parameters. The DTPA software also allows direct computation and analysis of time dependent texture parameters without first having to compute the TPMs on a fixed pixel grid. In this figure the time dependent TPs of EL versus NAWM tissue are displayed. Note that each TP has its own SNR. doi:10.1371/journal.pone.0067610.g007 cohort is determined by the epidemiology in MS. The incidence of MS in women is doubling males, and there is still a ongoing disproportional increase [36]. We have investigated different types of lesions in a small and female dominated patient collective diagnosed with different subtypes, which is a drawback. We conclude that large differences in TPs between different tissue types can be detected by DTPA. As a next step, we aim to perform a randomized clinical cohort study in RR-MS to analyze the differences between early and late ELs, NELS and NAWM. SNR of the TPM-SKEW and TPM-KURTs, depends on both SNR of the EPI input images (TPM-ORIG) series and kernel size. In this study we choose a kernel size of 5|5 voxels to compute TPMs (see Figure 3 on which the low SNR of the TPM-SKEW and TPM-KURT can be seen). Improved SNR could be obtained by increasing the kernel sizes to be used (e.g. 8|8), at the cost of spatial resolution however. We identified this SNR/spatial resolution problem during this study and extended the software to allow direct computation of texture parameters from user defined TOIs without the detour via TPMs computed using a fixed kernel. Figure 7 shows the time dependency for an EL and NAWM and using the direct computation method. Note that each TP has its own SNR. Whether the direct computation of texture parameters without first computing the texture parameter maps over a fixed kernel is better than direct computation from the defined TOI is still an open issue and has to be addressed in a separate study.
Finally, this study only examined a limited number of texture parameters and focused on histogram based texture parameters. In a follow up study we will also incorporate gradient-based, cooccurrence-matrix based and run length matrix based TPMs to explore the differences in TPMs of various subtypes of MS lesions.

Conclusion
This paper investigated the dynamics of texture parameter evolution in time series of DSCE-images in 18 patients with CIS or MS according to the 2010 McDonalds criteria. A novel software program is introduced here to investigate time dependent texture parameter maps derived from DSCE (difference) images. In order to enable comparison of texture parameters between two patients, a twofold image-normalization was performed. The first normalization, to compensate for variations in coil-load, is obtained by setting the mean NAWM voxel value during BLP to a specific value in the ORIG-TPM. The second normalization, to compensate for differences in cardiac output, is obtained by normalization of the time integral of the TPM-DIFF over the IFP and OFP for the same NAWM. After normalization, comparisons between patients by comparing the time averaged texture parameter values during the BLP, IFP, OFP and RPP are feasible. We detected dynamic texture features that revealed highly statistically significant differences between ELs, NELs, and NAWM. Based on these dynamic TPs, novel grading parameters for MS lesions can be introduced allowing for grading of MS-lesions on a numerical scale instead of an ordinal scale as is the case with pre/post contrast T1w image analysis. Our data support the hypothesis that, dependent on the tissue type, subtle differences in microcirculation are present in enhancing and non-enhancing MS lesions.