An Analytically Solvable Model for Rapid Evolution of Modular Structure

Biological systems often display modularity, in the sense that they can be decomposed into nearly independent subsystems. Recent studies have suggested that modular structure can spontaneously emerge if goals (environments) change over time, such that each new goal shares the same set of sub-problems with previous goals. Such modularly varying goals can also dramatically speed up evolution, relative to evolution under a constant goal. These studies were based on simulations of model systems, such as logic circuits and RNA structure, which are generally not easy to treat analytically. We present, here, a simple model for evolution under modularly varying goals that can be solved analytically. This model helps to understand some of the fundamental mechanisms that lead to rapid emergence of modular structure under modularly varying goals. In particular, the model suggests a mechanism for the dramatic speedup in evolution observed under such temporally varying goals.


Modularity measure
In this section, we define the modularity measure used to quantify the modularity of the interaction network defined by the matrix A. An NxN matrix A can be mapped to a 2N-node network. In this mapping, the matrix element A(i,j) represents the weighted directed interaction between input-node j to output-node (N+i) (as described in Fig. S1).
To quantify network modularity we employed the measure based on the approach of Newman and Girvan [1,2], refined to weighted directed networks.
Briefly, the Newman and Girvan algorithm finds the division of the nodes into modules that maximizes a measure Q. This measure is defined by the fraction of the edges in the network that connect between nodes in a module (or 'community' as termed by the authors) minus the expected value of the same quantity in a network with the same assignment of nodes into modules but random connections between the nodes [3]. where K is the number of modules, W is the total sum of the interaction weights in the network, w s is the total sum of the interaction weights between nodes in module s, d s in is the sum of the weighted in-degrees of the nodes in module s and d s out is the sum of the weighted out-degrees of the nodes in module s. We assume all weights are nonnegative, otherwise absolute values are taken.
We further refined this measure, as described in [4], by normalizing it with respect to randomized networks, (randomized matrices A of the same dimension N).
This is based on the observation that randomized networks typically do not show Q=0, due to the fact the even a random network has at least one partitioning of its nodes that yields a Q value well above zero. To address this, we used a normalized measure Q m : where Q real is the Q value of the network, Q rand is defined as the average Q value of randomized networks, and Q max is defined as the maximal possible Q value of the network.
To compute the modularity measure of a network, Q m , we first calculated its Q real . To measure Q rand we used 10 4 random matrices. To estimate Q max we used an upper bound Q max =1-1/K where K is the maximal number of modules, which can be considered here as N (the dimension of the matrix A).

Genetic Algorithm
The main text employed Hill-climbing dynamic to solve the evolution over time. Here, we describe results using a different evolutionary dynamics, employing a standard genetic algorithm (GA) [5][6][7]

Speedup Calculation
The speedup (S) is defined as the ratio of the time, starting from random initial conditions, to reach the solution in a fixed goal problem (T FG ) to the time the system reaches a solution in a modularly varying goal problem (T MVG ) : Note that the solution in an MVG problem is not a point but rather a limit cycle whose amplitude depends on the epoch time (E) and on the difficulty of the fixed goals composing the problem (ε). Accordingly, T MVG was calculated as the time at which the system reaches its limit cycle. The amplitude of the limit cycle was also used to determine the threshold for the convergence of the system to its fixed goals (see Fig.   S5). T FG is defined as the time it takes the system to reach the same distance from the optimal solution.
In case of nearly modularly varying goals, the definition above was modified to take into account the parameter characterizing the distance between the goals η (see Fig. S6). The limit cycle of the system was calculated, and the minimal distance between the solution and the corresponding modular solution was found in each epoch. The threshold for calculation of T FG was then taken as the mean of the two distances: R m = Mean(R 1 ,R 2 ).

The block structure of the evolved matrix is determined by the correlation between the goal input and output.
Here we show the relation between the block structure of the evolved matrix (A) and the covariance structure of the goal input and output vectors (V and U). Here we will treat the inputs and outputs as independent random variables ( ) By definition the goal G is modular if there exists a block diagonal matrix M such that MV=U (up to permutations of the columns of V and U). In terms of the v's and u's this reads: Since each variable is sampled k times, the means are: Here we used the linearity of the mean operator. Similarly, the second moments are: The covariance matrix of the goals input and output pairs can now be calculated: . Substituting this into Eq. (S5), doing the summation, and solving for M ij we find: Writing the covariance matrix in terms of the correlation matrix This can be further simplified by calculating ( ) i u σ : That is the matrix M normalized so that the norm of the rows is 1, we finally obtain (S11) Note that the structures of M and that ij M are the same. Now based on appendix A (of the main text) we know that under MVG, the matrix A evolves towards the matrix M.
And so it must evolve towards the structure of the correlation matrix between the input and output pairs.

Calculation of the critical epoch time for speedup.
In this section, we comment on the range of goal switching times where speedup can be observed in MVG relative to fixed goal problems. Following Eq. (5) in the main text, the evolution of the components a ij t ( ) of the mapping a ij * matrix A within a single epoch can be described by the following equation: (S12) a ij t () = K n e −λ n t n N ∑ + a ij * Where a ij * is the optimal solution and the prefactors K n { } are determined by the eigenvectors corresponding to λ n { } and the initial conditions in that particular epoch.
Taking the case of N=2 we can write: Where λ 1 is the large eigenvalue and λ 2 is the small one. The evolution can thus be Formally, for rapid evolution to occur, the rate of the rapid decay must be larger than the rate of the slow decay in each epoch. Hence: Thus a general condition for speedup is that the epoch time will be less than some a epoch time: (S16) E c = 1 Where F ( ) is a function of the initial conditions (and the eigenvalues) upon which T c depends weakly.
For epoch times obeying T < T c the dynamics is governed by the rapid decay component only (see Fig. S4 b). Hence, as long as T < T c the evolution time is essentially the same. However, if the epoch time increases beyond T c the added time will be spent along the slow component, and total evolution time will increase with the added time being linearly proportional to T − T c .

Other cost functions
As mentioned in the main text, one can test various cost functions (the benefit function is maintained as the Frobenius norm since we measure the Euclidian distance from the goal). The cost function we presented in the main text has the advantage of a fully solvable linear dynamics. To study the dynamics governed by other cost we used two different methods: numerical analysis and simulations using genetic algorithms.
We find qualitatively similar results and dynamics to those found in the main text, under a wide range of convex cost functions. Concave cost functions can lead to scenarios in which the modular solution is in fact the optimal solution (either local or global). Thus it can be found in FG problems as well as in MVG problems. However, it is not guaranteed that the system will actually reach the modular solution since it may be flanked by other local optima. This is avoided in MVG problems. The nonconvex cost |a| is of special interest since it is characterized by a true fitness plateau.
To study a broader range of cost functions, we examined a family of cost functions such as: (S17) , Accordingly in the steady state solution In order to proceed we make use of the fact that Accordingly, the equilibrium solution averaged over all initial conditions is: Thus we finally get which is the optimal solution with nonvanishing cost.

Modularity declines if goals become constant:
What happens to modularity under a constant goal if one begins with a modular solution as an initial condition? We find that modularity decays over time (Fig 8b). Generally, this decay corresponds to motion along the low gradient valley towards the optimal, non-modular fixed point. , and so 0 Q 2 D → = also in this limit.
Furthermore, Q D=2 can be calculated analytically by using the solution to the equations of motions as given by eq. (5) and using the fact that the modular initial conditions puts the system inside the valley, which means that fast terms in A can be neglected (expressions that depends on the large eigenvalue λ). Substituting into Q D=2 we get:  where w 1 and w 2 are constants. The result is that Q D=2 vanishes exponentially slowly (~t E ε − ). Note that since Q m is a linear function of Q (see supporting information), these results still hold for Q m .
The conclusion is that if left alone, the modularity of the system decays slowly to the modularity of the optimal solution. Thus, goals need to keep varying over time in order to maintain the modular structure.