Re: Haldane's Dilemma - clarifications - and Felsenstein [LONG]

"Perplexed in Peoria" <jimmenegay@xxxxxxxxxxxxx> wrote in message
news:e7eftb$1ai$1@xxxxxxxxxxxxxxxxxxxxxx

"Joe Felsenstein" <joe@xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx> wrote in message
news:e7cbcg$234r$1@xxxxxxxxxxxxxxxxxxxxxx
I'm not Walter, but do want to comment on the model. It's
a straightforward computer simulation of natural selection of
mutations, with a few wrinkles. I wonder if you don't want to
modify those wrinkles:
1. The alleles are either A or a. In this case the population will,
after enough generations, end with with mostly A's but some a's, as
there is back mutation at a rate equal to forward mutation. So it
will be approaching an asymptotic amount of adaptation. A useful
modification would be to have the alleles be integers which are not
just characters 'A' and 'a' but can be incremented. Thus we would
have
alleles 0, 1, 2, ... and this would allow new alleles to be replaced
by
even better alleles. Mutation would be incrementation of the integer.
Or if you wanted to have back mutation as well mutation could
sometimes
increment and sometimes decrement. The resulting model would gain
favorable mutations (with some but not much loss as well) and this
would go on indefinitely, or at least until you got near MAXINT.

A great idea. If you decide to implement this, Malcolm, don't forget
to introduce a factor of sizeof(int) into your memcpy's in the
recombination code in the makechild function.


Here we go.

Programmers will note a few technical problems introduced by Prof
Felsenstein's changes. These are to do with the fact that computers cannot
represent real numbers exactly, and when you ahve lots of genes the "w" or
fitness measure gets very large.

One interesting thing to note is that, in this model, when the number of
genes is large - I origininally set it to 30000 to represent the number
thought to exist in the human genome, then the number of mutations that make
it into the population is approximately 2 * SELECTIONCOEFFICIENT. When I
tried to speed the program up by a factor of 300 by reducing the numebr of
genes to 100, this relationship disappeared.

Maybe this is worth pursuing.

The number 250000 represents the numebr of generations since the human /
chimp split, or as ReMine would say, the alleged human/chimp split.

However this is only a model. it can be made better, for instance at the
moment males and females are absolutely symmetrial, whilst we know that in
human populations males have more variable fertility. So let's hope Walter
ReMine can, as Joe Felsenstein has done, suggest improvements to show the
effect he wants to show.

/*
Model of favourable mutations spreading through a sexually-
reproducing population.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define NGENES 1000 /* number of genes each creature has */
#define POPSIZE 1000 /* the size of the population */
#define NGENERATIONS 250000 /* number of generations to run simualtion
for */
#define RECOMBINATIONRATE 1.2 /* average number of recombinations per
chromosome */
#define MUTATIONRATE 0.001 /* chance of a mutation, per chromosome */
#define SELECTIONCOEFFICIENT 0.005
typedef struct
{
int matchromosome[NGENES]; /* maternal chromosome */
int patchromosome[NGENES]; /* paternal chromosome */
double fitness; /* fitness score (total of gene scores) */
double cumfitness; /* cumlative fitness for efficient choice */
} CREATURE;

typedef struct
{
CREATURE males[POPSIZE/2]; /* males in population */
CREATURE females[POPSIZE/2]; /* females in population */
} GENERATION;

typedef struct
{
GENERATION gen1; /* generation 1 */
GENERATION gen2; /* generation 2 */
int toggle; /* flag to switch generations between child
and parent */
} WORKSPACE;

/*
generate a binary random number (to fix up technical deficiencies in
rand())
*/
int coinflip(void)
{
return rand() > RAND_MAX/2 ? 1 : 0;
}

/*
generate a random number on 0-1
*/
double uniform()
{
return rand() / (RAND_MAX + 1.0) + rand() / (1.0 + (double) RAND_MAX *
RAND_MAX);
}

int poisson(double v)
{
int answer = 0;
double tot;
double product;

tot = exp(-v);
product = uniform();
while(product > tot)
{
product *= uniform();
answer++;
}

return answer;
}

/*
choose a parent
Params: candidates - list of potential parents
N - number in list
Returns: the selected parent
We take two at random, and return the one with the highest fitness,
or the second to be examined if equal.
*/
int choose(CREATURE *candidates, int N)
{
double target;
int top;
int bottom;
int mid;

target = uniform() * candidates[N-1].cumfitness;

bottom = 0;
top = N-1;
while(bottom != top && bottom != top-1)
{
mid = (top + bottom)/2;
if(candidates[mid].cumfitness < target)
bottom = mid;
else
top = mid;
}
if(candidates[mid].cumfitness < target)
mid++;

if(mid > 0 && candidates[mid-1].cumfitness > target)
mid--;

return mid;
}

/*
generate a child from two parents
Params: child - return for child
mother - parent 1
father - parent 2
The child has a maternal and a paternal chromosome, formed by
random recombination of the parental chromosomes.
This is as in most diploid organisms.
The child's fitness is his total of 'A' genes.
*/
void makechild(CREATURE *child, CREATURE *mother, CREATURE *father)
{
int fitness = 0;
int *chromosome1;
int *chromosome2;
int *temp;
int i;
unsigned char breakpoint[NGENES];
int N;

memset(breakpoint, 0, NGENES);
chromosome1 = mother->matchromosome;
chromosome2 = mother->patchromosome;
if(coinflip())
{
temp = chromosome1;
chromosome1 = chromosome2;
chromosome2 = temp;
}
N = poisson(RECOMBINATIONRATE);
for(i=0;i<N;i++)
breakpoint[rand() % NGENES] = 1;
for(i=0;i<NGENES;i++)
{
child->matchromosome[i] = chromosome1[i];
if(breakpoint[i])
{
temp = chromosome1;
chromosome1 = chromosome2;
chromosome2 = temp;
}
}

memset(breakpoint, 0, NGENES);
chromosome1 = father->matchromosome;
chromosome2 = father->patchromosome;
if(coinflip())
{
temp = chromosome1;
chromosome1 = chromosome2;
chromosome2 = temp;
}
N = poisson(RECOMBINATIONRATE);
for(i=0;i<N;i++)
breakpoint[rand() % NGENES] = 1;
for(i=0;i<NGENES;i++)
{
child->patchromosome[i] = chromosome1[i];
if(breakpoint[i])
{
temp = chromosome1;
chromosome1 = chromosome2;
chromosome2 = temp;
}
}

for(i=0;i<NGENES;i++)
{
fitness += child->matchromosome[i];
fitness += child->patchromosome[i];
}

child->fitness = pow(1 + SELECTIONCOEFFICIENT, fitness);
}

/*
generate a mutation.
A random target gene is flipped from 'a' to 'A' or back
*/
void mutate(GENERATION *gen)
{
int sex;
int chromosome;
int gene;
int id;
int fitness = 0;
CREATURE *creature;
int *target;
int i;


sex = coinflip();
id = rand() % (POPSIZE/2);
chromosome = coinflip();
gene = rand() % NGENES;

if(sex)
creature = &gen->males[id];
else
creature = &gen->females[id];

if(chromosome)
target = &creature->patchromosome[gene];
else
target = &creature->matchromosome[gene];

if( coinflip() )
*target += 1;
else
*target -= 1;


for(i=0;i<NGENES;i++)
{
fitness += creature->matchromosome[i];
fitness += creature->patchromosome[i];
}
creature->fitness = pow(1 + SELECTIONCOEFFICIENT, fitness);
}

/*
generate the next generation
*/
void nextgeneration(GENERATION *children, GENERATION *parents)
{
int i;
int mother;
int father;
int Nmutations;

for(i=0;i<POPSIZE/2;i++)
{
mother = choose(parents->females, POPSIZE/2);
father = choose(parents->males, POPSIZE/2);
makechild(&children->males[i], &parents->females[mother],
&parents->males[father]);
mother = choose(parents->females, POPSIZE/2);
father = choose(parents->males, POPSIZE/2);
makechild(&children->females[i], &parents->females[mother],
&parents->males[father]);

}

Nmutations = poisson(POPSIZE * 2 * MUTATIONRATE);
for(i=0;i<Nmutations;i++)
mutate(children);

children->females[0].cumfitness = children->females[0].fitness;
children->males[0].cumfitness = children->males[0].fitness;
for(i=1;i<POPSIZE/2;i++)
{
children->females[i].cumfitness = children->females[i-1].cumfitness +
children->females[i].fitness;
children->males[i].cumfitness = children->males[i-1].cumfitness +
children->males[i].fitness;
}
}

/*
run the simualtion for one step
(we toggle generations to avoid unnecessary copying of data)
*/
void step(WORKSPACE *wk)
{
if(wk->toggle)
nextgeneration(&wk->gen1, &wk->gen2);
else
nextgeneration(&wk->gen2, &wk->gen1);

wk->toggle = 1 - wk->toggle;
}

/*
set up intial creature with all 'a' alleles
*/
void initialisecreature(CREATURE *creature)
{
int i;

for(i=0;i<NGENES;i++)
{
creature->matchromosome[i] = 0;
creature->patchromosome[i] = 0;
}
creature->fitness = 1.0;
}

/*
set up the initial generations
*/
void initialisegen(GENERATION *gen)
{
int i;

for(i=0;i<POPSIZE/2;i++)
{
initialisecreature(&gen->males[i]);
initialisecreature(&gen->females[i]);
}
for(i=0;i<POPSIZE/2;i++)
{
gen->males[i].cumfitness = i + 1;
gen->females[i].cumfitness = i + 1;
}
}

/*
initialise the workspace
*/
void intialiseworkspace(WORKSPACE *wk)
{
initialisegen(&wk->gen1);
initialisegen(&wk->gen2);
wk->toggle = 0;
}

int main(void)
{
static WORKSPACE workspace;
int i;
int ii;
int score;

intialiseworkspace(&workspace);
for(i=0;i<NGENERATIONS;i++)
{
step(&workspace);
score = 0;
for(ii=0;ii<NGENES;ii++)
{
score += workspace.gen1.males[0].patchromosome[ii];
score += workspace.gen2.males[0].matchromosome[ii];
}
printf("%d: %d\n", i, score);
}

for(i=0;i<POPSIZE/2;i++)
printf("%g\n", workspace.gen1.males[i].fitness);

return 0;
}

--
Buy my book 12 Common Atheist Arguments (refuted)
$1.25 download or $7.20 paper, available www.lulu.com/bgy1mm




.

Relevant Pages