;;; Carlo C. Maley December 6, 1995 ;;; This program is supposed to replicate the Sepkoski & Kendrick (1993) model ;;; with the addition of a cladistic form of identifying taxa. Our intention ;;; is to remove the confounding effect of loss of data when comparing ;;; monophyletic and paraphyletic taxa. ;;; This is a sort of "null" model for the diversification of life. ;;; Perhaps there has been no change in conditions of speciation or ;;; extinction over time, but simply the probability of speciation has ;;; been higher than the probability of extinction. Thus resulting in ;;; an exponential (geometric?) curve for diversity over time. ;;; This model is an attempt to explore the characteristics of such a ;;; system. ;;; Sepkoski & Kendrick's concept of speciation has been layered on top of a more ;;; systematist concept by adding the "leaf-id" fields. When a speciation event ;;; occurs, the left leaf gets the leaf-id of the parent node, and the right branch ;;; gets a new leaf-id. ;;; Things to do:*** ;;; Done: 2. Extinction as a rise in extinction rate. This will be handled by ;;; a list that encodes pairs of the degree of extinction and when the number ;;; of nodes proportional to the maximum (*next-id*) that triggers it. Dp ;;; the list ((30 34) (60 83) (90 91)) would encode mass extinctions where ;;; about 34%, 83% and 91% of species went extinct after 30%, 60%, and 90% of ;;; the run was completed. ;;; Done: 3. Identifying paraphyletic taxa (random ancester plus all ;;; non-ancestor decendants) for one tenth of all nodes. ;;; Only accept one of these taxa with probability = 1 - 1/(n + k) ;;; where n is the smaller of the number of species in the taxon or ;;; the ancestral taxon (which is being reduced by having this new ;;; taxon removed from it). And K = 0.05, 0.5 and 100. ;;; Done: 4. Identify monophyletic groups (soft and hard cull). ;;; Done: 5. Implement sampling (constant rate every time step of producing ;;; fossil, and also ??? exponential distribution of probabilities which ;;; is used to determine a given time step's sampling rate) ;;; No. ??? 6. Do analysis for accumulations of 6 time steps. ;;; No.??? 7. Program partial coefficients (oy!) ;;; 8. Gather diversity statistics ;;; At first blush the div. stats. seem to be working for the paraphyletic ;;; groups and lineages. Next step is to gather all the different taxa and ;;; then all the forms of taxa under the sampling stats. (in-package 'USER) (setf *print-array* t) (format t "1 ") (defstruct snode ;species node (id 0) (leaf-id 0) (phenotype 0) ; Bit patterns representing morphology and characters. (prey 0) (generalism 0) (left ()) ; Left child snode. (right ()) ; Right child snode. (parent ()) ; Parental snode. (next ()) ; Next snode in a linked list. (age 0) ; The time when the snode went extinct or speciated. ;(depth 0) ; Depth of the sub-tree below this node. ;(tag ()) ; A tag for analyzing tree structure. ;(mutations ()) ; A list of which bits were flipped since speciation. ;(dominant ()) ; A list indicating membership in various depth ; dominant clades. (num-leaves 0) ; the number of leaves below this node. (first-sample -1) ; time when this lineage was first sampled in the strata. (last-sample -1) ; time when this lineage was last sampled in the strata. (paraphyle 0) ; default to the root taxon. (syst-monophyle -1) ; default to the root taxon. (hard-monophyle -1) (soft-monophyle -1) (full-monophyle -1) ; same as syst-monophyle + monotypic taxa. (rand-taxon -1) ; randomly assigned taxa (match-rand -1) (distrib-rand -1) (distrib-soft -1) (nested-monophyle -1) ) (format t "2 ") ;----------------------Global Variables------------------------------ (defvar *save-info* ()) ; Flag for saving log and tree info. (defvar *save-trees* ()) ; Flag for saving the classification trees. (defvar *root* ()) ; The root of the tree. (defvar *num-species* 0) ; number of species alive. (defvar *speciation* 0) ; Probability of speciation (out of PROBSCALE) (format t "2.1 ") (defvar *extinction* 0) ; Probability of extinction (out of PROBSCALE) (format t "2.2 ") (defvar *sample* 0) ; Probability of sampling a lineage in the strata (out of PROBSCALE) (defvar *k* 100) ; parameter for paraphyly. (defvar *paraphyly-prob* 0) ; Probability that a given node defines a taxon. (defvar *mutation* 0) ; Probability of macro mutation (out of PROBSCALE) (defvar *max-diversity* 0) ; Trigger for halting diversification. (format t "2.3 ") (defvar *first-epoch* ()) ; The extant list from the first epoch. (defvar *next-id* 0) ; Number of lineages in the model. (defvar *next-leaf-id* 0) ; Number of leaves in the model. (defvar *extinct-sched* ()) ; pairs of sizes and times of mass extinctions. (defvar *clock* 0) ; The keeper of Time. (format t "3 ") ;==============================Main================================== ;;; The amount-extinction is a proportion (0 to <1). The model will ;;; kill off between amount-extinction and 100% species depending on ;;; the bits that are selected and the population. Selection is ;;; complete -> no species with the susceptible bits will survive, but ;;; the exact percentage is not guaranteed. ;;; The Sepkoski and Kendric model flows in the following stages ;;; 1. Generate a tree (perhaps with extinctions) ;;; 2. Classify it into taxa ;;; 3. Add sampling data (first and last sample) ;;; 4. Compile diversity data ;;; Classifiying the tree is done in multiple ways ;;; 1. Label it with paraphyletic groups ;;; 2. Label it with monophyletic groups via a hard cull ;;; 3. Label it with monophyletic groups via a soft cull ;;; 4. Label it with a better monophyletic classification (?) ;;; 5. Label it with nested monophyletic groups ; (diversify 1 .4 t 1000 .2 .05 .1 .05 '((50 60)) t) ; (diversify 1 1.0 t 20 .5 .05 .4 .05 '((50 60)) t) ; (diversify 1 .8 t 20 .4 .05 .4 .05 '((50 60)) t) ; (diversify 1 .4 t 20 .2 .05 .2 .05 '((50 60)) t) ; (diversify 1 .4 t 1000 .2 .2 .1 .05 '((20 70) (45 95) (80 80)) t) ; (diversify 48 .4 t 1000 .2 .05 .1 .05 '((20 70) (45 95) (80 80)) t) ; (diversify 30 .4 t 1000 .2 .05 .1 .05 '((20 35) (45 48) (80 40)) t) run 1006 ;(progn ; (diversify 80 .4 t 1000 .2 .1 .1 .05 '((20 70) (45 95) (80 80)) t) ; (diversify 80 .4 t 1000 .2 .15 .1 .05 '((20 70) (45 95) (80 80)) t) ; (diversify 30 .4 t 1000 .2 .2 .1 .05 '((20 70) (45 95) (80 80)) t) ; (diversify 80 .4 t 1000 .2 .25 .1 .05 '((20 70) (45 95) (80 80)) t)) (defun diversify (num-trials ; number of runs at these parameters spec ; speciation rate &optional (save ()) ; save log files flag (max-div -1) ; number of leaves (species) (extinc -1) ; extinction rate (samp -1) ; sampling rate (para -1) ; proportion of taxa to species (k -1) ; parameter for taxa size (mass-ext ()) ; list of mass extinction events (save-trees ()) ; save trees flag (run-num -1)) ; repeat run number (when (>= run-num 0) (setf num-trials 1)) (dotimes (trial num-trials) (loop ; loops until a successful run. ;(close-files) (format t "Run number parameter: ~a ~%" run-num) (init-random run-num) (open-files) (initialize-tree spec extinc samp para k mass-ext max-div save save-trees) (if (run-trial) (progn (note-extant-sampled) (return *num-species*)) (undo-run *run-number*))) (save-first-epoch) ; Save the state of the first epoch (classify) (record-taxa-diversity mass-ext spec extinc) ;(record-results) (close-files))) (format t "4 ") ;---------------------------Initialization--------------------------- ;;; The mutation rate was set so that about half of the bits in the ;;; genome should flip by the end of the first epoch. (If it is a ;;; fully bushy tree, so probably more will flip.) ;;; Now it just defaults to twice the speciation rate, so a species on ;;; average will collect two character changes before splitting. ;;; I do some fancy manipulations to make sure that a species either speciates or ;;; goes extinct but not both. This is done by first checking for extinction. If ;;; the species survives it has a chance to speciate, but this chance has to be ;;; scaled up for the fact that the population has already been reduced by extinction. (defun initialize-tree (spec extinc samp para k mass-ext max-div save save-trees) (setf *save-info* save) (setf *save-trees* save-trees) (setf *extinction* (if (< extinc 0) (round (/ (* spec PROBSCALE) 2)) (round (* extinc PROBSCALE)))) (setf *speciation* (round (* PROBSCALE (/ (* spec PROBSCALE) (- PROBSCALE *extinction*))))) (setf *sample* (if (< samp 0) DEFAULTSAMPLE (round (* samp PROBSCALE)))) (setf *k* (if (< k 0) DEFAULTK k)) (setf *paraphyly-prob* (if (< para 0) DEFAULTPARAPROB (round (* para PROBSCALE)))) (setf *max-diversity* (if (< max-div 0) DEFAULTMAXDIV (1- max-div))) ; (setf *mutation* (if (< mut 0) ; (min PROBSCALE ; (round (* *speciation* 2))) ; (/ GENOMESIZE ; (* (log *max-diversity* 2) 2))))) ; (round (* mut PROBSCALE)))) ; (setf *min-death-count* (round (* amount-extinction ; *max-diversity*))) (setf *next-id* 0) (setf *next-leaf-id* 0) (setf *root* (create-snode ())) (iter-make-iter-list) (iter-insert-node *root*) (setf *num-species* 1) (setf *clock* 0) ; (setf *death-count* 0) (setf *extinct-sched* (process-extinctions mass-ext))) ;(save-header)) (defun process-extinctions (ext-list) (if (null ext-list) () (cons (process-an-extinction (first ext-list)) (process-extinctions (rest ext-list))))) (format t "5 ") ;-------------------------------Macros------------------------------- ;;; All macros are found in tree-macros.lisp... ;;; Create a new node, integrate it into the extant list, ;;; and increment the *num-species* count. Finally, return the new ;;; node so that the parent can set its appropriate child pointer. ;(defmacro add-snode (parent) ;;; Remove a node. You must have a pointer to the previous node in ;;; the *extant* list. ;(defmacro remove-next-snode (prev-snode) ;(defmacro speciate (prev-snode) ;(defmacro flip-bit (locus gene) ;(defmacro mutation-at (nde loc) ;(defmacro mutate (ptr-to-node) ;-----------------------------Running the Model---------------------- ;(defmacro step-trial () ;;; I use a dummy snode to head the list so that I can assume the ;;; first node won't be deleted. (defun run-trial (&optional (type ())) (do ((originations 0 (setf originations 0)) (extinctions 0 (setf extinctions 0)) (current-extinction 0) ; local variable for step-trial. (current-speciation 0)) ; local variable for step-trial. ((or (<= *num-species* 0) ; (>= *num-species* *max-diversity*)) (>= *next-leaf-id* *max-diversity*)) ;(if (= *num-species* *max-diversity*) t ())) (if (>= *next-leaf-id* *max-diversity*) (progn (iter-reset-list) ; consolodate the new species into the extants. t) ())) (step-trial) (when (= 0 (mod *clock* 10)) (format t "~a ~a ~a ~a~%" *clock* *num-species* originations extinctions)))) ;(when *save-info* ; (record-diversity *num-species* originations extinctions type)))) (format t "6 ") ;;; This makes a list equivalent to the *extant* list. It also saves ;;; the clock time. (defun save-first-epoch () (setf *first-epoch* (cons *clock* (copy-extant-list (iter-get-list) ())))) (defun copy-extant-list (nodes node-list) (if (null nodes) node-list (copy-extant-list (rest nodes) (cons (first nodes) node-list)))) ;;; *** This should probably also restore next-id, but that creates problems. (defun rewind-to-first-epoch () (setf *clock* (first *first-epoch*)) (setf *num-species* (1- (length *first-epoch*))) (iter-make-iter-list) (restore-extant (rest *first-epoch*))) ;(defmacro restore-node (node) (defun restore-extant (elist) (unless (null elist) (restore-node (first elist)) (restore-extant (rest elist)))) (defun unmutate (node mut-list) (unless (or (null mut-list) (<= (first-time mut-list) *clock*)) (mutation-at node (first-locus mut-list)) (unmutate node (rest mut-list)))) (defun note-extant-sampled () (dolist (node (iter-get-list)) (setf (snode-last-sample node) (1+ *clock*)) (when (< (snode-first-sample node) 0) (setf (snode-first-sample node) (1+ *clock*))))) (format t "done~%")