Morphological and Fractal Properties of Brain Tumors

Tumor interface dynamics is a complex process determined by cell proliferation and invasion to neighboring tissues. Parameters extracted from the tumor interface fluctuations allow for the characterization of the particular growth model, which could be relevant for an appropriate diagnosis and the correspondent therapeutic strategy. Previous work, based on scaling analysis of the tumor interface, demonstrated that gliomas strictly behave as it is proposed by the Family-Vicsek ansatz, which corresponds to a proliferative-invasive growth model, while for meningiomas and acoustic schwannomas, a proliferative growth model is more suitable. In the present work, other morphological and dynamical descriptors are used as a complementary view, such as surface regularity, one-dimensional fluctuations represented as ordered series and bi-dimensional fluctuations of the tumor interface. These fluctuations were analyzed by Detrended Fluctuation Analysis to determine generalized fractal dimensions. Results indicate that tumor interface fractal dimension, local roughness exponent and surface regularity are parameters that discriminate between gliomas and meningiomas/schwannomas.


INTRODUCTION
Tumor interface exhibits a complex and irregular geometry due to the dynamics involved in the tumor growth process, which in general takes into account tumor cell proliferation and invasion into the surrounding tissue. To characterize its complexity, fractal analysis has been used routinely for tumor detection (Iftekharuddin et al., 2009) and therapy monitoring (Di Ieva et al., 2012). In the case of brain tumors, magnetic resonance imaging techniques give detailed geometrical information with excellent spatial resolution and quality for the evaluation of the tumor interface. Parameters extracted from the tumor interface by scaling and fractal analysis have given relevant clues about the complex tumor growth dynamics (Brú et al., 2012) (Brú et al., 2003) (Brú et al., 2008) (Torres Hoyos and Martín-Landrove, 2012) which in turn can be used to further validate tumor growth models (Brú et al., 2014) for therapy simulation and prognosis. In a previous work (Martín-Landrove et al., 2020) it was demonstrated that scaling analysis provided a clear difference in the tumor growth model for gliomas, which follow a ballistic growth model in completely agreement with the Family-Vicsek ansatz (Family and Vicsek, 1991) (Barabasi and Stanley, 1995), compared to meningiomas/schwannomas. A different approach that includes fractal properties of the tumor interface or surface has been proposed by defining surface regularity measures (Pérez-Beteta et al., 2018) (Popadic et al., 2021). So far, scaling analysis parameters, such as fractal dimension and local roughness exponent, and surface regularity measures give a global picture of the fractal properties of the tumor interface. Since the tumor growth process is heterogeneous it is expected that its fractal properties should be heterogeneous as well and a more general approach is needed and general multifractal analysis methods (Kantelhardt et al., 2002) (Lopes and Betrouni, 2009) (Gu and Zhou, 2006) applied to the fluctuations over the interface have to be used. In the present work, an extended image database is used to determine an extended group of parameters that characterize the tumor interface, including those previously determined (Martín-Landrove et al., 2020).
Related to the assessment of the tumor interface, several methods have been proposed for the segmentation of brain tumors (Işın et al., 2016) (Wadhwa et al., 2019) (Zhao et al., 2020) from magnetic resonance images, which include conventional methods such as thresholding and region growing, supervised methods mainly represented by support vectors machines and artificial neural networks and unsupervised methods which include clustering methods and deformable models. Even though supervised methods perform better than unsupervised ones (Rao et al., 2018), for the purpose of the present work, unsupervised methods, such as K-means or Fuzzy C-means are preferable since these methods do not need any training set and in many cases are simpler in its numerical implementation. In the present work, an unsupervised method based on dynamic quantum clustering (Weinstein and Horn, 2009), (Horn and Gottlieb, 2001) is used for image segmentation to determine the tumor interface.
The article is organized as follows, in Section 2 it is described the selection of images for this study, the segmentation method employed for determination of the tumor interface and the different morphological parameters that describe the tumor interface such as fractal dimension and lacunarity, growth dynamics exponents, regularity measures, complex visibility graphs and parameters derived from multifractal analysis. The results are discussed in Section 3 and the conclusions are presented in Section 4.

Image Selection
Images for high grade gliomas were extracted from different collections in The Cancer Imaging Archive (Chang et al., 2011), (Clark et al., 2013); The Cancer Genome Atlas Low Grade Glioma (TCGA-LGG) data collection (Pedano et al., 2020), the Repository of Molecular Brain Neoplasia Data (REMBRANDT) (Scarpace et al., 2019) for astrocytomas and oligodendrogliomas of grades 2 and 3, and The Cancer Genome Atlas Glioblastoma multiforme [TCGA-GBM] collection (Scarpace et al., 2016) for glioblastoma multiforme. Also, data coming from the RSNA-ASNR-MICCAI Brain Tumor Segmentation (BraTS) Challenge 2021 (Baid et al., 2021) (Menze et al., 2015) (Bakas et al., 2017). For other brain neoplasias, such as meningiomas and acoustic schwannomas, local image datasets were used. Among these collections, only contrast enhanced T1-weighted images, with tumor lesions clearly identified as such and separated from anatomical structures, were selected for analysis.

Clustering of Data and Image Segmentation
Image digital levels were clustered using the Dynamic Quantum Clustering algorithm (DQC) (Weinstein and Horn, 2009) (Horn and Gottlieb, 2001) which assumes that data are described by a set of points, each one defined with some uncertainty, σ and the distribution for all points in data space is given by a Parzen estimator, φ, which satisfies the time independent Schrödinger equation for its ground state, which leads to the evaluation of an energy potential The number of potential minima (Horn and Gottlieb, 2001) was previously used (Martín-Landrove et al., 2020) to determine the number of classes in a K-means algorithm alone. In the present work, the full dynamic quantum clustering algorithm is used and by Ehrenfest theorem, digital levels evolve in time toward the potential minima according to the following equation of motion (Lafata et al., 2018), which is a second order Langevin equation with a dissipative term γ. The clustering process is performed with a suitable selection of parameters such as σ 2 , which is an equivalent to a "mass", and determines the number of classes among the digital levels (Horn and Gottlieb, 2001), dissipation, γ and time interval. These parameters are estimated according to the magnitude of the digital levels in the image and assess the evolution of all the data to the potential minima. An example of clustering of digital levels is shown in Figure 1. For image segmentation the dynamics is performed as a two step process (Sánchez and Martín-Landrove, 2021): Dynamic A: First application of the dynamics upon the original image using the potential calculated from the original image histogram. At the end of this step, a Parzen estimator is evaluated using the clustered image histogram, allowing for the calculation of a "trap" potential, which defines the number of classes if a further classification algorithm, such as K-means, is to be used. It has been proved that the dynamic quantum clustering algorithm provides the same set of centroids as the K-means algorithm (Sánchez and Martín-Landrove, 2021).
Dynamic B: Second application of the dynamics upon the original image using the "trap" potential.
In Figure 2 it is shown how the selection of the appropriate σ 2 determines the number of classes and therefore the image segmentation. It is important to note that the "trap" potential compared to the original Schrödinger potential, exhibit a well defined minima, which allows for an unsupervised definition of the number of clusters and an improvement in the assessment of the tumor interface.

Fractal Dimension and Lacunarity of the Tumor Interface
Fractal dimension (Iftekharuddin et al., 2003) and lacunarity (Plotnick et al., 1993) have been used as morphological parameters to characterize tumor interface and grading of brain tumors (Smitha et al., 2015) (Park et al., 2020). Both quantities can be calculated using a box counting algorithm and are defined as with N (ϵ), the number of boxes containing the fractal object and ϵ the size of the box. In a similar way, lacunarity is defined as FIGURE 2 | Dependence of the segmentation procedure on the value of the 'mass' of the particle, σ 2 . Arrows indicate how the potential energy minima, in the 'trap' potential, collapse as σ 2 is increased from top to bottom by a 3-fold factor. On the right, the corresponding segmentation.
where μ (ϵ) is the mean of point density within the box of size ϵ and σ (ϵ), the standard deviation.

Scaling Analysis and Tumor Growth Dynamics
Tumor interface, in both resected and in vitro samples (Brú et al., 2012) (Brú et al., 2003), and in vivo (Torres Hoyos and Martín-Landrove, 2012) (Martín-Landrove et al., 2016) (Martín-Landrove et al., 2020) has been characterized using scaling analysis techniques. These studies have shown that tumor contours exhibit super-rough scaling dynamics described by the Family-Vicsek ansatz (Family and Vicsek, 1991) (Barabasi and Stanley, 1995) that corresponds to a ballistic growth process or a proliferative-invasive tumor growth model. The roughness of the tumor interface can be parameterized at the global level, with an exponent α, and at the local level by a roughness exponent, α loc . In three dimensions, the local roughness exponent relates the scale-averaged width of the interface between tumor and host to the scale of growth s, exhibiting a power-law behavior for small s W s ( )~s α loc (6) with W given by (Brú et al., 2008), where < r i > s represents the average of the radius, measured from the tumor center, over a patch of scale s located at the tumor interface, and * { } Σ represents the average over all realizations (all possible patches of scale s) over the interface surface Σ. In order for the growing process to follow the Family-Vicsek ansatz (Family and Vicsek, 1991), fractal dimension and local roughness exponent are related in a general way (Family and Vicsek, 1991) (Barabasi and Stanley, 1995), i.e., their sum is equal to the embedding dimension of the shape, or Euclidean dimension, d E , Also, the saturation value of the interface width, W sat , scales with average of the tumor size, 〈R〉 as where α is the global roughness exponent.

Surface Regularity Measures
A surface regularity measure has been proposed to characterize glioblastoma multiforme (Pérez-Beteta et al., 2018) by the following equation, where TV is the total volume of the tumor and TV eq is the volume of a sphere which has the same surface area as the tumor, TS. Thus In a similar way, a surface factor has been proposed (Popadic et al., 2021) for grading meningioma tumors. the surface factor is defined as, where TS eq is the surface area of an sphere with a volume equal to the tumor volume, TV and TS is the surface area of the tumor. In terms of TV and TS, S F can be written, So any of the proposed regularity measures S R or S F are equivalent for the description of the surface regularity and therefore S R will be used in the present work. According to (Pérez-Beteta et al., 2018), S R is related to the fractality or roughness of the tumor surface, i.e., if S R ≪ 1, tumor exhibits a distinct fractal dimension with a rough surface, while otherwise its fractal dimension is close to the Euclidean one and the tumor surface is smooth. Nevertheless, the actual value of S R must be corrected by a shape factor due to the fact that the distribution of points over the tumor surface could be elongated along certain directions, departing from the roughed sphere condition, f shape is calculated in a likewise way as S R , where V ellipsoid is the volume of an ellipsoid obtained by principal component analysis of the point distribution over the tumor surface and V sphere is the volume of a sphere with the same surface area as the ellipsoid. For the tumor interface, the Fractional Anisotropy is defined as, where λ i , i = 1, 2, 3 are the ellipsoid axes. In Figure 3 it is shown how FA determines the tumor interface correction that leads to Equation 14. Another factor can be defined which takes into account the distribution of contrast inside the tumor, where TV C is the volume of the region with contrast and TS C is its total surface area, including inner and outer surfaces.

Ordered Series and Visibility Graphs
Ordered series and the associated visibility graphs (Lacasa et al., 2008) (Lacasa et al., 2009) that can be extracted from the tumor interface has been used to further discriminate between the dynamical models that describe tumor growth (Brú et al., 2014). In a previous work (Martín-Landrove et al., 2020), it has been demonstrated that the associated visibility graph, through its connectivity distribution function P (k) discriminates between gliomas and meningiomas/ schwannomas, in the exponent of its power law behavior, In the present work, the analysis of the ordered series and its associated graph is extended one step further to determine its scaling properties (Pérez-Beteta et al., 2018). Similarly to Equation 7, the local standard deviation of the vertex degree for a subset ϕ on the ordered series can be written as (Brú et al., 2014) (Estrada, 2010), It exhibits a power-law behaviour for small ϕ (Brú et al., 2014),  where a represents the local variance exponent and for ϕ = 2π, W k is related to the global variance or heterogeneity of the associated visibility graph, as shown in Figure 4.

Multifractal Analysis
In order to determine the multifractal scaling exponents for the one-dimensional ordered series extracted from the tumor interface, a general procedure of fluctuation analysis is used (Kantelhardt et al., 2002) (Lopes and Betrouni, 2009). First, the profile of the ordered series is determined by a cumulative sum, where 〈r〉 represents the mean radius of the tumor interface. The profile series is then partitioned into N s ≡int (N/s) segments of equal length s, the box probability p s (v), which is the sum of the values r k within each segment v of size s, is defined as,    The scaling properties and exponents can be obtained through the partition function, For large values of s, a power law behavior is obtained for Z q (s), allowing for a definition of the scaling exponent τ 1 (q) that characterizes the one-dimensional fluctuations, The generalized fractal dimensions are then defined as, and the generalized Hurst exponents can be obtained from the relation, The ordered series that can be extracted from any slice represent a one dimensional sampling of the tumor interface and as a consequence an incomplete picture of the tumor interface fluctuations. A more general approach is possible if a two dimensional detrended fluctuation analysis (Gu and Zhou,     2006) is performed. In this case, the tumor interface is parameterized as a two dimensional array with elements r (n ϕ , n Z ), the radii of the tumor interface, and is partitioned in two dimensional segments of size s. The detrended fluctuation in each segment is given by, where ϵ v,w (i, j) is the difference between the cumulative sum of r (i, j) and its trend over the segment (v, w). The average of the detrended fluctuation over all the segments is, for q ≠ 0 and for q = 0, For large values of s, F q behaves as a power law,   The multifractal nature of the fluctuation is characterized by the scaling exponents τ 2 (q) and related to the exponents h 2 (q) by, and Equation 25 holds for the generalized fractal dimensions D 2 (q).

RESULTS AND DISCUSSION
A total of 609 tumor interfaces were analyzed, discriminated as follows, 176 benign tumors including 99 meningiomas and 77 acoustic schwannomas, all of them coming from local databases, 46 Grade II and Grade III astrocytomas and oligodendrogliomas (Pedano et al., 2020) (Scarpace et al., 2019) and 387 glioblastoma multiforme Grade IV tumors (Scarpace et al., 2016) (Baid et al., 2021). The tumor interfaces were selected by its size, i.e., the number of points in the tumor interface must be greater than a certain threshold, allowing for adequate statistics in the evaluation of morphological parameters such as fractal dimension d F or local roughness exponent, α loc . Also, the distribution of points must be as isotropic as possible, with values of the fractional anisotropy, equation (Işın et al., 2016), closest to zero. This is performed by the procedure described in Section 2.3.3 and shown in Figure 3.

Scaling Analysis Results
Results are shown in Figure 5 which exhibits a high dispersion of the data. Average values are indicated by large circles and correspond from left to right to meningioma, acoustic schwannoma, Grade II and Grade III glioma TCGA-LGG and REMBRANDT databases), glioblastoma multiforme and high grade glioma (BraTS 2021 database) and glioblastoma multiforme (TCGA-GBM database). These average values are  6 | Generalized fractal dimensions obtained by Detrended Fluctuation Analysis on one dimensional ordered series, D 1 (1) and D 1 (2), and two dimensional interface data space r (n ϕ , n Z ).  , 2016). This fact is also shown in Figure 6 in the comparison of the corresponding histograms for d F and α loc . The values of the local roughness exponent α loc and the fractal dimension d F give important information about what proliferative-invasive process is taking place in the dynamics of tumor growth, if the sum of these parameters is close to the Euclidean dimension, d E , as reflected by Equation 8, the tumor growth dynamics corresponds to a ballistic growth model, characterized by the ansatz of Family-Vicsek (Family and Vicsek, 1991) ( Barabasi and Stanley, 1995 Figure 7 which depicts the dependence of W sat with respect to the average size of the tumor lesion, 〈R〉. Trend of the data conforms to a power law according to Equation 9 and the results are summarized in Table 2. There is a clear difference between tumor groups, for gliomas, α = 0.948 ± 0.038, and for meningiomas and acoustic schwannomas, α = 0.730 ± 0.087, which is what has to be expected if the growth dynamics governing the tumor interface corresponds to a more invasive process as it happens to be the case for malignant tumors.

Regularity Measures Results
The relationship between S C and S R are shown in Figure 8. Data points are dispersed below the diagonal since S C takes into account both inner and outer surfaces while S R only takes into account the outer surface (Pérez-Beteta et al., 2018). Data points along the diagonal correspond to tumors that lack the presence of either contrast free or necrotic volumes, a condition that occurs more frequently for meningiomas and acoustic schwannomas than for gliomas as seen in Figure 8. Large circles in Figure 8 correspond to average values that are summarized in Table 3. The average values of S R for meningiomas and acoustic schwannomas, 0.64 ± 0.16 and gliomas, 0.57 ± 0.23, do not differ significantly.
On the other hand, average values of S C are 0.61 ± 0.17 for meningiomas and acoustic schwannomas, and 0.31 ± 0.17 for gliomas, have a difference that clearly discriminates between these groups. In order to enhance the difference among tumor groups, the ratio S C /S R is used, as seen from Table 3. Figure 9 shows frequency distributions for meningiomas and acoustic schwannomas, with an average value 〈S C /S R 〉 = 0.94 ± 0.11, and gliomas with 〈S C /S R 〉 = 0.58 ± 0.23. Since regularity measures are related to the fractal properties of the tumor interface (Pérez-Beteta et al., 2018) there must be a certain correlation to the scaling parameters d F and α loc . For tumors with S C /S R close to 1, fractal dimension, d F , should be close to 2 and α loc should have its lowest value, i.e., tumor surface is regular and smooth. As S C /S R decreases it is expected an increase in the scaling parameters. Figure 10 shows this trend for the dependence of d F and α loc on S C /S R . The slopes that characterize the linear trend are summarized in Table 4.

Ordered Series, Visibility Graphs and Multifractal Analysis Results
Ordered series were extracted from the tumor interface for those slices that contain the maximum number of interface points, visibility graphs were generated and the degree distribution functions were obtained. Results for the average of P (k) distributions obtained for each of the different tumor types are shown in Figure 11A. For values of the connectivity index k > 20, P (k) decays abruptly due to the fact that the ordered series has a finite size which limits the probability P (k). For values of the connectivity index k < 20 all distributions exhibit a power law behavior as shown in Figure 12A with slopes that are summarized in Table 5 and shown in Figure 12B. If the tumor types are discriminated only into two classes: one including meningiomas and acoustic schwannomas and the other including high grade gliomas and glioblastoma multiforme, the slope distributions are clearly distinct for each tumor class as shown in Figure 11B. This result suggest the possibility of using γ as a possible parameter that characterizes tumor interface dynamics supported by the result shown in Figure 12B. Multifractal analysis results are summarized in Figure 13. As expected, the evaluation of generalized fractal dimensions for the one dimensional sampling of the tumor interface, Figure 13A, does not provide with enough information to discriminate between these two classes as seen in Figure 13B and summarized in Table 6. On the other hand, FIGURE 14 | Relevant morphological and scaling parameters used to discriminate between meningiomas and acoustic schwannomas (green) and gliomas (red). Average values are represented as large dots.
Frontiers in Physiology | www.frontiersin.org June 2022 | Volume 13 | Article 878391 two dimensional detrended fluctuation analysis of the tumor interface yield some differences in the average values of D (q), as shown in Table 6 and Figure 13D, that possibly could be improved by an adequate sampling of the tumor interface, i.e., higher image longitudinal and transverse resolution. In any case, generalized fractal dimensions do not discriminate appropriately between tumor types.

CONCLUSION
A method based on dynamic quantum clustering is used to perform contrast enhanced MRI of brain tumors. Tumor interfaces can be classified according to scaling analysis parameters such as the fractal dimension, d F and the local roughness exponent, α loc which clearly differentiate between the growth dynamics of different tumor types adding support of a ballistic growth model for gliomas and glioblastomas, following the Family-Vicsek ansatz and a non-ballistic growth model for other neoplasias such as meningiomas and acoustic schwannomas. Among the regularity measures, the ratio S C /S R exhibit some correlation with the scaling parameters and clearly discriminates between gliomas and meningiomas or acoustic schwannomas. The relation between d F , α loc and S C /S R is shown in Figure 14. Parameters obtained in series extracted from the tumor interface are size sensitive but nevertheless exhibit differences that could be used for tumor classification and in particular its growth dynamics, through the exponent γ.
Generalized fractal dimensions obtained by two dimensional detrended fluctuation analysis could possibly give significant differences if the tumor interface could be sampled with higher resolution. Further research should take into account combination of different MRI modalities.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.