BACKGROUND
the field of population genetics strives to determine, how elementary evolutionary forces  shape the genetic landscape within population of a species.
as for related areas in evolutionary biology, rapid advances in next-generation sequencing technology  have transformed the field with the completion of the first phase of the  <dig> genome project
<cit>  representing a recent highlight.
thus, population genetics is rapidly transforming into a data-driven science
<cit> , while many methods and tools are not up to the challenge yet
<cit>  with respect to scalability and efficiency.
because of the complex processes in population genetics, the inference machinery for analyzing the properties of a population is often limited and needs to be complemented by simulations.
some of the most common use cases of simulation software in population genetics are  the verification of novel methods/models
<cit>  and  the generation of datasets for assessing how using a different model affects the properties of of a population.
a promising, but computationally particularly expensive approach is to use simulations for inferring specific properties  of a given data sample, as implemented, for instance by using approximate bayesian computation techniques
<cit> .
for an overview of available software and applications, see
<cit> .
forward-in-time simulation represents a powerful approach to simulating the evolutionary processes that act on genomic regions .
essentially, sequences are represented in silico and the basic evolutionary operations/events  are explicitly applied to each simulated individual, on a generation-by-generation basis.
the inferred, simulated sequences represent the exact outcome of the underlying/assumed evolutionary process.
however, this high accuracy comes at a high computational cost.
execution times are quadratic with respect to the number of individuals, because the number of simulated generations needs to be increased with growing population size.
to obtain accurate results, the size of the population that needs to be simulated must be chosen realistically.
however, realistic population sizes range between a computationally already challenging  <dig>  individuals in humans to an entirely prohibitive number of  <dig> , <dig> individuals in some species of drosophila melanogaster or  <dig> , <dig> in escherichia coli <cit> .
these excessive computational requirements explain the popularity of approximate coalescent  simulations
<cit>  that are substantially faster.
the coalescent-based approach is faster because only genomic samples that survived until the present are simulated backward in time.
a substantial drawback of coalescent simulations is that, natural selection can only be incorporated into the model by means of a single mutation
<cit> .
to date, only few attempts have been made to improve the performance of forward simulators using optimization techniques from computer science and via low-level technical improvements.
chadeau-hyam et al.
developed a strategy to re-scale simulation parameters such as to decrease the effective population size that needs to be simulated
<cit> .
however, this shortcut induces a decrease in accuracy .
forwsim <cit>   implements a forward-algorithm that regularly simulates a small user-defined set of mutations under selection, forward in time.
the neutral part of each haplotype is simulated separately and executed with a delay of  <dig> generations.
thereby, one can circumvent simulating haplotypes, if no ancestral material has survived until the present.
distinguishing between neutral and non-neutral mutations for accelerating forward simulations represents a natural implementation choice, specifically since it is assumed that the majority of mutations is effectively neutral
<cit> .
finally, sfs_code <cit>  has become a state-of-the-art forward simulator, because of an efficient implementation in c coupled with a plethora of simulation features and parameters.
here, we describe design as well as optimization techniques that allow for substantially more efficient forward-in-time simulations.
we introduce a novel algorithm, for the a posteriori simulation of neutral mutations.
in addition to runtime gains, our algorithm generates a data structure that depicts the entire observable history of each individual sequence, a feature that –to the best of our knowledge– is unique among simulation software.
we make available a corresponding software package, the ancestry-aware forward-in-time simulator , a c++ implementation that can handle challenging post-genomic datasets for population genetics analyses.
methods
in the following, we describe our algorithm for ancestry-based simulation of neutral mutations.
subsequently, we discuss implementation and optimization issues.
algorithm for simulating neutral mutations
for forward simulation of a single generation of a population with effective population size ne, a scaled mutation rate of θ=4neμ , the scaled recombination rate ρ=4ner , the basic simulation steps are:  sample diploid individuals by their fitness,  determine the number of recombination events ) and create recombinants,  determine the number of non-neutral mutations ) and mutate sequences,  recalculate the fitness values for each individual.
since neutral mutations, by definition, do not affect the survival probability of an individual , we do not need to simulate neutral mutations forward in time.
this observation was first explored in forwsim <cit> , where the simulation of the neutral part of a sequence was delayed for a small number of generations.
thereby, the actual calculations required for this part of the simulation can be discarded, if the haplotype under consideration did not survive  until the present.
the authors demonstrated that, the simulation of neutral sequences can be accelerated by a factor of two to five by optimally choosing this so-called look-ahead parameter, that is, the number of generations, by which the simulation of neutral mutations is delayed.
in ana-fits, we deploy a different strategy: neutral mutations are not simulated forward in time at all.
instead, we keep track of the entire ancestry  of all surviving individuals , create a graph structure , and only extract the neutral part of all sequences from this graph structure once the forward simulation has been completed.
in other words, we retain all information accumulated during the forward simulation phase, exactly determine all parts of sequences, where neutral mutations could have occurred and finally insert neutral mutations into the sequence a posteriori.
note that, in analogy to forwsim and in contrast to coalescent simulations, our algorithm is equivalent to simulating all neutral mutations forward in time.
thereby, we can guarantee simulation accuracy, which is one of the key advantages of forward simulations.
simulating surviving neutral mutations
during l forward-simulation generations , ana-fits annotates an ancestry
al={a <dig> a <dig> …al}, where ai is a surjective function that maps the individual haplotypes in generation i to the haplotypes of the preceding generation i− <dig>  in the following, we outline this procedure for a single locus/chromosome.
the extension to multiple loci is described in the implementation section.
the variable ai is defined as ai:hi,j→hi− <dig> l, where hi,j is the haplotype with index j in generation i, that is mapped to its parental haplotype with index l in generation i− <dig>  haplotype indices are ordered, such that adjacent pairs that start with an even index represent the two homologous versions of each locus for diploid organisms  or  represent an individual).
as mentioned before, the algorithm also keeps track of all recombination events  that occurred during the forward phase.
for each event, we store the absolute sequence position, the haplotype index, and the generation number.
for a sequence of length l, the information in
al and the corresponding set of recombination events can be used to determine  the regions of ancestral haplotypes that survived until the present.
only mutations that occurred within these regions are observable in the present.
in the first step, our algorithm determines haplotypes in ancestral generations  that contributed material to the present generation .
for these regions, we simulate neutral mutations and identify those recombination events that contributed to mutations in the present haplotype instances.
the algorithm starts at the present generation and initializes surviving regions for each haplotype with a maximum interval of , where l is the sequence length.
regions for ancestral haplotypes are initialized by the tuple .
note that, the surviving regions of a haplotype instance can be fragmented/scattered into many segments by recombinations.
to economize on runtime and memory, we only keep track of the start position of the very first segment rs and the end position of the very last segment re.
by using this approximate surviving ’super-’region r=, we overestimate the size of the actual surviving region.
we correct for this in a later stage of the algorithm, when sequences are extracted from the graph structure.
finally, to determine the surviving haplotype region in generation i− <dig>  we iterate over all survivors
si of generation i starting with haplotypes contained in the present generation.
for each surviving haplotype
hi,l∈si with region , we propagate the surviving regions as described in the following algorithm to the members of the preceding generation.
let hi− <dig> k with region
be the parent haplotype and hi− <dig> k+ <dig> with region
be the homologous haplotype in the corresponding parent individual:
for simplicity, the algorithm described above assumes that the recombinant emerged from at most one recombination event.
for multiple recombination events per recombinant, the algorithm has to be adapted, such that surviving regions of parental haplotypes are extended until the next recombination breakpoint.
this has also been implemented in the ana-fits software.
during the procedure described above, we also determine the set of survivors
si− <dig> from the preceding generation.
in the case that a recombination event does not split the surviving region , we can simply ignore it, if the following property does not hold: assume, haplotype hi,l is mapped to hi− <dig> k with surviving region  and hi,l is a recombinant with a breakpoint b<rs.
this means that, hi,l does not inherit mutations from hi− <dig> k, but instead from the homologous haplotype hi− <dig> k+ <dig> it recombined with.
if we can not ignore the recombination, we create a node instance for the graph that will be constructed in the next step.
the survivors of each generation are stored for further processing in subsequent steps.
once survivors and surviving regions for a generation i have been determined, we can then simulate the neutral mutations that occurred in this generation and that are still observable in the present generation.
thus, the expected number of neutral mutations is reduced by a factor of
f=|si|2·ne.
furthermore, if a neutral mutation occurs outside the estimated surviving region, it can simply be ignored.
as shown in figure
<dig>  depending on the recombination rate r, the number of haplotypes that contribute genetic material to the present generation quickly converges and may cause values of f to become small.
note that, we expect the algorithm to become inefficient for exceptionally high recombination rates, such as r:= <dig> .
graph construction
in the algorithm step described above, we determined survivors
si for each generation i, the set of relevant recombinations , and the set of relevant neutral mutations .
for all nodes, we store the generation of origin and the corresponding haplotype index hi,l.
the goal of the procedure we describe here, is to construct a graph .
the beg is a directed acyclic graph with two types of nodes:  m-nodes  with out-degree  <dig> and  r-nodes  with out-degree  <dig>  thus, each node in the graph represents the state of an ancestral haplotype that reflects the history of all mutation or recombination events experienced by its ancestors.
a directed edge connecting node n <dig> to node n <dig> indicates that n <dig> emerged from n <dig> after mutation or recombination.
given this information, the graph can be constructed forward in time  using an array ai for the current generation and an array ai− <dig> for the preceding generation i− <dig> that contains node references/pointers.
array ai is populated by carrying over/propagating references from ai− <dig> according to the survival information stored in
si− <dig>  if a surviving haplotype of generation i underwent recombination or mutation, the respective r- or m-node nnew is used to replace the node nold that was propagated and an edge is created between nnew and nold.
we also create an additional edge to the haplotype which is homologous to nold for r-nodes .
figure
<dig> provides a complex example, where a recombinant emerges from three recombination events that originated at two homologous ancestral haplotypes  <dig> and  <dig>
extraction of sequences
the algorithm described above constructs the beg and an array of haplotype states ni , that survived until the present.
in the final step, we extract the exact sequence for each node ni from the beg by recursively traversing the graph starting at ni, while keeping track of the region of interest which is segmented by each recombination event which will be encountered.
since we execute the backtracking procedure for each node ni individually,  the exact borders of all segments of the surviving region for this node are known at any given point in time.
thus, we can ignore all mutation events that did not occur within the surviving region, but were inserted into the graph because of overestimation or because this mutation forms part of the surviving region of another haplotype.
in practice, the above extraction procedure is inefficient.
as a consequence, the sequence reconstruction step dominates execution times.
consider two haplotype states n <dig> and n <dig> that have survived, where n <dig> is an m-node with ancestor n <dig>  for both nodes, we have to traverse all ancestral nodes of n <dig> and n <dig> to obtain the two corresponding sequences s <dig> and s <dig>  to accelerate this process, it suffices to create sequence s <dig> via backtracing starting at n <dig> and to generate s <dig> as a copy of s <dig> and then simply add the mutation events induced by node n <dig> in the end.
in the following, we explain, how the computational cost of sequence extraction can be reduced.
prior to the actual sequence extraction step, the graph is traversed twice, starting with each haplotype ni that survived until the present.
in the first traversal, we approximate the surviving region for each node as described in the algorithm.
thus, if we decide to explicitly determine and store the sequence of a node, we already know the maximum size of the region of interest beforehand.
in a second traversal that starts with each haplotype in the present, we determine for which nodes intermediate sequences should be created.
in other words, the surviving mutations are represented explicitly in this node.
during these traversals, we actually stop a traversal, if the traversed node has already been visited by another traversal.
for each node, we also store the number v of times the node has been visited.
we refer to nodes with v≥ <dig> as coalescent nodes, since ancestral material  of two lineages coalesces in these nodes.
finally, we determine the distribution of v and the 5% nodes with highest v are represented explicitly, while all remaining nodes will be traversed several times during sequence extraction.
then, the sequences can be extracted via a final full traversal.
the 5% cut-off  was determined empirically.
it yields good results with respect to balancing runtime versus memory requirements.
note that, if every coalescent node  is represented explicitly, memory requirements would become prohibitive for whole-chromosome simulations with  <dig>  individuals.
in general, the 5% cut-off proves to be sufficient to attain a substantial decrease of the traversal cost for sequence extraction.
implementation and optimization
the ancestry-based algorithm for simulating neutral mutations as discussed in the previous section has specifically been developed to improve runtimes of forward-simulations.
in the following, we discuss low-level implementation issues and optimization techniques of ana-fits.
memory requirements of ancestry and beg
for the majority of possible ana-fits invocations with respect to user parameter combinations, storing the ancestry
al requires the largest amount of memory.
depending on the effective population size ne, we need to store ancestry information for 10·ne <dig> haplotypes.
if this exceeds the amount of available memory, the user can set a command line flag such that the memory utilization of the ancestry does not exceed a user-specified soft limit .
ana-fits then tries to split up the forward simulation into sections  and updates the graph after each section leading to increased runtimes.
in some instances the graph representation itself can become at least as memory intensive as the ancestry.
furthermore, for large graphs, the explicit bit vector representation of intermediate sequences may further increase memory requirements.
therefore, currently no hard upper memory utilization limit can be imposed on ana-fits runs.
the ancestry
al only works correctly for a single locus.
for multiple unlinked loci , we map each individual to its two parent individuals instead of mapping it to individual haplotypes.
for this, we use a per-chromosome bit vector to keep track of which of the two possible homologous haplotypes was inherited from the respective parent.
apart from their high memory-efficiency, bit vectors can also be efficiently initialized with random bits using a random number generator .
to avoid unnecessary and inefficient frequent memory allocations, the entire ancestry  and the array of survivors are allocated as single, monolithic blocks.
for the ancestry, we dynamically determine the minimal number of bytes that is required for storing a generation as a function of the population size of the preceding generation .
since, under most simulation scenarios, ne does not exceed  <dig>  individuals, it is mostly sufficient to use the unsigned short integer  datatype for saving memory.
note that, beside all efforts to keep the memory consumption of ana-fits as low as possible, the entire ancestry  is kept in main memory during simulation.
this necessarily means that ana-fits exhibits higher memory requirements than comparable programs.
random number generation
high quality random number generators are essential for simulation tools.
to date, the mersenne twister algorithm
<cit>  is considered as state-of-the-art in pseudo-random number generation .
a simd-based version of the mersenne twister
<cit>  employs streaming single instruction, multiple data extensions  instructions to generate random numbers at almost twice the speed as the original mersenne twister.
ana-fits comes with a distribution of randomlib
<cit> , a c++ implementation of sfmt.
for efficiency, randomlib initializes an array with random bytes and then transforms the bytes of this array  into the required primitive data type  for the requested range.
for forward simulations, we can exploit two application-specific properties to more efficiently use prngs:  usually for the ancestry, integer values of less than  <dig> bytes are required, and  uniform random numbers are frequently drawn from the same interval , where n is the effective population size of the previous generation).
when a large amount of random numbers with these two properties  is required, we can directly initialize a target array from the internal array of the prng and transform the minimal number of random bytes for the respective integer type into uniform random numbers for the required range using sse instructions.
for instance, if simulations require sampling  <dig> integers ∈[  <dig> ), we can copy  <dig> bytes  from the prng’s internal array into a location of the ancestry memory block and use vector instructions to transform the  <dig> bytes into  <dig> integers of the specified range.
without this optimization, randomlib would produce  <dig> bytes of randomness and transform  <dig> bytes into one integer separately.
non-neutral sequence representation
for forward simulators that work with a given number and pre-defined locations of polymorphic sites , it is straight-forward to represent haplotypes as bit arrays.
if we want to simulate a sequence with a huge number of base-pairs , bit arrays are not the data structure of choice any more because they become too sparse.
thus, in other forward simulator implementations, a haplotype is stored as a sequence of mutations with respect to an initial reference sequence.
for instance, sfs_code <cit>  uses self-balancing binary trees to represent the polymorphisms of a haplotype.
sorted binary trees allow for rapid insertion of novel mutations and fast extraction of sub-sequences which are required for creating recombinants.
an alternative data structure is used in forwsim that maintains sorted arrays of sequence positions.
since for ana-fits, we want to store additional information about each mutation , we use a data structure that is based on references/pointers to mutation object instances.
when sorted arrays are used, sub-sequence replication for creating recombinants can be efficiently implemented using the memcpy() system call.
this is because finding coordinates within the sequence via a binary search has the same time complexity of o) as searches in a sorted tree.
additionally, changing the array size induces a smaller memory allocation overhead than allocating additional nodes in sorted trees.
the main drawback of sorted arrays is the insertion cost of
o. in practice however, obtaining the correct relative position via a binary search is usually more expensive than executing an insertion as long as insertions are implemented via the memcpy() system call.
note that, the memcpy() system call is also vectorized via sse intrinsics.
finally, the number of insertion operations is relatively small in ana-fits, since we only simulate non-neutral mutations forward in time.
we experimented with various sorted tree implementations, for instance, the c++ stl std::set implementation or judy arrays by
<cit> , which are based on cache-optimized 256-ary trees.
for ana-fits, we found that the tree data structures did not outperform our sorted array implementation.
memory management
in forward simulations, the number of distinct haplotype instances containing the sorted array of mutation references and further meta-information and the number of mutation instances changes dynamically at a rapid pace.
thus, dynamic memory allocations  and free()) that are invoked in the inner parts of nested loops quickly dominate runtimes.
as mentioned in the previous section, ana-fits uses block allocation, instead of dynamic allocation whenever possible.
apart from being faster this also prevents heap fragmentation.
one approach to handle the rapid turnover of haplotype and mutation instances more efficiently is to use free lists in conjunction with reference counting.
the underlying idea is to keep track for each allocated instance, how often it is referenced in the current generation  or in present haplotypes .
instead of deallocating unreferenced instances, instances are appended to a list, the free list, from which new instances can then be allocated again, when required.
however, this allocation strategy still uses more than 10% of total runtime.
instead, we use two task-specific allocation schemes that keep memory allocation overhead to an absolute minimum.
for haplotype instances, we deploy a free ring scheme .
as the name suggests, a free ring contains references to allocated haplotype instances in a similar way as a free list.
whenever a haplotype is carried over into the next generation, we update its generation number .
thus, when haplotype instance memory is requested from the free ring structure, an iterator traverses the elements in the free ring and returns the first instance that is neither used in the current nor the previous generation.
if no free instance is available, the free ring capacity is doubled by allocating new haplotype instances.
since it is unlikely that all references of a haplotype will be carried over to the next generation in the sampling phase without being subject to recombination or mutation, the iterator typically finds an unused instance after a small number of search steps.
in contrast to this, mutation instances need to be allocated in blocks, since references to instances  represent an inefficient mechanism for this task, because mutation instances are comparably small in size.
moreover, for large numbers of mutations, a corresponding large number of memory allocation invocations would be required, with a negative impact on performance.
allocating monolithic blocks with dynamic re-sizing is also not possible, because addresses of mutation instances can change which would invalidate the references to mutations that are stored in the haplotype instances.
instead, we allocate multiple blocks  of instances and implement an iterator that searches for an unused instance, such that it jumps to the next block when the end of the current block has been reached.
when a free instance has been found, it is flagged as in use.
this flag is updated in a cleanup step that is applied after a fixed number of generations.
in this cleanup step, mutation instances that occur in all haplotypes  are removed from all haplotypes.
thus, this phase can be used to efficiently update the in use flag of each mutation without inducing a large overhead.
since, under this scheme, the addresses of instances do not change, we refer to this allocation structure as address-preserving free ring .
implementation of neutral sequence extraction
the main challenge in implementing the beg algorithm is to design an efficient method for extracting sequences from the graph.
for this task, we can use the following property: during the sequence extraction phase, all neutral mutations that occurred during simulation are known.
thus, bit arrays are the most adequate data structure for storing neutral sequences.
after sorting mutations by their location, each mutation is assigned its corresponding index in the bit array.
as a result, the cost of adding a mutation that is encountered during sequence extraction to a given sequence is reduced to merely setting a bit.
handling recombinations is computationally more expensive, since we have to determine the bit array indices of the start and end points for the surviving regions of the parent haplotype via a binary search.
carrying over the relevant mutations from parent haplotypes is implemented as a bit-wise or-operation.
this implementation is particularly efficient because the sequence output of the simulation can directly be written as bit array to a binary file.
thus, output files do not require a lot of disk space and programs that directly evaluate the output of a large simulation can efficiently parse and post-process this output data.
RESULTS
in this section, we evaluate the runtime performance of ana-fits.
in general, conducting a fair comparison between forward simulation tools represents a difficult task, since each simulator uses different assumptions .
for a meaningful technical analysis, we separately evaluate the two key aspects of ana-fits: efficient implementation of non-neutral forward simulation and ancestry-based simulation of neutral mutations.
for realistic use cases, both features independently contribute a certain amount of run time improvement.
for all runs discussed throughout this section, we set mutation and recombination frequencies to μ per base=rper base= <dig> ·10− <dig> which roughly correspond to the empirical rates for humans
<cit> .
for fixed μ and r, the runtime depends on the locus length l and the effective population size ne.
all runs were executed on a 4× <dig> core machine  with a total of  <dig> gb ram.
the execution times provided here are the result of  <dig> independent runs with different random number seeds.
we discard the slowest runtime in order to correct for the impact of seed values that generate exceptionally expensive simulations and average over the remaining three.
moreover, we validate ana-fits by comparing the output of our software to datasets generated with established simulation codes.
finally, we describe runtime gains for the simulation of partially neutral sequences, specifically in a whole-genome setting for human sequences.
runtime improvement: forward simulation only
here, we evaluate the efficiency of the pure forward simulation part of our implementation and compare it to sfs_code.
ana-fits can simulate each mutation forward in time by choosing a constant small selection coefficient s≠ <dig> for each mutation that is to be simulated.
if s is small enough, it is effectively neutral.
this way the beg algorithm is not applied.
simulating close to neutral mutations , forward in time, represents the most expensive case, since for large values of s the mutation quickly undergoes selection and for small negative values it is removed from the gene pool by selection.
in both cases the mutation quickly vanishes which in turn reduces runtimes.
in runs with sfs_code, we simulate neutral mutations forward in time.
in general, it is not possible to efficiently simulate one contiguous locus of the length of a chromosome with sfs_code .
instead, for simulating a locus of  <dig> , <dig> base pairs , the user has to break up the locus into a number of linked loci .
for obtaining sfs_code runtimes, we did not optimize the break-up strategy exhaustively, but followed the general suggestions by the author and simulated chromosome lengths listed in figure
<dig> as 50· <dig> kbp, 50· <dig> kbp and 500· <dig> kbp.
figure
<dig> depicts the relative runtimes of both programs.
overall, ana-fits is between  <dig>  and  <dig>  times faster than sfs_code.
in all except  <dig> runs, ana-fits is at least one order of magnitude faster than sfs_code.
the runtime improvements do not show a clear dependency on one of the two input dimensions .
we assume that this is because, for certain steps, sfs_code switches between alternative implementations depending on the number of recombination/mutation events that are expected per generation.
further parameters that may influence the runtime differences between the codes are:  frequency of the cleanup steps in either code  and  number of linked loci the sequence was split into .
runtime improvement: ancestry-based simulation
as already mentioned, the beg algorithm is inspired by the postponed simulation of neutral mutations as implemented in forwsim.
thus, in order to evaluate the stand-alone efficiency of the beg algorithm, we compare the runtime for simulating entirely neutral sequences between ana-fits and forwsim.
note that, forwsim simplifies the fisher-wright model in two ways:  while it assumes finite sites, at most one mutation per site is allowed and  haplotypes undergo at most one recombination event per generation.
we omit a comparison between forwsim and ana-fits that includes forward simulation of non-neutral mutations.
the reason for this is that, forwsim assumes a fitness model that differs from the model implemented in ana-fits and sfs_code: for each non-neutral site, the user has to specify selection coefficients for the two homozygous cases and the heterozygous case.
thus, forwsim needs to iterate over all non-neutral mutations of each individual and recompute the fitness for each individual as opposed to ana-fits and sfs_code, where one simply multiplies the per-haplotype coefficients.
in the final analysis, the non-neutral mutation model implemented in forwsim, is more powerful , but slower and has a hard limit with respect to the number of non-neutral mutations that can be simulated.
in ana-fits, the fitness of an individual is the product of the fitness effects f  of all mutations on either haplotype.
figure
<dig> depicts the relative runtimes of forwsim and ana-fits.
the beg algorithm allows for simulation of neutral sequences that is between  <dig>  and  <dig>  faster than the delayed simulation in forwsim.
in less than 25% of all cases , the runtime improvement is below a factor of  <dig>  since for these cases only a small number of recombination/mutation events occur along short sequences and a small effective population size of  <dig> is used, the runtime improvement in ana-fits is mainly due to the more efficient random number generation.
for longer sequences with a higher number of individuals, our beg algorithm is substantially more efficient than the postponing algorithm of forwsim.
note that, these runtimes already include the overhead that is associated with generating and writing the beg to file.
validation
we performed a validation analogous to the validation of sfs_code <cit> .
since under neutrality, summary statistics of datasets that are simulated using ana-fits can easily be compared with the equivalent summary statistics generated by the popular coalescent software ms <cit> , we mainly validate datasets that have been simulated under neutrality.
we also compare simulations with ana-fits under selection to datasets produced by sfs_code.
under neutrality
to ensure that population genetic simulation software produces correct results, it is common to compare summary statistics for datasets produced by the software that needs to be validated to the respective summary statistics of widely-used and well-tested software packages.
to date, ms is among the most cited population genetic simulation tools.
slight deviations in summary statistics may be expected, since for instance ms implements an infinite-sites model of mutation, whereas ana-fits implements a full finite-sites model.
we used ana-fits to simulate chromosome-sized sequences.
thus, the effect of non-single nucleotide polymorphisms or multiple mutation events that give rise to the same genotype should not influence the respective summary statistics.
in addition to this standard validation setup, for ana-fits we also have to show that  the beg-algorithm yields correct results  and that  the algorithm is correctly implemented.
thus, irrespective of the respective simulation scenarios we show, that all three simulation modes produce the same summary statistics as ms:  all mutations are simulated forward in time ,  a part of the mutations is simulated forward in time , while the rest is obtained via an ancestry-based simulation ,  all mutations are simulated via ancestry-based simulation ,
we test the correctness of our algorithm and implementation for  <dig> specific scenarios employing the number of haplotypes, the number of segregating sites, the nucleotide diversity, and the site frequency spectrum as summary statistics.
the scenarios are :  mutation only ,  mutation and recombination,  mutation, recombination, a doubling in population size taking place  <dig> ·np generations ago,  mutation, recombination, a with a factor-2-bottleneck taking place  <dig> ·np generations ago,  mutation, recombination, with an exponential population growth setting in  <dig> ·np generations ago.
for all scenarios, we assume an initial population size of ne:= <dig>  a sequence length l:= <dig>  a mutation rate m:=10− <dig> and a recombination rate r:=10− <dig> .
for the forward simulation scenario we calculated 10·ne generations in total.
for each scenario, the sampling size was set to  <dig> sequences, except in the exponential scenario, where we sampled  <dig> sequences.
figure
<dig> shows the distribution of summary statistics after  <dig>  simulations for each of the three ana-fits parameters and for corresponding ms runs.
to visualize the number of segregating sites and the nucleotide density, we used the freedman-diaconis’ rule
<cit>  for a binning that retains a high level of resolution.
even using such a high resolution, the distributions in figure
<dig> almost correspond perfectly, with the exception of minor range deviations around the mode of the distribution.
under selection
to validate our software for simulations with selection, we compared the same summary statistics  as obtained by ana-fits to summary statistics for datasets simulated with sfs_code.
to keep the computational effort within reasonable limits, we reduced the number of simulations to  <dig> , the sequence length to  <dig> bp, the effective population size to ne= <dig>  and set the mutation and recombination rate per sequence and per base pair to m:=r:=10− <dig> .
figure
<dig> shows a comparison of summary statistics for two scenarios.
in scenario one, we chose a fixed selection coefficient s:= <dig>  using either a positive  or a deleterious effect  for all mutations.
in scenario  <dig>  selection coefficients are drawn from a normal distribution with mean μ=− <dig>   and standard deviation σ= <dig> .
this particularly high value for the mean and the fixed value of selection coefficients respectively was chosen to ensure that selection does indeed affect the summary statistics that are being assessed.
for sfs_code simulations, we simulated the  <dig> kbp sequence as  <dig> linked loci of length  <dig> kbp.
as shown in the figure
<dig>  the distributions of summary statistics of ana-fits show a high level of agreement to the respective distributions obtained with sfs_code.
large-scale simulation of partially neutral sequences
our algorithmic techniques and low-level technical optimizations improve ana-fits execution times by one to two orders of magnitude compared to sfs_code for plain forward simulation of sequences.
if the user specifies a fraction of mutations to be neutral, ana-fits can become faster by more than one additional order of magnitude .
ana-fits  runtimes  for simulation of neutral sequences with a varying proportion f of mutations simulated forward in time ; sfs_code runtime for comparison.
simulations for 10·ne generations, with per-base recombination and mutation frequency of  <dig> ·10− <dig>
given these runtime improvements, coupled with the continuous advances in hardware, large-scale simulations of genomic regions at the chromosome level are becoming feasible.
ana-fits also allows for simulating multiple chromosomes .
if we reduce the fraction of non-neutral mutations to 10%, the entire  human genome  can be simulated for  <dig> individuals within  <dig>  seconds, resulting in a total of  <dig>  segregating sites.
for  <dig> individuals , our implementation requires  <dig>  seconds, yielding a dataset with  <dig>  segregating sites.
CONCLUSIONS
we discussed the design and optimization of a highly efficient simulator for the forward-in-time simulation of population genetic datasets.
our software is up to two orders of magnitude faster than competing state-of-the art software.
moreover, we described the beg algorithm, that may further accelerate simulation for the common scenario, when a fraction of mutations does not have an impact on the fitness of individuals.
our algorithm outperforms a similar exploratory approach by up to one order of magnitude, while it is at the same time capable of retaining the entire observable history for each haplotype.
we demonstrated that our software allows for forward-in-time simulation at the whole-chromosome and even whole-genome level within acceptable execution times.
our novel graph structure  may become useful for population geneticists.
graph metrics can be computed on this structure, and the surviving material of any ancestrial haplotype can be reconstructed for any given point in time.
alternatively, local trees could be extracted from the beg for each sequence segment that did not undergo recombination.
we believe that the availability of ana-fits will lead to an increased usage of forward simulation in population genetics, given its importance for research on natural selection, and in the context of approximate bayesian computation.
thus, our tool can handle the analytical challenges of the post-genomic era.
abbreviations
prng: pseudo random number generator.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
aja and as designed the study.
aja developed the algorithm and implemented/optimized the simulator.
aja and as wrote the paper.
both authors read and approved the final manuscript.

 the field of population genetics strives to determine, how elementary evolutionary forces shape the genetic landscape within population of a species. as for related areas in evolutionary biology, rapid advances in next-generation sequencing technology have transformed the field with the completion of the first phase of the genome project representing a recent highlight. thus, population genetics is rapidly transforming into a data-driven science, while many methods and tools are not up to the challenge yet with respect to scalability and efficiency. because of the complex processes in population genetics, the inference machinery for analyzing the properties of a population is often limited and needs to be complemented by simulations. some of the most common use cases of simulation software in population genetics are the verification of novel methods/models and the generation of datasets for assessing how using a different model affects the properties of of a population. a promising, but computationally particularly expensive approach is to use simulations for inferring specific properties of a given data sample, as implemented, for instance by using approximate bayesian computation techniques. for an overview of available software and applications, see. forward-in-time simulation represents a powerful approach to simulating the evolutionary processes that act on genomic regions. essentially, sequences are represented in silico and the basic evolutionary operations/events are explicitly applied to each simulated individual, on a generation-by-generation basis. the inferred, simulated sequences represent the exact outcome of the underlying/assumed evolutionary process. however, this high accuracy comes at a high computational cost. execution times are quadratic with respect to the number of individuals, because the number of simulated generations needs to be increased with growing population size. to obtain accurate results, the size of the population that needs to be simulated must be chosen realistically. however, realistic population sizes range between a computationally already challenging individuals in humans to an entirely prohibitive number of, individuals in some species of drosophila melanogaster or, in escherichia coli. these excessive computational requirements explain the popularity of approximate coalescent simulations that are substantially faster. the coalescent-based approach is faster because only genomic samples that survived until the present are simulated backward in time. a substantial drawback of coalescent simulations is that, natural selection can only be incorporated into the model by means of a single mutation. to date, only few attempts have been made to improve the performance of forward simulators using optimization techniques from computer science and via low-level technical improvements. chadeau-hyam et al. developed a strategy to re-scale simulation parameters such as to decrease the effective population size that needs to be simulated. however, this shortcut induces a decrease in accuracy. forwsim implements a forward-algorithm that regularly simulates a small user-defined set of mutations under selection, forward in time. the neutral part of each haplotype is simulated separately and executed with a delay of generations. thereby, one can circumvent simulating haplotypes, if no ancestral material has survived until the present. distinguishing between neutral and non-neutral mutations for accelerating forward simulations represents a natural implementation choice, specifically since it is assumed that the majority of mutations is effectively neutral. finally, sfs_code has become a state-of-the-art forward simulator, because of an efficient implementation in c coupled with a plethora of simulation features and parameters. here, we describe design as well as optimization techniques that allow for substantially more efficient forward-in-time simulations. we introduce a novel algorithm, for the a posteriori simulation of neutral mutations. in addition to runtime gains, our algorithm generates a data structure that depicts the entire observable history of each individual sequence, a feature that –to the best of our knowledge– is unique among simulation software. we make available a corresponding software package, the ancestry-aware forward-in-time simulator, a c++ implementation that can handle challenging post-genomic datasets for population genetics analyses. methods in the following, we describe our algorithm for ancestry-based simulation of neutral mutations. subsequently, we discuss implementation and optimization issues. algorithm for simulating neutral mutations for forward simulation of a single generation of a population with effective population size ne, a scaled mutation rate of θ=4neμ, the scaled recombination rate ρ=4ner, the basic simulation steps are: sample diploid individuals by their fitness, determine the number of recombination events ) and create recombinants, determine the number of non-neutral mutations ) and mutate sequences, recalculate the fitness values for each individual. since neutral mutations, by definition, do not affect the survival probability of an individual, we do not need to simulate neutral mutations forward in time. this observation was first explored in forwsim, where the simulation of the neutral part of a sequence was delayed for a small number of generations. thereby, the actual calculations required for this part of the simulation can be discarded, if the haplotype under consideration did not survive until the present. the authors demonstrated that, the simulation of neutral sequences can be accelerated by a factor of two to five by optimally choosing this so-called look-ahead parameter, that is, the number of generations, by which the simulation of neutral mutations is delayed. in ana-fits, we deploy a different strategy: neutral mutations are not simulated forward in time at all. instead, we keep track of the entire ancestry of all surviving individuals, create a graph structure, and only extract the neutral part of all sequences from this graph structure once the forward simulation has been completed. in other words, we retain all information accumulated during the forward simulation phase, exactly determine all parts of sequences, where neutral mutations could have occurred and finally insert neutral mutations into the sequence a posteriori. note that, in analogy to forwsim and in contrast to coalescent simulations, our algorithm is equivalent to simulating all neutral mutations forward in time. thereby, we can guarantee simulation accuracy, which is one of the key advantages of forward simulations. simulating surviving neutral mutations during l forward-simulation generations, ana-fits annotates an ancestry al={a a …al}, where ai is a surjective function that maps the individual haplotypes in generation i to the haplotypes of the preceding generation i− in the following, we outline this procedure for a single locus/chromosome. the extension to multiple loci is described in the implementation section. the variable ai is defined as ai:hi,j→hi− l, where hi,j is the haplotype with index j in generation i, that is mapped to its parental haplotype with index l in generation i− haplotype indices are ordered, such that adjacent pairs that start with an even index represent the two homologous versions of each locus for diploid organisms or represent an individual). as mentioned before, the algorithm also keeps track of all recombination events that occurred during the forward phase. for each event, we store the absolute sequence position, the haplotype index, and the generation number. for a sequence of length l, the information in al and the corresponding set of recombination events can be used to determine the regions of ancestral haplotypes that survived until the present. only mutations that occurred within these regions are observable in the present. in the first step, our algorithm determines haplotypes in ancestral generations that contributed material to the present generation. for these regions, we simulate neutral mutations and identify those recombination events that contributed to mutations in the present haplotype instances. the algorithm starts at the present generation and initializes surviving regions for each haplotype with a maximum interval of, where l is the sequence length. regions for ancestral haplotypes are initialized by the tuple. note that, the surviving regions of a haplotype instance can be fragmented/scattered into many segments by recombinations. to economize on runtime and memory, we only keep track of the start position of the very first segment rs and the end position of the very last segment re. by using this approximate surviving ’super-’region r=, we overestimate the size of the actual surviving region. we correct for this in a later stage of the algorithm, when sequences are extracted from the graph structure. finally, to determine the surviving haplotype region in generation i− we iterate over all survivors si of generation i starting with haplotypes contained in the present generation. for each surviving haplotype hi,l∈si with region, we propagate the surviving regions as described in the following algorithm to the members of the preceding generation. let hi− k with region be the parent haplotype and hi− k+ with region be the homologous haplotype in the corresponding parent individual: for simplicity, the algorithm described above assumes that the recombinant emerged from at most one recombination event. for multiple recombination events per recombinant, the algorithm has to be adapted, such that surviving regions of parental haplotypes are extended until the next recombination breakpoint. this has also been implemented in the ana-fits software. during the procedure described above, we also determine the set of survivors si− from the preceding generation. in the case that a recombination event does not split the surviving region, we can simply ignore it, if the following property does not hold: assume, haplotype hi,l is mapped to hi− k with surviving region and hi,l is a recombinant with a breakpoint b<rs. this means that, hi,l does not inherit mutations from hi− k, but instead from the homologous haplotype hi− k+ it recombined with. if we can not ignore the recombination, we create a node instance for the graph that will be constructed in the next step. the survivors of each generation are stored for further processing in subsequent steps. once survivors and surviving regions for a generation i have been determined, we can then simulate the neutral mutations that occurred in this generation and that are still observable in the present generation. thus, the expected number of neutral mutations is reduced by a factor of f=|si|2·ne. furthermore, if a neutral mutation occurs outside the estimated surviving region, it can simply be ignored. as shown in figure depending on the recombination rate r, the number of haplotypes that contribute genetic material to the present generation quickly converges and may cause values of f to become small. note that, we expect the algorithm to become inefficient for exceptionally high recombination rates, such as r:=. graph construction in the algorithm step described above, we determined survivors si for each generation i, the set of relevant recombinations, and the set of relevant neutral mutations. for all nodes, we store the generation of origin and the corresponding haplotype index hi,l. the goal of the procedure we describe here, is to construct a graph. the beg is a directed acyclic graph with two types of nodes: m-nodes with out-degree and r-nodes with out-degree thus, each node in the graph represents the state of an ancestral haplotype that reflects the history of all mutation or recombination events experienced by its ancestors. a directed edge connecting node n to node n indicates that n emerged from n after mutation or recombination. given this information, the graph can be constructed forward in time using an array ai for the current generation and an array ai− for the preceding generation i− that contains node references/pointers. array ai is populated by carrying over/propagating references from ai− according to the survival information stored in si− if a surviving haplotype of generation i underwent recombination or mutation, the respective r- or m-node nnew is used to replace the node nold that was propagated and an edge is created between nnew and nold. we also create an additional edge to the haplotype which is homologous to nold for r-nodes. figure provides a complex example, where a recombinant emerges from three recombination events that originated at two homologous ancestral haplotypes and extraction of sequences the algorithm described above constructs the beg and an array of haplotype states ni, that survived until the present. in the final step, we extract the exact sequence for each node ni from the beg by recursively traversing the graph starting at ni, while keeping track of the region of interest which is segmented by each recombination event which will be encountered. since we execute the backtracking procedure for each node ni individually, the exact borders of all segments of the surviving region for this node are known at any given point in time. thus, we can ignore all mutation events that did not occur within the surviving region, but were inserted into the graph because of overestimation or because this mutation forms part of the surviving region of another haplotype. in practice, the above extraction procedure is inefficient. as a consequence, the sequence reconstruction step dominates execution times. consider two haplotype states n and n that have survived, where n is an m-node with ancestor n for both nodes, we have to traverse all ancestral nodes of n and n to obtain the two corresponding sequences s and s to accelerate this process, it suffices to create sequence s via backtracing starting at n and to generate s as a copy of s and then simply add the mutation events induced by node n in the end. in the following, we explain, how the computational cost of sequence extraction can be reduced. prior to the actual sequence extraction step, the graph is traversed twice, starting with each haplotype ni that survived until the present. in the first traversal, we approximate the surviving region for each node as described in the algorithm. thus, if we decide to explicitly determine and store the sequence of a node, we already know the maximum size of the region of interest beforehand. in a second traversal that starts with each haplotype in the present, we determine for which nodes intermediate sequences should be created. in other words, the surviving mutations are represented explicitly in this node. during these traversals, we actually stop a traversal, if the traversed node has already been visited by another traversal. for each node, we also store the number v of times the node has been visited. we refer to nodes with v≥ as coalescent nodes, since ancestral material of two lineages coalesces in these nodes. finally, we determine the distribution of v and the 5% nodes with highest v are represented explicitly, while all remaining nodes will be traversed several times during sequence extraction. then, the sequences can be extracted via a final full traversal. the 5% cut-off was determined empirically. it yields good results with respect to balancing runtime versus memory requirements. note that, if every coalescent node is represented explicitly, memory requirements would become prohibitive for whole-chromosome simulations with individuals. in general, the 5% cut-off proves to be sufficient to attain a substantial decrease of the traversal cost for sequence extraction. implementation and optimization the ancestry-based algorithm for simulating neutral mutations as discussed in the previous section has specifically been developed to improve runtimes of forward-simulations. in the following, we discuss low-level implementation issues and optimization techniques of ana-fits. memory requirements of ancestry and beg for the majority of possible ana-fits invocations with respect to user parameter combinations, storing the ancestry al requires the largest amount of memory. depending on the effective population size ne, we need to store ancestry information for 10·ne haplotypes. if this exceeds the amount of available memory, the user can set a command line flag such that the memory utilization of the ancestry does not exceed a user-specified soft limit. ana-fits then tries to split up the forward simulation into sections and updates the graph after each section leading to increased runtimes. in some instances the graph representation itself can become at least as memory intensive as the ancestry. furthermore, for large graphs, the explicit bit vector representation of intermediate sequences may further increase memory requirements. therefore, currently no hard upper memory utilization limit can be imposed on ana-fits runs. the ancestry al only works correctly for a single locus. for multiple unlinked loci, we map each individual to its two parent individuals instead of mapping it to individual haplotypes. for this, we use a per-chromosome bit vector to keep track of which of the two possible homologous haplotypes was inherited from the respective parent. apart from their high memory-efficiency, bit vectors can also be efficiently initialized with random bits using a random number generator. to avoid unnecessary and inefficient frequent memory allocations, the entire ancestry and the array of survivors are allocated as single, monolithic blocks. for the ancestry, we dynamically determine the minimal number of bytes that is required for storing a generation as a function of the population size of the preceding generation. since, under most simulation scenarios, ne does not exceed individuals, it is mostly sufficient to use the unsigned short integer datatype for saving memory. note that, beside all efforts to keep the memory consumption of ana-fits as low as possible, the entire ancestry is kept in main memory during simulation. this necessarily means that ana-fits exhibits higher memory requirements than comparable programs. random number generation high quality random number generators are essential for simulation tools. to date, the mersenne twister algorithm is considered as state-of-the-art in pseudo-random number generation. a simd-based version of the mersenne twister employs streaming single instruction, multiple data extensions instructions to generate random numbers at almost twice the speed as the original mersenne twister. ana-fits comes with a distribution of randomlib, a c++ implementation of sfmt. for efficiency, randomlib initializes an array with random bytes and then transforms the bytes of this array into the required primitive data type for the requested range. for forward simulations, we can exploit two application-specific properties to more efficiently use prngs: usually for the ancestry, integer values of less than bytes are required, and uniform random numbers are frequently drawn from the same interval, where n is the effective population size of the previous generation). when a large amount of random numbers with these two properties is required, we can directly initialize a target array from the internal array of the prng and transform the minimal number of random bytes for the respective integer type into uniform random numbers for the required range using sse instructions. for instance, if simulations require sampling integers ∈[ ), we can copy bytes from the prng’s internal array into a location of the ancestry memory block and use vector instructions to transform the bytes into integers of the specified range. without this optimization, randomlib would produce bytes of randomness and transform bytes into one integer separately. non-neutral sequence representation for forward simulators that work with a given number and pre-defined locations of polymorphic sites, it is straight-forward to represent haplotypes as bit arrays. if we want to simulate a sequence with a huge number of base-pairs, bit arrays are not the data structure of choice any more because they become too sparse. thus, in other forward simulator implementations, a haplotype is stored as a sequence of mutations with respect to an initial reference sequence. for instance, sfs_code uses self-balancing binary trees to represent the polymorphisms of a haplotype. sorted binary trees allow for rapid insertion of novel mutations and fast extraction of sub-sequences which are required for creating recombinants. an alternative data structure is used in forwsim that maintains sorted arrays of sequence positions. since for ana-fits, we want to store additional information about each mutation, we use a data structure that is based on references/pointers to mutation object instances. when sorted arrays are used, sub-sequence replication for creating recombinants can be efficiently implemented using the memcpy() system call. this is because finding coordinates within the sequence via a binary search has the same time complexity of o) as searches in a sorted tree. additionally, changing the array size induces a smaller memory allocation overhead than allocating additional nodes in sorted trees. the main drawback of sorted arrays is the insertion cost of o. in practice however, obtaining the correct relative position via a binary search is usually more expensive than executing an insertion as long as insertions are implemented via the memcpy() system call. note that, the memcpy() system call is also vectorized via sse intrinsics. finally, the number of insertion operations is relatively small in ana-fits, since we only simulate non-neutral mutations forward in time. we experimented with various sorted tree implementations, for instance, the c++ stl std::set implementation or judy arrays by, which are based on cache-optimized 256-ary trees. for ana-fits, we found that the tree data structures did not outperform our sorted array implementation. memory management in forward simulations, the number of distinct haplotype instances containing the sorted array of mutation references and further meta-information and the number of mutation instances changes dynamically at a rapid pace. thus, dynamic memory allocations and free()) that are invoked in the inner parts of nested loops quickly dominate runtimes. as mentioned in the previous section, ana-fits uses block allocation, instead of dynamic allocation whenever possible. apart from being faster this also prevents heap fragmentation. one approach to handle the rapid turnover of haplotype and mutation instances more efficiently is to use free lists in conjunction with reference counting. the underlying idea is to keep track for each allocated instance, how often it is referenced in the current generation or in present haplotypes. instead of deallocating unreferenced instances, instances are appended to a list, the free list, from which new instances can then be allocated again, when required. however, this allocation strategy still uses more than 10% of total runtime. instead, we use two task-specific allocation schemes that keep memory allocation overhead to an absolute minimum. for haplotype instances, we deploy a free ring scheme. as the name suggests, a free ring contains references to allocated haplotype instances in a similar way as a free list. whenever a haplotype is carried over into the next generation, we update its generation number. thus, when haplotype instance memory is requested from the free ring structure, an iterator traverses the elements in the free ring and returns the first instance that is neither used in the current nor the previous generation. if no free instance is available, the free ring capacity is doubled by allocating new haplotype instances. since it is unlikely that all references of a haplotype will be carried over to the next generation in the sampling phase without being subject to recombination or mutation, the iterator typically finds an unused instance after a small number of search steps. in contrast to this, mutation instances need to be allocated in blocks, since references to instances represent an inefficient mechanism for this task, because mutation instances are comparably small in size. moreover, for large numbers of mutations, a corresponding large number of memory allocation invocations would be required, with a negative impact on performance. allocating monolithic blocks with dynamic re-sizing is also not possible, because addresses of mutation instances can change which would invalidate the references to mutations that are stored in the haplotype instances. instead, we allocate multiple blocks of instances and implement an iterator that searches for an unused instance, such that it jumps to the next block when the end of the current block has been reached. when a free instance has been found, it is flagged as in use. this flag is updated in a cleanup step that is applied after a fixed number of generations. in this cleanup step, mutation instances that occur in all haplotypes are removed from all haplotypes. thus, this phase can be used to efficiently update the in use flag of each mutation without inducing a large overhead. since, under this scheme, the addresses of instances do not change, we refer to this allocation structure as address-preserving free ring. implementation of neutral sequence extraction the main challenge in implementing the beg algorithm is to design an efficient method for extracting sequences from the graph. for this task, we can use the following property: during the sequence extraction phase, all neutral mutations that occurred during simulation are known. thus, bit arrays are the most adequate data structure for storing neutral sequences. after sorting mutations by their location, each mutation is assigned its corresponding index in the bit array. as a result, the cost of adding a mutation that is encountered during sequence extraction to a given sequence is reduced to merely setting a bit. handling recombinations is computationally more expensive, since we have to determine the bit array indices of the start and end points for the surviving regions of the parent haplotype via a binary search. carrying over the relevant mutations from parent haplotypes is implemented as a bit-wise or-operation. this implementation is particularly efficient because the sequence output of the simulation can directly be written as bit array to a binary file. thus, output files do not require a lot of disk space and programs that directly evaluate the output of a large simulation can efficiently parse and post-process this output data.  in this section, we evaluate the runtime performance of ana-fits. in general, conducting a fair comparison between forward simulation tools represents a difficult task, since each simulator uses different assumptions. for a meaningful technical analysis, we separately evaluate the two key aspects of ana-fits: efficient implementation of non-neutral forward simulation and ancestry-based simulation of neutral mutations. for realistic use cases, both features independently contribute a certain amount of run time improvement. for all runs discussed throughout this section, we set mutation and recombination frequencies to μ per base=rper base= ·10− which roughly correspond to the empirical rates for humans. for fixed μ and r, the runtime depends on the locus length l and the effective population size ne. all runs were executed on a 4× core machine with a total of gb ram. the execution times provided here are the result of independent runs with different random number seeds. we discard the slowest runtime in order to correct for the impact of seed values that generate exceptionally expensive simulations and average over the remaining three. moreover, we validate ana-fits by comparing the output of our software to datasets generated with established simulation codes. finally, we describe runtime gains for the simulation of partially neutral sequences, specifically in a whole-genome setting for human sequences. runtime improvement: forward simulation only here, we evaluate the efficiency of the pure forward simulation part of our implementation and compare it to sfs_code. ana-fits can simulate each mutation forward in time by choosing a constant small selection coefficient s≠ for each mutation that is to be simulated. if s is small enough, it is effectively neutral. this way the beg algorithm is not applied. simulating close to neutral mutations, forward in time, represents the most expensive case, since for large values of s the mutation quickly undergoes selection and for small negative values it is removed from the gene pool by selection. in both cases the mutation quickly vanishes which in turn reduces runtimes. in runs with sfs_code, we simulate neutral mutations forward in time. in general, it is not possible to efficiently simulate one contiguous locus of the length of a chromosome with sfs_code. instead, for simulating a locus of, base pairs, the user has to break up the locus into a number of linked loci. for obtaining sfs_code runtimes, we did not optimize the break-up strategy exhaustively, but followed the general suggestions by the author and simulated chromosome lengths listed in figure as 50· kbp, 50· kbp and 500· kbp. figure depicts the relative runtimes of both programs. overall, ana-fits is between and times faster than sfs_code. in all except runs, ana-fits is at least one order of magnitude faster than sfs_code. the runtime improvements do not show a clear dependency on one of the two input dimensions. we assume that this is because, for certain steps, sfs_code switches between alternative implementations depending on the number of recombination/mutation events that are expected per generation. further parameters that may influence the runtime differences between the codes are: frequency of the cleanup steps in either code and number of linked loci the sequence was split into. runtime improvement: ancestry-based simulation as already mentioned, the beg algorithm is inspired by the postponed simulation of neutral mutations as implemented in forwsim. thus, in order to evaluate the stand-alone efficiency of the beg algorithm, we compare the runtime for simulating entirely neutral sequences between ana-fits and forwsim. note that, forwsim simplifies the fisher-wright model in two ways: while it assumes finite sites, at most one mutation per site is allowed and haplotypes undergo at most one recombination event per generation. we omit a comparison between forwsim and ana-fits that includes forward simulation of non-neutral mutations. the reason for this is that, forwsim assumes a fitness model that differs from the model implemented in ana-fits and sfs_code: for each non-neutral site, the user has to specify selection coefficients for the two homozygous cases and the heterozygous case. thus, forwsim needs to iterate over all non-neutral mutations of each individual and recompute the fitness for each individual as opposed to ana-fits and sfs_code, where one simply multiplies the per-haplotype coefficients. in the final analysis, the non-neutral mutation model implemented in forwsim, is more powerful, but slower and has a hard limit with respect to the number of non-neutral mutations that can be simulated. in ana-fits, the fitness of an individual is the product of the fitness effects f of all mutations on either haplotype. figure depicts the relative runtimes of forwsim and ana-fits. the beg algorithm allows for simulation of neutral sequences that is between and faster than the delayed simulation in forwsim. in less than 25% of all cases, the runtime improvement is below a factor of since for these cases only a small number of recombination/mutation events occur along short sequences and a small effective population size of is used, the runtime improvement in ana-fits is mainly due to the more efficient random number generation. for longer sequences with a higher number of individuals, our beg algorithm is substantially more efficient than the postponing algorithm of forwsim. note that, these runtimes already include the overhead that is associated with generating and writing the beg to file. validation we performed a validation analogous to the validation of sfs_code. since under neutrality, summary statistics of datasets that are simulated using ana-fits can easily be compared with the equivalent summary statistics generated by the popular coalescent software ms, we mainly validate datasets that have been simulated under neutrality. we also compare simulations with ana-fits under selection to datasets produced by sfs_code. under neutrality to ensure that population genetic simulation software produces correct results, it is common to compare summary statistics for datasets produced by the software that needs to be validated to the respective summary statistics of widely-used and well-tested software packages. to date, ms is among the most cited population genetic simulation tools. slight deviations in summary statistics may be expected, since for instance ms implements an infinite-sites model of mutation, whereas ana-fits implements a full finite-sites model. we used ana-fits to simulate chromosome-sized sequences. thus, the effect of non-single nucleotide polymorphisms or multiple mutation events that give rise to the same genotype should not influence the respective summary statistics. in addition to this standard validation setup, for ana-fits we also have to show that the beg-algorithm yields correct results and that the algorithm is correctly implemented. thus, irrespective of the respective simulation scenarios we show, that all three simulation modes produce the same summary statistics as ms: all mutations are simulated forward in time, a part of the mutations is simulated forward in time, while the rest is obtained via an ancestry-based simulation, all mutations are simulated via ancestry-based simulation, we test the correctness of our algorithm and implementation for specific scenarios employing the number of haplotypes, the number of segregating sites, the nucleotide diversity, and the site frequency spectrum as summary statistics. the scenarios are: mutation only, mutation and recombination, mutation, recombination, a doubling in population size taking place ·np generations ago, mutation, recombination, a with a factor-2-bottleneck taking place ·np generations ago, mutation, recombination, with an exponential population growth setting in ·np generations ago. for all scenarios, we assume an initial population size of ne:= a sequence length l:= a mutation rate m:=10− and a recombination rate r:=10−. for the forward simulation scenario we calculated 10·ne generations in total. for each scenario, the sampling size was set to sequences, except in the exponential scenario, where we sampled sequences. figure shows the distribution of summary statistics after simulations for each of the three ana-fits parameters and for corresponding ms runs. to visualize the number of segregating sites and the nucleotide density, we used the freedman-diaconis’ rule for a binning that retains a high level of resolution. even using such a high resolution, the distributions in figure almost correspond perfectly, with the exception of minor range deviations around the mode of the distribution. under selection to validate our software for simulations with selection, we compared the same summary statistics as obtained by ana-fits to summary statistics for datasets simulated with sfs_code. to keep the computational effort within reasonable limits, we reduced the number of simulations to, the sequence length to bp, the effective population size to ne= and set the mutation and recombination rate per sequence and per base pair to m:=r:=10−. figure shows a comparison of summary statistics for two scenarios. in scenario one, we chose a fixed selection coefficient s:= using either a positive or a deleterious effect for all mutations. in scenario selection coefficients are drawn from a normal distribution with mean μ=− and standard deviation σ=. this particularly high value for the mean and the fixed value of selection coefficients respectively was chosen to ensure that selection does indeed affect the summary statistics that are being assessed. for sfs_code simulations, we simulated the kbp sequence as linked loci of length kbp. as shown in the figure the distributions of summary statistics of ana-fits show a high level of agreement to the respective distributions obtained with sfs_code. large-scale simulation of partially neutral sequences our algorithmic techniques and low-level technical optimizations improve ana-fits execution times by one to two orders of magnitude compared to sfs_code for plain forward simulation of sequences. if the user specifies a fraction of mutations to be neutral, ana-fits can become faster by more than one additional order of magnitude. ana-fits runtimes for simulation of neutral sequences with a varying proportion f of mutations simulated forward in time; sfs_code runtime for comparison. simulations for 10·ne generations, with per-base recombination and mutation frequency of ·10− given these runtime improvements, coupled with the continuous advances in hardware, large-scale simulations of genomic regions at the chromosome level are becoming feasible. ana-fits also allows for simulating multiple chromosomes. if we reduce the fraction of non-neutral mutations to 10%, the entire human genome can be simulated for individuals within seconds, resulting in a total of segregating sites. for individuals, our implementation requires seconds, yielding a dataset with segregating sites.  we discussed the design and optimization of a highly efficient simulator for the forward-in-time simulation of population genetic datasets. our software is up to two orders of magnitude faster than competing state-of-the art software. moreover, we described the beg algorithm, that may further accelerate simulation for the common scenario, when a fraction of mutations does not have an impact on the fitness of individuals. our algorithm outperforms a similar exploratory approach by up to one order of magnitude, while it is at the same time capable of retaining the entire observable history for each haplotype. we demonstrated that our software allows for forward-in-time simulation at the whole-chromosome and even whole-genome level within acceptable execution times. our novel graph structure may become useful for population geneticists. graph metrics can be computed on this structure, and the surviving material of any ancestrial haplotype can be reconstructed for any given point in time. alternatively, local trees could be extracted from the beg for each sequence segment that did not undergo recombination. we believe that the availability of ana-fits will lead to an increased usage of forward simulation in population genetics, given its importance for research on natural selection, and in the context of approximate bayesian computation. thus, our tool can handle the analytical challenges of the post-genomic era. abbreviations prng: pseudo random number generator. competing interests the authors declare that they have no competing interests. authors’ contributions aja and as designed the study. aja developed the algorithm and implemented/optimized the simulator. aja and as wrote the paper. both authors read and approved the final manuscript. 
in addition, we present a novel algorithm that further improves runtimes by up to an additional order of magnitude, for simulations where a fraction of the mutations is neutral. apart from simulated sequences, our tool also generates a graph structure that depicts the complete observable history of neutral mutations. the graph structure generated by our algorithm can give rise to novel approaches for visualizing and analyzing the output of forward simulations. in population genetics, simulation is a fundamental tool for analyzing how basic evolutionary forces such as natural selection, recombination, and mutation shape the genetic landscape of a population. we introduce ana-fits, a highly optimized forward simulation software, that is up to two orders of magnitude faster than current state-of-the-art software.
 in population genetics, simulation is a fundamental tool for analyzing how basic evolutionary forces such as natural selection, recombination, and mutation shape the genetic landscape of a population. forward simulation represents the most powerful, but, at the same time, most compute-intensive approach for simulating the genetic material of a population.  we introduce ana-fits, a highly optimized forward simulation software, that is up to two orders of magnitude faster than current state-of-the-art software. in addition, we present a novel algorithm that further improves runtimes by up to an additional order of magnitude, for simulations where a fraction of the mutations is neutral. apart from simulated sequences, our tool also generates a graph structure that depicts the complete observable history of neutral mutations.  the substantial performance improvements allow for conducting forward simulations at the chromosome and genome level. the graph structure generated by our algorithm can give rise to novel approaches for visualizing and analyzing the output of forward simulations. population geneticsforward-in-time simulationfisher-wright modelalgorithmsoftwarenatural selection 