/* * A C program to simulate the population dynamics of three morphs in Papilio dardanus written by Alex Freeman * * */ #include #include /* Declare subroutines */ void PrintHeadings(); void Display(int GEN, const char *comment, double Phh, double Phc, double Pcc, double Pht, double Pct, double Ptt); double sigmoid(double x, double a, double lxmin, double lxmax, double lymin, double lymax); double predation(double Pmodel, double Pmimic, double n_memory); void Normalise(double *a, double *b, double *c, double *d, double *e, double *f); void solve_GF(double H,double C,double T, double *h, double *c, double *t); void solve_MF(double *H, double *C, double *T, double h, double c, double t); double solve_quad(double a, double b, double c); void main(void) { /* declare variables */ double PH, PC, PT; /* frequency hippo, cenea, lamborni - calculated */ double Phh, Phc, Pcc, Pht, Pct, Ptt; /* frequency of each genotype in 1st generation - calculated */ double Ihh, Ihc, Icc, Iht, Ict, Itt; /* frequency of each genotype in next generation before predation */ double Ph, Pc, Pt; /* frequency of each gene in population - defined */ double Mhh, Mhc, Mcc, Mht, Mct, Mtt; /* probability of mating with each genotype - calculated*/ double M1hh, M1hc, M1cc, M1ht, M1ct, M1tt; /* probability of mating with each genotype on first mating - calculated*/ double CH, CC, CT; /* probability of choosing each morph - defined*/ double F, N; /* probability of follow-up mating or not a follow-up mating - defined*/ double S, D; /* probability of choosing same or different morph on follow-up - defined*/ int GEN; /* generation counter */ double N_MEM, Pencounter; /* Predator Memory function and encounter likelihood */ double RmodelH,RmodelC,RmodelT; /* Ratios of model to mimics [at start] */ double AmodelH,AmodelC,AmodelT; /* Abundance of models relative to dardanus morphs - calculated*/ GEN=0; /* probabilities for follow-up matings */ F = 0.59; //number of females mating more than once N = 1-F; S = 0.75; //probability of choosing same morph twice D = 1-S; /* morph proportions */ PH = 70./(70+25+8); // These are Ford's data for race meseres PC = 25./(70+25+8); PT = 8./(70+25+8); /* probabilities of choosing each morph on first mating - can be varied*/ CH=CC=CT=1; // Relative abundance of mimics & models => predation Pencounter=0.95; // Chance of encountering a predator before mating (VARIABLE) N_MEM=40; // Memory function of predators (VARIABLE) RmodelH=104/70.; // Model abundance (initially) - data from Ford for meseres RmodelC=232/8.; RmodelT=781/25.; /* initial gene proportions */ solve_GF(CH,CC,CT,&Ph,&Pc,&Pt); /* calculate genotype proportions from gene proportions */ Phh = Ph*Ph; Phc = 2*Ph*Pc; Pcc = Pc*Pc; Pht = 2*Ph*Pt; Pct = 2*Pc*Pt; Ptt = Pt*Pt; /* print out initial proportions */ PrintHeadings(); Display(0,"Initial proprtions",Phh,Phc,Pcc,Pht,Pct,Ptt); // // Generation LOOP // for (GEN=1;GEN<100;GEN++) { //PH,PC and PT are current morph proportions //Compute these from the genotypes Ph=Phh; Pc=Phc+Pcc; Pt=Pht+Pct+Ptt; solve_MF(&PH,&PC,&PT,Ph,Pc,Pt); /* calculate probabilities of mating with each genotype */ /* 1st mating */ M1hh = Phh*PH; M1hc = Phc*PC; M1cc = Pcc*PC; M1ht = Pht*PT; M1ct = Pct*PT; M1tt = Ptt*PT; /* Normalise */ Normalise(&M1hh,&M1hc,&M1cc,&M1ht,&M1ct,&M1tt); /* 2nd mating */ Mhh = Phh*(S*M1hh + D*(M1hc+M1cc) + D*(M1tt+M1ht+M1ct)); Mhc = Phc*(S*(M1hc+M1cc) + D*M1hh + D*(M1tt+M1ht+M1ct)); Mcc = Pcc*(S*(M1hc+M1cc) + D*M1hh + D*(M1tt+M1ht+M1ct)); Mht = Pht*(S*(M1tt+M1ht+M1ct) + D*M1hh + D*(M1hc+M1cc)); Mct = Pct*(S*(M1tt+M1ht+M1ct) + D*M1hh + D*(M1hc+M1cc)); Mtt = Ptt*(S*(M1tt+M1ht+M1ct) + D*M1hh + D*(M1hc+M1cc)); /* Normalise */ Normalise(&Mhh,&Mhc,&Mcc,&Mht,&Mct,&Mtt); /* Both matings together */ Mhh=N*M1hh + F*Mhh; Mhc=N*M1hc + F*Mhc; Mcc=N*M1cc + F*Mcc; Mct=N*M1ct + F*Mct; Mht=N*M1ht + F*Mht; Mtt=N*M1tt + F*Mtt; /* // print out probabilities of mating with each genotype if required Display(GEN,"Mating probability",Mhh,Mhc,Mcc,Mht,Mct,Mtt); */ /* calculate the proportion of each genotype in next generation */ /* based on gene dominance hierarchy */ Ihh = Mhh*(Phh + 0.5*Phc + 0.5*Pht) + 0.5*Mhc*(Phh + 0.5*Phc + 0.5*Pht) + 0.5*Mht*(Phh + 0.5*Phc + 0.5*Pht); Ihc = Mhh*(Pcc + 0.5*Phc + 0.5*Pct) + 0.5*Mhc*(Phh + Phc + Pcc + 0.5*Pct + 0.5*Pht) + Mcc*(Phh + 0.5*Phc + 0.5*Pht) + 0.5*Mht*(0.5*Phc + Pcc + 0.5*Pct) + 0.5*Mct*(Phh + 0.5*Phc + 0.5*Pht); Icc = 0.5*Mhc*(0.5*Phc + Pcc + 0.5*Pct) + Mcc*(0.5*Phc + Pcc + 0.5*Pct) + 0.5*Mct*(0.5*Phc + Pcc + 0.5*Pct); Iht = Mhh*(0.5*Pht + 0.5*Pct + Ptt) + 0.5*Mhc*(0.5*Pht + 0.5*Pct + Ptt) + 0.5*Mht*(Phh + 0.5*Phc + Pht + 0.5*Pct + Ptt) + 0.5*Mct*(Phh + 0.5*Phc + 0.5*Pht) + Mtt*(Phh + 0.5*Phc + 0.5*Pht); Ict = 0.5*Mhc*(0.5*Pht + 0.5*Pct + Ptt) + Mcc*(0.5*Pht + 0.5*Pct + Ptt) + 0.5*Mht*(0.5*Phc + Pcc + 0.5*Pct) + Mtt*(0.5*Phc + Pcc + 0.5*Pct) + 0.5*Mct*(Ptt + Pct + 0.5*Pht + Pcc + 0.5*Phc); Itt = 0.5*Mht*(0.5*Pht + 0.5*Pct + Ptt) + 0.5*Mct*(0.5*Pht + 0.5*Pct + Ptt) + Mtt*(0.5*Pht + 0.5*Pct + Ptt); /* Normalise */ Normalise(&Ihh,&Ihc,&Icc,&Iht,&Ict,&Itt); /* PREDATION */ /* // Sigmoidal predation: not used Phh= Ihh*sigmoid(Ihh/2,a,xmin,xmax,ymin,ymax); Phc= Ihc*sigmoid((Ihc+Icc)/2,a,xmin,xmax,ymin,ymax); Pcc= Icc*sigmoid((Ihc+Icc)/2,a,xmin,xmax,ymin,ymax); Pht= Iht*sigmoid((Iht+Ict+Itt)/2,a,xmin,xmax,ymin,ymax); Pct= Ict*sigmoid((Iht+Ict+Itt)/2,a,xmin,xmax,ymin,ymax); Ptt= Itt*sigmoid((Iht+Ict+Itt)/2,a,xmin,xmax,ymin,ymax); */ /* predation function for Batesian mimics - Huheey's function used*/ AmodelH=RmodelH*CH/(Ihh); AmodelC=RmodelC*CC/(Ihc+Icc); AmodelT=RmodelT*CT/(Iht+Ict+Itt); Phh= Ihh*((1-Pencounter) + Pencounter*predation(AmodelH,1-AmodelH,N_MEM)); Phc= Ihc*((1-Pencounter) + Pencounter*predation(AmodelC,1-AmodelC,N_MEM)); Pcc= Icc*((1-Pencounter) + Pencounter*predation(AmodelC,1-AmodelC,N_MEM)); Pht= Iht*((1-Pencounter) + Pencounter*predation(AmodelT,1-AmodelT,N_MEM)); Pct= Ict*((1-Pencounter) + Pencounter*predation(AmodelT,1-AmodelT,N_MEM)); Ptt= Itt*((1-Pencounter) + Pencounter*predation(AmodelT,1-AmodelT,N_MEM)); /* predation of males */ Phh+=Ihh*((1-Pencounter) + Pencounter*predation(0,1,N_MEM)); Phc+=Ihc*((1-Pencounter) + Pencounter*predation(0,1,N_MEM)); Pcc+=Icc*((1-Pencounter) + Pencounter*predation(0,1,N_MEM)); Pht+=Iht*((1-Pencounter) + Pencounter*predation(0,1,N_MEM)); Pct+=Ict*((1-Pencounter) + Pencounter*predation(0,1,N_MEM)); Ptt+=Itt*((1-Pencounter) + Pencounter*predation(0,1,N_MEM)); /* Normalise */ Normalise(&Phh,&Phc,&Pcc,&Pht,&Pct,&Ptt); /* print out frequencies from subsequent generations */ Display(GEN,"Gene frequencies",Phh,Phc,Pcc,Pht,Pct,Ptt); }//-----Loop Ends } // // Subroutines // void PrintHeadings() { /* print out column headings */ printf("Generation\t\t"); printf("hh\t"); printf("hc\t"); printf("cc\t"); printf("ht\t"); printf("ct\t"); printf("tt\t"); printf("\t"); printf("hippo gene\t"); printf("cenea gene\t"); printf("trophonius gene\t"); printf("\t"); printf("hippo morph\t"); printf("cenea morph\t"); printf("trophonius morph\n"); } void Display(int GEN, const char *comment, double Phh, double Phc, double Pcc, double Pht, double Pct, double Ptt) { printf("%d\t", GEN); printf("%s\t", comment); printf("%f\t", Phh); printf("%f\t", Phc); printf("%f\t", Pcc); printf("%f\t", Pht); printf("%f\t", Pct); printf("%f\t", Ptt); printf("\t"); printf("%f\t", (2*Phh+Phc+Pht)/2); printf("%f\t", (2*Pcc+Phc+Pct)/2); printf("%f\t", (2*Ptt+Pht+Pct)/2); printf("\t"); printf("%f\t", Phh); printf("%f\t", Pcc+Phc); printf("%f\n", Ptt+Pht+Pct); } /* Subroutine to normalise morph frequencies */ void Normalise(double *a, double *b, double *c, double *d, double *e, double *f) { double sum = *a + *b + *c + *d + *e + *f; *a/=sum; *b/=sum; *c/=sum; *d/=sum; *e/=sum; *f/=sum; } /* Batesian mimicry predation function */ double predation(double Pmodel, double Pmimic, double n_memory) { return 1-(1/(n_memory*Pmodel+Pmimic)); } /* Subroutine to calculate a sigmoidal predation function - not eventually used in this model */ double sigmoid(double xdash, double a, double lxmin, double lxmax, double lymin, double lymax) { double fmax, fmin, x; /* function for predator selection */ fmin=exp(a*lxmin)/(1+exp(a*lxmin)); fmax=exp(a*lxmax)/(1+exp(a*lxmax)); x=xdash*(lxmax-lxmin)+lxmin; return (1- (lymin+ (lymax-lymin)*(exp(a*x)/(1+exp(a*x)) -fmin)/(fmax-fmin))); } /* Subroutine to solve gene frequencies from morph frequencies */ void solve_GF(double H,double C,double T, double *h, double *c, double *t) { *h=sqrt(H); *c=solve_quad(1,2 * *h,-C); *t=solve_quad(1,2 * (*h+*c),-T); } /* Subroutine to solve morph frequencies from gene frequencies */ void solve_MF(double *H, double *C, double *T, double h, double c, double t) { *H=h*h; *C=c*c + 2*h*c; *T=t*t + 2*h*t + 2*c*t; } double solve_quad(double a, double b, double c) { //return the positive root of the quadratic return (sqrt(b*b-4*a*c)-b)/(2*a) ; }