Model

We simulate the evolution of a haploid asexual population of N binary sequences. In an individual genome, each locus (site, nucleotide position, amino acid position) numbered i = 1,2,…,L is occupied by one of two alleles, either the wild-type allele, denoted a_i = 0, or the mutant allele, a_i = 1. We use a discrete generation scheme in the absence of generation overlap (Wright-Fisher model). The evolutionary factors included in the model are random mutation with rate μL per genome, constant directional selection, and random genetic drift due to random sampling of progeny. Selection includes an epistatic network with a set strength and topology. Recombination is absent. A previous modeling study shows that moderate levels of recombination can enhance epistatic detection [40]. We use the standard model of fitness landscape with pairwise interaction. The logarithm of the average progeny number of an individual genome, W, depends on sequence [a_i], as given by

where T_ij = 0 or 1 is the binary matrix that shows interacting pairs. Here the selection coefficients s_i and s_j denote the individual fitness costs of two deleterious mutations that are partially compensated by each other. By the definition, E_ij is the degree of compensation of deleterious alleles at sites i and j. Values E = 0 and 1 represent no epistasis and full compensation, respectively.

The formalism applies at negative values of E as well, but we focus on positive E, which case is termed "diminishing returns epistasis". The reason for this choice is strong effects of epistasis and strong indirect interactions in this region. The analytic derivation in the Methods applies at any E < 1, i.e., below the full compensation point, which covers most basic types of epistasis.

In our simulation example in Fig 1, we consider a haploid population with the initial frequency of deleterious alleles, f₀ = 0.45, beneficial mutation rate, U = 0.07, fixed selection coefficient, s = 0.1, population of N = 1000 individuals, L = 40 sites, and fixed epistatic strength E = 0.75. The core Monte-Carlo simulation code for Fig 1 is written in MATLAB and deposited at site https://github.com/rbatorsky/hiv-recombination. It can be modified for different types of epistasis and recombination. The code for Fig 2 (data analysis) has been uploaded to https://github.com/irouzine/Pedruzzi.