Post date: Feb 17, 2017 5:1:45 AM
The purpose of this model is to ask how connectivity affects the extent to which demes differentiate (adapt) in a deme-by-deme vs. all demes at once (system) manner. Here is a short description of the model.
General conditions. We model a 10x10 matrix of demes. Each deme has a set population size (currently 100, this can be changed). Two habitat types exist. These can currently be arranged in one of three spatial configurations: sharp ecotone, a checkerboard (patchwork) or in a random arrangement. We monitor individuals. Each individual has a genome with a set of L loci that affect fitness. These are either unlinked or equally spaced on a single chromosome (so far I have been using the latter option). Each locus is bi-allelic. One allele is adaptive in one habitat type, the other in the other.
Starting conditions. We currently begin the simulations with all demes fixed for alleles adapted to one habitat type (no standing variation, this isn't the only option).
Simulations. Each generation consists of the following steps. Dispersal occurs probabilistically based on the migration rate (which we set and vary as a key parameter). The migration rate gives the probability that an individual disperses and the direction (up, down, left or right) is chosen at random. Dispersal out of the metapopulation is allowed and N can exceed 100 (carrying capacity) temporarily. We then have viability selection. Fitness is multiplicative, such that w = (1-s)^n, where n is the number of maladaptive alleles (across loci) and s is the selection coefficient (which is the same for all loci). w is the survival probability. Survival is then determine randomly. This is followed by reproduction for all remaining individuals in each deme. Pairs of parents are chosen at random to produce an offspring until the population is back to carrying capacity (N = 100 so far). An initial chromosome (gene copy) is chosen for each parent to give to the offspring for the first locus. We then walk along the chromosome allowing for mutations with probability mu (1e-6 so far) and recombination (assuming a chromosome size of 1 M). Double (and more) crossovers are allowed.
Data and stopping conditions. Simulations are run for a fixed amount of time. We record several key parameters. First, the mean allele frequency of alleles adaptive in the novel patch type for each deme. Note that this is a mean across loci (thus 0.5 could mean 0.5 for all loci or 0 for half 1 for half). This is our general measures of adaptedness with the expectation of 0 and 1 in the ancestral and novel patch type (respectively) for perfect adaptation. We also record the mean and s.d. in this mean allele frequency for all demes of a given patch type. This gets at overall adaptedness of demes of each type (mean) and the variability among demes in their degree of adaptedness for demes from each patch type (s.d.). The latter is the deme-by-demeness.
Here is the source code as of Feb. 16, 2017:
main.C
// file: main.C
// metapop - runs forward-time, individual-based simulations of
// metapopulations with dispersal, selection and recombination
// Time-stamp: <Friday, 27 January 2017, 20:50 MST -- zgompert>
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <time.h>
#include <getopt.h>
#include "metapop.H"
using namespace std;
gsl_rng * r; /* global state variable for random number generator */
/* ----------------- */
/* beginning of main */
/* ----------------- */
int main(int argc, char *argv[]) {
time_t start = time(NULL);
time_t end;
int rng_seed = 0, ch;
int gen;
ofstream os1; // outfile streams
ofstream os2;
// get command line arguments
if (argc < 2) {
usage(argv[0]);
}
// define main variables
parameters parms;
population MetaPop[DIM][DIM];
while ((ch = getopt(argc, argv, "g:k:p:l:m:s:r:f:w:v:u:")) != -1){
switch(ch){
case 'g':
parms.Ngen = atoi(optarg);
break;
case 'k':
parms.K = atoi(optarg);
break;
case 'p':
parms.patch = atoi(optarg);
break;
case 'l':
parms.L = atoi(optarg);
break;
case 'm':
parms.m = atof(optarg);
break;
case 's':
parms.s = atof(optarg);
break;
case 'r':
parms.linked = atoi(optarg);
break;
case 'f':
parms.ofile = optarg;
break;
case 'w':
parms.fitmod = atoi(optarg);
break;
case 'v':
parms.sgv = atoi(optarg);
break;
case 'u':
parms.mu = atof(optarg);
break;
case '?':
default:
usage(argv[0]);
}
}
// open outfiles
parms.o1 = string("summary_") + parms.ofile;
os1.open(parms.o1.c_str());
parms.o2 = string("demedata_") + parms.ofile;
os2.open(parms.o2.c_str());
// set up gsl random number generation
gsl_rng_env_setup();
r = gsl_rng_alloc(gsl_rng_default);
srand(time(NULL));
rng_seed = rand();
gsl_rng_set(r, rng_seed); /* seed gsl_rng with output of rand, which
was seeded with result of time(NULL) */
// create the metapopulation
makepop(parms, MetaPop);
// main loop over generations
for(gen=0; gen<parms.Ngen; gen++){
updatepop(parms, MetaPop);
// write results (not every generation)
if((gen % 10) == 0) // 10 gen write frequency
write(parms, MetaPop, gen, &os1, &os2);
}
// close outfiles
os1.close();
os2.close();
// prints run time
end = time(NULL);
cerr << "Runtime: " << (end-start)/3600 << " hr " << (end-start)%3600/60 << " min ";
cerr << (end-start)%60 << " sec" << endl;
return 0;
}
func.C
#include <iostream>
#include <sstream>
#include <fstream>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics_double.h>
#include <float.h>
#include <math.h>
#include "metapop.H"
using namespace std;
void usage(char * name){
fprintf(stderr,"\n%s version %s\n\n", name, VERSION);
fprintf(stderr, "Usage: metapop [options]\n");
fprintf(stderr, "-g Number of generations\n");
fprintf(stderr, "-k Carrying capacity or constant population size\n");
fprintf(stderr, "-p Patch organization: 0 = coarse grain, 1 = fine grain, 2 = random\n");
fprintf(stderr, "-l Number of causal loci\n");
fprintf(stderr, "-m Migration rate between neighboring demes\n");
fprintf(stderr, "-s Per allele selection coefficient\n");
fprintf(stderr, "-r Genetic arch.: 1 = SNPs liked, 0 = free recombination\n");
fprintf(stderr, "-v Standing genetic variation: 1 = true, 0 = false\n");
fprintf(stderr, "-u Mutation rate per locus\n");
fprintf(stderr, "-w Fitness function: 0 = additive, 1 = multiplicitive\n");
fprintf(stderr, "-f Outfile prefix\n");
exit(1);
}
// creates the metapopulation
void makepop(parameters parms, population MetaPop[][DIM]){
int i, j, n, l;
int a;
double x;
int g;
// new metapopulation
for(i=0; i<DIM; i++){
a = i % 2; // used to alternate patch types
for(j=0; j<DIM; j++){
// set popsize to carrying capacity
MetaPop[i][j].Nind = parms.K;
// set patch types
if(parms.patch==0){ // two big patches
if(j < ((double) (DIM-1)/2.0))
MetaPop[i][j].type = 0;
else
MetaPop[i][j].type = 1;
}
else if (parms.patch==1){
MetaPop[i][j].type = a;
a = abs(a - 1);
}
else if (parms.patch==2){
x = gsl_ran_flat(r, 0, 1);
if (x < 0.5)
MetaPop[i][j].type = 0;
else
MetaPop[i][j].type = 1;
}
// create genotype matrix
MetaPop[i][j].G1 = gsl_matrix_int_calloc(parms.K, parms.L);
MetaPop[i][j].G2 = gsl_matrix_int_calloc(parms.K, parms.L);
for(n=0; n<parms.K; n++){
for(l=0; l<parms.L; l++){
if(parms.sgv == 1){ // standing variation p = 0.5
x = gsl_ran_flat(r, 0,1);
g = roundZeroOne(x);
gsl_matrix_int_set(MetaPop[i][j].G1, n, l, g);
x = gsl_ran_flat(r, 0,1);
g = roundZeroOne(x);
gsl_matrix_int_set(MetaPop[i][j].G2, n, l, g);
}
else{
gsl_matrix_int_set(MetaPop[i][j].G1, n, l, 0);
gsl_matrix_int_set(MetaPop[i][j].G2, n, l, 0);
}
}
}
// allocate memory for additional variables in population
MetaPop[i][j].G1temp = gsl_matrix_int_calloc(parms.K * BUFFER, parms.L);
MetaPop[i][j].G2temp = gsl_matrix_int_calloc(parms.K * BUFFER, parms.L);
MetaPop[i][j].gtemp = gsl_vector_int_calloc(parms.L);
MetaPop[i][j].fitness = gsl_vector_calloc(parms.K * BUFFER);
}
}
}
// wrapper function for one generation of dispersal, selection and reproduction
void updatepop(parameters parms, population MetaPop[][DIM]){
// dispersal
dispersal(parms, MetaPop);
// selection
selection(parms, MetaPop);
// reproduction
reproduction(parms, MetaPop);
// mutation
if(parms.mu > 0) // allow mutations at this rate
mutation(parms, MetaPop);
}
// disperal function, probalistic dispersal between neighbors, moving off edges means death
void dispersal(parameters parms, population MetaPop[][DIM]){
int i, j, n;
int ii, jj;
double p[4] = {0.25, 0.25, 0.25, 0.25};
unsigned int dirv[4];
int dir;
// reset counters
setZero(MetaPop);
for(i=0; i<DIM; i++){
for(j=0; j<DIM; j++){
for(n=0; n<MetaPop[i][j].Nind; n++){
// probabilstic migration, directions determine by multinomial draw
if(gsl_ran_flat(r,0,1) < parms.m){
gsl_ran_multinomial(r, 4, 1, p, dirv);
dir = whichOne(4, dirv);
// get new coordinates
if(dir == 0){
ii = i;
jj = j - 1;
}
else if(dir == 1){
ii = i - 1;
jj = j;
}
else if(dir == 2){
ii = i;
jj = j + 1;
}
else {
ii = i + 1;
jj = j;
}
// assign genotype
if(ii >= 0 && ii < DIM && jj >= 0 && jj < DIM){ // didn't fall off the world
gsl_matrix_int_get_row(MetaPop[i][j].gtemp, MetaPop[i][j].G1, n);
gsl_matrix_int_set_row(MetaPop[ii][jj].G1temp, MetaPop[ii][jj].cntr, MetaPop[i][j].gtemp);
gsl_matrix_int_get_row(MetaPop[i][j].gtemp, MetaPop[i][j].G2, n);
gsl_matrix_int_set_row(MetaPop[ii][jj].G2temp, MetaPop[ii][jj].cntr, MetaPop[i][j].gtemp)
;
MetaPop[ii][jj].cntr++;
}
}
// non-migrant, move to temp matrix for the current pop.
else{
gsl_matrix_int_get_row(MetaPop[i][j].gtemp, MetaPop[i][j].G1, n);
gsl_matrix_int_set_row(MetaPop[i][j].G1temp, MetaPop[i][j].cntr, MetaPop[i][j].gtemp);
gsl_matrix_int_get_row(MetaPop[i][j].gtemp, MetaPop[i][j].G2, n);
gsl_matrix_int_set_row(MetaPop[i][j].G2temp, MetaPop[i][j].cntr, MetaPop[i][j].gtemp);
MetaPop[i][j].cntr++;
}
}
}
}
}
// hard viability selection, this is after migration
void selection(parameters parms, population MetaPop[][DIM]){
int i, j, n;
for(i=0; i<DIM; i++){
for(j=0; j<DIM; j++){
MetaPop[i][j].scntr = 0;
for(n=0; n<MetaPop[i][j].cntr; n++){
// determine absolute fitness
if(parms.fitmod == 0) // additive
calcfitness(parms, MetaPop, i, j, n);
else if(parms.fitmod == 1) // multiplicitive
calcmultfitness(parms, MetaPop, i, j, n);
// lives or dies
death(parms, MetaPop, i, j, n);
}
//cout << "pop = " << i << "," << j << "; start = " << MetaPop[i][j].cntr << "; end = " << MetaPop[i][j].scntr << endl;
}
}
}
// calculate fitness with additive directional selection
void calcfitness(parameters parms, population MetaPop[][DIM], int i, int j, int n){
int l;
double x = 0;
double w;
for(l=0; l<parms.L; l++){
x += gsl_matrix_int_get(MetaPop[i][j].G1temp, n, l);
x += gsl_matrix_int_get(MetaPop[i][j].G2temp, n, l);
}
if(MetaPop[i][j].type == 0) // 0 alleles more fit, w = 1 - sum(g=1) s
w = 1.0 - (parms.s * x);
else
w = 1.0 - (parms.s * (2.0 * parms.L - x));
if(w < 0)
w = 0; // ensure fitness is non-negative
gsl_vector_set(MetaPop[i][j].fitness, n, w);
}
// calculate fitness with multiplicitive directional selection
void calcmultfitness(parameters parms, population MetaPop[][DIM], int i, int j, int n){
int l;
int x = 0;
double w;
for(l=0; l<parms.L; l++){
x += gsl_matrix_int_get(MetaPop[i][j].G1temp, n, l);
x += gsl_matrix_int_get(MetaPop[i][j].G2temp, n, l);
}
if(MetaPop[i][j].type == 0) // 0 alleles more fit
w = gsl_pow_int(1.0 - parms.s, x);
else
w = gsl_pow_int(1.0 - parms.s, 2 * (int) parms.L - x);
gsl_vector_set(MetaPop[i][j].fitness, n, w);
}
// die, probablistically, and adjust temp genotype objects accordingly
void death(parameters parms, population MetaPop[][DIM], int i, int j, int n){
double w;
w = gsl_vector_get(MetaPop[i][j].fitness, n);
if(gsl_ran_flat(r, 0, 1) < w){ // survive
gsl_matrix_int_swap_rows(MetaPop[i][j].G1temp, n, MetaPop[i][j].scntr);
gsl_matrix_int_swap_rows(MetaPop[i][j].G2temp, n, MetaPop[i][j].scntr);
MetaPop[i][j].scntr++;
}
// die, nothing to do, will eventually be erased
}
// reproduce by chosing pairs of parents from survivors followed by recombination and gamete sampling
void reproduction(parameters parms, population MetaPop[][DIM]){
int i, j, n;
int p1, p2;
for(i=0; i<DIM; i++){
for(j=0; j<DIM; j++){
for(n=0; n<parms.K; n++){
p1 = gsl_rng_uniform_int(r, MetaPop[i][j].scntr); // index for parent 1
p2 = gsl_rng_uniform_int(r, MetaPop[i][j].scntr); // index for parent 2
if(parms.linked==1){
sampleLinked(parms, MetaPop, i, j, n, p1, p2);
}
else if(parms.linked==0){
sampleUnlinked(parms, MetaPop, i, j, n, p1, p2);
}
}
}
}
}
// mutations convert alleles from 0 -> 1 or 1 -> 0
void mutation(parameters parms, population MetaPop[][DIM]){
int l, g, i, j, n;
double x;
for(i=0; i<DIM; i++){
for(j=0; j<DIM; j++){
for(n=0; n<parms.K; n++){
for(l=0; l<parms.L; l++){
x = gsl_ran_flat(r, 0,1);
if(x < parms.mu){
g = abs(1 - gsl_matrix_int_get(MetaPop[i][j].G1, n, l));
gsl_matrix_int_set(MetaPop[i][j].G1, n, l, g);
}
x = gsl_ran_flat(r, 0,1);
if(x < parms.mu){
g = abs(1 - gsl_matrix_int_get(MetaPop[i][j].G2, n, l));
gsl_matrix_int_set(MetaPop[i][j].G2, n, l, g);
}
}
}
}
}
}
// sample gene copies and create genetic data when markers are unlinked (free recombination)
void sampleLinked(parameters parms, population MetaPop[][DIM], int i, int j, int n, int p1, int p2){
int l, g;
int brk, chrom;
double x;
// p1, sample initial gene copy and breakpoint
brk = gsl_rng_uniform_int(r, (parms.L-1)) + 1; // random int between 1 and L-2
x = gsl_ran_flat(r, 0,1);
chrom = roundZeroOne(x);
for(l=0; l<parms.L; l++){
if(chrom==0){
g = gsl_matrix_int_get(MetaPop[i][j].G1temp, p1, l);
gsl_matrix_int_set(MetaPop[i][j].G1, n, l, g);
}
else {
g = gsl_matrix_int_get(MetaPop[i][j].G2temp, p1, l);
gsl_matrix_int_set(MetaPop[i][j].G1, n, l, g);
}
if(l==brk)
chrom = abs(chrom-1);
}
// p2, sample initial gene copy and breakpoint
brk = gsl_rng_uniform_int(r, (parms.L-1)) + 1; // random int between 1 and L-2
x = gsl_ran_flat(r, 0,1);
chrom = roundZeroOne(x);
for(l=0; l<parms.L; l++){
if(chrom==0){
g = gsl_matrix_int_get(MetaPop[i][j].G1temp, p2, l);
gsl_matrix_int_set(MetaPop[i][j].G2, n, l, g);
}
else {
g = gsl_matrix_int_get(MetaPop[i][j].G2temp, p2, l);
gsl_matrix_int_set(MetaPop[i][j].G2, n, l, g);
}
if(l==brk)
chrom = abs(chrom-1);
}
}
// sample gene copies and create genetic data when markers are unlinked (free recombination)
void sampleUnlinked(parameters parms, population MetaPop[][DIM], int i, int j, int n, int p1, int p2){
int l;
int g;
double x;
for(l=0; l<parms.L; l++){
// chose parent 1 allele
x = gsl_ran_flat(r, 0,1);
if(x < 0.5){ // first gene copy
g = gsl_matrix_int_get(MetaPop[i][j].G1temp, p1, l);
gsl_matrix_int_set(MetaPop[i][j].G1, n, l, g);
}
else {
g = gsl_matrix_int_get(MetaPop[i][j].G2temp, p1, l);
gsl_matrix_int_set(MetaPop[i][j].G1, n, l, g);
}
// chose parent 2 allele
x = gsl_ran_flat(r, 0,1);
if(x < 0.5){ // first gene copy
g = gsl_matrix_int_get(MetaPop[i][j].G1temp, p2, l);
gsl_matrix_int_set(MetaPop[i][j].G2, n, l, g);
}
else {
g = gsl_matrix_int_get(MetaPop[i][j].G2temp, p2, l);
gsl_matrix_int_set(MetaPop[i][j].G2, n, l, g);
}
}
}
// write population allele frequencies
void write(parameters parms, population MetaPop[][DIM], int gen, ofstream * os1, ofstream * os2){
int i, j, n, l, a;
double x, sum, mn, sig, min, max;
double y[2], ysum[2];
int cnts[2];
gsl_vector * pvec;
gsl_vector * pvec1;
gsl_vector * pvec2;
int cnt = 0;
for(n=0; n<2; n++){
cnts[n] = 0;
y[n] = 0;
ysum[n] = 0;
}
// number of each type
for(i=0; i<DIM; i++){
for(j=0; j<DIM; j++){
if(MetaPop[i][j].type == 0)
cnts[0]++;
else
cnts[1]++;
}
}
pvec1 = gsl_vector_calloc(cnts[0]);
pvec2 = gsl_vector_calloc(cnts[1]);
pvec = gsl_vector_calloc(DIM * DIM);
for(n=0; n<2; n++){
cnts[n] = 0;
}
*os2 << gen;
for(i=0; i<DIM; i++){
for(j=0; j<DIM; j++){
x = 0;
sum = 0;
a = MetaPop[i][j].type; // 0 or 1
for(l=0; l<parms.L; l++){
for(n=0; n<MetaPop[i][j].Nind; n++){
sum += 2;
x += gsl_matrix_int_get(MetaPop[i][j].G1, n, l);
x += gsl_matrix_int_get(MetaPop[i][j].G2, n, l);
// type specific sums
ysum[a] += 2;
y[a] += gsl_matrix_int_get(MetaPop[i][j].G1, n, l);
y[a] += gsl_matrix_int_get(MetaPop[i][j].G2, n, l);
}
}
x = x/sum;
gsl_vector_set(pvec, cnt, x);
*os2 << " " << x;
cnt++;
// type specific sums
y[a] = y[a]/ysum[a];
if(a==0)
gsl_vector_set(pvec1, cnts[a], y[a]);
else
gsl_vector_set(pvec2, cnts[a], y[a]);
cnts[a]++;
y[a] = 0;
ysum[a] = 0;
}
}
*os2 << endl;
// summary stats
// mn = gsl_stats_mean(pvec->data, 1, DIM * DIM);
// sig = gsl_stats_sd(pvec->data, 1, DIM * DIM);
// min = gsl_stats_min(pvec->data, 1, DIM * DIM);
// max = gsl_stats_max(pvec->data, 1, DIM * DIM);
// *os1 << mn << " " << sig << " " << min << " " << max << endl;
// typs summary stats
mn = gsl_stats_mean(pvec1->data, 1, cnts[0]);
sig = gsl_stats_sd(pvec1->data, 1, cnts[0]);
min = gsl_stats_min(pvec1->data, 1, cnts[0]);
max = gsl_stats_max(pvec1->data, 1, cnts[0]);
*os1 << mn << " " << sig << " " << min << " " << max << " ";
mn = gsl_stats_mean(pvec2->data, 1, cnts[1]);
sig = gsl_stats_sd(pvec2->data, 1, cnts[1]);
min = gsl_stats_min(pvec2->data, 1, cnts[1]);
max = gsl_stats_max(pvec2->data, 1, cnts[1]);
*os1 << mn << " " << sig << " " << min << " " << max << endl;
// free memory
gsl_vector_free(pvec);
gsl_vector_free(pvec1);
gsl_vector_free(pvec2);
}
// ************ little helper functions *********** //
// function to round to 0 or 1
int roundZeroOne(double x){
int g;
if (x > 0.5)
g = 1;
else
g = 0;
return g;
}
// reset population counters to zeros
void setZero(population MetaPop[][DIM]){
int i, j;
for(i=0; i<DIM; i++){
for(j=0; j<DIM; j++){
MetaPop[i][j].cntr = 0;
}
}
}
// finds the array element with value 1
int whichOne(int N, unsigned int x[]){
int n = 0;
N = -9;
while(N < 0){
if(x[n] == 1)
N = n;
n++;
}
return N;
}
metapop.H
#include <iostream>
#include <sstream>
#include <fstream>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
// macro definitions
#define VERSION "1 -- Dec. 2016"
#define DIM 10
#define BUFFER 3
extern gsl_rng * r;
// function declarations
using namespace std;
// struct definitions
struct parameters{
int Ngen; // number of generations
int K; // carrying capacity, or consant N with soft selection (the only option for now)
int L; // number of SNP loci
int linked; // 1 = SNPs evenly spaced on one chromosome, 0 = free recombination
int patch; // defines organization of habitat pathes within metapopulation
// 0 = type 0 for columns 0:(DIM-1)/2, type 1 for others
// 1 = every other row and column has a different patch type
double m; // migration rate between neighboring demes/populaitons on the lattice
int fitmod; // 0 = additive, 1 = multiplicitive
double s; // +- additive/multiplictive fitness effect (per allele)
int sgv; // start with standing genetic variation, 1 = true, 0 = false
double mu; // mutation rate
char * ofile; // base for outfile names
string o1; // full outfile names
string o2;
};
struct population{
int Nind; // number of individuals in a population
gsl_matrix_int * G1; // genotype matrix 1 of 2, rows = inds, cols = snps
gsl_matrix_int * G2; // genotype matrix 2 of 2, rows = inds, cols = snps
// temps are larger, values before selection
gsl_matrix_int * G1temp; // genotype matrix 1 of 2, rows = inds, cols = snps
gsl_matrix_int * G2temp; // genotype matrix 2 of 2, rows = inds, cols = snps
gsl_vector_int * gtemp; // tempoary genotype vector
gsl_vector * fitness; // vector of fitness values for each individual
int type; // habitat type, defines whether the 0 or 1 allele has higher fitness (0 = 0, 1 = 1)
int cntr;
int scntr;
};
// function definitions
void usage(char *);
void makepop(parameters parms, population MetaPop[][DIM]);
void updatepop(parameters parms, population MetaPop[][DIM]);
void dispersal(parameters parms, population MetaPop[][DIM]);
void selection(parameters parms, population MetaPop[][DIM]);
void calcfitness(parameters parms, population MetaPop[][DIM], int i, int j, int n);
void calcmultfitness(parameters parms, population MetaPop[][DIM], int i, int j, int n);
void death(parameters parms, population MetaPop[][DIM], int i, int j, int n);
void reproduction(parameters parms, population MetaPop[][DIM]);
void mutation(parameters parms, population MetaPop[][DIM]);
void sampleLinked(parameters parms, population MetaPop[][DIM], int i, int j, int n, int p1, int p2);
void sampleUnlinked(parameters parms, population MetaPop[][DIM], int i, int j, int n, int p1, int p2);
void write(parameters parms, population MetaPop[][DIM], int gen, ofstream * os1, ofstream * os2);
int roundZeroOne(double x);
void setZero(population MetaPop[][DIM]);
int whichOne(int N, unsigned int x[]);