The role of modularity in self-organisation dynamics in biological networks
Bram A. Siebert Cameron L. Hall James P. Gleeson Malbor Asllani
nlin.PS nlin.AO physics.soc-ph





^1MACSI, Department of Mathematics and Statistics, University of Limerick, Limerick V94 T9PX, Ireland
^2Department of Engineering Mathematics, University of Bristol, Bristol Woodland Road, Clifton BS8 1UB, UK



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 propose 

a self-organising mechanism which explains the formation of macroscopic spatial patterns, which are homogeneous within modules. This may explain how spontaneous order in biological networks follows their modular structural organisation. 
We test our results on real-world networks to confirm the important role of modularity in creating macro-scale patterns. 


The role of modularity in self-organisation dynamics in biological networks
    Bram A. Siebert^1, Cameron L. Hall^1,2, James P. Gleeson^1, Malbor Asllani^1 

    December 30, 2023
==================================================================================





§ INTRODUCTION




Patterns are macroscopic structures that are the distinctive mark of the self-organisation in a system of microscopic interacting entities <cit.>. 
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 <cit.>. 
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 <cit.>. 
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. 

It can be 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 <cit.>. 

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 activator-inhibitor systems in discrete systems such as regular lattices, and they have been used in this way to describe pattern formations in cellular tissues <cit.>.
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 <cit.>, ecological meta-populations <cit.> or coupled chemical reactors <cit.>, researchers have extended the reaction-diffusion formalism to complex biological networks <cit.>. 
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. 


Hütt et al. <cit.> recently argued that the formalism of activator-inhibitor systems is relevant to  the dynamical processes evolving in the brain <cit.> .

The implementation of network tools for analysing the brain's structure has been used since the first years of network science <cit.>. 


In their seminal work, Watts and Strogatz <cit.> studied the topology of the neuronal network of the nematode C. elegans and discovered that these networks possess a “small-world” property. 
In the literature, it has also been argued that many brain networks might be small-world networks <cit.>. 
It is 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 <cit.>.        

A further property of brain networks is that they are often modular <cit.> 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 functional role that the modularity of brain connections has been discussed from several perspectives. 
For example, due to the increased structural stability <cit.>, the modularity might have been crucial in the evolution and development of the brain. 

According to <cit.> 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 <cit.>. 
To ensure both a low shortest path length, and a high clustering coefficient, brain networks are organised in a strict hierarchical manner <cit.> 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 <cit.>. 
There are many examples of this:
in protein interaction networks, the proteins that share  similar functions are grouped together in modules <cit.>; in metabolic networks, there are structural/functional communities corresponding to cycles or pathways <cit.>; and in citation networks, scientific papers are clustered according to their research topic <cit.>.

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 <cit.>.
We focus on modular networks, which in contrast to many other random networks, are characterised by a small spectral gap, i.e., a small distance of the second largest eigenvalue [The definition of the spectral gap depends on the way one defines the Laplacian matrix <cit.>. In our case the spectrum of the Laplacian is non positive.] of the Laplacian from the origin. Let us notice here that a small spectral gap is a characteristic also of large (dense) regular graphs, however, here we focus on random graphs.
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. <ref>, 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 <cit.> for networked biological systems where Turing patterns are known to exist.
However, it is important to note that this method partitions the network in a similar fashion to the Fiedler partitioning. Therefore, it is possible to underestimate the total number of communities. Additionally, using pattern formation for community detection does not distinguish between functional communities and structural communities.

In this paper we begin in Sec. <ref> 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. <ref>.
We describe the different types of patterns which form in Sec. <ref>, and show how increasing the modularity helps in the formation of patterns.
Finally in Sec. <ref> we look for Turing patterns in some real world networks.





§ PATTERN FORMATION ON A NETWORKED SYSTEM



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)) <cit.>. 
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 <cit.>.
In general, a two-species reaction-diffusion model on a network of N nodes will take the form, 

    d u_i/d t   =    f(u_i,v_i) + D_u ∑_j L_iju_j, ∀ i=1,…, N 
    d v_i/d t   =     g(u_i,v_i) + D_v ∑_j L_ijv_j, ∀ i=1,…, N,

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 𝐋 is the graph Laplacian operator. The entries L_ij of the graph Laplacian are defined by L_ij = A_ij - k_iδ_ij, where 𝐀 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:

    d (δx)/d t = (Ĵ + DL̂)δx,

where δ𝐱=(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 

    D = 
        [ D_uI_N      0;      0 D_v𝐈_N ]

is the diffusion constant matrix.
Note that 𝐈_N represents the N by N identity matrix, so that 𝐃 is 2N by 2N.
The Jacobian matrix and the extended Laplacian are correspondingly 

    Ĵ=
    [ f_uI_N f_vI_N; g_uI_N g_vI_N ],    L̂=
    [ L 0; 0 L ].









Note here that the notation J will be reserved to identify the Jacobian of the 2 × 2 reactions matrix:

    J=
    [ f_u f_v; g_u g_v ].

















We then look for solutions to Eq. (<ref>) of the form 

    δu = ∑_α = 1^N b_α e^σ(Λ_α) tΦ^α,
    δv = ∑_α = 1^N c_α e^σ(Λ_α) tΦ^α,

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.




As will be seen in the following, the description of the linear solution through the eigenvectors of the Laplacian matrix will be essential in our analysis for the prediction of the modularity of final nonlinear patterns. In fact, depending on which eigenvalues are positive, we can control the final shape of the pattern, as in Fig. <ref>.

 


Following the standard approach described by <cit.>, we substitute the expansion of the perturbations into Eq. (<ref>). This decomposes the extended Jacobian to a 2× 2 matrix (for each index α) for which the eigenvalue problem needs to be solved,

    J_α =
    [ f_u + D_uΛ_α          f_v;          g_u g_v + D_vΛ_α ],

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 - 4 detJ_α))/2 is positive. 
relation between the eigenvalues of the extended Jacobian and the eigenvalues of the Laplacian, σ(Λ_α), is known in the  literature as the dispersion relation <cit.>, for the continuous version see appendix <ref>.
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 concerns detJ_α=detJ + (f_uD_v + g_vD_u) Λ_α + D_uD_vΛ^2_α.
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.
Recalling the previous definition of an activator-inhibitor system, u (respectively,  v) increases (respectively, decreases), the production of both species u and 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 <cit.>, 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  (J_α)<0.




        





 §.§ The case for D_v ≳ D_u




From experimental observations <cit.> 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 <cit.> 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 network 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

    Λ_min = -f_uρ + g_v/2D_v.

From relation (<ref>) 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.










§ TURING PATTERNS ON MODULAR NETWORKS



It has been argued that that the existence of particular topological features in many types of networks, including brain networks, are of crucial importance in several important processes from neuronal communication <cit.> to structural robustness <cit.>.
Such functional properties are based on the short average path length that characterises this family of networks. 
We emphasised in the preceding subsection that the spectral gap is an important ingredient for the Turing instability. 
In this section, we further illustrate  this fact by taking into account a special family of networks, the modular ones, that are known for for their lack of spectral gap. As a comparison we contrast the process of pattern formation in a non-modular network such as a Newman–Watts (NW) network (a particular case of a small-world network) with the pattern formation on a modular network generated using the Stochastic Block Model (SBM).

As described in <cit.>, modular structure has been identified in many brain networks.
Since the FitzHugh–Nagumo model <cit.> is both useful for modelling neuronal dynamics  <cit.>, and since it can exhibit spatial pattern formation <cit.>, 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. (<ref>) where a, b, and c are constants.
The parameters of the model are always chosen such that we have a stable fixed point.


In Fig. <ref> 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. <ref> b), the distribution of the eigenvalues of the Laplacian matrix for the NW network, which 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 ρ, or optimise the rewiring to minimise the diameter.

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 non-modular network described above, this is shown in Fig. <ref> 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 <cit.>.
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 of the Laplacian matrix as the modular eigenvalues (for example in Fig. <ref> d) the first four non-zero eigenvalues). 
In Fig. <ref> d), we observe that the modular eigenvalues  are sufficiently close to 0 and in the interval of possible values of the spectrum where the instability can develop; in Fig. <ref> 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 small-world 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 <cit.>, Bojan shows that the absolute value of the second largest Laplacian eigenvalue |Λ_2| is bounded below by |4Nd|, 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 lower bound of the spectral gap (equivalently, |Λ_2|) is larger when the diameter d is smaller; impeding this way a non-modular network like the NW under consideration, having a smaller spectral gap than a modular network. [blueWe want to emphasise that regular networks (e.g., rings) have a large diameter, too, having this way a small spectral gap. However, our focus here is on random networks which, apart from the modular ones, are characterised by a small diameter.]




To further investigate how the spectral gap changes for different network topologies, we look at three different networks in Fig. <ref>.
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 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 <cit.>.
We first look at an ER graph, as shown in Fig. <ref> b).  
Notice that there is a large spectral gap in the corresponding dispersion relation, as shown in Fig. <ref> 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. <ref> c) and d).
Finally in Fig. <ref> f) we have reduced the number of inter-edges such that patterns form, and the spectral gap is greatly decreased, as in Fig. <ref> e).
Notice again that there are two set of eigenvalues, the first four non-zero eigenvalues (which we refer to as modular eigenvalues) and the remaining non-modular eigenvalues.
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 algebraic connectivity theorem <cit.> 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. [ Notice also that due to the algebraic connectivity, a network with a Laplacian spectral gap will always be modular.] 
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 <cit.>.   



We notice from Figure <ref> 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 <cit.> and in particular in dynamics of the brain <cit.>. While <cit.> is mainly an experimental paper, and first highlights the observation of spatial patterns on brain networks, we have laid down a rigorous mathematical foundation that explores the importance of modularity to the formation of Turing patterns. Additionally, to the best of our knowledge, we here propose the first self-organising 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 v observed at equilibrium by constructing and analysing 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 <cit.>.

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. <ref> a), b), and c). 
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. <ref> d), e), and f).
A hybrid state is obtained instead when both sets of eigenvalues contribute to the Turing instability as in Fig. <ref> g), h), and i). 
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 eigenvalues. 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. <ref> h) 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. <ref> i) 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. <ref>.

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. <ref> a); in particular, the fact that the components of the modular eigenvectors are segregated accordingly to the respective modules [ Notice here that it may be, as in the case of Fig. <ref>, that different modules might share by chance the same level of components. However, this should not be understood as these entries belonging to the same module.].  
To shed light on this peculiarity we will resort again to  spectral graph theory. 

As anticipated earlier, the smallest non zero eigenvalue of the Laplacian Λ_2 defines the spectral gap known also in the literature as the Fiedler eigenvalue and defines the algebraic connectivity <cit.>.

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 <cit.>.
However it should be noted that the Fiedler partitioning can underestimate the total number of modules as we show in Fig. <ref>.
The other modular eigenvectors also behave in a similar manner to the Fiedler eigenvector; their entries are segregated by module <cit.>. 

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. <ref> 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 <cit.>.
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.
Note that in order for a Turing pattern to form, we must begin from a stable fixed point.


In Fig. <ref> 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).
One could find patterns in the other regions by tuning the diffusion parameters, except in the red region due to the absence of a stable fixed point.


The result that brain networks have optimised their spatial interaction matrix in order to allow pattern formation has been already claimed by experimental observers <cit.>; to the best of our knowledge we present the first mechanism that explains the role of modularity in achieving this pattern formation.   
        
   
  
        


§ 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 <cit.>.
In Fig. <ref> a) we show the final modular pattern of the nematode P. pacificus <cit.>.
This follows from the theoretical prediction of the unstable Fiedler eigenvector, shown in Fig. <ref>b)
Here we have used the Fiedler eigenvector to identify the communities of neurons <cit.>.
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 <cit.> 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 <cit.>. 
Such modular contact networks have also been shown to be crucial for the pattern of disease spreading <cit.>.
In Fig. <ref> 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 <cit.> 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 <cit.>, which are decisive for the spreading of an epidemics in the entire urban area <cit.>.

These examples all show agreement with the mathematical analysis we have shown so far.     




§ 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 self-organisation 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 <cit.> in several experimental observations <cit.>.
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 <cit.>. 

Indeed, the (self-)segregation of the network structure in modules <cit.> influences also the shape of the dynamical pattern on it. 
Based on the fact (see <cit.> and Fig.<ref>) 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 <cit.> 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 <cit.> or the spread of an epidemic in a group of individuals <cit.>. 
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 <cit.>.



§ ACKNOWLEDGEMENTS

B. A. Siebert acknowledges funding from the Irish Research Council under grant GOIPG/2018/3026. The  work of J. P. Gleeson and M. Asllani is partly funded by Science Foundation Ireland (grant numbers 16/IA/4470, 16/RC/3918, 12/RC/2289 P2, 18/CRT/6049) and co-funded under the European Regional Development Fund.



§ APPENDIX



 §.§ The FitzHugh–Nagumo model


We have used the Fitzhugh–Nagumo model throughout this paper <cit.>, which is one of the first and best-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

    d u/d t   =    u - u^3 - v
    d v/d t   =    c(u - av - b)

where u is the membrane potential and v the recovery variable. The model itself was first introduced by FitzHugh <cit.> 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. <cit.> 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 <cit.>. In recent years, with the rapid development of network science, the FitzHugh–Nagumo model has been extended to diffusively coupled networks <cit.>.



 §.§ Continuous Formulation


The original continuous framework for pattern formation, in one dimension, is

    ∂ u/∂ t   =    f(u,v) + D_u ∂^2u/∂ x^2,
    ∂ v/∂ t   =    g(u,v) + D_v ∂^2v/∂ x^2,

where notation is as in Sec. <ref>.
The derivation of Turing patterns follows the same process as we describe in the main text, except in a continuous form.
Of note, the extended Jacobian is now
    
    J_k = J - Dk^2 =
            [ J_11 - k^2 D_u           J_12;           J_21 J_22 - k^2 D_v ]

where k is the wave number.
Then, the continuous dispersion relation, σ(k), is plotted against the wave number, k^2 instead of the eigenvalues of the Laplacian.


 §.§ 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 topological features, 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 <cit.>. 
In fact, most of the connectomes studied are organised in a modular structure, however each module is further organised in a small-world fashion. This is another amazing observation how nature tends to self-organize in order to better optimise the benefit from the both structural features, the modularity and the small-worldness.
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. <ref> 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. <ref>.


 
Thus, in the presence of hierarchy, the cyan and the yellow regions in Fig. <ref> 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 the modular networks studied in the main text.
   







