Maximum mutational robustness in genotype-phenotype maps follows a self-similar blancmange-like curve

Phenotype robustness, defined as the average mutational robustness of all the genotypes that map to a given phenotype, plays a key role in facilitating neutral exploration of novel phenotypic variation by an evolving population. By applying results from coding theory, we prove that the maximum phenotype robustness occurs when genotypes are organised as bricklayer’s graphs, so called because they resemble the way in which a bricklayer would fill in a Hamming graph. The value of the maximal robustness is given by a fractal continuous everywhere but differentiable nowhere sums-of-digits function from number theory. Interestingly, genotype-phenotype (GP) maps for RNA secondary structure and the HP model for protein folding can exhibit phenotype robustness that exactly attains this upper bound. By exploiting properties of the sums-of-digits function, we prove a lower bound on the deviation of the maximum robustness of phenotypes with multiple neutral components from the bricklayer’s graph bound, and show that RNA secondary structure phenotypes obey this bound. Finally, we show how robustness changes when phenotypes are coarse-grained and derive a formula and associated bounds for the transition probabilities between such phenotypes.

Methods
This data contains the data and code required to generate presented in "Maximum Mutational Robustness in Genotype-Phenotype Maps Follows a Self-similar Blancmange-like Curve" by Mohanty et al., published in Journal of the Royal Society Interface.
The exact maximum robustness curve corresponding to the robustness of the bricklayer's graphs, as well as the interpolated curve, can be generated from the RoBound Calculator, available free of charge and open source on GitHub (https://github.com/vaibhav-mohanty/RoBound-Calculator).
All bounds (e.g. Figure 1, 9, 10, and 11 as well as the bounds shown in Figure 3 or 7) can be calculated using the RoBound Calculator (https://github.com/vaibhav-mohanty/RoBound-Calculator).
In Figure 3, the RNA and HP model neutral component sizes and robustness values are provided in the files hp5x5_components.csv, hp24_components.csv, rna12_components.csv, and rna15_components.csv. These results obtainedd from Greenbury et al., "Genetic Correlations Greatly Increase Mutational Robustness and Can Both Reduce and Enhance Evolvability," PLOS Computational Biology, 2016.
In Figure 4, we show the deviation of RNA12 and RNA15 neutral networks from the maximum (bricklayer's) robustness. This data is obtained from rna12.csv (alternatively rna12.mat) and rna15.csv (alternatively rna15.mat), which were also results obtained from Greenbury et al., "Genetic Correlations Greatly Increase Mutational Robustness and Can Both Reduce and Enhance Evolvability," PLOS Computational Biology, 2016. The code to produce the figure is found in nc_err_corr.m
In Figure 7, the bounds can be calculated using the RoBound Calculator (https://github.com/vaibhav-mohanty/RoBound-Calculator). The raw data is obtained from using ViennaRNA (https://www.tbi.univie.ac.at/RNA/) to calculate the dot-bracket structures. These structures are then fed into the RNASHAPES software (https://bibiserv.cebitec.uni-bielefeld.de/rnashapes) to generate the coarse-grained data. The data files rna12abstract.mat and rna15abstract.mat contain the frequency and robustness values for the neutral networks at various levels of coarse-graining.
In Figure 8, the RNA12 transition probabilites between phenotypes can be obtained from rna12_theta.csv, which provides the phi_pq matrix when each column's sum is normalized to 1. The script phi_critical_ranges.m produces Figure 8.
Figures 2, 5, and 6 are schematics and have no associated data.