Small-Angle Scattering and Multifractal Analysis of DNA Sequences

The arrangement of A, C, G and T nucleotides in large DNA sequences of many prokaryotic and eukaryotic cells exhibit long-range correlations with fractal properties. Chaos game representation (CGR) of such DNA sequences, followed by a multifractal analysis, is a useful way to analyze the corresponding scaling properties. This approach provides a powerful visualization method to characterize their spatial inhomogeneity, and allows discrimination between mono- and multifractal distributions. However, in some cases, two different arbitrary point distributions, may generate indistinguishable multifractal spectra. By using a new model based on multiplicative deterministic cascades, here it is shown that small-angle scattering (SAS) formalism can be used to address such issue, and to extract additional structural information. It is shown that the box-counting dimension given by multifractal spectra can be recovered from the scattering exponent of SAS intensity in the fractal region. This approach is illustrated for point distributions of CGR data corresponding to Escherichia coli, Phospholamban and Mouse mitochondrial DNA, and it is shown that for the latter two cases, SAS allows extraction of the fractal iteration number and the scaling factor corresponding to “ACGT” square, or to recover the number of bases. The results are compared with a model based on multiplicative deterministic cascades, and respectively with one which takes into account the existence of forbidden sequences in DNA. This allows a classification of the DNA sequences in terms of random and deterministic fractals structures emerging in CGR.


Introduction
Understanding the correlations between DNA structure and its functions is one of the fundamental challenges in modern biology, with important implications in various biological processes, such as in replication or transcription [1]. Basically, eukaryotic DNA has a hierarchical organization in which the primary structure consists of a sequence of four nucleotides, i.e., adenine (A), thymine (T), guanine (G) and cytosine (C), arranged in two complimentary polynucleotide strands in the form of a double helix. In this configuration, each type of nucleotide on one strands, bonds with just one type of nucleotide on the other strand, according to base pairing rules: A with T, and C with G. At a higher level, DNA is tightly packed into nucleosomes (almost two turns of DNA wrapped around histone protein), separated by additional DNA fragments. Furthermore, the nucleosomes are arranged into a chromatin fiber, which forms loops, then chromatin domains, and finally chromosomes [2]. As such, the number of base pairs involved increases from about 10 2 bp (for short regions of DNA double helix), up to about 10 2 Mbp (for chromosomes).
It is argued that the primary structure contributes, through gene positions and transcriptional activity to the organization of chromatin at all scales, although the precise influence remains controversial [1]. Up to some extent, this may be linked to the existence of long-range correlations of nucleotide distributions [3][4][5][6][7]. The nature and origin of such correlations is intensively studied in the literature by performing various types of statistical analysis on DNA sequences [8][9][10][11][12][13]. One of the most accurate method in revealing the existence of power-law correlations with specific scale-invariance properties is the wavelet transform modulus maxima (WTMM) [13]. This method considers analyzing wavelets that make the wavelet transform microscope blind to low-frequency trends, and can reveal the hierarchy that governs the spatial distribution of multifractal measures [1].
More recently, with the advent of bioinformatics, numerical representation of DNA sequences has become an important approach in analyzing long-range correlations in big data sets as generated from high-throughput methods for sequencing. To this aim, a popular method, which is gaining an increased interest is the chaos game representation (CGR), introduced by Jeffrey in Ref. [14]. This is an iterative mapping technique which assigns to each nucleotide unique coordinates in 2D space, and it allows a visual representation of both local and global patterns in DNA sequences, revealing previously unknown structures [14]. One of the most important property, which makes CGR very useful in numerical encoding is that it gives a one-to-one representation, i.e., given a CGR point in the plane, one can trace back to the origin of the sequence, and therefore the original DNA sequence can always be reconstructed [15]. In addition, the 2D point distributions generated by CGR can be further analyzed using different methods, to extract additional structural information.
To this aim, one of the most common methods, used either standalone or in combination with other methods, is multifractal analysis, since it allows distinguishing between coding and non-coding sequences, or it can be used in phylogeny reconstruction and in investigating the clustering of protein structures. In particular, by using a CGR in which a DNA symbolic sequence is mapped onto a singular measure on the attractor of a particular iterated function system (IFS), the multifractal spectrum of the resulting measure is shown to be more sensitive for detecting dependence structures within DNA than the averaged contribution given by redundancy [16]. WTMM applied to multifractal analysis of CGR images show that the scale-invariance range of CGR edge can be extended to three orders of magnitude, and complete singularity spectra can be calculated [17]. By using CGR of protein sequences based on the detailed HP model, together with their multifractal and correlation analyses, a more precise phylogenetic tree of bacteria has been proposed [18]. CGR of randomly linked functional protein sequences used together with recurrent IFS helps to extract some biological functions of these proteins [19]. Multifractal detrended cross-correlation analysis of genome sequences using CGR shows the existence of multifractal nature and power-law correlation behavior between any pair of genome sequences [20]. A version of CGR, modified to visualize heteroplasmic mutations, used together with lacunarity analysis, reveals fractal properties of mitochondrial DNA sequences, and can quantitatively characterize Parkinson's disease [21].
However, as we shall see below, the multifractal spectra of different point distributions can be very similar, and additional analysis is required to discriminate between such arrangements, and which can provide new insights into their structure. To address this issue, a useful approach is to use the formalism underlying analysis of small-angle scattering (SAS) [22], largely used in molecular biology for investigating the structure of complex macromolecules. Physically, when X-rays are used, SAS is based on the electrostatic interaction of the electromagnetic wave with electrons, while in the case of neutrons, it is based on their interaction with the atomic nuclei (nuclear scattering), and in ferromagnetic materials, on interaction of neutron's magnetic moment with that of the atom (magnetic scattering). Theoretically, depending on the behavior of the SAS intensity, i.e., if we have a simple or a generalized power-law decay (i.e., a succession of maxima and minima superimposed on a simple power-law decay), one can use models specific to random (statistically self-similar) [23,24], and respectively to deterministic (exact self-similar) fractals [25].
A first step in this direction has been recently performed in Ref. [26], where SAS technique has been applied to study the structural properties of various fractal patterns generated by CGR, and an analytic expression of the corresponding SAS intensity has been derived for Sierpiński triangles. The SAS intensity from CGR has been calculated numerically by using a simplified version of the Debye formula [27]. The good agreement between numerical data and the theoretical model indicates that SAS technique can be successfully used to obtain the main structural characteristics, including the overall fractal size, scaling factor, fractal dimension, or the number of units composing the fractal, corresponding to CGRs of DNA sequences, and for which analytic expressions of intensity are not available or are hard to be derived. Such information may be useful to distinguish functional regions of DNA sequences or solving issues related to the classification of organisms. In particular, knowledge of the fractal dimension may provide information about the percentage content of coding regions, i.e., the higher the fractal dimension, the higher the coding percent [28]. In addition, the scaling factor may provide information about the relative sizes of coding regions, and which, in turn, describes the degree of self-similarity for different organisms, since it is related to the value of the fractal dimension [26].
Here, the results obtained in Ref. [26] are extended by performing a combined analysis of CGR using SAS and multifractal analysis of Escherichia coli, Phospholamban and Mouse mitochondrial DNA experimental data. The complementarity of the two methods is illustrated first on a model based on multiplicative deterministic cascades. A model in which both approaches can be hardly used to distinguish between various patterns is also provided, and is based on the existence of forbidden sequences in DNA. Furthermore, the two models are used as benchmarks against which the structures of the experimental data are compared in terms of the interplay between mass and surface-like fractals structures emerging in CGR. As such, it is shown that for Escherichia coli the structure resemble closer a random fractal, while in the case of Phospholamban and Mouse mitochondrial DNA, the corresponding structure is closer to deterministic ones.

Iterated Function Systems and Chaos Game Representation of DNA Sequences
The theoretical framework provided by iterated function systems (IFS) is very useful for classification and description of fractals. Mathematically, a (hyperbolic) IFS is given by a complete metric space (X, d) together with a finite set of contraction mappings w n : X → X, with respective contractivity factors s n , n = 1, 2, · · · , N [29]. Using a shorthand notation, an IFS is {X; w n , n = 1, 2, · · · , N} and s = max {s n , n = 1, 2, · · · , N}. Please note that a transformation f : X → X on a metric space (X, d) is a contraction mapping if there is a constant (contractivity factor) 0 ≤ s < 1 such that d ( f (x), f (y)) ≤ s · d(x, y) ∀x, y ∈ X.
If one considers an IFS with contractivity factor s, and (H(X), h(d)) the space of nonempty compact subsets with the Hausdorff metric h(d), then the transformation W : H(X) → H(X) defined by is a contraction mapping on the complete metric space (H(X), h(d)) with contractivity factor s [29], i.e., h(W(B), W(C)) ≤ s · h(B, C) ∀B, C ∈ H(X).
Its unique fixed point, A ∈ H(X) obeys A = ∪ N n=1 w n (A), is given by A = lim m→∞ W •m (B) for any B ∈ H(X), and is called the attractor of the IFS [29].
Here, for rendering pictures of attractors we play chaos game using the random iteration algorithm. As such, one start by assigning the probability p n > 0 to w n for n = 1, 2, · · · , N, where ∑ N n=1 p n = 1. Then, a point x 0 ∈ X is chosen, and then recursively, next points are obtained according to: The probability of the event x k = w n (x k−1 ) is p n , and k = 1, 2, · · · . This generates the sequence {x k : k = 0, 1, · · · } ⊂ X which converges to the attractor of IFS.
When the chaos game is played with 4 points (i.e., N = 4), the matrix representation of the IFS of affine maps is: where the coefficients a i , b i , c i , d i , e i , f i with i = 1, 2, 3, 4. For a uniformly filled square, they are given in Table 1.
To display visually the underlying structure of the attractor resulting from a string of four letters, Jeffrey proposed [14] to construct a square with vertices "A", "C", "G" and "T", and to control the chaos game, not with the probabilities p i (see Table 1), but with the DNA sequence. In particular, to choose the next point, one could use the next base in the DNA sequence. For a sequence of three letters "CGT", the algorithm translates as follows: first, plot "C" at half distance between the center of the ACGT square and "C" vertex, then, the next base "G" is plotted at half distance between the previous point and "G" vertex, and finally, the last bases "T" is plotted at half distance between the previous point and the "T" vertex.

Fractals and Multifractals
Mathematically, description of fractals is based on concepts from measures theory [30], and a rigorous definition has been introduced by Hausdorff in Ref. [31]. This involves a subset S of the n-dimensional Euclidean space, and let us consider that {C i } is a cover of S with c i = diam (C i ) ≤ s and s ∈ S. Then, the Hausdorff measure m α (S) is given by taking the infimum over all possible coverings, i.e., and thus, the fractal dimension D can be written as: This corresponds to the value of α for which the Hausdorff measure changes from zero to infinity. When α = D, m α (S) can have arbitrary values within this range.
However, for most practical purposes, it is very difficult to use Equation (6) for determination of the fractal dimension. To avoid this issue, another approach is to determine the variation of the fractal measure M, inside a sphere of dimension n and radius r centered on the fractal. For fractal systems, one can write a mass-radius relation of the type [32]: where lim r→∞ log A (r) / log r → 0. The above equation plays a central role in development of the concepts used thereafter, since it allows us to describe mass and surface fractals, and provide analytical expression for fractal dimensions for a large class of fractals.
Thus, let us consider further a fractal of size L composed of balls of size a. Then, the number of balls enclosed by the imaginary sphere of radius r with a ball in the center, is given by [32]: with l r L. If the fractal is a line, one has D = 1, for a smooth surface, D = 2, while for a regular Euclidean 3D object, D = 3. For fractals with a scaling factor β s , one can use the property that at first iteration the fractal consists of k copies of itself, each of size β s L, and write that [32]: Then, by using Equation (7), one obtains: which can be used to obtain the fractal dimension D. For fractals with multiple scaling factors β si and k i copies with i = 1, · · · , n, Equation (9) is rewritten as: Similarly, for surface fractals one can write: where S(r) represents the area between the boundary of the (rough) surface and the envelope of all spheres of radius r centered on the boundary. Here, S 0 is a constant, which is the surface area itself for a smooth surface, i.e., when D s = 2. The fractal dimensions D which appear in the above relations, are equivalent to the box-counting dimension and describe fractals with a single scaling factor. However, when several scaling factors are present, a more detailed description is required, and can be achieved by using the multifractal formalism [33,34], where one considers an object S covered by a grid of boxes B i (l) of size l. By considering that the measure determined by the probability of hitting the object in the box B i is µ(B), the number of covered boxes N at resolution l is N ∝ 1/l 2 , and thus one can write [35]: where Z is called the "partition function" and has a power-law behavior when l → 0 and N → ∞, so that Z s ∝ l D s (s−1) . Here, i denotes each individual box, and are the hitting probabilities. Then, the generalized dimension spectrum is given by [35]: By considering the ratio p i ≡ N i (l)/N, which gives the relative weight of the i-th box, one can write: The generalized dimension spectrum is a monotonically decreasing function, with horizontal asymptotes at α max = lim q→−∞ D s and α min = lim q→∞ D s . Their values can be used to describe the heterogeneity, i.e., if α max = α min the fractal is heterogeneous (multifractal), and homogeneous otherwise.
An equivalent representation of generalized dimension spectrum is provided by f (α) spectrum, which gives a mathematically precise and intuitive description of the multifractal measure in terms of interwoven sets, with singularity strength α, and Hausdorff dimension f (α). Without going into the details of derivation here, it can be shown that the generalized dimension and f (α) spectra are related through a Legendre transform [36], i.e., where α(s) = dτ(s)/ds and τ(s) = (s − 1)D s . In the remaining of the paper, Equations (15) and (16) are used to characterize the multifractal properties of point distributions generated by CGR described in the previous section. To illustrate the method, analytically solvable models given by multiplicative deterministic cascades are first analyzed in Section 3.1.1.

Small-Angle Scattering
One considers here a two-phase approximation in which microscopic scattering objects with scattering length b j have the scattering length density (SLD) ρ s (r) = ∑ j b j δ(r − c j ), where c j are the position vectors of the objects. The differential elastic cross section is defined by dσ/dΩ = |A t (q)| 2 , where q is the scattering vector. For a three-dimensional object A t (q) = V ρ s (r) exp(iq · r)dr is the total scattering amplitude and V is the irradiated volume. When the objects are embedded in a solid matrix of SLD ρ 0 , then the scattering contrast is defined by ∆ρ = ρ − ρ 0 , and the total scattering intensity can be written as [22]: where c is the concentration of objects, V is their volume, and F(q) ≡ (1/V) V exp(−iq · r)dr is the normalized form factor, with F(0) = 1, and the symbol · · · denotes the ensemble averaging over all orientations. In what follows, since the object are two-dimensional ACGT squares, the volume V in Equation (17) shall be replaced by the corresponding surface area. Please note that the above averaging procedure allows the rotation of the squares in three-dimensional space, with equal probability. Since the positions of the points are generated by CGR, we can start with the Debye formula to calculate the scattering intensity, and thus Equation (17) can be written as [27]: where I s (q) = 1 is the intensity scattered by each point, q is the magnitude of the scattering vector q, and r ij is the distance between arbitrary points i and j. When the number of points exceeds a few thousand, the computation of the term sin(qr ij )/(qr ij ) is very time consuming, and thus it is handled via a pair-distance histogram g(r), with a bin-width commensurate with the experimental resolution [37]. Therefore Equation (18) becomes: where M bins is the number of bins, and g(r i ) is the pair-distance histogram at pair-distance r i . The latter quantity is calculated from the positions of points inside the ACGT square.

Analysis of Theoretical Models
To illustrate the complementarity between SAS and multifractal techniques in analyzing point distributions, one considers a model of multiplicative deterministic cascades. This is followed by a second model, based on forbidden sequences in DNA, and which shows that both the SAS and multifractal spectra can hardly be used to differentiate between various structures.

Multiplicative Deterministic Cascades
The multiplicative deterministic cascade model is constructed by dividing a square into for equal squares. To each of the subsquare, one assigns the probabilities p i ∈ [0, 1], with i = 1, 2, 3, 4. This is called the first iteration, and it is denoted n = 1 (Figure 1a). Then, at second iteration (n = 2), each of the four subsquares is further divided in four squares, and to each of them are assigned probabilities given by the Kronecker product of the matrices with the elements p i placed as in Figure 1a. The results are shown in Figure 1b. At third iteration, one perform a similar division of squares, and to each of them one assigns the probabilities given by the Kronecker product of the matrices with as elements the probabilities from iteration n = 1, and respectively from n = 2. The multiplicative cascade model is obtained in the limit of high number of iteration, and the distribution of square values therefore depends on the initial choices of probabilities p i . This model is similar to the model for the displacement of a viscous fluid by a nonviscous fluid in a porous medium, developed in Ref. [38]. However, in this reference, at n = 2 the probabilities p i associated with each such division are multiplied in random order by probabilities p i . The same random assignment of probabilities is kept for subsequent iterations, thus leading to a multiplicative random cascade model.   For these models of multiplicative cascades, it can be shown that the generalized dimension spectrum has an analytic expression given by [39]: where f j = p j / ∑ 4 i=1 p i , and j = 1, 2, 3, 4. The corresponding dimension spectra are presented in Figure 3a. The results clearly show that the spectra of models M1, M2, M3, M5 and M6 have a decreasing behavior, indicating that the corresponding structures are multifractals. However, for model M4 the spectrum is constant and indicates a simple fractal structure (Sierpiński gasket) with fractal dimension about 1.58, as expected. Also, in the case of model M7 (uniform square) the fractal dimension 2 is recovered. Please note that for models M5 and M6 the dimension spectra coincide for all s. They are also the most heterogeneous structures, since the difference α min − α max 3.45 − 1.48 = 1.97 is maximum, as compared to the other models.
The f (α) spectra are calculated by using Equation (16), and the results are presented in Figure 3b. Generally, these spectra are convex functions with a single maximum at α = α 0 . This gives the most frequent value of the scaling indices α, and appears at s = 0, where f (α) = D 0 [39]. Also, the scaling indices take vales in the finite range [α min , α max ]. The smaller the length of this interval, the more homogeneous the structure. In particular, for models M4 and M7, the spectra degenerate into a single point at f (α) 1.58 and respectively f (α) = 2. Please note that for models M1, M3, M5 and M6, the f (α) spectra are asymmetric with respect to their maxima. Lower values of f (α) for s < 0 (models M1, M5 and M6) indicate a higher influence of the lowest values of the fractal measure in the spectral complexity. Similarly, the appearance of lower values of f (α) for q > 0 (model M3) indicate a higher influence of the lowest values of the fractal measure [40].  Table 1 for the corresponding probabilities p 1 , p 2 , p 3 and p 4 .
In the following one obtains the SAS from models M5 and M6. To this aim, the coordinates of the points needs to be extracted from Figure 2 and used as input in Debye equation (Equation (19)). However, a general property of such images is that the distribution of grey levels is concentrated within a short range, and thus image binarization at various thresholds becomes impracticable. Figure 4 shows this situation for models M5 and M6, where the number of pixels are mostly concentrated in the 0.7 ÷ 0.95 range. Despite this local pixel concentration, their distribution is different, and this indicate that SAS technique can distinguish between the two structures.  Therefore, an equalization of histograms such that they span the entire range from 0 (black) to 1 (white) needs first to be performed. This is done in Figure 5, which shows both the transformed image and the corresponding pixel distribution (upper part-model M5, and lower part-model M6).
The histogram distribution of the models show a periodicity specific to multiplicative deterministic cascade structures, since the grey levels alternate in a regular fashion. The heights of the single maxima (for model M5) and of the groups of maxima (for model M6) follow a parabola-like behavior over the entire range of pixel values. These features now allows the analysis of the transformed images at various thresholds t. Figure 6 shows the corresponding binarized images for models M5 (upper part) and M6 (lower part) at thresholds 0.2, 0.4, 0.6 and 0.8 (from left to right).
The structure factor of the point distributions in Figure 6 is calculated by using Equation (19) and the results are shown in Figure 7. The main feature of the scattering curves in all cases is the presence of three structural regimes: S(q) ∝ q 0 at ql 2π (Guinier region), S(q) ∝ q −D , with D 2 at 2π ql 2πl 0 (fractal region), and S(q) 1/k at 2πl 0 ql (asymptotic region). Here, l is the overall fractal size (in pixels), k is the number of points composing the fractal, and l 0 is the pixel size. The length of the Guinier region provides information about the overall fractal size and they are very similar at each threshold, since the dimensions of squares are the same. In the fractal region one can see a generalized power-law decay (GPLD), i.e., the presence of a succession of maxima and minima superimposed on a simple power-law decay. A GPLD is characteristic to deterministic fractals [25] and can be used, as we shall see below, in extracting the fractal iteration number and the scaling factor.
The differences between scattering curves for models M5 and M6 are clearly visible at smaller values of the thresholds, and less visible at higher values, in particular, at t = 0.8 the differences are practically indistinguishable. This is because by increasing t the distribution of points becomes more uniform (see last column in Figure 6), and the deterministic character cease to dominate. Consequently, the maxima and minima are smoothed to the extent that the resulting curves coincide. Also, the scattering exponent D 2 is kept in all cases ant it coincides with the numerical values of D 0 given by multifractal spectra (see Figure 3). Finally, in the asymptotic region one can see that the total number of pixels k increases with t (since the ratio 1/k decreases), and this can be used to obtain the total number of points present in the fractal.

Missing Sequences Models
Once a complete genome is available, a natural question which arise is whether there are subsequences of short strings of letters "A", "C", "G" or "T" that are missing, since this may reveal some biological meanings such as evolutionary relatedness of species [41]. To this aim, important steps have been performed in Ref. [42], where sequences are transformed into portraits, and thus the missing sequences can be visualized. Although this method has also the advantage of providing analytic expressions of the fractal dimensions of the patterns emerging from this visualization scheme [43], here CGR shall be used since in the case of missing of certain types of string composed on 2-contiguous letters, the resulting pattern has an obvious fractal structure (see Figure 8), while the fractal dimension can be extracted from either SAS or multifractal spectra, as discussed in the previous section. Moreover, by establishing relations between DNA sequences with missing subsequences, and the generalized Cantor sets, such as those presented in Ref. [44], the possibility of using SAS technique in analyzing structures with missing subsequences can be greatly extended, since theoretical expressions for SAS intensities are already available for some classes of generalized Cantor fractals [45].
In the following, one considers a random sequence of length 8 × 10 4 consisting from the four letters "A", "C", "G" or "T", but that has never "A" followed by "T" (model TR12), "A" followed by "G" (model TR13), or "G" followed by "G" (model TR33). The corresponding structures obtained from CGR are shown in Figure 8. The resulting patterns has fractal properties in which the complementary structure, i.e., white regions, have various geometrical shapes (model dependent) and with sizes following a power-law distribution. The simplest complementary structure is shown in Figure 8 (middle part), and which consists from one square of edge length 1/8, 4 1 squares of edge lengths 1/16, 4 2 squares of edge lengths 1/32 and so on. Such a distribution is specific to surface-like fractals and has important implications in behavior of the corresponding scattering curves. In particular, when the associated measure is multifractal, then the SAS intensity shows a not only a single mass fractal region, as in Figure 7, but a succession of mass fractal regime followed by a surface fractal one [46].
Since for all three models, both the number of complementary regions forming the surface fractal, as well as their scaling factor are unchanged, one may expect that the corresponding SAS curves from these models to be quite similar. This is confirmed numerically in Figure 9a, which shows that the SAS curves decay ∝ q −D , with D 2 in the fractal regime, and they are hardly distinguishable. For comparison, the SAS intensity of a square, i.e., D = 2 is also included (black curve). The absence of maxima and minima in the fractal regime reflect a nearly uniform arrangement of the points, resulting from the random process used in generating their positions. The f (α)-spectra for these models are shown in Figure 9b. They are confined in a narrow range of α values, i.e., 1.62 α 1.95 for models TR12 and TR13, and 1.72 α 1.95 for model TR33, thus indicating that the multifractal nature is not so pronounced as in the case of multiplicative cascade models (see also Figure 3b). Also, for comparison, Figure 9b includes the f (α)-spectrum of a uniform distribution of points in the "ACGT" square, which is an Euclidean object. This is represented as a single point (black) of coordinates (2,2).
The results in Figure 9 illustrate that generally, SAS and f (α) can hardly be used to perform a clear differentiation between such structures. This issue may be eventually addressed by establishing first the relationship with the corresponding generalized Cantor set, as performed in Ref. [44], followed then by calculating the corresponding SAS intensities.   9. Scattering intensities (a) and f (α) spectra (b) for models TR12, TR13 and TR33. l is the overall size of the "ACGT" square. The data for a uniform square are included for comparison.
The protein encoded by Phospholamban is found as a pentamer, is a major substrate for the cAMP-dependent protein kinase in cardiac muscle, and is an inhibitor of cardiac muscle sarcoplasmic reticulum Ca(2 + )-ATPase in the unphosphorylated state [47]. However, inhibition is relieved upon phosphorylation of the protein, and subsequent activation of the Ca(2 + ) pump leads to enhanced muscle relaxation rates, thereby contributing to the inotropic response elicited in heart by beta-agonists. The encoded protein is a key regulator of cardiac diastolic function. Mutations in this gene are a cause of inherited human dilated cardiomyopathy with refractory congestive heart failure, and familial hypertrophic cardiomyopathy [47].
In the case of Mouse mitochondrion, the genome displays exceptional economy of organization, with tRNA genes interspersed between rRNA and protein-coding genes with zero or few non-coding nucleotides between coding sequences [48]. The genome is highly homologous in overall sequence as well as the organization to human mitochondrial DNA, and an important feature in that the translational start codon is AUN, with any of the four nucleotides in the third position, whereas the only translational stop codon is the orthodox UAA [48].
O145:H28 is one of the major non-O157 Shiga toxin (Stx)-producing Escherichia coli (STEC) lineages that causes severe diseases and is frequently isolated from humans, animals and foods [49]. Highly dynamic features of prophages and plasmids in the diversification of the STEC lineage, and the differential dynamics and impacts of these mobile genetic elements on the pangenome and virulence factor repertoire among O145:H28 strains [50]. The CGR of Phospholmban and Mouse mitochondrion in Figure 10 resemble more closer a deterministic-like fractal structure, in which a pattern repeats itself at increasingly smaller scales. In particular, the structure of Phospholamban resemble, in a first approximation, the structure of model TR12 shown in Figure 8.
Note that in "CG" sub-quadrant of CGR from Phospholamban, i.e., in the upper-left region of the G-quadrant in Figure 10, the number of points are relatively much smaller as compared to other sub-quadrants. This indicates a paucity of subsequences ending in "CG" [14]. Similarly, the high number of points in "TT", "AT", "TA" and "AA" sub-quadrants, indicates a higher number of subsequences ending with these suffixes.
However, for Mouse mitochondrion, the white regions specific to model TR12, are partially populated with points, while their overall density is higher at the bottom. The paucity of subsequences ending in "CG" is still preserved, but however, the "AA" sub-sequence seems to be predominant. As it shall be seen below, this has important consequences on the behavior of SAS and f (α)-spectra. A somehow similar arrangement can be seen also for Mouse mitochondrion, where a higher density of points is formed within the "ATC" triangle. The structure of Escherichia coli is characterized by the coexistence of regions with different densities, distributed more uniformly, as compared to Phospholamban and Mouse mitochondrion. Figure 11a shows the SAS intensities from Phospholamban, Mouse mitochondrion and Escherichia coli calculated with Equation (19). For comparison, the scattering curve of a uniformly filled square of the same size as the "ACGT" square is included. The results show the presence of the fractal region at π ql π/r min , where l is the overall size of the square (360 × 360 pixels 2 ), and r min is related to the smallest size of the pixels in the image. The common feature is the power-law decay with scattering exponent D 2, revealing a mass fractal structure in all three cases. However, for Phospholamban and Mouse mitochondrion, the scattering curve has a generalized power-law decay (see Section 1). To investigate this behavior in more details, Figure 11b shows the curves q D I(q)/I(0), which reveal more clearly the periodicity of Phospholamban and Mouse mitochondrion in the fractal regime. This is a signature of their deterministic nature, as indicated also by their structures shown in Figure 10a,b.
Then, the periodicity and the number of minima can be related to the scaling factor β and to the number of iterations [25]. In particular, the curve q D I(q)/I(0) is approximately log-periodic with the period equal to the inverse scaling factor, and thus β 1/2, while the number of iterations is equal to three. Please note that the value of the scaling factor is related to the CGR algorithm, and thus if one plot a given point, not at half distance between the previous point and the corresponding vertex, but at another fraction, or even using a 3D CGR as described in Ref. [51], one expects to recover also the corresponding scaling factor. Therefore, a SAS analysis could be employed to check the nature of the algorithm used to generate the points distribution. Symmetry properties occurring in the "ACGT" square can also be revealed by performing a structural analysis in real space. Thus, by taking the Fourier transform of Equation (19), one obtains the corresponding pair-distance distribution function (pddf) p(r), which gives the probability density of finding the distance r between two arbitrary points. For fractal structures, pddf can be used to reveal the fractal scaling factors, since it is characterized by a succession of groups of maxima and minima distributed periodically on a double logarithmic scale [45,52]. Figure 12a shows the pddf of Phospholamban, Mouse mitochondrion and Escherichia coli. For comparison, the pddf of a uniform distribution of points is included. The common feature is that all curves follow the same trend as that given by the uniform distribution, thus reflecting the symmetry properties of the square. In addition, the most common distances occur at r/l 0.5, i.e., at half the edge length of "ACGT" square, while the maximum distance between two arbitrary points occurs at r/l 1.4. However, the pddf of Phospholamban and Mouse mitochondrion are characterized by small maxima and minima in the form of wriggles, in the range 0.15 r/r 0.8. This is the result of their deterministic nature and correspond to the generalized power-law decays occurring in Figure 11a,b in the fractal region. Finally, the pddf of Escherichia coli shows slight variations in the range 0.2 r/l 1. In particular, for 0.2 r/l 0.7 the pddf has smaller values as compared with the pddf of the uniform distribution of points, thus indicating the presence of more rarefied regions. However, for 0.7 r/l 1, the pddf is slightly higher as compared to the pddf of the uniform distribution, thus indicating the presence of more dense regions (see also Figure 10c).
The non-uniform points density distribution in the CGR for these genomes is also reflected in Figure 12b, which shows that all structures are characterized by a multifractal structure. In particular, the Phospholamban appears to have the most heterogeneous distributions, reflected by the highest range of α values spanned by f (α) curve. Please note that the maxima of all spectra occurs at f (α) 2, in agreement with the values of fractal dimensions obtained from the scattering exponents in the fractal regime, shown in Figure 11a.

Conclusions
In this work SAS and multifractal techniques are employed to perform a structural analysis of CGR of DNA sequences. To address the complementarity of the two techniques, they are applied first to two theoretical models, i.e., to a deterministic multiplicative cascades and to missing sequences, and then to genome data of Phospholamban, Mouse mitochondrion and Escherichia coli.
The deterministic multiplicative cascades model may serve as a benchmark against which other structures can be compared, since it provides an analytical expression for the generalized dimension spectra. It has been shown that within this model, multifractal analysis may not always differentiate between the generalized dimension, and consequently between f (α) spectra (Figure 3), but however, when an appropriate threshold is used, SAS may distinguish between such structures. The differentiation is based on the behavior of the scattering curves in the fractal regimes and on the specific way, i.e., amplitude, periodicity and number of most pronounced maxima and minima, change with the magnitude of the scattering vector ( Figure 7). The missing sequences model addresses the question of whether there are subsequences of short strings of letters ("A", "C", "G" and "T") that are missing. Here, a detailed analysis is performed on structures in which "A" is never followed by "T" or "G", and when "G" is never followed by "G". It is shown that SAS and multifractal techniques can hardly differentiate between such structures, since their corresponding spectra coincide (Figure 9). This may be attributed to the interplay between mass and surface-like fractals appearing in the "ACGT" square. While the mass fractals consist from the points resulted from CGR, surface-like fractals are formed by the power-law distribution of unoccupied regions (Figure 8).
Analysis of DNA sequences of Phospholamban, Mouse mitochondrion and Escherichia coli reveals that both SAS and multifractal techniques can be used to distinguish between structures (Figures 11 and 12), as well as to extract additional structural information. SAS shows that Phospholamban and Mouse mitochondrion DNA generate a more ordered arrangement in the CGR, as compared to Escherichia coli. Figure 10 shows the corresponding arrangements. In the former case, these resemble closer deterministic structures, and in addition to the fractal dimension, it allows us to obtain the scaling factor (from the log-periodicity of minima in the curve q D I(q)/I(0)) and the fractal iteration number (from the number of these minima; Figure 9b). Escherichia coli generates a more uniform distribution but still with some density fluctuations (Figure 10c). For such a structure, SAS is characterized by the absence of clear periodicity of maxima and minima in the fractal regions (Figure 11b). The heterogeneity of such structures is clearly revealed also in the f (α) spectra (Figure 12b), and which show that Phospholamban and Mouse mitochondrion have the highest structural heterogeneity.
The analysis performed in this work provides good insights on the discriminating power of both SAS and multifractal formalism, and gives useful structural information which may help to distinguish between coding and non-coding sequences in DNA, phylogeny reconstruction or clustering of protein structures. In particular, the fractal dimension D of coding segments has smaller values as compared to non-coding segments [53]. This can be used to quantify the excess of short runs and the deficit of long runs of weak and of strong hydrogen bases in coding sequences. The dimensions D −2 , D −1 and D 1 from the multifractal spectra (Equation (15)) can be used to perform a classification of bacteria by assigning to each sequence a point in 2D and 3D spaces (D −1 , D 1 ), and respectively (D 1 , D 1 , D −2 ). In this representation, bacteria that are close phylogenetically, are almost close in these spaces [54].