Controlling the Multifractal Generating Measures of Complex Networks

Mathematical modelling of real complex networks aims to characterize their architecture and decipher their underlying principles. Self-repeating patterns and multifractality exist in many real-world complex systems such as brain, genetic, geoscience, and social networks. To better comprehend the multifractal behavior in the real networks, we propose the weighted multifractal graph model to characterize the spatiotemporal complexity and heterogeneity encoded in the interaction weights. We provide analytical tools to verify the multifractal properties of the proposed model. By varying the parameters in the initial unit square, the model can reproduce a diverse range of multifractal spectrums with different degrees of symmetry, locations, support and shapes. We estimate and investigate the weighted multifractal graph model corresponding to two real-world complex systems, namely (i) the chromosome interactions of yeast cells in quiescence and in exponential growth, and (ii) the brain networks of cognitively healthy people and patients exhibiting late mild cognitive impairment leading to Alzheimer disease. The analysis of recovered models show that the proposed random graph model provides a novel way to understand the self-similar structure of complex networks and to discriminate different network structures. Additionally, by mapping real complex networks onto multifractal generating measures, it allows us to develop new network design and control strategies, such as the minimal control of multifractal measures of real systems under different functioning conditions or states.


Results
Weighted multifractal graph model. In order to generate weighted networks with multifractal characteristics and to minimally control complex networks, we introduce the weighted multifractal graph (WMG) model to generate and capture random weighted graphs with multifractal topology. By choosing a few parameters, the model can generate a wide variety of weighted multifractal topologies with arbitrary statistical properties, such as generalized degree distribution and clustering coefficient. This multifractal model can be used to fit diverse realworld datasets in fields such as biological systems, geoscience and financial markets. Also, it enables to control multifractality of networks with mild modifications and adjustments.
Along the same lines as in the multifractal network generator 32 , we generate the WMG model recursively from simple square structures. First, we define simple generating measures on R-layer unit squares. The x-axis and y-axis of each layer are identically divided into M intervals, and each interval length is denoted as l i , i = 1…M. The unit square is divided into M 2 rectangles, each rectangle (i, j) is assigned with a set of probabilities p ij (r), i, j = 1…M, r = 1…R. The WMG model is developed from this initial generating measure. The self-similarity appears when we divide each rectangle in the unit square into M 2 sub-rectangles, which has the same structure as the unit one. The self-similar WMG model is formed after repeating the operation K times. Meanwhile, the assigned probabilities p r ( ) When generating a random multifractal graph, we first select a discrete weight set {w(r)} r=1:R . The r and w(r) are the index of the discrete weight level and the corresponding weight value, respectively. p r ( ) can be any desired positive and real values according with specific distributions for specific (desired) network features. Here R is the www.nature.com/scientificreports www.nature.com/scientificreports/ resolution level of the discrete weights. Choosing larger R leads to higher graph resolution and more parameters will be needed to be introduced. Like De Wijs model 24,25 , Kronecker random graph [29][30][31] , stochastic block model 34 and multifractal network generator 32 , a weighted random graph with N nodes is generated from the linking probability measures using the following strategy: (1) Generate a series of uniformly distributed random variables U(n) for each node n = 1…N independently. (2) Find the category i(n) of each node n such that . In the recursive model with K iterations, this suggest that the node n falls into sub-interval i(n) on the unit edge. Since l i and l i K ( ) are identical for all layers, the categories i(n) assigned for the node n in every weight layer are also identical. (3) For each pair of nodes (u, v), generate an uniformly distributed random variable L(u, v), and find index r such that it satisfies ∑ < ≤∑ . Next, we assign the weight value w(r) to the link between nodes u and v. If there is no r satisfying this condition, then there will be no edge between this pair of nodes (u, v). It is noteworthy that when M K ≈ N and l { } i i M 1: = is equally sized, if each node n corresponds with one category i(n) without repetition, this special case of WMG model would reduce to Kronecker graph [29][30][31] . The iterative WMG model building procedure is illustrated in Fig. 1(a). The presented model example has 3 layers. In this example, l { } [0 027 0 063 0 063 0 147 0 063 0 147 0 147 0 343] . . Two networks in Fig. 1(b,c) are generated under the same model parameters M = 2, R = 2, 1 2] . The color of nodes represents generating category and the stroke width of edges displays the weight. We highlight one node in each graph, its neighbours and the edges connecting them. As shown in Fig. 1 , the edges connecting them have weight w(2) = 2. By contrast, in Fig. 1(c), the categories of nodes are more evenly distributed because the interval lengths are = . for i = 1. . . 8. The highlighted node in Fig. 1(c) has category i = 1. In this example, for i = 1. . . 8. Therefore, in this generated graph, all the links containing the highlighted node have weight w(1) = 1. Fig. 1(d,e) show two networks generated by different weight distributions w(r) = r and = + w r r ( ) log (1 ). While all other parameters are kept the same, from Fig. 1(d,e), we can see that a change in w(r) from r to + r log ( 1) leads to networks with different weight features. In Fig. 1(f) we show the normalized histogram of the weights in two networks (d,e). Next, we present the multifractal analysis of this WMG model.

Multifractal analysis.
Traditional approaches for quantifying the multifractal properties of networks rely on box-counting and box-covering method 18,20,37 . The fractal dimension can also be obtained with the help of closeness centrality 38 . However, finding the minimum number of boxes of a given radius required to cover the graph is as hard as the graph coloring problem 20 . Introducing randomness in the graph model brings uncertainty to the box-counting problem. Although there exist few analytical results that determine inequality constraints on the chromatic number for the graph coloring problem in Erdős-Rényi model [39][40][41] , these inequalities cannot be analytically extended to networks with higher-order correlations and so determine their box covering.
To elucidate the multifractal properties of the above-mentioned weighted network construction strategy, we describe next an analytical framework capable of estimating the partition function, the mass exponent, the generalized fractal dimension, the Lipschitz-Hölder exponent and the multifractal spectrum. Since the link generating process contains both node attributes and edge formation characteristics, the proposed mathematical formalism estimates the number of edges which are generated under the same rule. While this formalism bears some similitudes with the multiplicative cascade model [26][27][28] , it counts the number of edges in the generated graph rather than counting sub-blocks with same structure. Our edge-based strategy enables us to analytically compute the partition function and derive the multifractal related metrics (e.g., generalized fractal dimension, multifractal spectrum). For simplicity, let The analytical derivation of the partition function allows us to calculate the mass exponent as τ , the Lipschitz-Hölder exponent (also refer to coarse Hölder exponent or singularity index in prior work) α = τ q ( ) d q dq ( ) and the multifractal spectrum from f(α) = α(q)q − τ(q).
The networks generated by the proposed WMG model inherit a rich variety of multifractality compared with the unweighted MFNG model 32 . For example, when R = 1 (corresponding to MFNG model), the partition function in equation (3) shows that the variation in K does not influence significantly the mass exponent τ(q), the generalized multifractal dimension D(q), the Lipschitz-Hölder exponent α(q) and the multifractal spectrum f(α) because of the log operation. However, as shown in Fig. 2(a), for M = 2, R = 3 and same {p ij (r)} and {l i }, the variation in the model iteration step K causes a shift of multifractal spectrum towards a support on higher Lipschitz-Hölder exponents. One can observe that the shift in the support of the multifractal spectrum is smaller with increasing changes ΔK in the model iteration steps K.
The parabola-like curve shows the multifractality of the WMG model. Fig. 2(b) presents more results of multifractal spectrum curves with four different sets of initial unit square model parameters Of note, for identical parameters (i.e. p ij (r) = (0.25, 0.25, 0.25, 0.25)), the WMG generates a monofractal network (represented in Fig. 2(b) by a blue diamond marker). More precisely, this implies that a more evenly assigned p ij (r) leads to more symmetric and narrower curve of the multifractal spectrum (which corresponds to a point for identical p ij (r)). A shift and reshape of the curve can also be caused by changing l i .
We observe that there is an interesting effect of the sub-interval length {l i } on the multi-fractal spectrum f(α). For example, as presented in Fig. 2(c), when node generating probabilities are identical ( = . . l { } (0 5, 0 5) i -yellow curve), the maximum value of multifractal spectrum f(α) places at leftmost compared with other curves given by .
. -green), the maximum value of the multifractal spectrum remains unchanged but the left and right hand side end points of the spectrum change. Similarly, the influence of {p ij (r)} can also be explored. retains more information of the network topology compared to the straightforward strength distribution, which helps us to better understand the network structure and its multifractality.
Via the extension of the multifractal network generator 32 , in the WMG model, the generalized degree distribution is expressed as Detailed derivation of Eq. (4) is provided in Methods section 'Generalized degree distributions' . Fig. 3(a) illustrates the comparison between the analytical and empirical generalized degree distribution. The comparison analysis considers 100 runs, where networks of size N = 5000 nodes are generated from a set of random parameters {p ij (r)} and {l i }, with M = 2, K = 5, and R = 4. As one can notice, the WMG model can generate a wide variety of generalized degree distributions.
The clustering coefficients measure the local cliquishness and neighborhood connectivity 5 . Concerned with the density of connected triangles, it reflects the extent to which the neighbors of a certain node are also connected. In weighted graphs, the clustering coefficient of a node can be computed as the mean weights of triangles containing the node, divided by the mean weights of triplets (link pairs) originating from that node 42 . Following the derivation in the work of Palla et al. 32 , the average local clustering coefficient of a node falling into interval i can be calculated as c w l l p r p r p r w l l p r p r = + +  are the arithmetic mean of triangles and triplets edge weights. Index i marks the interval which the node falls into when generating the graph. i = 1…M K . Comparison between empirical and analytical clustering coefficients is shown in Fig. 3(b). For most of the intervals, the empirical result fits well with the analytical one. However, smaller clustering coefficient corresponds with sparser local graph structure, thus standard deviation tends to be larger than the ones with comparative larger clustering coefficients.
The proposed WMG model can reproduce the statistical properties of classic theoretical models. We test our model on a weighted graph generated by the Erdős-Rényi model. The weighted version of Erdős-Rényi model is is the corresponding linking probability with the discrete weight w(r). We apply the WMG model to a network with 500 nodes generated by a 4-layer Erdős-Rényi model. Fig. 3(c) illustrates the generalized degree distribution of the simulated network and the recovered WMG model. The blue asterisks denote the target distribution from the simulated ER graph, while the red circles represent the generalized degree distribution obtained from the identified WMG model. It shows that the proposed WMG model can capture the statistical properties of the classic random graph model.

Multifractality in biological systems.
In two recent works 18,43 , the authors show that the interactions among different chromosomal regions and brain networks display a non-trivial multifractal property. We also applied the box-covering method proposed in 18 to show the existence of multifractality in real networks. The box covering algorithm returns the mass of each boxes with radius ϵ to fully cover the network with N nodes. As illustrated in Fig. 4(a,b), the plots of partition functions and the scales reveal the multifractality behavior in the brain networks and the chromosome networks.
However, as discussed in section 'Multifractal Analysis' , finding boxes to cover the graph is as hard (from a computational complexity perspective) as graph coloring problem and analytical expressions cannot be introduced. By mapping the real-world networks to multifractal generating measures that allow to construct weighted graph models, it allows us to characterize the complexity of real networks in function level and also makes it possible to better understand and control the topology structure. In this work, instead of counting boxes and calculating the multifractality in a given graph, we consider the multifractality in the recovered generating measures and its corresponding graph model. Therefore, we exploit the proposed WMG model for investigating two biological systems, namely the chromosome interactions within the yeast genome and brain functionality networks during various stages of Alzheimer's disease. Our goal is not only to characterize the generalized degree distributions in these real biological networks, but also to exemplify how a new network wide control strategy can be employed for ensuring specific mathematical network characteristics.
Yeast genome. Capturing the chromatin interactions via the Hi-C technique allows us to study the three-dimensional structure of the genome and understanding their observed behavior [44][45][46] . For instance, significant topological reorganization of yeast cell in the quiescence state has been observed 47 . Starting from these premises, we consider the yeast chromosome interaction data from Rutledge et al. 47 , where the Hi-C experiments are conducted on yeast cells in exponential growth on glucose-containing medium and on yeast cells in quiescence induced by glucose starvation. The data of Yeast genome are from the library of Juicebox software 44,45 . By interpreting the Hi-C matrix as an adjacency matrix corresponding to a complex weighted network, we estimate the WMG model from the yeast chromosome interaction data for both the exponential growth and quiescence states. We then regulate the multifractal spectrum of yeast cell in quiescence and transit it to the state of exponential growth using the proposed multifractality control strategy (see Methods section on 'Multifractal control'). Fig. 5(a,b) illustrate the generalized degree distribution for the yeast cells in the exponential growth state and quiescence state, respectively. The blue asterisks denote the target distribution from chromosome interaction matrix, while the red circles represent the generalized degree distribution obtained from the identified WMG model. The simulated annealing algorithm can reproduce the model with similar statistical properties (in terms of shape and trend) to the real-world data (see Methods section on 'Reconstructing weighted multifractal model'). Eq. (16) can also be redefined with any other desired network structure properties such as strength distribution and average nearest neighbour strength, but this is left for future work.
As discussed in Methods section 'Multifractal control' , by solving Eq. (15), we modulate quiescent yeast cell to the desired state which yeast cell is growing exponentially. In Fig. 5(c), the red line represents the desired multifractal spectrum from actively growing cell. The curve with blue dots is the multifractal spectrum of quiescent www.nature.com/scientificreports www.nature.com/scientificreports/ yeast cell. The magenta asterisks are the optimal modulated solution given by simulated annealing algorithm. By controlling real-world network multifractality, we could enforce or track cell evolution procedure and therefore regulate bio-feature mechanisms yeast growth.
Alzheimer's disease. Alzheimer's disease (AD) is a neurodegenerative disease that leads to progressive cognitive decline. While it's widely known that no cure exists for AD or terminating the neurodegeneration, early medical treatment might help to relieve the symptoms and slow the deterioration. To investigate and exemplify the benefit of proposed formalism within the context of this problem, we use data from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). The rfMRI subject data are processed by BrainSuite fMRI Pipeline (BFP) and grayordinate representations are generated 48,49 .
The cerebral cortex is modeled as a surface mesh and globular subcortical nuclei are modeled as volume parcels 50 . The grayordinate data were downsampled to 445 nodes. The Pearson correlation matrix was computed using fMRI time-series 51 .
To control the brain interactions matrix from late mild cognitive impairment (LMCI) patients and improve their brain structure, we apply our model and modulate it's multifractality to fit the one from cognitively normal (CN) people. We follow the same procedure as discussed in Hi-C case study. Fig. 6(a,b) illustrate the generalized degree distributions of CN matrix and LMCI matrix, respectively. The blue asterisks represent the target distribution from CN/LMCI matrix and the red circles are the distribution of recovered model via the simulated annealing-based reconstruction algorithm (see section 'Reconstructing weighted multifractal model' in  www.nature.com/scientificreports www.nature.com/scientificreports/ Methods). Simulated annealing (SA) algorithm provides a probabilistic approach for searching a large discrete space and finding approximate global optimum (e.g., parameters of a model) for a nonconvex optimization problem within a limited amount of time. Although it does not return the global optimum and it might not return identical results when running the SA algorithm multiple times, we find that one can distinguish between different networks from the retrieved WMG parameters. We individually run the simulated annealing algorithm on CN and LMCI matrix 200 times. We use 150 results of CN and 150 results of LMCI to train an SVM classifier with linear kernel, and using the remaining 50 of each class to test the dataset. The classification accuracy on the test dataset is 94%. It shows that the WMG model characterizes the network patterns and it could also discriminate between different traits of these complex networks. Fig. 6(c) shows the multifractal spectrum of LMCI matrix before and after the control. The red curve is the desired multifractal spectrum from CN data. The blue dots and magenta asterisks represent the multifractal spectrum curves of LMCI and optimal regulated solution obtained by adapting the LMCI network structure. The controlled multifractal spectrum also fits well with the target curve.

Discussion
We propose the weighted multifractal graph (WMG) model to capture the nodes attributes and interactions in weighted complex networks. We show that it fits a wide variety of multifractal spectrum and statistical properties of real-world networks including generalized degree distribution, clustering coefficients and joint degree distribution. Beyond the self-similar constructing procedure of the recursive generating model, we present the multifractal analysis of the random graph model and show that the WMG model could generate a variety of multifractal spectrum curves with different shapes and varying degree of asymmetry. More importantly, the introduction of weights in the WMG model brings flexibility in multifractal spectrum. The proposed model has potential applications in many disciplines including biological systems, geoscience, financial markets and social networks. In order to ground this extended model into real-world applications, one needs to develop rigorous, efficient and accurate model identification techniques that need to take into account the relationship between the complexity of real networks and the generating measures, assess the impact of model parameters on network properties (e.g., information flow performance, robustness), all of which are not trivial mathematical tasks.
Apart from modelling real networks, an important implication of the WMG model is to minimally control the multifractality of complex networks. The experiments on real-world datasets of yeast genome and Alzheimer's disease reveal that it is possible to regulate the multifractal spectrum of complex networks with minimal adjustments in the WMG model parameters. The proposed model is a step towards exploring the underlying growing mechanism of evolving networks such as brain structure of Alzheimer's disease patients, or actively growing yeast cells. Moreover, the control of multifractality can provide a novel way to treat brain diseases, control growing states or recover from potential attacks in biological and social systems.
Future extensions of the WMG model can consider the generalization of l i . In the proposed WMG model, we consider that the values of the interval length l i are identical for all R layers. However, one can consider that l i depends on r, but this will introduce more complexity and model identification strategies. Moreover, each node can have different category i r in different weight layer r and the linking probability would become p r ( ) Another important addition would be to have temporal dimension into the proposed model. This corresponds to modelling and controlling the time varying complex networks by regulating the multifractality of graphs. It will help to answer research questions such as how does network structure evolve over time from one state to another. For example, in decoding human behavior, what are the causes of changing genome or brain structure and what does it lead to. It is also crucial to understand the physical meaning of network multifractality and deciphering the hidden information of real-world complex networks with self-similar patterns.

Methods
Multifractal analysis. We study the multifractal analysis of the WMG model based on the partition function. Let can be defined as  Clustering coefficients. Clustering coefficient is a measure of graph cliquishness and neighborhood connectivity 5 . Following the definition of clustering coefficients in weighted graphs 42 and the multifractal network generator 32 , the clustering coefficient of a node in category i is written as c w l l p r p r p r w l l p r p r are the arithmetic mean of triangles and triplets edge weights. The numerator is the mean weights of triangles containing a node falling into sub-interval i, i = 1…M K , the denominator is the mean weights of triplets (link pairs) originating from that node. Joint degree distribution. The joint degree distribution of graphs with discrete weight set w(r) is defined as www.nature.com/scientificreports www.nature.com/scientificreports/ Four multifractal control problem simulations are shown in Fig. 8. The red lines are the desired multifractal spectrums, blue dot lines are the starting multifractal spectrums given by random parameters. Magenta asterisks are the optimal solutions given by simulated annealing algorithm. Different deviations of shifting and scaling are shown here. In each case, the multifractal spectrums are almost identical to the desired ones. It shows that by minimizing the cost function in equation (15), we could control network multifractality with minimal changes. Reconstructing weighted multifractal model. We use simulated annealing algorithm 53  where p * (d, r) and p(d, r) are the generalized degree distributions generated from the real-world network and the proposed WMG model by equation (4), ⟨ ⟩ c i is the average clustering coefficient in equation (5), c * is the empirical clustering coefficient. Equation (16) could also be redefined by other desired statistical properties of the networks. Minimizing Eq. (16) is a non-convex problem. And as the weight resolution in the WMG model increases, more parameters need to be estimated. Other optimization methods may also be applied.