Spatiotemporal Heterogeneity in Multiparametric Physiologic MRI Is Associated with Patient Outcomes in IDH-Wildtype Glioblastoma

Intratumoral heterogeneity in glioblastoma represents a biologically complex spatial variation in gene expression, histopathology, and macroscopic structure (1). Importantly, spatial heterogeneity frequently indicates a poor prognosis (2), because treatment response is not uniform across the tumor and resistance to therapy occurs in different tumor regions. The combination of genetic heterogeneity with metabolic differences are known to result in spatial habitats, which comprise subpopulations of cells preserving distinct characteristics (3–5). Tumor regions showing differences in metabolism, vascularity, and cellularity can be captured with physiologic MRI, including cerebral blood volume (CBV) and apparent diffusion coefficient (ADC) mapping, which provide data directly reflecting biologic information (6, 7).

Current tumor analysis using histograms (8) or radiomics analysis (9) focuses on quantifying the heterogeneity and complexity by calculating the relationship between voxels (5). However, the spatial distribution of voxels within the tumor are not considered, making biologic validation difficult (10). Identification of similar voxels differs in the underlying assumptions and methodology from that of radiomics, which may define multiple subregions at the voxel level with a common tumor biology (3, 5). This shows spatial correspondence with genomics and histopathology and provides clinically useful information for identifying subregions of tumor progression that may require biopsy, are resistant to current therapy, and may become a therapeutic target in the future.

Single-parametric changes, such as parametric response maps, have been used for predicting the treatment response of glioma (11–13). However, binary thresholds are not absolute and depend on the acquisition and analysis of images (5). In addition, tracking individual voxels is both extremely difficult and cannot directly identify tumor subregions based on multiparametric imaging parameters (5). A data-driven approach using voxel-wise clustering may identify functionally coherent tumor subregions based on multiparametric physiologic MRI (5, 14).

The pathophysiologic correlates of ADC and CBV are cell density and necrosis (15, 16), and vessel density (17, 18), respectively, but the measurement of a single quantitative parameter is limiting owing to tumor heterogeneity and treatment-related changes (19). We hypothesized that the temporal changes in spatial habits derived from multiparametric physiologic MRIs potentially provide information on both spatial and temporal heterogeneity in posttreatment glioblastoma. The technique does not require voxel-by-voxel matching for serial follow-ups and accounts for newly appeared or disappeared voxels by creating spatial habitats with clustering of similar ADC and CBV values. The purpose of this study was to develop multiparametric physiologic MRI-based spatiotemporal habitat analysis and to validate the association between temporal changes in these habitats with progression-free survival (PFS) after concurrent chemoradiotherapy (CCRT) in patients with glioblastoma.

This retrospective clinical study was approved by the institutional review board (IRB) of Asan Medical Center (local approval number: 2020-0363) and was conducted in compliance with the U.S. Health Insurance Portability and Accountability Act regulations and the Declaration of Helsinki. For all the patients in the validation set, written and signed informed consent was obtained as per IRB-approval protocol. Figure 1 shows the inclusion process. A total of 150 patients with a histopathologic diagnosis of isocitrate dehydrogenase (IDH)-wildtype glioblastoma according to the World Health Organization classification 2016 (1), and a contrast-enhancing lesion (CEL) on MRI after CCRT were identified from the neuro-oncologic database of our institution between April 2016 and August 2018. Patients were enrolled according to the following inclusion criteria: (i) minimum age of 18 years; (ii) underwent standard treatment including surgery, CCRT, and adjuvant temozolomide; (iii) underwent conventional and advanced MRI, including diffusion-weighted imaging (DWI) and dynamic susceptibility contrast (DSC) imaging; (iv) had a measurable CEL of more than 1 × 1 cm² on the first postchemoradiation examination; (v) clinicoradiologic follow-up with at least three successive imaging examinations or pathologic confirmation to monitor treatment response; and (vi) image quality with no artifacts affecting the analysis. Fifty-three patients did not meet the inclusion criteria and were excluded: 15 patients did not undergo DSC imaging, 3 had poor-quality imaging/perfusion data, and 35 did not undergo 1-year follow-up or more than three successive MRIs. Therefore, 97 consecutive patients [mean age, 58.1 years (range, 21–86); 48 (50.5%) women] were included in the study to determine the physiologic MRI-based spatial habitats.

An independent temporal validation set of IDH-wildtype glioblastoma was used to evaluate the true clinical performance. This independent dataset was from a prospective brain glioma registry (ClinicalTrials.gov identifier: NCT02619890) and included patients diagnosed and followed up between September 2018 and May 2020. Informed consent was obtained from all patients in a prospective brain glioma registry. The final validation set consisted of 53 patients [mean age, 57.8 ± 13.4 years (range, 32–74); 28 (52.8%) women].

The training and validation set shared the same imaging acquisition protocol. The brain tumor imaging protocol was acquired on a 3-T scanner (Ingenia 3.0 CX, Philips Healthcare) and included conventional and advanced sequences, including T2-weighted imaging (T2WI), T2-weighted fluid-attenuated inversion recovery (FLAIR), and precontrast and postcontrast T1-weighted imaging (T1WI). The FLAIR and contrast-enhanced T1WI imaging parameters are shown in Supplementary Table S1.

The DWI parameters were as follows: repetition time (TR)/echo time (TE), 3,000/56 ms; diffusion gradient encoding, b = 0 and 1,000 seconds/mm²; field of view (FOV), 250 × 250 mm; matrix, 256 × 256; and slice thickness/gap, 5/2 mm. ADC images were calculated from b = 1,000 and b = 0 second/mm² DWI images.

DSC imaging was acquired using a gradient-echo echo-planar imaging protocol. A preload of 0.01 mmol/kg gadoterate meglumine (Dotarem; Guerbet) was administered, followed by a dynamic bolus of a standard dose of 0.1 mmol/kg gadoterate meglumine delivered at a rate of 4 mL/second by an MRI-compatible power injector (Spectris; Medrad). The contrast material bolus was followed by injecting 20 mL saline at the same injection rate. The DSC imaging parameters included the following: TR/TE, 1,808/40 ms; flip angle, 35°; FOV, 24 × 24 cm; slice thickness/gap, 5/2 mm; matrix, 128 × 128; total acquisition time, 1 minute and 54 seconds. The dynamic acquisition was performed with a temporal resolution of 1.5 seconds, and 60 dynamics were acquired. DSC imaging was acquired with complete tumor volume coverage and the same section orientation as that used in conventional MRI. The MRI examinations performed at the first visit after CCRT (examination 1, Post-CCRT #1) and the second visit after CCRT (examination 2, Post-CCRT #2) were used in the analysis.

For the three-dimensional contrast-enhanced T1WI (3D CE-T1WI) and FLAIR data, skull stripping was performed with an algorithm (https://github.com/MIC-DKFZ/HD-BET) optimized for processing heterogeneous MRI data with varying degrees of pathology or posttreatment alterations. For all patients, lesion segmentation masks were generated on the 3D CE-T1WI and FLAIR data using 3D UNet-based method (https://github.com/MIC-DKFZ/nnUNet; ref. 20) from the PyTorch package version 1.1 in Python 3.7 (www.python.org).

For the DSC analysis, pharmacokinetic map calculation was performed using Nordic ICE (NordicNeuroLab). Using the integrated DSC module, which incorporates a relative CBV (rCBV) leakage-correction algorithm and manual noise thresholding, the amount of blood in a given volume of tissue was obtained and expressed as mL per 100 mL tissue. The calculations using the Weisskoff–Boxerman method based on time-dependent deviation of the pixel-wise concentration–time curve from a reference curve were assumed to be unaffected by leakage (21). The rCBV maps were normalized according to the normal-appearing white matter, to create normalized CBV (nCBV) maps.

The 3D CE-T1WI images of each patient in the training dataset were coregistered and resampled into isometric-voxel sizes to evaluate changes across consecutive scans. Then, the FLAIR, nCBV, and ADC images were coregistered and resampled to the iso-voxel CE-T1WI images using rigid transformations with six degrees of freedom in SPM package (version 12, www.fil.ion.ucl.ac.uk/spm/). This step allowed continuous slices without gap and enabled tracking of habitats by voxel-wise analysis.

The final voxel classifications based on nCBV and ADC values were implemented using a k-means clustering module in the scikit-learn python package.

All voxels from all segmented masks in the training set were first aggregated; then, all voxels in the ADC and nCBV maps were split into three clusters, which were functionally coherent subregions within a CEL. The individual voxels in each cluster were grouped on the basis of their similarities and differences using the k-means clustering algorithm with squared Euclidean distances between voxel intensities as the similarity metric. Because the optimal number of clusters in a dataset is an important issue in k-means clustering, 2, 3, 4, and 5 clusters were initially set. The aggregation of all voxels from all enhancing lesions in the primary cohort and the three clusters are demonstrated in Supplementary Fig. S1. The three clusters demonstrated both ADC and nCBV effects, while the fourth cluster was dominated by nCBV and the fifth cluster by ADC and displayed a narrow range of clusters. Finally, three clusters were chosen based on ADC and nCBV differences and to avoid overparameterized models (22).

Three clusters were set using two distinct feature maps: cluster 1 represented “hypervascular cellular tumor” with high CBV value and low ADC value; cluster 2 represented “hypovascular cellular tumor” with low CBV value and low ADC value; and cluster 3 represented “nonviable tissue” with low CBV value and high ADC value. All voxels were allocated into one of the three clusters and were displayed as spatial habitats in the original image space. The range for boundary of spatial habitats was 4.37–4.44 for nCBV and 150–187 (×10⁻⁶ mm²/s) for ADC.

The two consecutive MR examinations acquired at the first visit after CCRT (examination 1, Post-CCRT #1) and the second visit after CCRT (examination 2, Post-CCRT#2) were used in the analysis. First, changes in the numbers of voxels in the entire enhancing lesion and in each habitat were calculated between the two examinations (examination 2 − examination 1). Second, percentage of each habitat to the CEL volume was calculated (number of voxels in habitat 1/number of voxels in CEL volume) and change in the percentages between the two examinations (examination 2 − examination 1) was calculated.

Clinical, radiologic, and pathologic data were retrieved from medical records. Baseline characteristics including sex, age, Karnofsky performance score (KPS; binary, score >70 or ≤70), initial tumor volume, O6-methylguanine-DNA-methyltransferase (MGMT) promoter status, and extent of surgery (biopsy, partial resection, or gross total resection) were collected.

The primary endpoint was PFS. Tumor progression on MRI was assessed according to the Response assessment in neuro-oncology (RANO) criteria (23), using serial measures of the product of the two largest cross-sectional diameters. An objective response was defined as a complete or partial response on two consecutive MRI acquisitions with reduced or stable doses of corticosteroids. The consecutive clinicoradiologic diagnoses were made by consensus between a neuro-oncologist (J.H. Kim, with 26 years of experience in neuro-oncology practice) and a neuroradiologist (H.S. Kim, with 21 years of experience in neuro-oncologic imaging) after complete imaging and medical chart review (23). The overall survival (OS) was calculated from the day of diagnosis until the day of death, using the information obtained from the national health care data linked to our hospital.

RANO assessment was separately performed by two central reviewers (Y.-H. Kim and J.E. Park with 10 years of experience in neuro-oncologic practice and imaging, respectively), who reached consensus and had not participated in any other imaging review. The RANO assessment was based on postoperative imaging and examination 2 (Post-CCRT #2), and the C-index was compared with that from spatiotemporal habitat.

The site of progression was analyzed using the follow-up examination at the time of progression. The CEL volume at the time of progression was matched to the habitats at examination 2 (Post-CCRT #2). The overlap of each spatiotemporal habitat with the CEL volume at the time of progression was calculated using the DICE similarity coefficient, |{\rm{DICE}} = {\frac{{2| {P \cap R} |}}{{| P | + | R |}}$| (24). P indicates each spatiotemporal habitat and R indicates the CEL volume at the time of progression. The DICE ranges between 0 (no overlap) and 1 (perfect agreement).

The methodology written here adheres to the REMARK (Reporting Recommendations for Tumor Marker Prognostic Studies) recommendations (25). The sample size of the validation set was determined according to power analysis for a Cox regression of the log HR, which was performed using Power Analysis Software (PASS) version 15.0.7.

Univariate analysis was performed for analyzing the association of spatiotemporal habitats with PFS using Cox regression or the Kaplan–Meier method (log-rank test). HRs indicate relative change in hazard incurred by 1 unit increase in each parameter, and 20,000 voxels (20 k voxels) were considered as 1 unit in this study. For the significant predictors identified in the univariate Cox regression, an optimal cutoff for stratifying high- and low-risk groups was estimated using maxstat in R with 10-fold cross-validation, which ensured unbiased prediction within the sample (26). The z-score is the ratio of each regression coefficient to its SE, a Wald statistic which is asymptotically standard normal under the hypothesis that the corresponding β is zero (27). The z-score was used to demonstrate statistical significance of each spatiotemporal habitat.

Using this cut-off from the exploratory analysis, a habitat risk score was developed to realize risk stratification for patients; each predictor was assigned a discrete score of 1 if it was larger than its cut-off, and 0, if it was lesser. For example, when two habitats relevant to PFS showed a discrete increase above the cut-off, the habitat risk score was 2.

Then, the C-index was calculated using independent predictors from multivariate Cox regression analyses to predict survival, by combining the habitat risk score and clinical predictors.

All results are reported as median with range or 95% confidence interval (CI) for continuous variables, and as frequencies or percentages for categorical variables. Differences in clinicopathologic factors between the training and validation cohorts were assessed using the χ² test for discrete variables and Student t test for continuous variables. All statistical tests were two-sided, and a P value < 0.05 was considered statistically significant. Statistical analyses were performed using R software version 3.4.3 (https://www.r-project.org).

The clinical characteristics of the patients in the training and validation sets are summarized in Table 1. No differences in sex, age, initial KPS, treatment regimen, and MGMT promoter methylation status were noted between the training and validation cohorts.

The data from 42 patients (28%) who were alive at the time of analysis were censored and included. The median OS time of the deceased patients was 15 months (interquartile range, 15–29 months) in the training set and 13 months (interquartile range, 6–18 months) in the validation set. The median follow-up periods were 39 months and 19 months for patients in the training set and validation set, respectively.

The mean interval between the two serial examinations (Post-CCRT#1 and #2) was 75 ± 22 days in the training set (mean ± SD) and 67 ± 22 days in the validation set (mean ± SD); no significant difference was observed between the two datasets.

The results of the univariate analysis to evaluate the change in three spatiotemporal habitats are summarized in Table 2. The distribution of the average percentage of the three spatiotemporal habitats (hypervascular cellular cluster, hypovascular cellular cluster, and nonviable tissue cluster) was similar between examination 1 and 2; for examination 1, the values were 25.0%, 50.1%, and 25.9%, respectively, and for examination 2, 23.5%, 52.1%, and 24.4%, respectively, for each habitat (Supplementary Fig. S2). The habitat characteristics are shown in Supplementary Table S2.

Among the various changes in spatial habitats, an increase in the hypovascular cellular habitat (20 k voxels) was most significantly associated with a relatively short PFS (HR, 3.78; 95% CI, 2.38–6.00; P < 0.001). An increase in the hypervascular cellular habitat (HR, 40.0; 95% CI, 4.71–339.2; P = 0.001; z-score = 3.38) and percentage of the hypovascular cellular habitat were also associated with a relatively short PFS (HR, 1.01; 95% CI, 1.00–1.02; P = 0.02; z-score = 2.31); however, the association was weaker than that observed for the increase in the hypovascular cellular habitat (z-score = 5.6). Figure 2 demonstrates the association of PFS with change in the hypovascular cellular habitat. An increase in the hypervascular cellular habitat or nonviable tissue habitat was not associated with PFS (P > 0.05).

The optimal cut-off of the spatiotemporal habitat needed for stratifying the longer and shorter PFS groups was an increase of more than 130 voxels in the hypovascular cellular habitat. The cut-off stratified the longer and shorter PFS groups with significant difference in the log-rank test (P = 0.016). Supplementary Figure S3 shows the Kaplan–Meier survival curve and the risk table generated on the basis of the hypovascular cellular habitat. The optimal cutoff for the hypervascular cellular habitat was an increase of more than 0 voxel. The habitat risk score, based on the discrete increase of the hypervascular and hypovascular cellular habitats was calculated for the patients. For example, when a patient showed an increase of both hypervascular cellular habitat (>0 voxel) and hypovascular cellular habitat (>130 voxels), the habitat risk score was 2.

In the training set, the habitat risk score stratified patients according to low risk (median PFS, 518 days; 95% CI, 441–672 days), intermediate risk (median PFS, 316 days; 95% CI, 180–527 days), and high risk (median PFS, 160 days; 95% CI, 149–279 days; log-rank test; P < 0.001).

A multivariate Cox analysis was conducted with the habitat risk score, age, KPS scores, extent of surgery, and the initial tumor volume as the factors. Only a high initial KPS score (HR, 0.43; 95% CI, 0.25–0.75; P = 0.003) and the habitat risk score (HR, 2.17; 95% CI, 1.59–2.96; P < 0.001) were identified as the independent predictors of PFS. The C-index values for PFS and OS were 0.79 (95% CI, 0.72–0.85) and 0.67 (95% CI, 0.64–0.71), respectively.

The survival prediction value for our model was compared with that of the null model. The hypovascular cellular habitat showed HR 3.78, and HR 3.00 was conservatively set as the meaningful difference. The event rate (progression) was 0.7 and SD of the hypovascular cellular habitat was 0.9 in the training set. A sample size of 49 would achieve 80% power to have meaningful survival prediction based on HR 3.00, with an estimated SD of differences of 1.0 with a significance level (alpha) of 0.05 using a Cox regression model. Thus, the sample size of 53 in the validation set was sufficient for predicting survival.

The three spatial habitats displayed similar patterns in the validation cohort in contrast to those displayed in the training cohort. The habitat risk score successfully stratified patients according to low risk (median PFS, 329 days; 95% CI, 229–518 days), intermediate risk (median PFS, 122 days; 95% CI, 86–180 days), and high risk (median PFS, 80 days; 95% CI, 78–82 days; log-rank test; P = 0.02; Fig. 3). In combination with initial KPS score and the habitat risk score, the C-index values for PFS and OS were 0.79 (95% CI, 0.72–0.86) and 0.74 (95% CI, 0.70–0.78), respectively.

The C-index for the RANO-based prediction was 0.58 (95% CI, 0.56–0.61) in the training set and 0.54 (95% CI, 0.50–0.58) in the validation set. The C-index from the habitat risk score only was 0.67 (95% CI, 0.64–0.71) in the training set and 0.71 (95% CI, 0.66–0.75) in the validation set, which was significantly better at predicting OS (P < 0.001) than that from RANO assessment.

The tumor progressed in 46 patients (86.7%) in the validation set. Figure 4 demonstrates the relationship between the temporal changes in the hypovascular cellular habitat and the site of progression. The DICE index for the site of progression was calculated between examination 2 and the time of recurrence (examination 3). The mean DICE index was 0.13 (range, 0–0.6; SD 0.10) for hypervascular cellular habitat, 0.34 (range, 0–0.72; SD 0.12) for hypovascular cellular habitat, and 0.07 for nonviable tissue (range, 0–0.18; SD, 0.07). The DICE index for the hypovascular cellular habitat was highest in 34 of the 46 cases (73.9%).

In this study, we demonstrated that PFS could be predicted in patients with IDH-wildtype glioblastoma using the spatiotemporal habitats, which were defined by multiparametric physiologic MRI values, quantified in the early stages of pathophysiologic alterations occurring after concurrent chemotherapy. An increase in the hypovascular cellular habitat with both low ADC and CBV values showed the most significant association with PFS. Furthermore, the increase in both the hypervascular and hypovascular cellular tumor was scored and stratified with a habitat risk score, and the patients in a prospective registry-based validation set were successfully stratified into groups with low, intermediate, and high risk of early progression.

The heterogeneous nature of glioblastomas after treatment (28) poses significant challenges in their evaluation using single quantitative parameters; such binary approaches are often affected by different image settings and algorithms (5). Furthermore, a threshold-based approach has an intrinsic limitation because it is difficult to label all tissue voxels when considerable heterogeneity is present, which results in tissue “gray-zones” (29). Alternatively, the clustering of voxels is advantageous because it requires no prior assumption (30) and is not affected by acquisitions that are normalized by the normal brain tissue. The centroids in k-means clustering can be stabilized without using a cut-off value (31); this approach has better implications for minimizing gray-zone tissue and applications in different imaging settings.

Temporal changes in physiologic MRI parameters when added to the information of spatial heterogeneity can enable clinical outcome prediction. Besides treatment-related changes such as low perfusion and high diffusion, the posttreatment glioblastoma tissues may also exhibit a spectrum of low-to-high perfusion and low diffusion (32). The hypervascular cellular habitat, with high CBV and low ADC values, represents the most malignant region of a tumor. However, the association between the change in the hypervascular cellular habitat with PFS was weaker than that of the hypovascular cellular habitat with PFS. Temporal changes in spatial habitats, such as “spatiotemporal habitats,” would be more helpful than a single snapshot of the posttreatment glioblastomas to predict PFS because of the inherent heterogeneity in the histopathology of posttreatment glioblastoma (33). Moreover, longitudinal analysis of an entire tumor with physiologic MRI might help to overcome this ambiguity and provide adequate information for subsequent treatment decisions.

The spatiotemporal habitat risk score was significantly better than the RANO-based assessment at predicting OS. This result is similar to that of a recent study based on the EORTC-26101 trial that showed the RANO assessment produced C-index of 0.57 (95% CI, 0.54–0.61) for predicting OS (34). Recent deep learning segmentation or T1-subtraction methods have emphasized the importance of objective tumor volume measurements for determining PFS (34, 35). As most glioblastomas recur locally (36, 37), a spatial habitat-driven approach may refine the volumetric analysis to predict their progression. From the results of our analysis, we speculate that ADC and CBV subdivides the contrast-enhancing tumor volume and this may provide more sensitive and specific information of tumor progression. Using the habitat risk score described herein, the change in spatiotemporal habitats can be easily determined, which might help to stratify patients with glioblastoma according to the risk of early progression.

Not only did the spatiotemporal habitats improve general survival performance, they also highlighted the importance of identifying subregions to determine intratumoral heterogeneity and identify regions with treatment resistance. We observed that the tumor recurrence foci matched the regions that showed short-term increase in the hypovascular cellular tumor habitat using quantitative analysis of the DICE index. The hypovascular cellular tumor habitat showed the highest DICE index with CEL at the time of progression. Parcellations ensure that the images are obtained at a single time point and that the functionally coherent regions can be followed up and matched to the tumor recurrence foci.

This study has several limitations. First, the study included a small cohort size, especially in the validation set. Second, strict pathologic correlations with image-based segmentations were lacking. However, such correlations would be very difficult to achieve, and unnecessary surgery should be avoided; therefore, we made assumptions regarding the clusters using a logical approach. Third, the three habitats defined here are relative indices of cellular habitats showing relatively low ADC values compared to those of other habitats. The b values selected in our study were only b = 0 and b = 1,000 seconds/mm², thus, there is a risk of perfusion contamination in ADC values. One strategy is to obtain images with b value = 150–200 and 400–600 seconds/mm² and use them to fit the data (38). This will provide more cellular components with less perfusion contamination, which is recommended in the future study. Fourth, although the study results were validated on a separate, prospectively enrolled cohort, this was a single-center study; thus, application of different protocols from other institutions with a relatively large number of patients is necessary. The normalization of quantitative values of ADC and CBV holds potential to increase reproducibility in a multicenter study. However, there is a big gap between taking the data as presented and turning them into a clinically useful tool that could guide patient therapy until it is proved in a multicenter setting, considering the cut-off for the habitat risk score depends on the training set. In the future, a prospective multicenter study using ADC and CBV in comparison with a RANO-based assessment to predict OS is further needed to confirm its utility. Finally, the technology used in this study is complex; the development of an easy-to-use software is warranted for widespread clinical adoption. A specialized postprocessing software for spatial habitat analysis and deep learning segmentation is being developed for clinical implementation.

In conclusion, the temporal change in spatial habitats derived from multiparametric physiologic MRIs may be a useful predictor of clinical outcome in patients with posttreatment glioblastoma. An increase in the hypovascular cellular habitat can be regarded as a strong imaging biomarker to predict early progression, as well the site of progression.

No disclosures were reported.

J.E. Park: Conceptualization, formal analysis, methodology, writing-original draft. H.S. Kim: Conceptualization, supervision, methodology, writing-review and editing. N. Kim: Software, formal analysis, investigation, visualization. S.Y. Park: Software, formal analysis. Y.-H. Kim: Data curation, project administration. J.H. Kim: Resources, project administration.

This research was supported by a National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIP; grant numbers: NRF-2020R1A2B5B01001707 and NRF-2020R1A2C4001748).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.