the field of population genetics endeavors to ascertain how basic evolutionary forces , such as natural selection , recombination , and mutation , shape the patterns of genetic variation within and between populations.1,2 computer simulation software has traditionally been used to explore analytically intractable genetic models.3,4 continuous advances in numerical simulation and the wide availability of computational resources allow researchers to use numerical simulation to test mathematical models in virtual populations , and even to analyze genetic data.5,6 as a direct consequence , there is a plethora of simulators available , each one tailored to a specific scenario . this forces geneticists to choose a simulator depending on the research being conducted.3 the challenge is to provide flexible simulation programs with the capacity for dealing with complex and realistic demographic / evolutionary models , using efficient algorithms to deal with genome - wide data sets processing complexity.4,7,8 the availability of cost - effective next - generation sequencing technologies has democratized whole - genome sequencing , making large genomic data sets available to most researchers.2,9,10 among the two approaches to simulation algorithms , forward - in - time and coalescent - based , the latter is most widely used because of its efficiency and flexibility.3,11,12 in particular , the standard coalescent approach was shown to be extremely efficient for short sequences ( less than few mega base pairs ) , but it becomes computationally demanding for simulating long genome regions with large recombination rate.13 new approaches have emerged to overcome this issue , such as the sequential markov coalescent and markov chain monte carlo methods5,8,1416 and alternative implementations of the exact coalescent , including algorithmic and data structure optimizations , such as scrm and msprime.17,18 in this article , we focus on coalescent simulators , which are shown to be computationally intractable when working at the genome scale,8,11,19 and more specifically on mspar,20 which is the parallel version of hudson s ms coalescent simulator,21 ie , the most classical and widely used coalescent simulator.3,7,22 mspar addresses both the problem of sampling long genomic regions with large recombination rates and running analysis requiring a large number of samples.20 we have previously shown that mspar execution time is orders of magnitude faster than hudson s ms when simulating large samples of long genomic sequences and large recombination rates using a high performance computing ( hpc ) cluster,20 but it is less efficient than the other approximated standard coalescent simulators such as macs and fastsimcoal.1315 our interest in further developing the standard coalescent simulators resides in the underlying side effects shown by approximated methods , like the one implemented by macs : 
 we can suffer loss of accuracy for certain evolutionary models as when sampling from populations separated by reduced gene flow.23ignoring type 2 recombination events ( from marjoram and wall s classification24 ) may impact on the mean and variance of most recent common ancestor times when long sequences are simulated.25 we can suffer loss of accuracy for certain evolutionary models as when sampling from populations separated by reduced gene flow.23 ignoring type 2 recombination events ( from marjoram and wall s classification24 ) may impact on the mean and variance of most recent common ancestor times when long sequences are simulated.25 regarding the implementation details , we have considerably improved mspar s memory management , improved its master  worker implementation , removed the data structure dependencies inherited from hudson s ms , and made use of the vader shared memory byte transport layer ( btl ) provided by openmpi v1.8.26 by applying a message - passing / shared - memory parallelization strategy , msparsm is able to achieve similar execution times to those of macs , simulating long stretches of dna sequences with large recombination rates . thus , we significantly reduce the computational burden of the standard coalescent method . for the purpose of comparison , we have used two case studies , with fixed mutation rates approximating the human and drosophila melanogaster average values . accuracy of simulations is maintained to be the same as the one from hudson s ms . it is ensured by the fact that the code in charge of sample generation was maintained unchanged from original hudson s ms code . it is also important to note that msparsm , as same as mspar and hudson s ms , only generates samples under a wright  fisher neutral model of genetic variation . in this section , we describe the parallelization approach taken over mspar , the experimental setup we designed for evidence gathering , and the hardware / software used to run the experiments . worker application design that includes the following : removing the data structure dependencies inherited from hudson s ms , refactoring mspar s memory management , and improving the communication patterns used in the mspar s master  worker implementation . first , we have taken the advantage of the openmpi s new transport layer ( vader shared - memory btl ) for transferring data between communication end points , which was introduced in version 1.8 . the message passing interface ( mpi ) transport layer used by mspar is known as sm btl ( shared - memory btl ) , which follows a copy - in / copy - out pattern : when an mpi process x sends a message to a process y , the message is first copied from the x s buffer to the shared memory , and then the receiver ( ie , the process y ) copies the message from the shared memory into a buffer.27 the vader shared memory btl provides support for xpmem linux kernel module and allows to directly transfer messages between sender and receiver buffers when both mpi processes are in the same node , thereby saving one copy operation compared against the previous transport layer used included in openmpi ( ie , sm btl ) . to take the advantage of vader btl , we have changed the master  worker topology used by mspar . given n as the total number of available processes for parallel computation , we define one single global master process coordinating the n1 remaining worker processes . in msparsm , we have one master process per compute node ( henceforth referred to as node master ) and one global master process that coordinates the m node masters ( where m is the number of compute nodes ) . each master node coordinates the worker processes located in its node , returning a consolidated message with all of the generated samples back to the global master when computation is performed . it is important to note that there is going to be one node hosting two master processes : the global master and the node master . by using non - blocking mpi collectives , we enabled node masters to also generate samples while waiting for local worker communications ; therefore , the loss of worker processes compared with mspar is compensated as much as possible . in a setup with m compute nodes and p cores per node , in mspar , we have a total number of m*p  1 dedicated worker processes , while in msparsm the number is m*(p  1 )  1 . for the use case of the scientist using a single node for the computation , eg , a single fat node or even a workstation , the application is not going to use a global master but a single master that is also going to generate replica samples . another improvement we made is to allow mpi processes located in the same compute node as the global master , to directly output the generated samples avoiding the use of mpi point - to - point communications with a master process . summarizing the major improvements , we 
 greatly decreased the number of memory allocation calls by favoring memory reusing and reallocation as much as possible;refactored the master  worker strategy to leverage the use of the vader btl;included the master process in the replica generation process;made worker processes more autonomous to reduce the number of mpi communications;removed all global structures inherited from hudson s ms , making it possible to explore a fine - grained parallelization approach ( eg , using openmp and/or cuda ) . greatly decreased the number of memory allocation calls by favoring memory reusing and reallocation as much as possible ; refactored the master worker strategy to leverage the use of the vader btl ; included the master process in the replica generation process ; made worker processes more autonomous to reduce the number of mpi communications ; removed all global structures inherited from hudson s ms , making it possible to explore a fine - grained parallelization approach ( eg , using openmp and/or cuda ) . we compared our msparsm with hudson s ms and macs in terms of execution time , using two case studies with fixed population mutation rate values ( ie ,  = 4n , where n is the current population size and  is the mutation rate per site ) , one of them approximating the estimated human average ( case 1 ) and the other one approximating the estimated d. melanogaster average ( case 2 ;  = 0.001 and 0.005 values , meaning that 1 and 5 are differences between two random individuals for each 1000 sites , respectively ) . the number of generated replicas was set at 300 , each one with 100 chromosomes , a genetic region of 1  10 bp , and population recombination per site ( 4n , where  is the mutation rate between two contiguous sites ) ranging from 0.0005 to 0.08 values . the parameters used in each of the simulations performed in this work are given in supplementary table 1 . we implemented the parallel application in ansi c and compiled it using the linux gnu compiler 4.9.1 . inter - process communication was implemented using openmpi 1.10.1.28 we ran the binaries on a cluster using up to eight nodes each with two intel xeon x5660 six - core processors using hyper - threading technology running at 2.80 ghz , 12 mb l3 cache , and 96 gb of double data rate synchronous dynamic - random access memory ( ddr - ram ) . worker application design that includes the following : removing the data structure dependencies inherited from hudson s ms , refactoring mspar s memory management , and improving the communication patterns used in the mspar s master  worker implementation . first , we have taken the advantage of the openmpi s new transport layer ( vader shared - memory btl ) for transferring data between communication end points , which was introduced in version 1.8 . the message passing interface ( mpi ) transport layer used by mspar is known as sm btl ( shared - memory btl ) , which follows a copy - in / copy - out pattern : when an mpi process x sends a message to a process y , the message is first copied from the x s buffer to the shared memory , and then the receiver ( ie , the process y ) copies the message from the shared memory into a buffer.27 the vader shared memory btl provides support for xpmem linux kernel module and allows to directly transfer messages between sender and receiver buffers when both mpi processes are in the same node , thereby saving one copy operation compared against the previous transport layer used included in openmpi ( ie , sm btl ) . to take the advantage of vader btl , we have changed the master  worker topology used by mspar . given n as the total number of available processes for parallel computation , we define one single global master process coordinating the n1 remaining worker processes . in msparsm , we have one master process per compute node ( henceforth referred to as node master ) and one global master process that coordinates the m node masters ( where m is the number of compute nodes ) . each master node coordinates the worker processes located in its node , returning a consolidated message with all of the generated samples back to the global master when computation is performed . it is important to note that there is going to be one node hosting two master processes : the global master and the node master . by using non - blocking mpi collectives , we enabled node masters to also generate samples while waiting for local worker communications ; therefore , the loss of worker processes compared with mspar is compensated as much as possible . in a setup with m compute nodes and p cores per node , in mspar , we have a total number of m*p  1 dedicated worker processes , while in msparsm the number is m*(p  1 )  1 . for the use case of the scientist using a single node for the computation , eg , a single fat node or even a workstation , the application is not going to use a global master but a single master that is also going to generate replica samples . another improvement we made is to allow mpi processes located in the same compute node as the global master , to directly output the generated samples avoiding the use of mpi point - to - point communications with a master process . summarizing the major improvements , we 
 greatly decreased the number of memory allocation calls by favoring memory reusing and reallocation as much as possible;refactored the master  worker strategy to leverage the use of the vader btl;included the master process in the replica generation process;made worker processes more autonomous to reduce the number of mpi communications;removed all global structures inherited from hudson s ms , making it possible to explore a fine - grained parallelization approach ( eg , using openmp and/or cuda ) . greatly decreased the number of memory allocation calls by favoring memory reusing and reallocation as much as possible ; refactored the master worker strategy to leverage the use of the vader btl ; included the master process in the replica generation process ; made worker processes more autonomous to reduce the number of mpi communications ; removed all global structures inherited from hudson s ms , making it possible to explore a fine - grained parallelization approach ( eg , using openmp and/or cuda ) . we compared our msparsm with hudson s ms and macs in terms of execution time , using two case studies with fixed population mutation rate values ( ie ,  = 4n , where n is the current population size and  is the mutation rate per site ) , one of them approximating the estimated human average ( case 1 ) and the other one approximating the estimated d. melanogaster average ( case 2 ;  = 0.001 and 0.005 values , meaning that 1 and 5 are differences between two random individuals for each 1000 sites , respectively ) . the number of generated replicas was set at 300 , each one with 100 chromosomes , a genetic region of 1  10 bp , and population recombination per site ( 4n , where  is the mutation rate between two contiguous sites ) ranging from 0.0005 to 0.08 values . the parameters used in each of the simulations performed in this work are given in supplementary table 1 . we implemented the parallel application in ansi c and compiled it using the linux gnu compiler 4.9.1 . inter - process communication was implemented using openmpi 1.10.1.28 we ran the binaries on a cluster using up to eight nodes each with two intel xeon x5660 six - core processors using hyper - threading technology running at 2.80 ghz , 12 mb l3 cache , and 96 gb of double data rate synchronous dynamic - random access memory ( ddr - ram ) . we evaluated the performance of msparsm by analyzing the speedup and efficiency as a function of increasing numbers of compute nodes . the speedup is defined as sp = t1/tp , where p is the number of processors , t1 is the execution time of the sequential application , and tp is the execution time of the parallel application with p processors . efficiency is defined as ep = sp / p . in figure 1 , the speedup of msparsm compared with the sequential hudson s ms version is shown , which was registered in case study 1 . we observed a high parallel efficiency of msparsm and a strong increase in speed in relation to sequential hudson s ms version . this happened because the node memory was mostly exhausted by the mpi processes ( 90 gb in use out of 96 gb available ) , implying a penalty because of the page swapping . in figure 2 , the same graphic data as before are shown , but for the case study 2 . compared with the case study 1 2a ) , but there was a similar performance for upper values ( ie , 2000 , 4000 , and 8000 ) . an explanation for this behavior is associated with the message payload exchanged between mpi processes . average message with sample data is 4.5 times higher for the case study 2 , suggesting that the time spent in communication is less relevant as long as the recombination rate increases . in figure 3 , the average consumed memory in function of the recombination rate and cores used for the computation is shown , which is registered in case study 1 . in figure 4 , the same exact data are shown , but for the case study 2 . from both figures , we can observe how the node s main memory is exhausted when the maximal recombination rate is reached ( fig . 4d ) , irrespective of the number of involved cores , explaining why the efficiency is impacted given the current experimental setup . the most important contribution for the geneticist is the reduction in the overall running time of the simulation process . in table 1 , the average time hudson s ms , macs , mspar , and msparsm spent running both case studies using different problem sizes ( ie , recombination rate values ) is given . we observe that how the execution time of msparsm is orders of magnitude faster than hudson s ms , faster by a factor than mspar , and significantly better than macs in most configurations , mainly when mutation rate is higher . although the parallel execution of msparsm is behind macs , the scalability of presented msparsm is generally good when adding more computational resources with the possible limitation of the available memory of the nodes . in this work , we developed a new application for parallelizing the hudson s ms , a standard coalescent simulator , coined as msparsm , addressing the sampling of long stretches of dna sequences with large recombination rates , and we showed that this new application has comparable performance with macs when working with large recombination rate , without the side effects of macs . interestingly , the application achieves high parallel efficiency figures ( > 70% ) when working with large recombination rate , suggesting that it could be run on more nodes than used for this work and get better execution times . in addition , our application can outperform macs in most of the cases , even when running the application in one single node . we believe that this new application will facilitate simulating coalescent processes with long genomic regions and large recombination rates in population genomics and evolutionary biology , without compromising the accuracy of the standard coalescence .
the msparsm application is an evolution of mspar , the parallel version of the coalescent simulation program ms , which removes the limitation for simulating long stretches of dna sequences with large recombination rates , without compromising the accuracy of the standard coalescence . 
 this work introduces msparsm , describes its significant performance improvements over mspar and its shared memory parallelization details , and shows how it can get better , if not similar , execution times than macs . two case studies with different mutation rates were analyzed , one approximating the human average and the other approximating the drosophila melanogaster average . 
 source code is available at https://github.com/cmontemuino/msparsm .