#include #include #include #include #include #define SEQUENCELENGTH 16 #define SPECIESPERCLADE 16777216 #define STARTINGSPECIES 4096 #define STARTINGCLADES 256 #define MAXLINELENGTH 66536 //Define lineage elements (species and higher taxa) typedef struct species { int index; //This is just an index int ancestor; //The index of the species this species is derived from int extant; //0 if the species is extinct, 1 if not extinct float extinction; //A variable measuring the probability of the species' extinction in a given timestep float speciation; //A variable measuring the probability of the species giving rise to new species in a given timestep float evolution; //The rate of evolution of the species' gene sequence over its history int speciesage; //The number of timesteps in which the species has existed int highertaxon; //The index of the higher taxon (clade) to which the species belongs int pollination; //0 if the species is wind-pollinated, 1 if insect-pollinated int dispersal; //0 if the species is wind-dispersed, 1 if vertebrate-dispersed float range; //Geographic range of a species, in arbitrary units char Sequence[SEQUENCELENGTH]; //Sequence of 512 nucleotides, to be modeled each timestep int couplet; //Couplet is used for tree building int *Descendents; //Descendents is for the finding of couplets, the list of descendent species int numdescendents; //numdescendents keeps track of the number of descendents. int currentlyancestral; //Used each timestep to see if it's currently ancestral to a living species. } species; typedef struct clade { int index; //Also just an index char Name[MAXLINELENGTH]; //Name of the clade - for the human to keep track of it with. int ancestor; //Index of the species this clade is derived from or sister to int ancestorclade; //Index of the clade to which the ancestral species belongs (yes, that "clade" is paraphyletic if this one is excluded from it - but this is how evolution happens) int specieslist[SPECIESPERCLADE]; //Array containing the indeces of all of the species in the clade; if it gets over 256, create a new clade int diversity; //Number of species in the clade, so we don't run off the end of specieslist. int standingdiversity; //Number of species still surviving in the clade. int extant; //0 if the entire clade is extinct, 1 if at least some of its species are still extant float extinction; //A variable measuring the rate of species extinction probability over the entire clade (used in calculating probabilities for new species) float speciation; //As above, but for extinction float evolution; //As above, but for sequence evolution float cladogenesis; //Probability of a clade giving rise to an entirely new higher taxon (yes, this makes the clade paraphyletic) int age; //Age of the clade int pollinationtendency; //Overall preferred pollination type of the clade, used for making new species (together with a global clade sensitivity constant). Numbers same as species. int dispersaltendency; //As above, but dispersal float rangetendency; //Tendency of the clade to have a large range or a small range. Species will gravitate toward occupying this entire range over a clade. } clade; typedef struct couplet { int index; int firstmember; //member is the index of the member of the couplet int firsttype; //type is 0 if species, 1 if another couplet int secondmember; //age is the age of the split int secondtype; //CladeName is the string name of the couplet for a NEXUS file. int age; char *CladeName; int cladenamelength; } couplet; //Define globals to be directly used species *Lineages = 0; //Array containing all of the species, to be reallocated as necessary clade *Classes = 0; //Array containing all of chte clades, also may be reallocated couplet *Pairings = 0; int Speciescount = 0; //Number of species int Speciesextinct = 0; //Number of species which have become extinct int Classcount = 0; //Number of clades int Cladesextinct = 0; //Number of clades which have become extinct int Timestep = 0; //Timestep int CurrentMaxSpecies = STARTINGSPECIES; int CurrentMaxClades = STARTINGCLADES; int ExtinctionPulse = 0; int ExtinctionPulseCount = 0; int CoupletCount = 0; //Constants to be read from file. float SpeciationDecayConstant; //Rate at which speciation rate decays with time float ExtinctionDecayConstant; //Rate at which extinction rate decays with time float EvolutionDecayConstant; //Rate at which sequence evolution rates decay with time float PollinationSpeciationMultiplier; //Arithmetic multiplier to speciation rate for animal pollination (lowers speciation rate if < 1.0) float PollinationExtinctionMultiplier; //As you would expect, but extinction float PollinationSpeciationDecayConstant; //Arithmetic multiplier to SpeciationDecayConstant for animal-pollinated things float PollinationExtinctionDecayConstant; //Extinction float PollinationEvolutionMultiplier; //Multiplies rate of molecular evolution for animal-pollinated things float PollinationRangeMultiplier; //Increases (if > 1.0) geographic range tendency for animal-pollinated (reduces if < 1.0) float DispersalSpeciationMultiplier; //These are the same things, but for dispersal float DispersalExtinctionMultiplier; float DispersalSpeciationDecayConstant; float DispersalExtinctionDecayConstant; float DispersalEvolutionMultiplier; float DispersalRangeMultiplier; float StochasticRangeMultiplier; float ExtinctionPulseFrequency; float ExtinctionPulseIntensity; float ExtinctionPulsePollinationSelectivity; float ExtinctionPulseDispersalSelectivity; float RangeImportance; float TransitionRate; //Intrinsic rate of transitions (C to T, T to C, A to G, G to A), as per a Kimura two-parameter evolutionary model float TransversionRate; //Intrinsic rate of transversions (A to C, C to G, G to T, T to A, C to A, G to C, T to G, A to T), as per a Kimura two-parameter evolutionary model float SpeciationStateChangeRates[2][2][2][2]; //Rates of character state transitions during speciation; first and second indicies represent the state changing from, the third and fourth the state changing to; float BackgroundStateChangeRates[2][2][2][2]; //Rates of character state transitions during background evolution; work the same way as speciation //first and third are pollination, second and fourth are dispersal; 0 in all cases is wind-, and 1 is animal-. float CladePollinationImportance[2]; //Multiplication constant on the rates of state changes TO wind (0) or animal (1) pollination. //Clade tendency has no effect if this constant is 1.0, but has a stronger positive effect if > 1.0. Clades will actually move AWAY from their "preferred" type if < 1.0. float CladeDispersalImportance[2]; //Exactly the same as above, but for dispersal int MaxSpecies; //If the simulation reaches this number of species, terminate. Ignored if 0. int MaxClades; //If the simulation reaches this number of clades, terminate. Ignored if 0. int MaxTimesteps; //When the simulation has run for this number of timesteps, terminate. Ignored if 0 (not recommended, as it may run on forever). float MaxSpeciesExtinction; //When this proportion of the total number of species that have become extinct reaches this value, terminate. float MaxCladeExtinction; //As above, but for clades. //File management FILE *model; //Input file that describes the model of evolution FILE *biota; //Input file that describes the starting clade(s) FILE *output; //Output file containing a table of the final parameters of each species FILE *nexus; //Output file of the surviving species in a phylogenetic tree FILE *csv; //Output file of the final character states of surviving species, for use with diversitree FILE *diversitycurve; //Output file that prints diversity over time. FILE *finalstats; //Output file that gives general statistics (used in the R code invoking both this and diversitree) char modelname[STARTINGCLADES]; char biotaname[STARTINGCLADES]; char outname[STARTINGCLADES]; char nexusname[STARTINGCLADES]; char csvname[STARTINGCLADES]; char curvename[STARTINGCLADES]; char statsname[STARTINGCLADES]; void CheckCladogenesis(int cladenumber); void CheckExtinction(int species); void CheckExtinctionPulse(); void CheckRangeDynamics(int species); void CheckSpeciation(int species); void ClearExtinctionPulse(); int DefineCouplet(int firstmember, int firstiscouplet, int secondmember, int secondiscouplet, int divergenceage); void DoBackgroundEvolution(int species); void ExitBadBiota(); void ExitBadModel(char *Bad, char *Explanation); void Extinction(int species); int FindExtantDescendentSpecies(int ancestorspecies); void GenerateDiversityCurves(int *truediversity); double GetRandomFloat(double randommax); int GetRandomInt(int randommax); void IdentifyExtinctClades(int cladenumber); void InitializeStartingClade(char *CharacterLine); float NthWordFloat(int wordnumber, char *Line); int NthWordInt(int wordnumber, char *Line); int RecursiveBuildTreeFromSpecies(int focalspecies); void SequenceEvolution(int species); void Speciate(int ancestor); int main (int argc, char **argv) { char line[MAXLINELENGTH]; int i, j, k, l; int *truediversity = 0; // int iscomma = 0; double numstates[2]; Lineages = (species*) malloc(STARTINGSPECIES * sizeof(species)); Classes = (clade*) malloc(STARTINGCLADES * sizeof(clade)); //Obsolete code for when this program was compiled on a different system /* if (argc != 8) { fprintf(stderr, "Speciation Extinction Model for Bisse data: wrong number of arguments used (needed 8, used %d). Remember, the number of arguments is 1 + the number of filenames given.\n", argc); fprintf(stderr, "First argument given should be the model, second argument the starting clade(s), and the final argument the name of the file to be output.\n"); exit(1); } */ sprintf(modelname, "%s/Written_C_Programs/src/SpeciationExtinctionModel/Data/Model.txt", getenv("HOME")); sprintf(biotaname, "%s/Written_C_Programs/src/SpeciationExtinctionModel/Data/Biota.txt", getenv("HOME")); sprintf(outname, "%s/Written_C_Programs/src/SpeciationExtinctionModel/Output/output.txt", getenv("HOME")); sprintf(nexusname, "%s/Written_C_Programs/src/SpeciationExtinctionModel/Output/TreeFile.txt", getenv("HOME")); sprintf(csvname, "%s/Written_C_Programs/src/SpeciationExtinctionModel/Output/CharStates.csv", getenv("HOME")); sprintf(curvename, "%s/Written_C_Programs/src/SpeciationExtinctionModel/Output/DiversityCurve.txt", getenv("HOME")); sprintf(statsname, "%s/Written_C_Programs/src/SpeciationExtinctionModel/Output/FinalStats.csv", getenv("HOME")); //Seed the random number generator srand48(time(0)); //Load Files model = fopen(modelname, "r"); biota = fopen(biotaname, "r"); output = fopen(outname, "w"); nexus = fopen(nexusname, "w"); csv = fopen(csvname, "w"); diversitycurve = fopen(curvename, "w"); finalstats = fopen(statsname, "w"); //Input model parameters //NthWordFloat needs to use strtod(), which has the same format as strtol(). if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "SpeciationDecayConstant"); SpeciationDecayConstant = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "ExtinctionDecayConstant"); ExtinctionDecayConstant = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "EvolutionDecayConstant"); EvolutionDecayConstant = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "PollinationSpeciationMultiplier"); PollinationSpeciationMultiplier = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "PollinationExtinctionMultiplier"); PollinationExtinctionMultiplier = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "PollinationSpeciationDecayConstant"); PollinationSpeciationDecayConstant = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "PollinationExtinctionDecayConstant"); PollinationExtinctionDecayConstant = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "PollinationEvolutionMultiplier"); PollinationEvolutionMultiplier = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "PollinationRangeMultiplier"); PollinationRangeMultiplier = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "DispersalSpeciationMultiplier"); DispersalSpeciationMultiplier = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "DispersalExtinctionMultiplier"); DispersalExtinctionMultiplier = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "DispersalSpeciationDecayConstant"); DispersalSpeciationDecayConstant = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "DispersalExtinctionDecayConstant"); DispersalExtinctionDecayConstant = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "DispersalEvolutionMultiplier"); DispersalEvolutionMultiplier = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "DispersalRangeMultiplier"); DispersalRangeMultiplier = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "StochasticRangeMultiplier"); StochasticRangeMultiplier = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "RangeImportance"); RangeImportance = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "ExtinctionPulseFrequency"); ExtinctionPulseFrequency = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "ExtinctionPulseIntensity"); ExtinctionPulseIntensity = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "ExtinctionPulsePollinationSelectivity"); ExtinctionPulsePollinationSelectivity = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "ExtinctionPulseDispersalSelectivity"); ExtinctionPulseDispersalSelectivity = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "TransitionRate"); TransitionRate = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "TransversionRate"); TransversionRate = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "SpeciationStateChangeRates"); for (i = 0; i < 2; i++) { for (j = 0; j < 2; j++) { for (k = 0; k < 2; k++) { for (l = 0; l < 2; l++) SpeciationStateChangeRates[i][j][k][l] = NthWordFloat(i + (2 * j) + (4 * k) + (8 * l), line); } } } if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "BackgroundStateChangeRates"); for (i = 0; i < 2; i++) { for (j = 0; j < 2; j++) { for (k = 0; k < 2; k++) { for (l = 0; l < 2; l++) BackgroundStateChangeRates[i][j][k][l] = NthWordFloat(i + (2 * j) + (4 * k) + (8 * l), line); } } } if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "CladePollinationImportance"); CladePollinationImportance[0] = NthWordFloat(0, line); CladePollinationImportance[1] = NthWordFloat(1, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "CladeDispersalImportance"); CladeDispersalImportance[0] = NthWordFloat(0, line); CladeDispersalImportance[1] = NthWordFloat(1, line); //Input termination conditions if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "MaxSpecies"); MaxSpecies = NthWordInt(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "MaxClades"); MaxClades = NthWordInt(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "MaxTimesteps"); MaxTimesteps = NthWordInt(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "MaxSpeciesExtinction"); MaxSpeciesExtinction = NthWordFloat(0, line); if (!fgets(line, MAXLINELENGTH, model)) ExitBadModel("", "MaxCladeExtinction"); MaxCladeExtinction = NthWordFloat(0, line); //Input starting clades for (i = 0;; i++) { if (!fgets(line, MAXLINELENGTH, biota)) break; if (strcmp(line, "\n") == 0) ExitBadBiota(); InitializeStartingClade(line); Classcount++; } if (i == 0) ExitBadBiota(); //Run the evolutionary model for (Timestep = 0; (Timestep < MaxTimesteps) || (MaxTimesteps == 0); Timestep++) { CheckExtinctionPulse(); numstates[0] = 0; numstates[1] = 0; if (Timestep == 128) i = 0; //Species evolution for (i = 0; i < Speciescount; i++) { Lineages[i].speciesage++; if (Lineages[i].extant) { DoBackgroundEvolution(i); CheckRangeDynamics(i); CheckSpeciation(i); CheckExtinction(i); } if (Lineages[i].extant) { SequenceEvolution(i); numstates[Lineages[i].pollination]++; } } //Anything happening at the class level for (i = 0; i < Classcount; i++) { Classes[i].age++; if (Classes[i].extant) IdentifyExtinctClades(i); if (Classes[i].extant) CheckCladogenesis(i); } printf("Timestep %d done.\n", Timestep); printf("Ancestral state count: %f. Derived state count: %f\n", numstates[0], numstates[1]); printf("Proportion with derived state %f\n", numstates[1] / (numstates[0] + numstates[1])); //break conditions if ((Speciescount >= MaxSpecies) && (MaxSpecies > 0)) break; else if ((Classcount >= MaxClades) && (MaxClades > 0)) break; else if (((float) Speciesextinct / (float) Speciescount) >= MaxSpeciesExtinction) break; else if (((float) Cladesextinct / (float) Classcount) >= MaxCladeExtinction) break; ClearExtinctionPulse(); //Need to initialize truediversity for the diversity curve. if (Timestep == 0) truediversity = calloc(1, sizeof(int)); else truediversity = realloc(truediversity, (1 + Timestep) * sizeof(int)); for (i = 0, truediversity[Timestep] = 0; i < Speciescount; i++) { truediversity[Timestep] += Lineages[i].extant; } } //Once broken from previous loop, we're done with the simulation. Make the general output table for (i = 0; i < Speciescount; i++) { fprintf(output, ">Species %d %d %d ", Lineages[i].index, Lineages[i].ancestor, Lineages[i].extant); fprintf(output, "%f %f %f ", Lineages[i].extinction, Lineages[i].speciation, Lineages[i].evolution); fprintf(output, "%d %d %d %d %f\n", Lineages[i].speciesage, Lineages[i].highertaxon, Lineages[i].pollination, Lineages[i].dispersal, Lineages[i].range); fprintf(output, "%s\n", Lineages[i].Sequence); } fprintf(output, "\n\n\n\n"); for (i = 0; i < Classcount; i++) { fprintf(output, "%d %s ", Classes[i].index, Classes[i].Name); fprintf(output, "%d %d %d %d ", Classes[i].ancestor, Classes[i].ancestorclade, Classes[i].diversity, Classes[i].standingdiversity); fprintf(output, "%f %f %f %f ", Classes[i].extinction, Classes[i].speciation, Classes[i].evolution, Classes[i].cladogenesis); fprintf(output, "%d %d %d:\n",Classes[i].age, Classes[i].pollinationtendency, Classes[i].dispersaltendency); fprintf(output, "Species: "); for (j = 0; j < Classes[i].diversity; j++) fprintf(output, "%d ", Classes[i].specieslist[j]); fprintf(output, "\nSurviving species: "); for (j = 0; j < Classes[i].diversity; j++) { if (Lineages[Classes[i].specieslist[j]].extant == 1) fprintf(output, "%d ", Classes[i].specieslist[j]); } fprintf(output, "\n\n"); } //Now, make the tree file! fprintf(nexus, "#NEXUS\n"); //Taxa first fprintf(nexus, "BEGIN TAXA;\n"); fprintf(nexus, "\tTITLE Taxa;\n"); fprintf(nexus, "\tDIMENSIONS NTAX=%d;\n", Speciescount - Speciesextinct); fprintf(nexus, "\tTAXLABELS\n"); fprintf(nexus, "\t\t"); for (i = 0; i < Speciescount; i++) { if (Lineages[i].extant) fprintf(nexus, "%d ", i); } fprintf(nexus, "\n\t;\n\nEND;\n\n"); //Characters next fprintf(nexus, "BEGIN CHARACTERS;\n"); fprintf(nexus, "\tTITLE Character_Matrix;\n"); fprintf(nexus, "\tDIMENSIONS NCHAR=2;\n"); fprintf(nexus, "\tFORMAT DATATYPE = STANDARD GAP = - MISSING = ? SYMBOLS = \"0 1\";\n"); fprintf(nexus, "\tMATRIX\n"); for (i = 0; i < Speciescount; i++) { if (Lineages[i].extant) fprintf(nexus, "\t%d %d%d\n", i, Lineages[i].pollination, Lineages[i].dispersal); } fprintf(nexus, "\n\n;\n\nEND;\n\n"); //Tree last fprintf(nexus, "BEGIN TREES;\n"); fprintf(nexus, "\tTitle 'Simulated taxon tree';\n"); fprintf(nexus, "\tLINK Taxa = Taxa;\n"); /* fprintf(nexus, "\tTRANSLATE"); for (i = 0; i < Speciescount; i++) { if (Lineages[i].extant) { if (iscomma) fprintf(nexus, ","); else iscomma = 1; fprintf(nexus, "\n\t\t%d Spc%d", i, i); } } fprintf(nexus, ";\n"); */ RecursiveBuildTreeFromSpecies(0); if (CoupletCount > 1) fprintf(nexus, "\tTREE simulatedtree = %s;\n", Pairings[CoupletCount - 1].CladeName); else { fprintf(stderr, "Warning: surviving diversity insufficient to generate tree (i.e. one or zero species surviving in total)!\n"); fprintf(finalstats, "finalspecies,finalsurvivors\n%d,%d\n", Speciescount, Speciescount - Speciesextinct); exit(1); } fprintf(nexus, "\nEND;\n\n"); //Now we make the csv text file because BiSSE can't seem to read characters out of a nexus file! fprintf(csv, "species,pollination,dispersal\n"); for (i = 0; i < Speciescount; i++) { if (Lineages[i].extant) fprintf(csv, "%d,%d,%d\n", i, Lineages[i].pollination, Lineages[i].dispersal); } printf("Total species: %d\nTotal surviving species: %d\nTotal number of extinction pulses: %d\n", Speciescount, Speciescount - Speciesextinct, ExtinctionPulseCount); //And print out the final stats! fprintf(finalstats, "finalspecies,finalsurvivors\n%d,%d\n", Speciescount, Speciescount - Speciesextinct); GenerateDiversityCurves(truediversity); exit(0); } void CheckCladogenesis(int cladenumber) { //No constants currently implemented. } void CheckExtinction(int species) { float actualextinction, survival, multiplier, timemeasure, randomint; actualextinction = Lineages[species].extinction; survival = 1 - actualextinction; timemeasure = Timestep; if (Lineages[species].pollination == 0) { if (Lineages[species].dispersal == 0) multiplier = pow(ExtinctionDecayConstant, timemeasure); else multiplier = pow(ExtinctionDecayConstant * DispersalExtinctionDecayConstant, timemeasure); } else { if (Lineages[species].dispersal == 0) multiplier = pow(ExtinctionDecayConstant * PollinationExtinctionDecayConstant, timemeasure); else multiplier = pow(ExtinctionDecayConstant * PollinationExtinctionDecayConstant * DispersalExtinctionDecayConstant, timemeasure); } multiplier *= pow(Lineages[species].range, RangeImportance); multiplier *= pow(ExtinctionPulseIntensity, ExtinctionPulse); if (Lineages[species].pollination) multiplier *= ExtinctionPulsePollinationSelectivity; //state 1 if (Lineages[species].dispersal) multiplier *= ExtinctionPulseDispersalSelectivity; actualextinction *= multiplier; if (Lineages[species].pollination == 1) actualextinction *= PollinationExtinctionMultiplier; if (Lineages[species].dispersal == 1) actualextinction *= DispersalExtinctionMultiplier; randomint = GetRandomFloat(actualextinction + survival); if (randomint < actualextinction) Extinction(species); } void CheckExtinctionPulse() { for (ExtinctionPulse = 0; GetRandomFloat(1) < ExtinctionPulseFrequency; ExtinctionPulse++); if (ExtinctionPulse) { //printf("Extinction pulse of intensity %f in timestep %d\n", pow(ExtinctionPulseIntensity, ExtinctionPulse), Timestep); ExtinctionPulseCount++; } } void CheckRangeDynamics(int species) { int highertaxon = Lineages[species].highertaxon; float targetrange = Classes[highertaxon].rangetendency / sqrt((float) Classes[highertaxon].standingdiversity); float stochasticvariation = pow(StochasticRangeMultiplier, 1 - GetRandomFloat(2)); if (Lineages[species].pollination == 1) targetrange *= PollinationRangeMultiplier; if (Lineages[species].dispersal == 1) targetrange *= DispersalRangeMultiplier; Lineages[species].range = sqrt(Lineages[species].range * targetrange); Lineages[species].range *= stochasticvariation; } void CheckSpeciation(int species) { float actualspeciation, nonspeciation, multiplier, timemeasure, randomint; actualspeciation = Lineages[species].speciation; nonspeciation = 1 - actualspeciation; timemeasure = Timestep; if (Lineages[species].pollination == 0) { if (Lineages[species].dispersal == 0) multiplier = pow(SpeciationDecayConstant, timemeasure); else multiplier = pow(SpeciationDecayConstant * DispersalSpeciationDecayConstant, timemeasure); } else { if (Lineages[species].dispersal == 0) multiplier = pow(SpeciationDecayConstant * PollinationSpeciationDecayConstant, timemeasure); else multiplier = pow(SpeciationDecayConstant * PollinationSpeciationDecayConstant * DispersalSpeciationDecayConstant, timemeasure); } actualspeciation *= multiplier; if (Lineages[species].pollination == 1) actualspeciation *= PollinationSpeciationMultiplier; if (Lineages[species].dispersal == 1) actualspeciation *= DispersalSpeciationMultiplier; randomint = GetRandomFloat(actualspeciation + nonspeciation); if (randomint < actualspeciation) { Speciate(species); CheckSpeciation(species); } } void Cladogenesis(int ancestorspecies, float statechangeprob, int extranewspecies) { int ancestorclade = Lineages[ancestorspecies].highertaxon; float nonectarnofruit, nonectarfruit, nectarnofruit, nectarfruit, randomstate, i; if ((CurrentMaxClades - Classcount) == 0) { CurrentMaxClades += STARTINGCLADES; Classes = (clade*) realloc(Classes, CurrentMaxClades * sizeof(clade)); } if ((CurrentMaxSpecies - Speciescount) == 0) { CurrentMaxSpecies += STARTINGSPECIES; Lineages = (species*) realloc(Lineages, CurrentMaxSpecies * sizeof(species)); } //Make the new species first. Lineages[Speciescount].index = Speciescount; Lineages[Speciescount].ancestor = ancestorspecies; Lineages[Speciescount].extant = 1; //Come back to evolution, extinction, and speciation; they're complicated. Lineages[Speciescount].speciesage = 0; Lineages[Speciescount].highertaxon = ancestorclade; //Now, character evolution nonectarnofruit = SpeciationStateChangeRates[Lineages[ancestorspecies].pollination][Lineages[ancestorspecies].dispersal][0][0]; nectarnofruit = SpeciationStateChangeRates[Lineages[ancestorspecies].pollination][Lineages[ancestorspecies].dispersal][1][0]; nonectarfruit = SpeciationStateChangeRates[Lineages[ancestorspecies].pollination][Lineages[ancestorspecies].dispersal][0][1]; nectarfruit = SpeciationStateChangeRates[Lineages[ancestorspecies].pollination][Lineages[ancestorspecies].dispersal][1][1]; if (Classes[ancestorclade].pollinationtendency == 0) { nonectarnofruit *= CladePollinationImportance[0]; nonectarfruit *= CladePollinationImportance[0]; } else { nectarnofruit *= CladePollinationImportance[1]; nectarfruit *= CladePollinationImportance[1]; } if (Classes[ancestorclade].dispersaltendency == 0) { nonectarnofruit *= CladeDispersalImportance[0]; nectarnofruit *= CladeDispersalImportance[0]; } else { nonectarfruit *= CladeDispersalImportance[1]; nectarfruit *= CladeDispersalImportance[1]; } randomstate = GetRandomFloat(nonectarnofruit + nonectarfruit + nectarnofruit + nectarfruit); if (randomstate < nonectarnofruit) { Lineages[Speciescount].pollination = 0; Lineages[Speciescount].dispersal = 0; } else if (randomstate < (nonectarnofruit + nonectarfruit)) { Lineages[Speciescount].pollination = 0; Lineages[Speciescount].dispersal = 1; } else if (randomstate < (nonectarnofruit + nonectarfruit + nectarnofruit)) { Lineages[Speciescount].pollination = 1; Lineages[Speciescount].dispersal = 0; } else { Lineages[Speciescount].pollination = 1; Lineages[Speciescount].dispersal = 1; } strcpy(Lineages[Speciescount].Sequence, Lineages[ancestorspecies].Sequence); //At separation, two lineages share identical sequence. Lineages[Speciescount].range = GetRandomFloat(Lineages[ancestorspecies].range); //New ranges start out small. Lineages[ancestorspecies].range -= GetRandomFloat(Lineages[Speciescount].range); //They also take away from existing ranges. Lineages[Speciescount].evolution = Lineages[ancestorspecies].evolution; Lineages[Speciescount].speciation = Lineages[ancestorspecies].speciation; Lineages[Speciescount].extinction = Lineages[ancestorspecies].extinction; //Cladogenesis-related changes in state if (GetRandomFloat(1 + statechangeprob) < statechangeprob) Lineages[Speciescount].pollination = GetRandomInt(2); if (GetRandomFloat(1 + statechangeprob) < statechangeprob) Lineages[Speciescount].dispersal = GetRandomInt(2); //Incomplete code from an earlier draft of the model, used to set up new radiations from within otherwise declining groups /* if (statechangeprob) { Lineages[Speciescount].evolution *= pow(1 / EvolutionDecayConstant, Classes[ancestorclade].age / 2); Lineages[Speciescount].extinction *= pow(1 / ExtinctionDecayConstant, Classes[ancestorclade].age / 2); Lineages[Speciescount].speciation *= pow(1 / SpeciationDecayConstant, Classes[ancestorclade].age / 2); } */ //Now, create the clade. Classes[Classcount].index = Classcount; sprintf(Classes[Classcount].Name, "%s.%d", Classes[ancestorclade].Name, Classcount); Classes[Classcount].ancestor = ancestorspecies; Classes[Classcount].ancestorclade = ancestorclade; Classes[Classcount].specieslist[0] = Speciescount; Classes[Classcount].diversity = 1; Classes[Classcount].extant = 1; Classes[Classcount].extinction = Lineages[Speciescount].extinction; Classes[Classcount].speciation = Lineages[Speciescount].speciation; Classes[Classcount].evolution = Lineages[Speciescount].evolution; Classes[Classcount].cladogenesis = Classes[ancestorclade].cladogenesis; Classes[Classcount].age = 0; Classes[Classcount].pollinationtendency = Lineages[Speciescount].pollination; Classes[Classcount].dispersaltendency = Lineages[Speciescount].dispersal; if (statechangeprob) { if (GetRandomInt(2)) Classes[Classcount].rangetendency = Classes[ancestorclade].rangetendency * (1 - statechangeprob); else Classes[Classcount].rangetendency = Classes[ancestorclade].rangetendency / (1 - statechangeprob); } else Classes[Classcount].rangetendency = Classes[ancestorclade].rangetendency; printf("New clade %s originated, founded by species %d\n", Classes[Classcount].Name, Speciescount); Speciescount++; for (i = 0; i < extranewspecies; i++) Speciate(Classes[Classcount].specieslist[GetRandomInt(Classes[Classcount].diversity)]); Classcount++; } void ClearExtinctionPulse() { ExtinctionPulse = 0; } int DefineCouplet(int firstmember, int firstiscouplet, int secondmember, int secondiscouplet, int divergenceage) { int firstbranch, secondbranch, neededlength, i; char firststring[40], secondstring[40]; //40 digits is way more than needed! char *firstpointer, *secondpointer; if (CoupletCount == 0) Pairings = calloc(Speciescount, sizeof(couplet)); for (i = 0; i < CoupletCount; i++) { if ((Pairings[i].firstmember == firstmember) && (Pairings[i].firsttype == firstiscouplet)) return i; if ((Pairings[i].secondmember == firstmember) && (Pairings[i].secondtype == firstiscouplet)) return i; if ((Pairings[i].firstmember == secondmember) && (Pairings[i].firsttype == secondiscouplet)) return i; if ((Pairings[i].secondmember == secondmember) && (Pairings[i].secondtype == secondiscouplet)) return i; } Pairings[CoupletCount].index = CoupletCount; Pairings[CoupletCount].firstmember = firstmember; Pairings[CoupletCount].firsttype = firstiscouplet; Pairings[CoupletCount].secondmember = secondmember; Pairings[CoupletCount].secondtype = secondiscouplet; Pairings[CoupletCount].age = divergenceage; if (firstiscouplet) { firstbranch = divergenceage - Pairings[firstmember].age; firstpointer = Pairings[firstmember].CladeName; neededlength = strlen(Pairings[firstmember].CladeName); } else { firstbranch = divergenceage; sprintf(firststring, "%d", firstmember); firstpointer = firststring; neededlength = strlen(firststring); Lineages[firstmember].couplet = CoupletCount; } if (secondiscouplet) { secondbranch = divergenceage - Pairings[secondmember].age; secondpointer = Pairings[secondmember].CladeName; neededlength += strlen(Pairings[secondmember].CladeName); } else { secondbranch = divergenceage; sprintf(secondstring, "%d", secondmember); secondpointer = secondstring; neededlength += strlen(secondstring); Lineages[secondmember].couplet = CoupletCount; } neededlength += 2 + log10(MAX(firstbranch,1)) + log10(MAX(secondbranch,1)); Pairings[CoupletCount].CladeName = calloc(10 + neededlength, sizeof(char)); //need only 6 characters (two colons, a comma, open and close parens, and a \0), but give a cushion sprintf(Pairings[CoupletCount].CladeName, "(%s:%d,%s:%d)", firstpointer, firstbranch, secondpointer, secondbranch); Pairings[CoupletCount].cladenamelength = strlen(Pairings[CoupletCount].CladeName); CoupletCount++; return Pairings[CoupletCount - 1].index; } void DoBackgroundEvolution(int species) { int highertaxon = Lineages[species].highertaxon; float nectarnofruit, nonectarnofruit, nonectarfruit, nectarfruit, randomstate; nonectarnofruit = BackgroundStateChangeRates[Lineages[species].pollination][Lineages[species].dispersal][0][0]; nectarnofruit = BackgroundStateChangeRates[Lineages[species].pollination][Lineages[species].dispersal][1][0]; nonectarfruit = BackgroundStateChangeRates[Lineages[species].pollination][Lineages[species].dispersal][0][1]; nectarfruit = BackgroundStateChangeRates[Lineages[species].pollination][Lineages[species].dispersal][1][1]; if (Classes[highertaxon].pollinationtendency == 0) { nonectarnofruit *= CladePollinationImportance[0]; nonectarfruit *= CladePollinationImportance[0]; } else { nectarnofruit *= CladePollinationImportance[1]; nectarfruit *= CladePollinationImportance[1]; } if (Classes[highertaxon].dispersaltendency == 0) { nonectarnofruit *= CladeDispersalImportance[0]; nectarnofruit *= CladeDispersalImportance[0]; } else { nonectarfruit *= CladeDispersalImportance[1]; nectarfruit *= CladeDispersalImportance[1]; } randomstate = GetRandomFloat(nonectarnofruit + nonectarfruit + nectarnofruit + nectarfruit); if (randomstate < nonectarnofruit) { Lineages[species].pollination = 0; Lineages[species].dispersal = 0; } else if (randomstate < (nonectarnofruit + nonectarfruit)) { Lineages[species].pollination = 0; Lineages[species].dispersal = 1; } else if (randomstate < (nonectarnofruit + nonectarfruit + nectarnofruit)) { Lineages[species].pollination = 1; Lineages[species].dispersal = 0; } else { Lineages[species].pollination = 1; Lineages[species].dispersal = 1; } } void ExitBadBiota() { fprintf(stderr, "Biota file not initialized!\n"); fprintf(stderr, "Copy the results file into the biota for a (simple) starting clade (note all parameters, including speciation and extinction, will be set to neutral values!).\n"); fprintf(output, "Amborellales 0.5 0.5 1.0 1.0 0 0 1.0 20"); exit(1); } void ExitBadModel(char *Bad, char *Explanation) { if (Bad[0] == '\0') fprintf(stderr, "Error in input of model parameter %s: either model file truncated or extra line input.\n", Explanation); else fprintf(stderr, "There is an error in the model input (1st input file):\n%s, bad value %s.\n", Explanation, Bad); fprintf(stderr, "A sample model (all parameters set to 1.0, meaning no effects) file that you could use is being printed to the output.\n"); fprintf(output, "1.0 --SpeciationDecayConstant; //Rate at which speciation rate decays with time (does not decay if 1.0, grows if > 1)\n"); fprintf(output, "1.0 --ExtinctionDecayConstant; //Rate at which extinction rate decays with time (ditto)\n"); fprintf(output, "1.0 --EvolutionDecayConstant; //Rate at which rates of sequence evolution decay with time\n"); fprintf(output, "1.0 --PollinationSpeciationMultiplier; //Arithmetic multiplier to speciation rate for animal pollination (lowers speciation rate if < 1.0)\n"); fprintf(output, "1.0 --PollinationExtinctionMultiplier; //As you would expect, but extinction\n"); fprintf(output, "1.0 --PollinationSpeciationDecayConstant; //Arithmetic multiplier to SpeciationDecayConstant for animal-pollinated things\n"); fprintf(output, "1.0 --PollinationExtinctionDecayConstant; //Extinction\n"); fprintf(output, "1.0 --PollinationEvolutionMultiplier; //Multiplies rate of molecular evolution for animal-pollinated things\n"); fprintf(output, "1.0 --PollinationRangeMultiplier; //Increases (if > 1.0) geographic range tendency for animal-pollinated (reduces if < 1.0)\n"); fprintf(output, "1.0 --DispersalSpeciationMultiplier; //These are the same things, but for dispersal\n"); fprintf(output, "1.0 --DispersalExtinctionMultiplier;\n"); fprintf(output, "1.0 --DispersalSpeciationDecayConstant;\n"); fprintf(output, "1.0 --DispersalExtinctionDecayConstant\n"); fprintf(output, "1.0 --DispersalEvolutionMultiplier\n"); fprintf(output, "1.0 --DispersalRangeMultiplier\n"); fprintf(output, "1.0 --StochasticRangeMultiplier //Inter-timestep random volatility in species ranges; a value of 1.0 means no variation.\n"); fprintf(output, "0.0 --RangeImportance //Exponent on the effect of geographic range on extinction. 0.0 = no effect.\n"); fprintf(output, "0.0 --ExtinctionPulseFreuency //Frequency of elevated extinction timesteps. 0.0 = no effect.\n"); fprintf(output, "1.0 --ExtinctionPulseIntensity //Intensity of extinction pulses. 1.0 = extinctions have no effect.\n"); fprintf(output, "1.0 --ExtinctionPulsePollinationSelectivity //Selectivity of extinction pulses. 1.0 = no effect; >1 means animal pollination suffers more severe extinction.\n"); fprintf(output, "1.0 --ExtinctionPulseDispersalSelectivity //Selectivity of extinction pulses. 1.0 = no effect; >1 means animal dispersal suffers more severe extinction.\n"); fprintf(output, "1.0 --TransitionRate; //Intrinsic rate of transitions (C to T, T to C, A to G, G to A), as per a Kimura two-parameter evolutionary model\n"); fprintf(output, "1.0 --TransversionRate; //Intrinsic rate of transversions (A to C, C to G, G to T, T to A, C to A, G to C, T to G, A to T), as per a Kimura two-parameter evolutionary model\n"); fprintf(output, "1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0"); fprintf(output, " --Rates of character state transitions during speciation ; first and second indicies represent the state changing from, the third and fourth the state changing to;"); fprintf(output, " first and third are pollination, second and fourth are dispersal; 0 in all cases is wind-, and 1 is animal-. This configuration generates no character evolution.\n"); fprintf(output, "1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0"); fprintf(output, " --Rates of character state transitions during background evolution; mechanics are the same as for speciation, but happen when a species already exists."); fprintf(output, "1.0 1.0 --CladePollinationImportance[2]; //Multiplication constant on the rates of state changes TO wind (0) or animal (1) pollination.\n"); fprintf(output, "1.0 1.0 --CladeDispersalImportance[2]; //Exactly the same as above, but for dispersal; "); fprintf(output, "clade tendency has no effect if this constant is 1.0, but has a stronger positive effect if > 1.0. Clades will actually move AWAY from their 'preferred' type if < 1.0.\n"); fprintf(output, "0 --MaxSpecies; //If the simulation reaches this number of species, terminate. Ignored if 0.\n"); fprintf(output, "0 --MaxClades; //If the simulation reaches this number of clades, terminate. Ignored if 0.\n"); fprintf(output, "1024 --MaxTimesteps; //When the simulation has run for this number of timesteps, terminate. Ignored if 0 (not recommended, as it may run on forever).\n"); fprintf(output, "1.0 --MaxSpeciesExtinction; //When this proportion of the total number of species that have become extinct reaches this value, terminate.\n"); fprintf(output, "1.0 --MaxCladeExtinction; //As above, but for clades.\n"); fprintf(output, "\n\nOnly the numbers at the beginning of each line before the hyphens are actually necessary; everything else is comments. Hyphens do not create comments, however;"); fprintf(output, " including hyphenated comments will cause the program to abort and print this file (the position of the numbers in each line is what matters). Similarly, the program does not"); fprintf(output, " read after the MaxCladeExtinction variable, so additional comments can be put here without aborting the program.\n"); fprintf(output, "\nBe advised, however: rates of speciation, extinction, and evolution input into this model are not the actual rates (the true rates are modified by range, pollination, dispersal, and time, as determined by other parameters above).\n"); fflush(output); exit(1); } void Extinction(int species) { Lineages[species].extant = 0; //printf("Species %d of clade %d (%s) has become extinct.\n", species, Lineages[species].highertaxon, Classes[Lineages[species].highertaxon].Name); Speciesextinct++; } int FindExtantDescendentSpecies(int ancestorspecies) { int i, targetspecies; if (Lineages[ancestorspecies].extant) return ancestorspecies; else { for (i = 0; i < Lineages[ancestorspecies].numdescendents; i++) { targetspecies = FindExtantDescendentSpecies(Lineages[ancestorspecies].Descendents[i]); if (targetspecies >= 0) return targetspecies; } return -1; } } void GenerateDiversityCurves(int *truediversity) { int i, j, k, origindate; int *treediversity; treediversity = calloc(Timestep, sizeof(int)); fprintf(diversitycurve, "Timestep\tTrue diversity\tTree diversity"); for (i = 0; i < Timestep; i++) //Cycle through each timestep to generate the diversity curve. { //Housekeeping: no species yet ancestral to anything in this timestep. for (j = 0; j < Speciescount; j++) Lineages[j].currentlyancestral = 0; for (j = 0; j < Speciescount; j++) //Cycle through each species to find who is extant. { if (Lineages[j].extant) //Species that are extant, loop back and find the ancestor. { for (k = j; k >= 0; k = Lineages[k].ancestor) { origindate = Timestep - Lineages[k].speciesage; //Origindate is the timestep in which the ancestor species speciated. //If i >= origindate, then the species was in existence by timestep i, so this was the species in our lineages through time plot. if (i >= origindate) { Lineages[k].currentlyancestral = 1; //Set it as an ancestor (it may already be from another living descendent). break; } } } } //Now, we generate the original diversity curve data for timestep i (need truediversity from outside); for (j = 0; j < Speciescount; j++) treediversity[i] += Lineages[j].currentlyancestral; //Print out the table fprintf(diversitycurve, "\n%d\t%d\t%d", i, truediversity[i], treediversity[i]); } fflush(diversitycurve); } double GetRandomFloat(double randommax) { double random = drand48() * randommax; return random; } int GetRandomInt(int randommax) { int random = GetRandomFloat(randommax); return random; } void IdentifyExtinctClades(int cladenumber) { int i, cladeextant = 0; if (Classes[cladenumber].extant == 0) //Already extinct return; for (i = 0; i < Classes[cladenumber].diversity; i++) cladeextant += Lineages[Classes[cladenumber].specieslist[i]].extant; if (cladeextant == 0) { //printf("Clade %d (%s) has become extinct in timestep %d!\n", cladenumber, Classes[cladenumber].Name, Timestep); Classes[cladenumber].extant = 0; } else { printf("In timestep %d, ", Timestep); printf("clade %d ", cladenumber); printf("(%s) ", Classes[cladenumber].Name); printf("comprises %d species.\n", cladeextant); } } void InitializeStartingClade(char *CharacterLine) { int i, randomint, diversity; if ((CurrentMaxClades - Classcount) == 0) { CurrentMaxClades += STARTINGCLADES; Classes = (clade*) realloc(Classes, CurrentMaxClades * sizeof(clade)); } Classes[Classcount].index = Classcount; Classes[Classcount].ancestor = -1; //Signifies that the clade is initial Classes[Classcount].ancestorclade = -1; //Also signifies that clade is initial Classes[Classcount].extant = 1; for (i = 0; CharacterLine[i] != ' '; i++) Classes[Classcount].Name[i] = CharacterLine[i]; Classes[Classcount].Name[i] = '\0'; Classes[Classcount].extinction = NthWordFloat(1, CharacterLine); Classes[Classcount].speciation = NthWordFloat(2, CharacterLine); //Note: speciation and extinction for initialized species does NOT account for dispersal or pollination! Classes[Classcount].evolution = NthWordFloat(3, CharacterLine); Classes[Classcount].cladogenesis = NthWordFloat(4, CharacterLine); Classes[Classcount].age = 0; Classes[Classcount].pollinationtendency = NthWordInt(5, CharacterLine); Classes[Classcount].dispersaltendency = NthWordInt(6, CharacterLine); Classes[Classcount].rangetendency = NthWordFloat(7, CharacterLine); //Now, add species diversity = NthWordInt(8, CharacterLine); if (diversity == 0) { fprintf(stderr, "Error: clade %d, %s, has been coded as extinct upon initialization! Not much evolution in an already extinct clade...\n", Classes[Classcount].index, Classes[Classcount].Name); Classes[Classcount].extant = 0; Cladesextinct++; return; } else { if ((CurrentMaxSpecies - Speciescount) == 0) { CurrentMaxSpecies += STARTINGSPECIES; Lineages = (species*) realloc(Lineages, CurrentMaxSpecies * sizeof(species)); } Classes[Classcount].specieslist[0] = Speciescount; Classes[Classcount].diversity++; Classes[Classcount].standingdiversity++; Lineages[Speciescount].index = Speciescount; Lineages[Speciescount].ancestor = -1; Lineages[Speciescount].extant = 1; Lineages[Speciescount].extinction = Classes[Classcount].extinction; Lineages[Speciescount].speciation = Classes[Classcount].speciation; Lineages[Speciescount].evolution = Classes[Classcount].evolution; Lineages[Speciescount].speciesage = 0; Lineages[Speciescount].highertaxon = Classcount; Lineages[Speciescount].pollination = Classes[Classcount].pollinationtendency; Lineages[Speciescount].dispersal = Classes[Classcount].dispersaltendency; Lineages[Speciescount].range = GetRandomFloat(Classes[Classcount].rangetendency); Lineages[Speciescount].couplet = -1; for (i = 0; i < SEQUENCELENGTH; i++) { randomint = GetRandomInt(4); if (randomint == 0) Lineages[Speciescount].Sequence[i] = 'A'; else if (randomint == 1) Lineages[Speciescount].Sequence[i] = 'C'; else if (randomint == 2) Lineages[Speciescount].Sequence[i] = 'G'; else Lineages[Speciescount].Sequence[i] = 'T'; } Speciescount++; for (i = 1; i < diversity; i++) Speciate(Classes[Classcount].specieslist[GetRandomInt(Classes[Classcount].diversity)]); } } float NthWordFloat(int wordnumber, char *Line) { int i, wordlength, countdown; char NewWord[MAXLINELENGTH]; char *bad = 0; float targetvalue; for (i = 0, wordlength = 0, countdown = wordnumber; countdown >= 0; i++) { if (Line[i] == ' ') countdown--; else if ((Line[i] == '\n') || (Line[i] == '\0')) { if (countdown == 0) break; else { sprintf(NewWord, "End of line reached finding Nth word in line, word %d", wordnumber); ExitBadModel(Line, NewWord); } } else if (countdown == 0) { NewWord[wordlength] = Line[i]; wordlength++; } } NewWord[wordlength] = '\0'; targetvalue = strtod(NewWord, &bad); if (*bad) ExitBadModel(NewWord, "Bad floating point number"); else return targetvalue; } int NthWordInt(int wordnumber, char *Line) { int i, wordlength, countdown, targetvalue; char NewWord[MAXLINELENGTH]; char *bad = 0; for (i = 0, wordlength = 0, countdown = wordnumber; countdown >= 0; i++) { if (Line[i] == ' ') countdown--; else if ((Line[i] == '\n') || (Line[i] == '\0')) { if (countdown == 0) break; else { sprintf(NewWord, "End of line reached finding Nth word in line, word %d", wordnumber); ExitBadModel(Line, NewWord); } } else if (countdown == 0) { NewWord[wordlength] = Line[i]; wordlength++; } } NewWord[wordlength] = '\0'; targetvalue = strtol(NewWord, &bad, 10); if (*bad) ExitBadModel(NewWord, "Bad integer"); else return targetvalue; } int RecursiveBuildTreeFromSpecies(int focalspecies) { int i, firstspecies, secondspecies, firstcouplet, secondcouplet, targetspecies; firstcouplet = -1; secondcouplet = -1; firstspecies = -1; secondspecies = -1; for (i = Lineages[focalspecies].numdescendents - 1; i >= 0; i--) { targetspecies = Lineages[focalspecies].Descendents[i]; if (firstcouplet < 0) firstcouplet = RecursiveBuildTreeFromSpecies(targetspecies); else secondcouplet = RecursiveBuildTreeFromSpecies(targetspecies); //There will be no uncoupled extant species from the descentent if there are more than one living descendents. //This will only return -1 if there are either no living descendents or if there is only one. //FindExtantDescendentSpecies will find the one survivor (it will find *a* survivor if there are more than one, so we don't want to use it unless we're sure there are only one or zero). if (firstcouplet >= 0) { if (Lineages[focalspecies].extant && (Lineages[focalspecies].couplet < 0)) firstcouplet = DefineCouplet(firstcouplet, 1, focalspecies, 0, Lineages[targetspecies].speciesage); else if (secondcouplet >= 0) { firstcouplet = DefineCouplet(firstcouplet, 1, secondcouplet, 1, Lineages[targetspecies].speciesage); secondcouplet = -1; } else if (firstspecies >= 0) { if (Lineages[firstspecies].couplet < 0) { firstcouplet = DefineCouplet(firstcouplet, 1, firstspecies, 0, Lineages[targetspecies].speciesage); firstspecies = -1; } } else { firstspecies = FindExtantDescendentSpecies(targetspecies); if (firstspecies >= 0) { if (Lineages[firstspecies].couplet < 0) { firstcouplet = DefineCouplet(firstcouplet, 1, firstspecies, 0, Lineages[targetspecies].speciesage); firstspecies = -1; } else firstspecies = -1; } } } else if (firstspecies < 0) firstspecies = FindExtantDescendentSpecies(targetspecies); else secondspecies = FindExtantDescendentSpecies(targetspecies); if (firstcouplet < 0) { if (Lineages[focalspecies].extant && (Lineages[focalspecies].couplet < 0) && (firstspecies >= 0)) { if (Lineages[firstspecies].couplet < 0) { firstcouplet = DefineCouplet(firstspecies, 0, focalspecies, 0, Lineages[firstspecies].speciesage); firstspecies = -1; } else firstspecies = -1; } if ((secondspecies >= 0) && (firstspecies >= 0)) { if ((Lineages[secondspecies].couplet < 0) && (Lineages[firstspecies].couplet < 0)) { firstcouplet = DefineCouplet(firstspecies, 0, secondspecies, 0, Lineages[targetspecies].speciesage); firstspecies = -1; secondspecies = -1; } else { if (Lineages[secondspecies].couplet >= 0) secondspecies = -1; else firstspecies = -1; } } } } return firstcouplet; } void SequenceEvolution(int species) { int i; for (i = 0; i < SEQUENCELENGTH; i++) { if (GetRandomFloat(1) < TransversionRate) { if ((Lineages[species].Sequence[i] == 'A') || (Lineages[species].Sequence[i] == 'G')) { if (GetRandomInt(2)) Lineages[species].Sequence[i] = 'C'; else Lineages[species].Sequence[i] = 'T'; } else { if (GetRandomInt(2)) Lineages[species].Sequence[i] = 'A'; else Lineages[species].Sequence[i] = 'G'; } } else if (GetRandomFloat(1) < TransitionRate) { if (Lineages[species].Sequence[i] == 'A') Lineages[species].Sequence[i] = 'G'; else if (Lineages[species].Sequence[i] == 'C') Lineages[species].Sequence[i] = 'T'; else if (Lineages[species].Sequence[i] == 'G') Lineages[species].Sequence[i] = 'A'; else Lineages[species].Sequence[i] = 'C'; } } } void Speciate(int ancestor) { int highertaxon = Lineages[ancestor].highertaxon; float nonectarnofruit, nonectarfruit, nectarnofruit, nectarfruit, randomstate; //Determine highertaxon if (Classes[highertaxon].diversity == SPECIESPERCLADE) { Cladogenesis(ancestor, 0, 0); //0 and 0 mean no more than usual chance of character state changes return; } if ((CurrentMaxSpecies - Speciescount) == 0) { CurrentMaxSpecies += STARTINGSPECIES; Lineages = (species*) realloc(Lineages, CurrentMaxSpecies * sizeof(species)); } //Note this function only speciates new species, not new clades. Use Cladogenesis for that. Lineages[Speciescount].index = Speciescount; Lineages[Speciescount].ancestor = ancestor; Lineages[Speciescount].extant = 1; //Come back to evolution, extinction, and speciation; they're complicated. Lineages[Speciescount].speciesage = 0; Lineages[Speciescount].highertaxon = highertaxon; Classes[highertaxon].specieslist[Classes[highertaxon].diversity] = Speciescount; Classes[highertaxon].diversity++; Classes[Classcount].standingdiversity++; if (Timestep > 0) { //printf("Species %d originated of clade %s, timestep %d.\n", Speciescount, Classes[highertaxon].Name, Timestep); Lineages[Speciescount].range = GetRandomFloat(Lineages[ancestor].range); //New species have small ranges Lineages[ancestor].range -= GetRandomFloat(Lineages[Speciescount].range); //But they take a chunk out of the original species range. } else Lineages[Speciescount].range = GetRandomFloat(Classes[highertaxon].rangetendency); //Now, character evolution nonectarnofruit = SpeciationStateChangeRates[Lineages[ancestor].pollination][Lineages[ancestor].dispersal][0][0]; nectarnofruit = SpeciationStateChangeRates[Lineages[ancestor].pollination][Lineages[ancestor].dispersal][1][0]; nonectarfruit = SpeciationStateChangeRates[Lineages[ancestor].pollination][Lineages[ancestor].dispersal][0][1]; nectarfruit = SpeciationStateChangeRates[Lineages[ancestor].pollination][Lineages[ancestor].dispersal][1][1]; if (Classes[highertaxon].pollinationtendency == 0) { nonectarnofruit *= CladePollinationImportance[0]; nonectarfruit *= CladePollinationImportance[0]; } else { nectarnofruit *= CladePollinationImportance[1]; nectarfruit *= CladePollinationImportance[1]; } if (Classes[highertaxon].dispersaltendency == 0) { nonectarnofruit *= CladeDispersalImportance[0]; nectarnofruit *= CladeDispersalImportance[0]; } else { nonectarfruit *= CladeDispersalImportance[1]; nectarfruit *= CladeDispersalImportance[1]; } randomstate = GetRandomFloat(nonectarnofruit + nonectarfruit + nectarnofruit + nectarfruit); if (randomstate < nonectarnofruit) { Lineages[Speciescount].pollination = 0; Lineages[Speciescount].dispersal = 0; } else if (randomstate < (nonectarnofruit + nonectarfruit)) { Lineages[Speciescount].pollination = 0; Lineages[Speciescount].dispersal = 1; } else if (randomstate < (nonectarnofruit + nonectarfruit + nectarnofruit)) { Lineages[Speciescount].pollination = 1; Lineages[Speciescount].dispersal = 0; } else { Lineages[Speciescount].pollination = 1; Lineages[Speciescount].dispersal = 1; } strcpy(Lineages[Speciescount].Sequence, Lineages[ancestor].Sequence); //At separation, two lineages share identical sequence. Lineages[Speciescount].evolution = Lineages[ancestor].evolution; Lineages[Speciescount].speciation = Lineages[ancestor].speciation; Lineages[Speciescount].extinction = Lineages[ancestor].extinction; //Manage clades if (Lineages[ancestor].numdescendents == 0) Lineages[ancestor].Descendents = calloc(256, sizeof(int)); else if ((Lineages[ancestor].numdescendents % 256) == 0) Lineages[ancestor].Descendents = realloc(Lineages[ancestor].Descendents, (Lineages[ancestor].numdescendents + 256) * sizeof(int)); Lineages[ancestor].Descendents[Lineages[ancestor].numdescendents] = Speciescount; Lineages[ancestor].numdescendents++; Lineages[Speciescount].numdescendents = 0; Lineages[Speciescount].couplet = -1; Speciescount++; }