To explore the interplay between developmental and functional coupling on the evolution of brain structure, we devised a model that simulates the evolution of a population of individuals in an environment, where fitness is determined by the size and costs of brain components.
Definitions
We employ terminology from the evolutionary neurobiology literature, but the debate could also be understood in terms of principles derived from quantitative genetics. ‘Concerted evolution’ is the result of correlated evolution among all brain components, which can be caused by two main processes. The first, referred to as ‘developmental coupling’, occurs through genetic correlations that act through development, such that variation at a particular locus affects the development of multiple, or all, brain components. ‘Developmental decoupling’ refers to the breakdown, or absence, of these genetic correlations. The second, referred to as ‘functional coupling’, instead occurs through correlated selection pressures, which either arise due to components contributing to shared behavioural or computational functions, or the environment selecting on multiple aspects of brain structure (e.g. dark environments selecting for increased olfactory processing and against visual processing; in this case, the selection correlation is negative). Both correlated selection and genetic correlations can lead to coevolution between traits [35], and our model is designed to explore this in the context of brain evolution. However, we note it is likely applicable to other anatomies.
Base model components
The model is characterised by three components: the population of individuals, the environment and an evolutionary algorithm. We provide a simplified description of the model in Additional file 1: Figure S1.
Each individual is characterised by a number (N) of brain components, each component having a specific size (S_{i,j}), where i denotes the individual and j denotes the jth of N components. For example, among individuals with three brain components (N = 3), the brain of ‘Individual A’ could have a first component of size S_{A,1} = 5, a second component of size S_{A,2} = 7 and a third component of size S_{A,3} = 2, which we denote as S_{A} = (5, 7, 2). A second individual, ‘Individual B’, might have component brain sizes of 3, 3 and 1, respectively, denoted S_{B} = (3, 3, 1). In this case, the total size of Individual A’s brain would be 14, while the size of Individual B’s brain would be 7. Our model assumes that all brain components have the same fitness cost per unit size (see below), so it is total brain size that determines the evolutionary cost of an individual’s brain, while the size of individual components contributes differently to fitness benefits. Whether evolution is mosaic or concerted depends upon the factors influencing changes in brain component sizes over time. The sizes of brain components are allowed to vary, through mutation, and this variation is influenced by developmental coupling (D), which takes values between 0.0 (no developmental coupling, i.e. a fully ‘mosaic’ brain) and 1.0 (complete developmental coupling, i.e. a brain structure fully determined by total brain size). D is analogous to the strength of genetic correlations between components. For example, a D of 0.5 would indicate that 50% of the variation in each component is determined by variation in total brain size, and 50% of the variation in each brain component is independent of variation in both total brain size and other components. When a mutation event occurs, the program generates N + 1 random mutation factors (m_{j}), for example between 0.5 and 1.5 for a 50% mutation step size, where there is one factor for the whole brain (m_{0}) and one for each component (m_{1} to m_{N}), each brain component is then scaled by these mutation factors, with the variation in mutations affecting particular brain components being flattened depending on D’s value. For example, when D is 0, m_{0} for the total brain size mutation will be multiplied by 0 but other mutation factors will vary independently according to a scaling factor (1 – D, i.e. unscaled when D is 0), whereas when D is 1 all mutation factors for individual brain components will be multiplied by 0 and only the mutation factor for total brain size will persist. Models with intermediate values of D fall between these extremes. This is determined according to the following formula:
$$ {S}_{i,j(new)}=\left[\left(D\times {m}_0\right)+\left(\left(1D\right)\times {m}_j\right)\right]\times {S}_{i,j(old)};j\in \left\{1,\dots, N\right\} $$
(1)
The Environment is characterised by three factors. First, an environment is characterised by Functional Coupling (F), between 0.0 (no functional links between two brain components) and 1.0 (complete functional interdependence between two brain components), which determines how similar brain component benefits are to each other (for F = 1, all benefits are identical) — see below for implementation details. Second, a set of benefits (B_{j}), one for each of the N brain components, which represents the fitness contribution of each size unit an individual gains from that particular brain component, at a particular size. For example, in an environment with benefits B = (1, 2, 3), Individual A, from above, will have total fitness contributions from brain component sizes of 5 × 1 + 7 × 2 + 2 × 3 = 25, while Individual B will have total fitness benefits of 3 × 1 + 3 × 2 + 1 × 3 = 12. This is implemented by varying B_{max}, the highest possible benefit a unit of size could give, with the benefit per component size B_{j} given as a random fraction (uniformly sampled from the range 0 to 1) of this maximum benefit. Where F = 1, the generated fraction is the same for all components; if F = 0, then the benefit per component is independent (i.e. three generated fractions), and for any intermediate value of F, the benefit provided by each component is determined by contributions from both the crossbrain fraction and the percomponent fractions, scaled so that the sum of fractions per component never exceed 1. With higher F values, the individual component benefits are constrained to be more similar. The average benefit per component is always half the maximum benefit permitted, and the degree of functional coupling (F) determines how correlated they are. Third, the environment imposes a fitness cost (C) per size unit of the brain, which is uniform for all brain components based on the assumption that there is a linear ‘per neuron’ energetic cost [61], and that units of ‘size’ in the model are analogous to neuron number. More specifically, the units of ‘size’ in the model are analogous to the ratio between the neuron number of the current organism and the neuron number for the common ancestor, with the common ancestor starting with equal sized brain components, S_{0} = (1,1,1), in all our simulations. Total fitness for an individual i with brain component sizes S_{i,j} in an environment defined by benefits B_{j} and cost C is thus given by \( \sum \limits_j\left({S}_{i,j}\times {B}_j\right)\left(C\times \sum \limits_j{S}_{i,j}\right) \). For example, if the environment described above had a cost C = 1, Individual A will experience a fitness cost of 14, giving a total fitness of 25–14 = 11, while Individual B will experience a fitness cost of 7, giving a total fitness of 5. To visualise the effects of varying C and B_{max}, we measure the ratio of the average B across components (annotated B̅), to C.
The Evolutionary Process progresses through the following steps:

1.
Determine the number of ‘offspring’ for each individual in a population, and the age of all individuals, measured in the number of generations (these are identical across the population).

2.
Initialise an environment and a population of individuals with identical, uniform brain component sizes (S_{i,j} = 1).

3.
For a given number of simulation steps

a.
Generate offspring for all individuals

b.
Mutate all offspring

c.
Increase the age of all individuals

d.
Remove individuals whose age exceeds the maximum age

e.
Rank individuals according to their total fitness in the environment (calculated as described above)

f.
Remove the lowest ranking individuals until the population size returns to the origin size (i.e. population size is stable over time)
Given the above model parameters, we can examine the effects of developmental coupling (D) and functional coupling (F) on the evolution of brain component sizes. We can also explore the evolution of brain structure in populations of individuals with an intermediate value of D (unless specified otherwise, we use D = 0.5), which are neither fully concerted nor fully mosaic. We call these ‘partially mosaic’ individuals, where some of the mutations affect total brain size, scaling each component equally, and some affect each component independently, as given by Eq. (1) above. We can then assess which mechanism, for example, a fully mosaic brain (D = 0), a fully concerted brain (D = 1) or a partially mosaic brain (D = 0.5), is most successful in different scenarios by measuring the frequency of individuals in a population with that D value, as a proportion of total population size, after n generations. The model was initially implemented with no upper ceiling on overall brain size, in which case it most accurately simulates periods reflecting directional increases/decreases in brain components, and by extension brain size, which is a common but not universal trend [70, 71]. This means that under static selection regimes the base model can lead to continuous directional changes in total brain size. We also allow the environment to change randomly over the course of a simulation to examine how temporal heterogeneity in selection regimes affects the longterm success of alternative brain models.
Introducing constraints on total brain size
As described above, the initial model imposes no upper or lower limit on total brain size. While this will reflect periods of directional changes in the brain, or brain component, size [70, 71], the close correlation between brain and body size, which may evolve under contrasting selection pressures [35], may impose limitations on how brains respond to selection that is not captured in the base model. However, these upper and lower boundaries can be envisaged in terms of nonlinear relationships between the size of brain components and their costs. As brain size approaches an upper/lower ceiling, the relative cost of increasing/decreasing each component is likely increased, resulting in increased disparity of B̅/C ratios for each component when selection is favouring increases in particular brain regions. To compare how such boundaries may impact the probability of different outcomes for D values, we implemented an extension to the model in which C is scaled according to its distance from the brain size of the common ancestor, with costs becoming increasingly prohibitive near an upper or lower boundary, set as 1.5 and 0.5 times the starting combined size of all three components, respectively. This scaling factor was implemented as:
$$ {C}_i=C\times 0.25\times {\left(\frac{\sum_j{S}_{0,j}}{\sum_j{S}_{i,j}}+\frac{\sum_j{S}_{i,j}}{\sum_j{S}_{0,j}}\right)}^2 $$
(2)
where C is the environmentaldetermined base cost, ∑_{j}S_{i, j} is the total size of all brain components for the current individual and ∑_{j}S_{0, j} is the total brain size of the common ancestor. Under these conditions, the model aims to examine situations where variation in total brain size is under stabilising selection or some form of constraint.
Comparing effects of parameter variation on probabilities of mosaicism
Using the base model, we first conducted a series of simulations to explore four key questions identified in the introduction:

1.
What mechanisms can produce concerted evolution? Here, we fixed C to equal 1 and fix B_{max} to be 1, 2 or 4, giving a range of average B̅/C conditions, while varying D and F. We then ran simulations to examine the degree of mosaicism observed under high, moderate and low levels of F and D.

2.
Can both mechanisms be adaptive? Here, we repeat the comparisons above, but in competitive environments to examine the probability of obtaining coordinated changes in brain components under high, moderate and low levels of F.

3.
Do the costs of neural tissue select against concertedness when selection acts on specific brain components? i.e. How does variation in fitness contributions from different components affect the way brains respond to selection? This can be addressed in two ways: first, by varying F, which determines how correlated B values of each structure are, or by varying B_{max}, to alter the B̅/C ratio. Here, our aim was to test whether different levels of variation in the fitness contribution of additional brain tissue alter the probability of obtaining a mosaic or concerted brain.

4.
Is developmental integration evolutionary labile? i.e. in a fluctuating environment do strong developmental constraints evolve and/or collapse? In this comparison, we took a different approach. We introduced a starting population of brains with a range of component sizes and D. We then allowed C and B_{max} to vary randomly every 2 generations to test what combination of factors persist over time. Here, costs were sampled uniformly in the range 0.5–5, and B_{max} was sampled uniformly in the range 1–10. As a result, the average benefit (B̅) is in the range 0.5–5, which is the same range as the cost, but the actual B̅/C ratio will vary widely between generations. F was also sampled uniformly in the range 0–1. We subsequently explored how varying key life history traits (numbers of offspring and maximum age) might buffer the effects of random environmental fluctuations.

5.
How do constraints on brain size interact with the probability of mosaicism? Finally, we subsequently examined how imposing upper and lower bounds on total brain size impact the probability of obtaining a mosaic or concerted brain under each of the conditions described above.
In all case, the simulations were run over 100 generations, with 1000 iterations, a fixed population size of 300 individuals, initiated with 100 individuals per D value. In the main text, we present results of simulations with a mutation step size of 5%, but runs using the base model were also repeated with a larger mutation step size of 50% to examine how effects were influenced by mutation size (full results are presented in the Supplementary Information). To explore the early stages of the simulations, we also present results from a subset of simulations with 10 generations for these base comparisons. In experiments 2 and 3, simulations were run with an initial population containing equal numbers of individuals with different D values (D = 0, 0.5 or 1), which were then evolved under different environment conditions (determined by F and the ratio of benefits to cost). In experiment 4, environmental conditions were randomly varied every 2 generations for 150 generations. When imposing the upper and lower bounds on brain size, we repeated experiments 1–4, as described above, with a mutation step size of 5%. The ‘success’ of a D value was determined by the proportion of individuals in a population with that D value at the end of the simulation run. The full output of all models are summarised in Additional file 1: Figure S2S23. The frequency of D values was compared using generalised linear models and the glm() function in R [72] across batches of nine, 1000iteration simulations where F, D and B_{max} and/or C varied (experiments 1–3, see for example Additional file 1: Figure S2,S6), or where F, D, maximum lifespan and offspring number varied (experiment 4, see, for example Additional file 1: Figure S14) (all experiments total n = 27,000 iterations). We estimated the effects of all parameters and interactions where indicated, with a Gaussian distribution when comparing ‘degrees of mosaicism’ (defined as the natural log of the ratio between the largest brain component and the smallest brain component in each individual, averaged across the population) and a quasibinomial distribution when comparing proportional frequencies.
The code files in which the model is implemented are openly available for readers to implement additional parameter settings, or to extend the model, and can be accessed from github.com/shaharavin/BrainEvolutionSimulator [73]. Biplots were made using PlotsOfData [74]. Many simulations produce a bimodal distribution of frequencies, we therefore display the mean of these iterations solely to illustrate the skew in the outcome.