Multiple Loci Selection with Multi-Way Epistasis in Coalescence with Recombination

Abstract

As studies move into deeper characterization of the impact of selection through non-neutral mutations in whole genome population genetics, modeling for selection becomes crucial. Moreover, epistasis has long been recognized as a significant component in understanding the evolution of complex genetic systems. We present a backward coalescent model, EpiSimRA, that accommodates multiple loci selection, with multi-way (k-way) epistasis for any arbitrary k. Starting from arbitrary extant populations with epistatic sites, we trace the Ancestral Recombination Graph (ARG), sampling relevant recombination and coalescent events. Our framework allows for studying different complex evolutionary scenarios in the presence of selective sweeps, positive and negative selection with multiway epistasis. We also present a forward counterpart of the coalescent model based on a Wright-Fisher (WF) process, which we use as a validation framework, comparing the hallmarks of the ARG between the two. We provide the first framework that allows a nose-to-nose comparison of multiway epistasis in a coalescent simulator with its forward counterpart with respect to the hallmarks of the ARG. We demonstrate, through extensive experiments, that EpiSimRA is consistently superior in terms of performance (seconds vs. hours) in comparison to the forward model without compromising on its accuracy.

1. Introduction

EpiSimRA (both backward and forward) source, executable, user manuals are available at: https://github.com/ComputationalGenomics/SimRA (accessed on 24 April 2021). Nothing in Biology Makes Sense Except in the Light of Evolution [1], and simulating the evolution process, whether of multi-cellular humans, unicellular micro-organisms or even cancer-tumors, continues to be an important device in understanding the observed molecular profiles of populations. These profiles are an attribute of the genetic variability due to mutations and the change in frequency of alleles within populations over time. The selectively neutral infinite-sites model [2] is often used to analyze this variation [3]. Simulating random populations plays a significant role in investigating the effects of complex evolutionary processes on genetic diversity [4]. There are mainly two types of simulation algorithms: backward-time or coalescent and forward-time. The coalescent simulation [5] allows for fast approximation of the neutral Wright-Fisher (WF) model with natural selection, shaping patterns of variation in populations. The ARG [6] is a variant of Kingman’s coalescent and is used to reconstruct the most recent common ancestor (MRCA), starting backwards from the extant populations or leaves, using coalescent and recombination events. Once it finds the MRCA or if it involves all the trees, the grand most recent common ancestor (GMRCA), it traverses the ARG to the extant populations and introduce mutations or other genetic information in the genealogy. Forward-time simulators are more precise than their backward (coalescent) counterparts in modeling selection along with other factors as it starts from an initial population and tracks its evolution under the influence of various factors such as recombination, mutations, varying effective population size, fitness effects, and so forth. It progresses over multiple sequential generations, usually drawing random samples from the last generation to construct an ARG and its hallmarks. However, coalescent processes are much faster than forward-time simulation algorithms [7].

The classical approach for coalescent simulation as defined by Hudson in the seminal ms tool [3], applied the effects of recombination and coalescence to the ancestors of the samples going back in time in the extant population. This was later computed more efficiently in msprime [8], which used a new encoding for correlated trees resulting from simulations of the coalescent with recombination. Some approximations to the coalescent algorithms, which are fast, also exist, such as SMC [9], MaCS [10] and fastsimcoal [11]. Many programs have been developed to simulate scenarios not captured by ms such as selection [12,13,14,15], demographic inference [12,16] and admixture [17] among others. Coalescent models tracking genealogies in the presence of selection can also build an Ancestral Selection Graph (ASG), which is a branching-coalescing random graph within which the genealogy of a sample is embedded [18] conditional on the frequencies of the selected allele of the sample [19]. However, none of these methods take into account epistasis, which has long been recognized as a significant component in understanding genealogies and the evolution of complex genetic systems [20]. Here, we present the first coalescent simulator EpiSimRA which captures multiway epistasis, that is, allowing for interaction between alleles in multiple loci under selection. EpiSimRA tracks the ARG from randomly sampled extant populations and unlike ASG, is not conditional on allele frequencies. It constructs the genealogy dependent on the time to the closest recombination and coalescent event going backwards. Along with this, we also present an alternative, simple forward-time algorithm fwd-EpiSimRA, which efficiently simulates epistatic scenarios, to provide a validation framework to the coalescent simulator.

Forward simulators usually track the complete ancestral information, that is, studying all the lineages that survived until the last generation as a result of recombination events. Although forward simulations have been around for decades, only a few forward simulators exist to provide a framework to model multi-way epistasis, such as SELAM [21], allowing for pairwise epistatic selection to model the process and consequences of admixture or SLiM [22,23], which constructs ecologically realistic scenarios while accounting for a host of complex biological processes beyond the WF framework. Specifically, its functionality of tree-sequence recording draws parallels to fwd-EpiSimRA in a WF framework, providing support for epistatic interactions. However, none of these packages [21,23] can be used to compare with EpiSimRA as it is not possible to reconstruct the ARG from random extant samples. As fwd-EpiSimRA traces the ARG to obtain the MRCA of the random extant samples and record its hallmarks, we use it for a nose-to-nose comparison with the coalescent simulator as a validation framework. In the remainder of the paper we introduce the coalescent simulator and explain how it tracks multiway epistasis in the presence of recombination, followed by an overview of the forward simulator and the ARG tracking algorithm. Thereafter, we show the concordance between the coalescent and forward models for complex evolutionary scenarios and finally conclude by discussing multiway epistasis in simulating real world scenarios of admixture, cryptic relatedness and viral phylodynamics.

2. Materials and Methods

2.1. The Coalescent Simulator

The algorithm works back-in-time starting from the present (time 0), moving back into the past. Here we focus on how EpiSimRA is able to simulate multiway epistasis (the interested reader is referred to [17] for the neutral scenario). Let the number of loci under selection be l, possibly with multiway epistasis. As an illustration let l be 3 with selection values $s_1$, $s_2$ and $s_3$. The algorithm will assign three random locations on the genetic segment, unless the locations are explicitly specified and we assume that one of the alleles (either major or minor) is under selection while the other is neutral. The possible multiway epistasis are $e_{12}$, $e_{13}$, $e_{23}$ and $e_{123}$. If no value is specified then the epistasis is assumed to be neutral. Given this, we get $2^l$ possible types of lineages, each of them denoted as $l_z$. Let $l_0$ be the lineage type with no selection. For the example we ran, the other lineage types are $l_1$, $l_2$, $l_3$, $l_{12}$, $l_{13}$, $l_{23}$ and $l_{123}$. For two lineage types $z_a$ and $z_b$, let

$$l_{z_a}\prec l_{z_b}\mathrm{when} z_a\supset z_b.$$

For example, $l_{12}\prec l_1$ and $l_{12}\prec l_2$. Also, $l_{123}\prec l_{12}$. For the lineage type z, let $N_z$ be the effective population size.

2.1.1. Selection Scenarios

Effective population size is the reciprocal of the probability that two individuals will have the same parent—or that two chromosomes in the next generation will share the same parent chromosome. Fitness, in this case, would just be the ratio of the probability that the parent chromosome lineage with the allele will have an offspring chromosome to the probability that a parent chromosome lineage without the allele (neutral) will have an offspring chromosome [24]. The fitness, $1+s$ is thus the expected number of copies that a copy of the allele gives rise to in generation $t+1$, relative to the expected number that a neutral allele will give rise to. Thus, in infinitely large populations, the proportion of alleles under selection in a generation is related to the effective population size. Let $N_s$ be the partially effective population size with the allele under selection and $N_{\tilde{s}}$ ($=N-N_s$) is the partially effective population size with the reference or ancestral allele which is not under selection, giving:

$$\frac{N_s}{N_s+N_{\tilde{s}}}:\frac{N_{\tilde{s}}}{N_s+N_{\tilde{s}}}=1+s:1\Rightarrow N_s=\frac{1+s}{2+s}N=f_{s}N.$$

(1)

Thus $-1<s$, extendable to multiple loci with or without epistasis and the fitness defined as

$$f_s=\frac{1+s}{2+s}.$$

(2)

The fitness coefficient is a representative average of the allele frequency of the selected alleles in a generation, p. With the allele frequency, the effective population size with selection at a single locus would be $\mathbb{E}\left[N_s\right]=\frac{(1+s)pN}{(1+ps)}$ where the fitness $f_s=\frac{p(1+s)}{(1-p)}$. The effective population size is defined as $N_z=2N{f}^{\prime }$ where $f^{\prime }$ is the fitness for $l_z$-coalescence in the coalescent simulator and $N_z$ is the effective population size for $l_z$ lineage coalescence with alleles under selection.

The $f^{\prime }$ varies with neutral or epistatic scenarios for the loci under selection. For a neutral scenario with no selection $f^{\prime }$ is defined as follows

$$f^{\prime }=1-\sum _i{f}_{s_i}+\sum _{i,j}{f}_{s_{ij}}-f_{s_{\alpha \beta \gamma }},$$

(3)

where we remove the fitness effects of odd sites under selection and add the effects from even sites in a simulation scenario with three Single Nucleotide Polymorphisms or SNPs ($\alpha$, $\beta$ and $\gamma$) are considered to be under selection. Alternatively, for a single locus ($\alpha$) under selection with no epistasis in effect $f^{\prime }$ will be defined as follows with the signs reversed for even and odd sites under selection

$$f^{\prime }=f_{s_{\alpha }}-\sum _i{f}_{s_{\alpha i}}+\sum _{i,j}{f}_{s_{\alpha ij}}.$$

(4)

For two or multiple loci under selection there can be two cases with differing $f_s$. We define it as follows:

$$f_s= \left\{\begin{array}{cc}{f}_{e_{s_{\alpha \beta ..\omega }}}\hfill & \mathrm{with}\mathrm{epistasis}\hfill \\ \prod _i^{\omega }{f}_{s_i}\hfill & \mathrm{without}\mathrm{epistasis},\hfill \end{array}\right.$$

where $e_s$ is the user defined epistatic coefficient when epistasis is in effect across all $\omega$ sites. For a scenario with all three sites are under selection $f^{\prime }=f_{s_{\alpha \beta \gamma }}$.

2.1.2. EpiSimRA: Multiple Loci Selection & Multiway Epistasis

If $s_i$ and $s_j$ are two locations with the minimum (or derived) allele under selection at locus i and j respectively, then $e_{ij}$ denotes the epistasis between the two. The fitness coefficients in Equation (2) takes into account the $e_{ij}$ with respect to the s. If it is not explicitly specified then a neutral case (without epistasis) is assumed. The algorithm randomly chooses the location of the SNPs on the genetic segment being simulated.

We assume that no more than one event, coalescent or recombination, occurs at a generation and there is no back mutation, that is, a base undergoes no more than one mutation in the entire ARG. The mutation and recombination rates are uniform over the segment being simulated. If there is recombination, the lineages are randomly assigned but if $r=0$, the lineages are so assigned that no pair of types of lineages straddle (either they are disjoint or one is contained in the other). Lineage $l_0$ corresponds to lineage with no alleles under selection. For each lineage $l_z$, the algorithm only appends each node to a list when a recombination occurs with time $t>T_z$ where $T_z$ keeps track of the time to GMRCA. The recombination rate for $l_z$, $r_l^{\prime }$ is defined as

$$r_l^{\prime }=N_{z}grI.$$

For nodes which are not leaf nodes the length of the genetic material, s is proportional to the recombination rate as the rate is governed by the effective population size $N_z$. The stochastic nature of the method allows for a loop which pools lineages together at each iteration to find the event closest to the time $t_z$, over all lineages $l_z$. For each lineage $t_z=N_z\times t$ is computed where t is the time to next event using

$$\begin{array}{ccc}\hfill t& =& \min \left(\overbrace{\underset{1\le a<b\le L_z}{\min }\left(t_{ab}^{\mathrm{coal}}\right)},\underbrace{\underset{1\le i\le L_z}{\min }\left(t_i^{\mathrm{rcmb}}\right)}\right)\hfill \\ & =& \mathrm{Exp}\left(\overbrace{1+1+...+1}+\underbrace{r_1^{\prime }+r_2^{\prime }+..+r_{L_z}^{\prime }}\right)\hfill \\ & =& \mathrm{Exp}\left(\overbrace{L_z}+\underbrace{r_1^{\prime }+r_2^{\prime }+..+r_{L_z}^{\prime }}\right).\hfill \end{array}$$

(5)

$t^{\mathrm{coal}}$ is the time to coalescence and similarly $t^{\mathrm{rcmb}}$ is the time to the next recombination event. Equation (5) computes the closest event to this time (coalescent or recombination) where the overbraces capture the $\left(\genfrac{}{}{0pt}{}{L_z}{2}\right)$ coalescent events and the underbraces capture the $L_Z$ recombination events. When there is only one lineage in the pool, only recombination event can occur. Otherwise, we use the three properties outlined in the Appendix A to find the next event closest to the time computed in Equation (5). The time T, aggregated over t is thus the time to GMRCA as outlined in Algorithm 1. The event is a coalescence is chosen with the probability

$$\frac{\left(\genfrac{}{}{0pt}{}{L_z}{2}\right)}{\left(\genfrac{}{}{0pt}{}{L_z}{2}\right)+{\sum }_l{r}_l^{\prime }}$$

(6)

and recombination at lineage $1\le k\le L_z$ with probability,

$$\frac{r_k^{\prime }}{\left(\genfrac{}{}{0pt}{}{L}{2}\right)+{\sum }_l{r}_l^{\prime }}.$$

(7)

Equations (6) and (7) are used in a single draw of a random number such as in unit interval [0,1] broken up into $1+L_z$ sub intervals with cumulative ratio. The first interval implies a coalescent event and $k_{th}$ interval ($k>1$) implies a recombination at the lineage $l_{k-1}$. Since the events are randomly selected, t is estimated first and then the lineages are picked at random from $L_z$ active lineages. If the $T_z$ falls in the respective interval for a coalescent event, then the next event is coalescence and otherwise, if $r>0$, then the next event is recombination.

Algorithm 1 EpiSimRA

Coalescence Event

In a coalescence event $L_z$ is decremented by 1 as two random lineages of type $l_z$ are coalesced into one at time $T_z$ and the outgoing edge of the coalesced node is labeled by lineage $l_z$. If $\left|L\right|=1$, z is a singleton label (such as $s_1$ but not $s_1{s}_2$ or $s_1{s}_2{s}_3$), and, there exist no active lineage $l_z^{\prime }$ such that $z^{\prime }\prec z$, then the mutation(s) corresponding to lineage $l_z$ is assigned to this edge (using an approach in [17]) and the label of the outgoing edge of the new node is changed to $l_0$ and L is incremented by 1. Next, L is set to 0 and thus the lineage $l_z$ is made inactive.

Recombination Event

In a recombination event a lineage of type $l_z$ is randomly picked and a node v is created at $T_z$. The label of z is randomly split it into two lineage labels that is compatible with the location of the SNPs on the genetic segment I carried by the node v. Thereafter, L is incremented by 1.

The algorithm for EpiSimRA is described in Algorithm 1 and see Appendix A for an illustrative example of the algorithm for a three-way epistatic scenario.

2.2. The Forward Simulator

The model simulates evolution for a full population, forward in time with each generation containing N individuals with equal number of males and females, each carrying two chromosomes (see Appendix A for a detailed discussion and extension to selection on multiple loci). The complex evolutionary relationships between generations yields a number of mutations, recombinations, selected allele inheritance, linkage disequilibrium, and so forth, along the length of chromosome for each individual. These data are recorded in a data structure, which we call the “book of populations”, keeping a record of the past genealogy of the population. We trace the lineage of each site along the chromosome while tracing the ‘book’ and constructing the ARG. Inheritance follows the convention of a standard WF model applied to diploid organisms [3], with children randomly picking their parents weighted by the fitness coefficients when selection is in effect.

2.2.1. Simulating the “Book of Populations”

Each chromosome is represented by the alleles at each locus $l\in [1,g]$, which is randomly assigned initially. We use same notations as defined in

Table 1. Input parameters of the coalescent simulator.

	Parameters	Example Values	User-Specified Units	Units in bp for the Algorithm	Scaling Factor
g	seqment length	25; 75	Kb	$\times {10}^3$ bp	$\times {10}^3$
m	extant units	10; 20; 30; 40	-	-	× 1
N	population size	100; 200; 500; 1000	-	-	× 1
I	length of genetic material	1000	bp	1 bp	× 1
rates/generation
r	recombination rate	1	bp/gen $\times {10}^{-7}$	bp/gen	$\times {10}^{-7}$
$\mu$	SNP mutation rate	1.5	mut/bp/gen $\times {10}^{-8}$	×1 mut/bp/gen	$\times {10}^{-8}$
selection, epistasis parameters
$s_i$	fitness	0.3	-	× 1
$e_{ij}$	epistasis	0.1, 0.15	-	-	× 1

to describe fwd-EpiSimRA. The model assumes that each locus l has a fitness function $s_l\left(a\right)\in \mathbb{R}$, where a is an allele comprising the genotype. An individual i with allele $a_{il}$ at locus l is assigned a selection coefficient $s_{il}=s\left(a_{il}\right)$ which is user-defined, similar to EpiSimRA. The function $s(.)$ denotes the selective pressure and can be varied by intentional specification of recessive, dominant, additive, and other configurations, including homozygous advantage. This function encompasses selection at both single and multiple loci allowing flexible user-defined variations. When selection is not present, we set $s_{il}=0$.

For an individual i, the probability that it has children is given by

$$p_i=\frac{{\prod }_l\left(1+s_{il}\right)}{{\sum }_i{\prod }_l\left(1+s_{il}\right)}.$$

(8)

(See Appendix A for derivation). In each new generation, as in the WF model, the N children pick their parents with replacement according to the parent probabilities $p_i$. The simulation is run for $t=\{0,1,\dots ,G\}$ discrete generations with the $t=0$ being the base generation, outlined in Figure 1.

In each new generation, as in the WF model, the N children pick their parents with replacement according to the parent probabilities $p_i$. The simulation is run for $t=\{0,1,\dots ,G\}$ discrete generations with the $t=0$ being the base generation, outlined in Figure 1.

2.2.2. Modeling Multiway Epistasis

Multiway epistasis requires multiple interacting loci with similar selection effects. We assign selection coefficients to interacting sites for k-way epistasis, where k is the maximum number of interacting sites. Let there be q groups of loci, each containing at most k elements and we re-compute Equation (8) accounting for fitness related to interacting sites as,

$$p_i=\frac{{\prod }_q\left(1+S_{iq}\right)}{{\sum }_i{\prod }_q\left(1+S_{iq}\right)}.$$

(9)

If a group only has one element, that is if the selected locus is non-interacting, then we allow $S=s$, the user defined selection coefficient. For all other cases, we select S from a matrix or tensor of all possible allele combinations with respect to the number of interacting sites. S, the combined fitness coefficient is calculated by taking the fitness product of each interacting site as,

$$S=\prod _j\left(1+\sum _i{s}_i^{\left(j\right)}\right)+e_q.$$

(10)

$e_q$ is the epistatic interaction coefficient for each combination of interacting sites and $s_i^{\left(j\right)}$ is the selection coefficient at allele j in individual i’s chromosome.

2.2.3. Tracing the ARG

Detecting the past recombination events from extant sequences and specifying the place of each recombination is well studied [25,26,27]. The ARGs define a genealogical graph for all of the chromosomes in a population. Recent advances in population genetics simulators have resulted in tree-sequence recordings, which obtains the genealogical history of all genomes in a simulated population [28]. However, no natural ARG is recorded for the interacting loci with epistasis in effect and randomly sampling populations from extant generation, in forward simulators. It is traced from the “book of populations” from a number of extant haplotypes. We start from m randomly selected extant populations and trace the recombination and coalescent events back each generation. We keep track of each lineage corresponding to every site along the chromosome and stop when we have found a convergence for all lineages. This final coalescent event along the entire “book” is known as MRCA and we output the corresponding ARG.

3. Results

3.1. Comparison Study

Comparing the two models under selection calls for an assessment of the values. In both the models, common phenomena such as faster coalescence, decreasing diversity, decreasing number of recombination events occur when we study the individuals under selection. Hence, we compare the H, the height of the ARG or the time to MRCA, as it is a significant hallmark of the common history of a sample.

We run simulations for different parameter set-ups for the forward and backward model by running each experiment 100 times. We demonstrate the accuracy of the two algorithms by comparing H under different simulation scenarios allowing at most three interacting loci. The simple scenarios in this case is when there is no selection in effect, that is, the neutral coalescent model and selection at a single locus. We show that the two proposed models EpiSimRA and fwd-EpiSimRA show agreement in all of the different epistatic scenarios including selection in single locus.

The results for the complex scenario in this setting, accounting for epistasis with three loci, are shown in Figure 2, where we show the concordance for the forward and backward simulation with box-whisker plots, QQ-plots, CDF plots and PP plots (Figure A5). To obtain further validation we observed similar agreement in the Kolmogorov-Smirnov (KS) test on the distributions of H as returned by fwd-EpiSimRA and EpiSimRA for all scenarios. We found that for each, the null hypothesis that the two samples are drawn from the same distribution is never rejected and the test statistic is very small (

Table A1. K-S test statistics with corresponding p-values showing that the probability distributions of H as returned by fwd-sSimRA and back-sSimRA abstracts each other very closely.

3 Interacting Loci			${\mathit{e}}_{\mathit{s}}$	m	p-Value	Test Statistic
${\mathbf{s}}_{\mathbf{1}}$	${\mathbf{s}}_{\mathbf{2}}$	${\mathbf{s}}_{\mathbf{3}}$
				10	0.1400	0.16
				20	0.4431	0.12
×	×	×	×	30	0.3439	0.13
				40	0.9995	0.05
				10	0.6766	0.08
				20	0.7942	0.08
$s_1$	×	×	×	30	0.6766	0.10
				40	0.5750	0.11
				10	0.9921	0.06
				20	0.5560	0.11
$s_1$	$s_2$	×	×	30	0.7942	0.09
				40	0.8938	0.08
				10	0.8938	0.08
				20	0.9995	0.05
$s_1$	$s_2$	×	$0.1$	30	0.9710	0.06
				40	0.7942	0.09
				10	0.3439	0.13
				20	0.7942	0.08
$s_1$	$s_2$	$s_3$	×	30	0.6766	0.10
				40	0.5576	0.11
				10	0.9610	0.07
				20	0.9610	0.07
$s_1$	$s_2$	$s_3$	$0.1$	30	0.3556	0.13
				40	0.6766	0.10

).

3.2. Evaluating Epistatic Scenarios

We compare the H under selection in EpiSimRA and show how different scenarios impact the height of the ARG (Figure A1). We find that positive selection affects the time to coalescence inversely with more selective pressure results in less time to coalescence when simulated with $N=1000$ samples and genome length of $g=250$ kbp. Epistasis in two and three interacting loci results in lower time to MRCA (MRCA) than single locus selection and the neutral case amounting to higher selective pressure. We further studied effects of epistasis by simulating populations of $N=$ 10,000 with three-way epistasis. We find epistasis leads to a more complex evolutionary history resulting in longer time to coalescence in MRCA (Figure A2 in Appendix A). When epistasis is not in effect, that is, when $f_{s^{*}}={\prod }_i{f}_{s_i}$, when $i\in [1,3]$ results in a shorter TMRCA with simpler evolutionary history. In addition, an exhaustive comparison between the two simulators for all scenarios with or without epistasis is included in the Appendix A.

4. Discussion

Selection in a diploid organism impacts heritability two ways: (1) heterozygosity can mediate the impact of selection on the transmission of a selected haploid lineage, and (2) recombination can hide the impact selection may have had on the ARG. This model focuses more on the impact of the latter. Selection in a diploid heterozygous sample can boost, for one generation, the non-selected chromosome. This can complicate the impact of selection on lineages in the diploid forward model, but not the haploid. We expected the impact of boosted preference to be minimal along any given lineage since such a boost only occurs for dominant or additive alleles, and then for only one generation, with combinations in a population over time, this effect could be more significant.

The coalescent model and its forward counterpart under epistatic selection scenarios were concordant in the simulation studies. Epistasis makes the evolutionary history of extant populations more complex, but with selective pressure on certain alleles, the TMRCA appears to be shorter than the single locus selection or the neutral scenarios, respectively. It is expected that the H would be shorter for any selection scenario when compared to the neutral case, however, it is particularly intriguing to observe how the epistatic two and three-way scenarios have more selective pressure with a cumulative effect resulting in a decreased TMRCA. Although, the mean of fwd-EpiSimRA distribution (Figure 2) is a bit lower than EpiSimRA, we see concordance in the overall distribution, including the QQ and CDF plots and as well as in the KS tests. We posit that the difference in mean may be due to underlying differences in model assumptions such as diploid mechanisms for fwd-EpiSimRA in comparison to a haploid structure in EpiSimRA.

Computational complexity of EpiSimRA is directly proportional to N, the number of individuals per generation; g, length of the genome under simulation and k for k-way epistatic interactions. As we increase these parameters, we obtain a more complex evolutionary history leading to longer running time due to complex interactions between the inherited loci from one generation to the other, for randomly sampled extant populations. As we observe concordance in the observed TMRCA for the coalescent simulator EpiSimRA as well as in fwd-EpiSimRA when multiway epistasis is in effect, we obtain validation about the empirical correctness of the coalescent simulator. As in the coalescent simulator we cannot assume correctness until after the ARG has been established, we used the forward model to show the correctness of EpiSimRA, under varying values of selection coefficients and epistatic scenarios. The coalescent simulator, EpiSimRA is extremely fast in finding approximations to TMRCA, in comparison to fwd-EpiSimRA, as the latter has to build the entire “book of populations” and trace it. This leads to a difference in running time of hours for the forward model vs. seconds for its coalescent counterpart with varying input dimensions.

5. Conclusions

We present an algorithm that builds multi-locus selection and multiway epistasis into the backward coalescent model with recombinations, as well as in a forward scheme. Moreover, to the best of our knowledge, this is the first model that has taken a backward simulator with multiway epistasis and compared it nose-to-nose with its forward counterpart. Through extensive empirical comparison studies, albeit for small populations due to the time constraint of the forward model, we show that for complex scenarios with selection and epistasis (or even under neutral scenarios) the hallmark values by the backward and the forward schemes approximate each other. As the distributions of both the schemes are concordant, we conclude that either of the simulators (EpiSimRA or fwd-EpiSimRA) can be used to understand the effects of negative and positive selection, with multiway epistasis, along with selective sweeps across generations. Due to the lack of similar assumptions, parameters and hallmarks of ARGs returned, we did not compare EpiSimRA with present coalescent simulators for selection at a single locus. As fwd-EpiSimRA is based on the Wright-Fisher model and allows for epistatic interactions, we used it as a validation framework.

Multiway interaction across multiple loci leads to complex population genetic history but has a shorter height of the ARG relative to non-epistatic interactions. EpiSimRA encompasses all such scenarios with the potential for further exploration for viral phylodynamics with random sampling of a bacteria or virus populations. The time to coalescence when reconstructing its phylogeny under selection and epistasis allows us to study important epidemiological, immunological and evolutionary processes of viruses [29] such as the recent SARS-CoV-2 or similar Coronaviridae. This allows a validation framework for including selection and epistasis into standard population genetic models where we can now study the different scenarios when all the diploids associated with mutated sites along the chromosome with differing fitness values corresponding to the alleles.

6. Patents

There is no patent resulting from this work.