For the simplest case \(K=0\) , with no epistatic interaction between neighbouring loci, we can verify in Fig REF , that the average fitness for point mutations and inversions are equal \(\langle f_{P} \rangle _{K=0}=\langle f_{I} \rangle _{K=0}\simeq 0.667\) .
In this case, the landscape is smooth with only a single peak.
Hence, mutations that increase fitness are not hard to find and \(\langle f_{\nu } \rangle _{K=0}\)  is independent of mutation types.
This result also agrees with the Kauffman's (analytical) result \(\langle f \rangle _{K=0} = \frac{2}{3}\) , which was calculated using order statistics arguments [1]}.
Then, for \(K>0\) , the average fitness \(\langle f_{\nu }\rangle _{K}\)  increases with \(K\)  until reaching a maximum value of fitness.
In relation to this maximum, in Fig REF  we can identify the following four cases:
i) point mutations with adjacent epistatic interactions, \(\langle f_{P}\rangle _{K}=0.707\)  for \(K=2\) ;
ii) point mutations with random epistatic interactions, \(\langle f_{P}\rangle _{K}=0.722\)  for \(K=5\) ;
iii) inversions with adjacent epistatic interactions, \(\langle f_{I}\rangle _{K}=0.747\)  for \(K=4\) ; and
iv) inversions with random epistatic interactions \(\langle f_{I}\rangle _{K}=0.732\)  for \(K=4\) .
After these maximum, the fitness values decrease as \(K\)  increases.
This trend is also consistent with the seminal simulations carried out by Kauffman and Weinberger [1]}, [3]}.
For example, when \(K\rightarrow N-1\) , the mean fitness converge to the same value, regardless of the type of epistatic interaction neighbourhood.
For point mutations with random and adjacent epistatic interactions, we obtained \(\langle f_{P} \rangle _{K=96}\simeq 0.580\)  and this agrees very well with the Kauffman's numerical outcomes for \(K=96\) , c.f. [1]}.
It is worth mentioning that these trends—lower fitness being associated with increasing epistatic interactions—correspond to the well-known “complexity catastrophe” described by Kauffman [1]} (see also Refs.[6]}, [7]}).
These numerical outcomes confirm that our numerical set-up reproduces, with point mutations, what is known about \(\langle f_P \rangle _{K}\)  vs K in the NK model [1]}, [3]}.
Now, what's new is that for inversions, the average fitness values are higher than those for point mutations.
Indeed, for \(K\rightarrow N-1\) , the average fitness trend is very different from that of point mutations.
For example, the complexity catastrophe estimates that as \(K\)  increases, the expected fitness of the local maximum (for point mutations) decreases toward \(1/2\)  [1]}, which is indeed verified here.
But for inversions, the evolutionary process reaches higher expected fitness values \(\langle f_{I}\rangle _{\mathrm {K=99}}=0.610>\langle f_{P}\rangle _{\mathrm {K=99}}=0.579\) , for both random and adjacent epistatic interactions.
To generalize these results, we reproduced this experiment with other, more restrictive, definition of inversion mutations. More specifically, we tested inversions on linear chromosomes (with boundary conditions) and circular inversions with upper size limit \(s\le N\) , ranging from \(1\%\) —a single locus (i.e. inversions are like point mutations)—up to \(100\%\)  of chromosome size.
Simulations with linear chromosomes show no significant difference from our reference circular model (see supplementary S1 Text, section 2 for detailed results).
Simulation with an upper size limit show that the final fitness values increases as the size limit increases. However, the gain is maximal for small \(s\)  values (typically up to \(s = 16\) ) showing that small and mid-sized inversions are sufficient to reach high fitness peaks.
