{'fraction_non_alphanumeric': 0.02900452170653422, 'fraction_numerical': 0.008926270412107865, 'mean_word_length': 4.7378679930098135, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The role of modularity in self-organisation dynamics in biological networks

Interconnected ensembles of biological entities are perhaps some of the most complex systems that modern science has encountered so far. In particular, scientists have concentrated on understanding how the complexity of the interacting structure between different neurons, proteins or species influences the functioning of their respective systems. It is well-established that many biological networks are constructed in a highly hierarchical way with two main properties: short average paths that join two apparently distant nodes (neuronal, species, or protein patches) and a high proportion of nodes in modular aggregations. Although several hypotheses have been proposed so far, still little is known about the relation of the modules with the dynamical activity in such biological systems. Here we show that network modularity is a key ingredient for the formation of self-organising patterns of functional activity, independently of the topological peculiarities of the structure of the modules. In particular, we show that macroscopic spatial patterns at the modular scale can develop in this case, which may explain how spontaneous order in biological networks follows their modular structural organisation. Our results also show that Turing patterns on biological complex networks can be a signature of the presence of modular structure and consequently a possible protocol for community detection. We test our results on real-world networks to confirm the important role of modularity in creating macro-scale patterns.


I. INTRODUCTION
Patterns are macroscopic structures that are the distinctive mark of the self-organisation in a system of microscopic interacting entities [1]. They are ubiquitous in nature and can be seen in the spots of a leopard's fur or the coloured scales of a butterfly's wing [2]. In 1952, Alan Turing published his seminal work on pattern formation, The Chemical Basis of Morphogenesis where he laid down an elegant and plausible theory that can be used to explain the formation of patterns [3]. Turing developed a simple model of pattern formation that established the minimal requirements for a biochemical system to self-organise. Turing's minimal system is composed of two "competing" chemicals, an activator and an inhibitor, which share the same spatial domain where they react and diffuse.
Based on a diffusion-driven instability mechanism, today known as Turing instability, Turing showed that it is possible to explain and predict the growth of spatially inhomogeneous perturbations away from a spatially homogeneous steady state. These perturbations in concentration are later stabilised by nonlinearities in the system, yielding the celebrated Turing patterns. The emergence of stable heterogeneous patterns, the result of Turing instabilities, was initially at odds with the general understanding that diffusion by itself is a smoothing process. Nevertheless, it has been shown that the right combination of short-range activation and long-range inhibition, caused by slowly diffusing activators and rapidly diffusing inhibitors, enables the pattern forming phenomenon [4].
Conventionally, an activator-inhibitor system is modelled using a set of reaction-diffusion equations that describe the evolution of the concentrations of activator and inhibitor throughout a continuous medium. These equations can readily be adapted to describe activatorinhibitor systems in discrete systems such as regular lattices, and they have been used in this way to describe pattern formations in cellular tissues [4,5]. However, biological tissue often takes more complex forms, and the spatial support cannot always be adequately formulated via regular lattices. Inspired by the network structures of early stages of embryogenesis [6], ecological meta-populations [7] or coupled chemical reactors [8], researchers have extended the reaction-diffusion formalism to complex biological networks [5,[9][10][11]. These discrete structures consist of graphs where the nodes usually represent the cells inside which reactions occur, and the edges usually represent the routes through which cells communicate by exchanging chemicals. Abundant studies in recent years have concentrated on the effect that different topologies of interactions (a peculiarity of the given cellular tissue) have on the formation or destruction of patterns [9,10,12,13] Hütt et al. [14] recently argued that the formalism of activator-inhibitor systems is relevant to the dynamical processes evolving in the brain [14] . The brain consists of billions of cells, the neurons, connected together into an extremely complex system [15,16]. The implementation of network tools for analysing the brain's structure has been used since the first years of network science [17]. In their seminal work, Watts and Strogatz [17] studied the topology of the neuronal network of the nematode C. elegans and discovered that these networks possess small average shortest paths, as well as a high clustering coefficient They described networks with these two properties as being "small-world" networks. Later, it was verified that many brain networks are small world networks [18][19][20]. It has been ar-gued and widely accepted that the small-world property of brain connectomes should help the communication between neurons inside the brain by integrating multiple segregated sources of information [21].
A further property of brain networks is that they are often modular [18] so that the neurons can be segregated into communities (referred to as modules) where two neurons chosen at random from the same module are much more likely to be connected than two neurons chosen at random from different modules. The effect of network modularity on pattern formation will be the main focus of this paper.
The functional role that the topology of brain connections has been discussed from several perspectives. For example, due to the increased structural stability [21,22], the modularity might have been crucial in the evolution and development of the brain. According to [18,22] modular topology can also optimise the wiring cost in the case of spatial networks. A small number of long range (and thus costly) connections reduces the diameter of the network, and allows the remaining nodes, now grouped into communities or modules, to form dense small world networks. Also, more compact segregation of neurons may contribute to the specialisation of the neurons in their functional duties [21]. To ensure both a low shortest path length, and a high clustering coefficient, brain networks are organised in a strict hierarchical manner [16,18,23] where at the first level of the hierarchy sets of nodes (the modules) are connected to mimic a small-world topology and the same happens at the second level of hierarchy and so on, until the single node level. For a more detailed discussion of the role of the hierarchy in the pattern formation process see the Appendix.
More generally, modularity is a common topological property that naturally emerges in biological, ecological, and social scenarios where the different communities are associated with different functions of the system represented by the network as a whole [24]. There are many examples of this: in protein interaction networks, the proteins that share similar functions are grouped together in modules [25]; in metabolic networks, there are structural/functional communities corresponding to cycles or pathways [26]; and in citation networks, scientific papers are clustered according to their research topic [27]. In addition to these properties, in this paper, we propose a new mathematical mechanism that highlights the role that modularity takes in self-organising processes in biological networks.
Using the Turing theory of pattern formation, we show that spatially extended patterns can be triggered by the segregation of the nodes (neurons) in distinguishable communities. To formally analyse the chances of such networks self-organising, we use a linear stability approach known in the literature as the dispersion relation [2]. We focus on modular networks, which (in contrast to other networks, e.g., small-world ones), are characterised by a small spectral gap, i.e., a small distance of the second largest eigenvalue [28] of the Laplacian from the origin. To anticipate some of the technical details, we discuss the key features of modular networks in the following paragraphs and outline how these affect pattern formation.
For modular networks, the Laplacian eigenvalues that may be responsible for the Turing instability can be split into two sets. In one set, we have the eigenvalues emerging due to the global modularity of the network, which we denote as "modular eigenvalues". In Sec. III, we will show that when only this part of the spectrum is responsible for the instability, then the shape of the associated pattern follows that of the network in the sense that nodes belonging to the same modules have very similar concentrations of the species among themselves but these concentrations are distinctly different from the concentrations in other modules. In contrast, if the instability is caused by the remaining set of eigenvalues, which correspond to the local connectivity of nodes, here denoted as "non-modular eigenvalues", then all the nodes have (in principle) different concentrations making the pattern globally heterogeneous. In this latter case, if the eigenvalues responsible for the instability are limited to the eigenvalues belonging to a single module, then the pattern will first emerge in that module.
We aim to create a bridge between the role of the structure in many biological networks with the dynamical activity therein. In particular, in our model, we explain how communities of biological entities (cells, individuals, etc.) can act as functional units in their corresponding biological systems. As a consequence, we argue that this approach can potentially be used in community detection methods [29][30][31] for networked biological systems where Turing patterns are known to exist.
In this paper we begin in Sec. II with a description of the mathematical background of Turing patterns. This will lead us into a discussion as to why modularity is critical to the formation of patterns in Sec. II A. We describe the different types of patterns which form in Sec. III, and show how increasing the modularity helps in the formation of patterns. Finally in Sec. IV we look for Turing patterns in some real world networks.

II. PATTERN FORMATION ON A NETWORKED SYSTEM
The process of the formation of Turing patterns was originally put forward by Turing [3] on continuous domains to explain the emergence of spontaneous order in biological contexts, but it was formulated on networked systems only in the 1970s by Othmer & Scriven [5,9]. In a continuous domain, the most simple Turing mechanism is given in terms of reaction-diffusion equations that describe the evolution through time and space of the concentrations of two competing chemical species, called the activator (with concentration denoted u(x, t)) and the inhibitor (with concentration denoted v(x, t)) [2,3]. In general, an activator increases production of both itself and the inhibitor. The inhibitor, in turn, slows down the growth in activator. When the spatial support is instead discrete, constituted by spatial patches (nodes) connected through communicating routes (links) the reaction-diffusion mechanism can be formulated using ODEs, instead of PDEs [5]. In general, a two-species reaction-diffusion model on a network of N nodes will take the form, where u i and v i represent the concentrations of activator and inhibitor respectively at node i, f and g are nonlinear functions that describe the net production rates of activator and inhibitor respectively, D u and D v are the diffusion coefficients of activator and inhibitor respectively, and L is the graph Laplacian operator. The entries L ij of the graph Laplacian are defined by where A is the adjacency matrix, k i is the degree of node i, δ is the Kronecker delta, and where we do not sum over repeated indices. In order to understand the development of spatial patterns, we analyse the linear stability of the system starting from a homogeneous steady state (u * , v * ) that is stable in the absence of diffusion. If the diffusion coefficients are nonzero and the ratio ρ = D v /D u is large enough, the steady state (u * , v * ) becomes unstable and small random perturbations of the previous steady state will grow. This growth is exponential in the initial linear regime, and may then be stabilised by the nonlinear terms of the functions f and g so that the system reaches a stable but spatially inhomogeneous steady state. Such a mechanism is responsible for the emergence of Turing patterns.
The linearised system in matrix form reads: where δx = (u − u * 1 N , v − v * 1 N ) is the perturbations vector of the activator u and inhibitor v species, 1 N is the all-ones N −dimensional vector, and is the diffusion constant matrix. Note that I N represents the N by N identity matrix, so that D is 2N by 2N . The Jacobian matrix and the extended Laplacian are corre-spondinglyĴ Note here that the notation J will be reserved to identify the Jacobian of the 2 × 2 reactions matrix: We then look for solutions to Eq. (2) of the form where Λ α , Φ α are respectively the eigenvalues and eigenvectors of the Laplacian L matrix, σ(Λ α ) are the eigenvalues of the extended Jacobian (Ĵ + DL), and α is the index term. Following the standard approach described by [2,5,9], we substitute the expansion of the perturbations into Eq. (2), which decomposes the extended Jacobian to a 2 × 2 eigenvalue problem for each index α: where subscripts on the activation function f (u, v) and the inhibition function g(u, v) represent partial derivatives evaluated at (u * , v * ) . To study the stability of the linear system we look for positive real parts of the eigenvalues of J α . Turing instability occurs when the real part of the larger of the two eigenvalues σ(Λ α ) = trJ α + (trJ α ) 2 − 4detJ α 2 is positive. The function σ(Λ α ), is known in the literature as the dispersion relation [2]. For an activator-inhibitor system the necessary conditions for stability are trJ α < 0 and detJ α > 0. The first condition is always true, since trJ α = trJ + (D u + D v )Λ α , and this is negative since the stability of the fixed point in the absence of diffusion implies that trJ < 0, while the non-positivity of the Laplacian spectrum implies Λ α < 0. We therefore turn our attention to the second condition for stability, which In order for a Turing instability to occur, we require detJ α < 0. Noting that the stability of the fixed point in the absence of diffusion implies that det(J) > 0 and noting that Λ α < 0, it is straightforward to conclude that the only way for detJ α to be negative is for (f u D v + g v D u ) to be positive. Without loss of generality we define u to be the activator and v to be the inhibitor. Since u corresponds to the activator, its increased presence increases the production of u; since v corresponds to the inhibitor, its increased presence decreases the production of v; as a result of this, the signs of the respective partial derivatives are f u > 0, and g v < 0. Therefore, we require ρ = D v /D u > 1 for instability, implying that the inhibitor should diffuse faster than the activator in order for Turing patterns to arise. In many practical cases, this difference needs to be very large in order to achieve det(J α ) > 0 . Notice here the presence of unstable eigenvalues (inset) and that the eigenvalues are separated in two sets by an important gap, between the first and second set of eigenvalues. The first four non-zero one are denoted as the modular eigenvalues and the rest non-negative ones as the non-modular eigenvalues. The parameters of the Fitzhugh-Nagumo model are in both cases Du = 1, ρ = 5.5, a = 0.7, b = 0.05, c = 1.7. Finally, note the different colormaps used between panels a) and c) to highlight the lack of patterns in the former.

A. The case for Dv Du
From experimental observations [32][33][34][35] it is rarely true that the inhibitor diffuses much faster than the activator, but instead the chemicals diffuse with similar rates. In the case where D v D u , it can be shown that the dispersion relation is positive only for values of the spectrum of the Laplacian very near to the origin. To prove this we analyse the behaviour of det(J α ) when considered as a function of Λ α ; more precisely, we focus on the value of Λ α corresponding to a minimum of det(J α ). It is known in literature [2] that for the continuous case, it will always exist a non-positive value of Λ α such that the det(J α ) < 0 or, in other words, that Turing instability can occur. In order to proceed with our analysis, in the following, we will consider that Λ α takes continuous values and will see that the spectrum of a (strongly) modular networks fits in the domain of the continuous dispersion relation for which the instability occurs for the particular case, D v D u . We start by differentiating with respect to Λ α and after some algebraic manipulation, we find that the minimum of detJ α is found at Λ α = Λ min where From relation (5) we note that if D v is kept fixed while ρ → 1 then Λ min → 0. To show this we set ρ = 1 + . Under the conditions of the Turing instability, Λ min is non-positive, so (1 + )f u + g v > 0. Rearranging, we can write (1 + )f u + g v = trJ + f u and, noting that trJ is necessarily negative, we conclude that the positive quantity trJ + f u can be at most of order , since f u > |trJ|. This shows that Λ min is of order . Therefore, as decreases, the value of Λ α for which detJ α is at its minimum tends towards the origin. Hence, the possible values of Λ α that may permit Turing instabilities tend towards zero as the ratio ρ of diffusivities tends to 1. In practice, this implies that the range of values of Λ α for which instabilities can occur decreases in size and is restricted to small values of Λ α . Therefore, a small spectral gap is needed to allow patterns to form. This is significant for the analysis of modular networks that follows since, as shown in the following section, modular networks are characterised by a small spectral gap |Λ 2 − Λ 1 |. Hence the Laplacian of a modular network will have eigenvalues close to the origin. Because of this, we are able to find modular networks where Turing instabilities, and thus pattern formation, may occur where otherwise (i.e., in non-modular networked systems) they would not. This modular pattern formation may even occur for values of ρ that are close to those observed in real systems.

III. TURING PATTERNS ON MODULAR NETWORKS
It has been argued that that the existence of a smallworld topology in many types of networks, including brain networks, is of a crucial importance in several important processes from neuronal communication [38] to structural robustness [39]. Such functional properties are based on the short average path length that characterises this family of networks. Nevertheless these advantages also present a counter effect, that of a large spectral gap. We emphasised in the preceding subsection that the spectral gap is an important ingredient for the Turing instability. In this section, we illustrate the importance of a small spectral gap to Turing instabilities by comparing pattern formation on a Newman-Watts (NW) network (as an example of a small-world network) with pattern formation on a modular network generated using the Stochastic Block Model (SBM).
As described in [18,21,40], modular structure has been identified in many brain networks. Since the FitzHugh-Nagumo model [41,42] is both useful for modelling neuronal dynamics [2], and since it can exhibit spatial pattern formation [2,11], we will use this model throughout this paper. In dimensionless form, FitzHugh-Nagumo dynamics correspond to using the functions f (u, v) = u − u 3 − v and g(u, v) = c(u − a + bv) to describe the net production of activator and inhibitor in Eq. (1) where a, b, and c are constants.
In Fig. 1 we compare the pattern on a single-module NW network (of 125 nodes and 660 edges) and a modular network with 5 communities, each with 25 nodes and a local Erdős-Rényi (ER) topology. As can be observed from the dispersion relation in Fig. 1 b), the distribution of the eigenvalues of the Laplacian matrix for the NW network shows a large spectral gap. This makes the Turing instability impossible for the given choice of parameters (including ρ = 5.5), since the instability (i.e., values of Λ corresponding to positive values of the continuous curve) is concentrated near the origin. We could potentially create an instability by significantly increasing ρ. As t → ∞, the Fitzhugh-Nagumo models considered in this paper will tend to an equilibrium. One way to depict these equilibria is to plot the concentration of the activator species at long times. For the NW network described above, this is shown in Fig. 1 a) and we see that the activator concentration is homogeneous across all nodes as expected.
In contrast to this, for a strongly modular topology the spectrum is divided into two distinct sets of eigenvalues. The first set is those nonzero eigenvalues near the origin (of which there are M − 1 where M is the number of the modules) and the second set is composed of all the remaining eigenvalues that are far from the origin [43]. We note that both the NW network and the modular network have the same number of nodes and edges, so the difference between the networks' spectra cannot be attributed to a difference in the number of nodes or in the average degree of these nodes. As already anticipated, we will refer to the first set of nonzero eigenvalues as the modular eigenvalues (for example in Fig. 1 d) the first four non-zero eigenvalues). In Fig. 1 d) we observe that the modular eigenvalues are sufficiently close to 0 that one of them corresponds to an instability; in Fig. 1  c) we see that this leads to a pattern in the activator concentrations at equilibrium.
To understand the reason why the spectrum of a modular network can be divided into two subsets we should first explain the reason behind the spectral gap in smallworld networks. As mentioned earlier, the denomination "small-world" refers to a certain class of networks, one feature of which is the small average distance between nodes. In [44], Bojan shows that the absolute value of the second largest Laplacian eigenvalue Λ 2 is bounded below by 4 N d , where N is the number of nodes in the network and d is the diameter. This means that for a fixed value of the size N of the network, the spectral gap (equivalently, |Λ 2 |) is larger when the diameter d is smaller; thus, a NW network will have a larger spectral gap than a modular network.
To further investigate how the spectral gap changes for different network topologies, we look at three different networks in Fig. 2. We create these networks in a simple way. First we divide our 125 nodes into five modules of nodes, and define the total number of intra-edges (connections within modules) and inter-edges (connections between modules). Then we allocate each module an equal number of intra-edges and inter-edges and randomly connect nodes within and between the modules, while avoiding double entries in both cases. If we define the number of intra-edges to be much larger than Also we used the algorithms described in refs. [36,37] with resolution parameter γ = 1 the modularity Q in each case.
the number of inter-edges, then this process will yield a network with as strong modular structure. We describe three examples of these networks with increasing "modularity", where modularity is defined by the Q function described in [36,37]. We first look at an ER graph, as shown in Fig. 2 b). Notice that there is a large spectral gap in the corresponding dispersion relation, as shown in Fig. 2 a). By simply modifying the ratio of inter-edges to intra-edges, we can then generate a new network which begins to close the spectral gap, as in Fig. 2 c) and d). Finally in Fig. 2 f ) we have reduced the number of inter-edges such that patterns form, and the spectral gap is greatly decreased, as in Fig. 2 g). This leads us to ask why a highly modular network closes the spectral gap so well. Note that in the Appendix, we consider the hierarchical case where each module is arranged in a small-world fashion.
To understand the small spectral gap of modular networks, we first imagine a scenario in which the modules are disconnected from each other. Individually, these modules are denser and smaller than the Newman-Watts network, therefore each of them is expected to have a relatively large spectral gap. From the Perron-Frobenius theorem [45] we know that the number M of the connected components (the modules in this case) corresponds with the number of zero eigenvalues of the Laplacian operator. However, once the modules are connected with a small number of links then M − 1 of these eigenvalues will move away from zero. From a spectrum perturbative analysis, we find that these become very small nonzero eigenvalues, with only one zero eigenvalue still remaining to signify that the whole modular network is connected. This explains the small size of the spectral gap in modular networks and consequently the emergence or not of Turing patterns respectively in small-world and modular networks [46,47].
We notice from Figure 1 c) that although the pattern is highly heterogeneous at a global level, the patterns on nodes within each single module are quite homogeneous, having almost the same concentration of the species for each node in the module. Such macroscopic spatially extended patterns where densely connected entities (e.g., of biological nature) show the same amount of activity have been observed in different biological contexts [24][25][26] and in particular in dynamics of the brain [48,49]. To the best of our knowledge, we here propose the first selforganising mechanism that explains the uniformity at the module level of Turing patterns in biological networks.
We can obtain insight into the patterns of u and the eigenvectors associated with the Turing instabilities.
From an initial condition close to the unstable homogeneous equilibrium, the rate of change in the concentrations u and v will initially be dominated by the eigenvector associated with the largest positive eigenvalue of the Jacobian. This initial growth will ultimately be stabilised by nonlinear terms, and we expect that the state equilibrium pattern of concentrations will be reminiscent of the eigenvectors associated with the instability [2,3]. To begin our analysis of the resultant patterns, we select parameters which lead to a single modular eigenvalue being positive, and observe the final "homogeneous by module" pattern as in Fig. 3 a), b). The situation changes when the instability is exclusively induced from the non modular eigenvalues. In this case the concentration is no longer uniform for each module as shown in Fig. 3 c), d). A hybrid state is obtained instead when both sets of eigenvalues contribute to the Turing instability as in Fig. 3 e), f ). These hybrid states can lead to patterns that are similar to either the modular patterns or the heterogeneous patterns. This is because the Turing instability in this case involves a competition between the eigenvectors associated with the unstable modular eigenvalues and the eigenvectors associated with the non-modular eigenval-ues. The dominant instability (and therefore the eigenvector that we expect to be most similar to the equilibrium pattern) will be the eigenvector associated with the largest eigenvalue of the Jacobian. In Fig. 3 f ) for example, we observe that the largest eigenvalue of the Jacobian is associated with one of the modular eigenvalues of the Laplacian, and this is associated with a pattern in Fig.  3 e) that could be described as almost being modular.
In the Supplementary Material (SM) we discuss several criteria to establish which eigenvalue is dominating over the others.
In order to understand why the final shape of the pattern can be modular we focus on the study of the eigenvectors as plotted in Fig. 4. From the stability analysis we know that initially the pattern is shaped according to the unstable eigenvectors and this form is largely retained in the final nonlinear regime. Nevertheless, what surprises is the particular form of the eigenvectors associated with the modular eigenvalues as in Fig. 4 a); in particular, the fact that the components of the modular eigenvectors are very small in all but one module. To shed light on this peculiarity we will resort again to spectral perturbation theory. , ρ = 6. The portion of the parameter space indicated in red represents the region where no Turing patterns are allowed, as the system is not in a steady state. The blue part is when the system is Turing stable, that is the system is at a steady state but no Turing patterns form. The rest of the region is when patterns may occur: in the magenta region patterns only form in the continuous domain case, in green we have "modular" patterns, Fig. 3 a), yellow "mixed" state patterns, Fig. 3 c), and cyan heterogeneous patterns, Fig. 3 b). Notice that as the ratio of diffusivites approaches 1, ρ → 1, the only patterns which form are the modular patterns, showing that in a real scenario modularity is the only way to induce pattern formation in networks with modular structures, e.g., brain networks.
of the Laplacian Λ 2 defines the spectral gap known also in the literature as the Fiedler eigenvalue and defines the algebraic connectivity [50,51]. Its corresponding eigenvector is known as the Fiedler eigenvector and has the property that the entries of the nodes corresponding to the same modules take very similar values. Because of this property, the Fiedler eigenvector has been extensively used as the basis of several community detection methods [36,46,47]. The other modular eigenvectors also behave in a similar manner to the Fiedler eigenvector; their entries are segregated by module [46,47]. Since the modular eigenvectors are often the fastest growing modes in the Turing instability, this means that the modular shape of the global pattern is a consequence of the modularity of the structure of the network itself.
On the other hand, when the instability is caused strictly by the non modular eigenvalues, another behaviour occurs during the pattern forming phenomenon. This is best considered by again considering a modular network to be a perturbation of a network with initially M disconnected components. In such a case, each nonzero eigenvalue of the Laplacian will correspond to an eigenvector whose components are all zero outside a single component. A modular network will be a small perturbation to this, and so the non-modular eigenvectors will also be close to zero except within a single component. If only one non-modular eigenvalue corresponds to a Turing instability, then only one module of the network will show pattern formation, as illustrated in Fig. 4  b). Thus, we can predict the module on which pattern formation will occur by looking at the components of the eigenvector whose eigenvalue corresponds to the fastest growing mode of the Turing instability.
So far we have considered the contribution in the formation of patterns of both modular and non modular eigenvalues, however when we deal with Turing patterns in real scenarios the ratio ρ = D v /D u is quite close to one [32][33][34][35]. To evaluate the conditions under which different patterns form in real conditions we now explore the parameter space of the Fitzhugh-Nagumo dynamics in more detail.
In Fig. 5 it can be observed that although different types of patterns can be found in the space of the parameters a and c, as the ratio of diffusivities gets closer to 1 the region where patterns can occur shrinks and, more importantly, the only possible Turing patterns are modular ones (indicated in green colour). The result that brain networks have optimised their spatial interaction matrix in order to allow pattern formation has been already claimed by experimental observers [48,49]; to the best of our knowledge we present the first mechanism that explains the role of modularity in achieving this pattern formation.

IV. SELF-ORGANISATION IN REAL MODULAR NETWORKS
Heretofore we have discussed the role of modularity in the formation of patterns only for synthetic networks. In this part we will illustrate our findings in real examples of biological or ecological networks. The neuronal networks of several primitive animals such as nematodes have been well characterised. Indeed, it was the study of nematode neuronal networks that first inspired the development of small-world network models [17]. In Fig. 6a) we show the final modular pattern of the nematode P. pacificus [52]. This follows from the theoretical prediction of the unstable Fiedler eigenvector, shown in Fig. 6b) Here we have used the Fiedler eigenvector to identify the communities of neurons [36]. In this particular case two modules are clearly distinguishable and the level of activity of the nodes inside the modules are quite homogeneous. Other examples of Turing patterns in neuronal networks are presented in the Supplementary Material. Although the modularity of brain networks has been well-studied [18,21,38,40] other types of natural networks manifest this property also. For instance, this is the case for ecological networks where the individuals are connected to each other through trophic relations [2,53]. Such modular contact networks have also been shown to be crucial for the pattern of disease spreading [54,55]. In Fig. 6 c) and d) we present respectively the equilibrium pattern of the FitzHugh-Nagumo equations and its comparison to the unstable eigenvector of the contact network of a zebra herd [53] where a community of 11 individuals out of a total of 23 is clearly visible. However, the formation of patterns of spreading are not limited only to contact networks, which in general can be small in size. Modularity is a common property in other types of networks which, although they are not directly related to biological systems, are still essential for biological phenomena occurring on them. This is for instance, the case for networks of human mobility, such the roads networks in the city of Chicago presented in the SM [56][57][58][59], which are decisive for the spreading of an epidemics in the entire urban area [54,55]. These examples all show agreement with the mathematical analysis we have shown so far.

V. DISCUSSION AND CONCLUSIONS
In this paper we have analytically and numerically explored pattern formation on modular networks. We have shown that modularity, a ubiquitous topological feature of many biological networks, is crucial for the selforganisation of the global dynamics on a network. To study this behaviour we have considered here the Turing instability as a paradigmatic mechanism for pattern formation in biology, ecology or neuroscience. The possibility of pattern formation via the Turing mechanism on non-modular networks is limited to unrealistically extreme ratios of the diffusion constants of the activator and inhibitor species making the small spectral gap of the Laplacian matrix a fundamental requirement for the Turing instability. This feature is a structural advantage of modular networks which follows from spectral perturbation theory. A strongly modular network can be considered as a set of connected components weakly attached with a small number of intermodule links. From spectral perturbation theory this yields a number -equal to one fewer than the number of modules -of non zero eigenvalues very near to the origin. This characterisation at the linear stability level influences the shape of the spatially extended patterns. Due to the segregation of the entries of the eigenvectors corresponding to the set of modular eigenvalues, we are able to explain why Turing patterns are homogeneous per module on these networks.
This result opens to an important aspect regarding the functional resolution of the brain modes which was hypothesised [14,18,60] in several experimental observations [48,49]. To the best of our knowledge, the model we present here constitutes the first self-organising mechanism where the modules are presented as functional blocks of biological networks. In this sense, we argue that the module is the smallest spatial unit to be taken into account from the functional point of view i.e. if we "zoom" out far enough from a modular network, the individual modules behave like individual supernodes. For the particular example of the brain the modules might be the super-nodes of the functional connectomes [40,48,60]. Indeed, the (self-)segregation of the network structure in modules [24] influences also the shape of the dynamical pattern on it. Based on the fact (see [14,18,40] and Fig.6) that in real scenarios Turing patterns should be exclusively modular, we believe that the results we have shown here can be potentially used to formulate a community detection protocol [24,31,36] in the case where patterns of self-organised activity are known to exist.
In the case when we relax Turing conditions to allow the instability for the non-modular part of the spectrum, then we can use the eigenvector corresponding to the largest eigenvalue to indicate the module in which the Turing pattern is first seeded before finally spreading to the rest of the network. This behaviour can potentially make the pattern formation process a powerful diagnostic tool for studying and eventually controlling the emergence of abnormal dynamics which characterise many neurological diseases [61] or the spread of an epidemic in a group of individuals [54,55]. We test our theoretical results on several real connection data sets of neuronal, ecological and infrastructure networks verifying the correctness of our findings, that modularity is crucial for the development of patterns, and that when the instability is derived from the first set of modular eigenvalues, that the resultant self-organisation follows the modular structure of the network.
The results we have presented here can extend also to more complicated scenarios. This is, for example, the case when the hierarchy of a network is considered as a complement to its modularity. In the Appendix we show that in a hierarchical modular network the modular eigenvalues are even more relevant for the Turing pattern forming process. Further extensions of our approach are also possible; for example to consider the effect of directed edges in a modular network. In this case we expect a richer dynamics where travelling Turing waves should emerge in a directed modular networks [11]. We have used the Fitzhugh-Nagumo model throughout this paper [41,42], which is one of the first and most well known mathematical models used to describe the spiking dynamics of neurons. In terms of mathematical equations the behaviour of a single neuron is described by where u is the membrane potential and v the recovery variable. The model itself was first introduced by FitzHugh [41] to explain the generation of spikes in excitable systems, i.e., neurons. A spike is a short-lasting elevation of the membrane voltage u diminished over time by a slower and linear recovery variable v once the system is periodically excited by an external current. The following year Nagumo et al. [42] developed the electric circuit which mimics such behaviour. However, although the model itself is mainly used to describe the oscillatory behavior of neurons, it also admits a stable fixed point, which is a necessary requirement for Turing instabilities. Once this model is equipped with a diffusion term, it turns out in a perfect candidate for pattern formation [2]. In recent years, with the rapid development of network science, the Fitzhugh-Nagumo model has been extended to diffusively coupled networks [11,62].
B. The role of hierarchy of the brain networks in the pattern formation We have discussed the role that modularity has on pattern formation, isolating it from other important features such as the small-world property, which is in fact an integral aspect of many networks, including brain networks. So then a question that arises naturally is, how does the brain cope with maintaining both features and their functional properties at the same time? We now are able to answer this question by recalling an important empirical results that characterises most real networks, their hierarchical structure [18,31,63]. In fact, most of the connectomes studied are organised in a modular structure, however each module is further organised in a small-world fashion. In a hierarchical modular network the entire network is organised in modules which are attached to each other so as to have a small diameter and at the same time the nodes in the modules are connected in such way to form sub-modules again minimising their diameter and this process goes on this way up to smallest building unity, the single nodes. A hierarchical structure stresses once more the necessity of modularity for the self-organising phenomena in the networks. In Fig. 7 we show that the difference of the smallest non modular eigenvalue from the origin is larger when the modules have a small-world topology compared to when they are organised at random (e.g. ER network) for the same number of nodes, edges and modules. The reason for this can be found once more by taking a perturbative approach. The spectral gap of an individual module (disconnected from the rest of the network) is larger when its diameter is smaller, as it is in the Newman-Watts network used in Fig. 7.
Thus, in the presence of hierarchy, the cyan and the yellow regions in Fig. 5 would be even smaller making the modularity region shown in green larger compared to the previous two. We notice, however, that the instability invariance is still valid for values of the diffusivites ratio ρ near to 1, that is when only the green region in the parameter space is available. In conclusion, a hierarchical arrangement where each module is arranged in a small world fashion, and these modules are again connected in a small world fashion, are even better candidates for forming modular patterns, than then modular networks studied in the main part of this paper.