## Compare PFC white matter growth between humans and chimpanzees. ## The data are from Sakai et al., 2011 and from this study. ## These data consider the correction to the article published in Current Biology ## I extract age of growth cessation in humans and chimpanzees, and subsample values to ensure they are not driven by outliers ## Comparative analysis of dynamic time warping and epochs to ensure that the timetable of prefrontal cortex white matter growth holds up regardless of statistics used #library packages library(easynls) library(dtw) ##open file and set the working directory Chimp_Human_PFC<-read.table("TableS3.txt", header=TRUE, stringsAsFactors=FALSE) #Subset data tables by species Human<-subset(Chimp_Human_PFC, Species=="Human") Chimp<-subset(Chimp_Human_PFC, Species=="Chimpanzee") ##Plot of Prefrontal cortex volume versus age in days after conception plot(Human$Human_Age_days_post_conception, Human$Human_PFC_white_matter, log="y", col="darkred", pch=16, cex=3.5, cex.axis=2.5, ylab="PFC white matter cm3", xlab="Age in days after conception", cex.lab=2.5, xlim=c(200, 12000)) points(Chimp$Human_Age_days_post_conception, Chimp$Human_PFC_white_matter, col="cornflowerblue", pch=16, cex=3.5) ########################################################################## ###Assess when the prefrontal cortex white matter ceases to growth in humans test<-cbind.data.frame(Human$Human_Age_days_post_conception, log10(Human$Human_PFC_white_matter)) nlsfit(test, model=3) ##fit the model to test for a linear plateau abline(v=911.64920731, col="darkred") ##estimated age in which the growth of the PFC ceases to grow abline(v=270) ##birth in humans x<-c(300, 911.64920731, 12000) ##fit the model to the data y<-c(10^(1.64987070+0.00223226*(300-911.64920731)), 10^1.64987070, 10^1.64987070) ##fit the model to the data lines(x, y, col="darkred", lwd=2) ##fit the model to the data ###Assess when the prefrontal cortex white matter ceases to grow in chimpanzees test<-cbind.data.frame(Chimp$Human_Age_days_post_conception, log10(Chimp$Human_PFC_white_matter)) test<-na.omit(test) nlsfit(test, model=3) ##fit the model to test for linear plateau abline(v=2694.45393971, col="cornflowerblue") ##fit the model to the data x<-c(300, 2694.45393971, 12000) ##fit the model to the data y<-c(10^(1.05290393+0.00025848*(300-2694.45393971)), 10^1.05290393, 10^1.05290393) ##fit the model to the data lines(x, y, col="cornflowerblue", lwd=2) ##fit the model to the data ########################################################################## ##Capture epochs in prefrontal cortex growth in humans ##Capture 100% of adult volume (1.64987070*1-1.64987070)/0.00223226+911.64920731 ##Capture 90% of adult volume (1.64987070*0.9-1.64987070)/0.00223226+911.64920731 ##Capture 80% of adult volume (1.64987070*0.8-1.64987070)/0.00223226+911.64920731 ##Capture 70% of adult volume (1.64987070*0.7-1.64987070)/0.00223226+911.64920731 ##Capture 60% of adult volume (1.64987070*0.6-1.64987070)/0.00223226+911.64920731 ##Capture epochs in prefrontal cortex growith in chimpanzees ##Capture 100% of adult volume (1.05290393*1-1.05290393)/0.00025848+2694.45393971 ##Capture 90% of adult volume (1.05290393*0.9-1.05290393)/0.00025848+2694.45393971 ##Capture 80% of adult volume (1.05290393*0.8-1.05290393)/0.00025848+2694.45393971 ##Capture 70% of adult volume (1.05290393*0.7-1.05290393)/0.00025848+2694.45393971 ##Capture 60% of adult volume (1.05290393*0.6-1.05290393)/0.00025848+2694.45393971 ##Add these values to graph in humans abline(v=911.6492, col="darkred", lwd=2) abline(v=837.7389, col="darkred", lwd=2) abline(v=763.8286, col="darkred", lwd=2) abline(v=689.9182, col="darkred", lwd=2) abline(v=616.0079, col="darkred", lwd=2) ##Add these values to graph in chimpanzees abline(v=2694.454, col="cornflowerblue", lwd=2) abline(v=2287.109, col="cornflowerblue", lwd=2) abline(v=1879.765, col="cornflowerblue", lwd=2) abline(v=1472.421, col="cornflowerblue", lwd=2) abline(v=1065.076, col="cornflowerblue", lwd=2) ##Compare extracted epochs across humans and chimpanzees Human_epochs<-c(616.0079, 689.9182, 763.8286, 837.7389, 911.6492) Chimp_epochs<-c(1065.076, 1472.421, 1879.765, 2287.109, 2694.454) ##Plot the data plot(Human_epochs, Chimp_epochs, xlim=c(500, 2800), ylim=c(500, 2800), log="xy", col="plum4", cex=4, cex.axis=3, cex.lab=3, ylab="Chimpanzee", xlab="Human", pch=16) x<-c(500, 2000) y<-c(500, 2000) lines(x, y) ########################################################################## ###Subsample data to see how cessation of corpus callosum growth varies with sample size ##humans test<-cbind.data.frame(Human$Human_Age_days_post_conception, log10(Human$Human_PFC_white_matter)) dim(test) #sample size out1<-NULL outi<-NULL outR<-NULL outP1<-NULL outP2<-NULL outP3_1<-NULL outP3_2<-NULL outP3_3<-NULL for (i in (32:20)){ tryCatch({ set.seed(15) da<-test[sample(nrow(test), i), ] da<-as.data.frame(da) a<-nlsfit(da, model=3) a<-as.data.frame(a) Critical_point<-a[12,2] R_squared<-a[7,2] p1<-a[4,2] p2<-a[5,2] p3_1<-a[4,2] p3_2<-a[5,2] p3_3<-a[6,2] out1<-rbind(out1, Critical_point) outR<-rbind(outR, R_squared) outi<-rbind(outi, i) outP1<-rbind(outP1, p1) outP2<-rbind(outP2, p2) outP3_1<-rbind(outP3_1, p3_1) outP3_2<-rbind(outP3_2, p3_2) outP3_3<-rbind(outP3_3, p3_3) }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) } ##Make graphs to assess timing of cessation in growth plot(outi, (out1-270)/365, cex=4, col="darkred", pch=16, ylim=c(0, 5), xlab="Sample size", ylab="Age of plateau in years", cex.axis=3, cex.lab=3) da<-smooth.spline(outi, (out1-270)/365, df=4) lines(da$x, da$y, col="darkred", lwd=2) abline(h=quantile((out1-270)/365, probs = c(0.05, 0.95)), lty=2, col="darkred") quantile((out1-270)/365, probs = c(0.05, 0.95)) ##95% CI Human_plateau<-cbind.data.frame(outi, (out1-270)/365) Human_plateau_line<-cbind.data.frame(da$x, da$y) Human_plateau_quantile<-quantile((out1-270)/365, probs = c(0.05, 0.95)) ##95% CI ##Make graphs to assess how sample size impacts variance plot(outi, 100*outR, ylim=c(10, 100), cex=4, pch=16, col="darkred", cex.axis=3) da<-smooth.spline(outi, 100*outR, df=5) lines(da$x, da$y, col="darkred", lwd=2) abline(h=quantile(100*outR, probs = c(0.05, 0.95)), lty=2, col="darkred") quantile(100*outR, probs = c(0.05, 0.95)) ##95% CI ##Make graphs to assess how sample size impacts test statistics plot(outi, -log10(outP3_2), cex=4, col="grey", cex.axis=3, pch=16, ylim=c(0, 10)) abline(h=-log10(0.05)) TOG<-cbind.data.frame(outi, outP3_2) subset<-subset(TOG, TOG$outP3_2<0.05) points(subset[ ,1], -log10(subset[ ,2]), col="darkred", cex=4, pch=16) ########################################################################## ###Subsample the data to see how subsampling impacts of corpus callosum growth cessation in chimpanzees test<-cbind.data.frame(Chimp$Human_Age_days_post_conception, log10(Chimp$Human_PFC_white_matter)) dim(test) #sample size out1<-NULL outi<-NULL outR<-NULL outP1<-NULL outP2<-NULL outP3_1<-NULL outP3_2<-NULL outP3_3<-NULL for (i in (31:20)){ tryCatch({ set.seed(15) da<-test[sample(nrow(test), i), ] da<-as.data.frame(da) a<-nlsfit(da, model=3) a<-as.data.frame(a) Critical_point<-a[12,2] R_squared<-a[7,2] p1<-a[4,2] p2<-a[5,2] p3_1<-a[4,2] p3_2<-a[5,2] p3_3<-a[6,2] out1<-rbind(out1, Critical_point) outR<-rbind(outR, R_squared) outi<-rbind(outi, i) outP1<-rbind(outP1, p1) outP2<-rbind(outP2, p2) outP3_1<-rbind(outP3_1, p3_1) outP3_2<-rbind(outP3_2, p3_2) outP3_3<-rbind(outP3_3, p3_3) }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) } ##Make graphs to assess timing of cessation in growth plot(outi, (out1-243)/365, cex=4, col="cornflowerblue", pch=16, ylim=c(0, 10), xlab="Sample size", ylab="Age of plateau in years", cex.axis=3, cex.lab=3) da<-smooth.spline(outi, (out1-270)/365, df=4) lines(da$x, da$y, col="cornflowerblue", lwd=2) points(Human_plateau[ ,1], Human_plateau[ ,2], col="darkred", cex=4, pch=16) lines(Human_plateau_line[ ,1], Human_plateau_line[ ,2], col="darkred") abline(h=Human_plateau_quantile, col="darkred", lty=2) abline(h=quantile((out1-270)/365, probs = c(0.05, 0.95)), lty=2, col="cornflowerblue") quantile((out1-270)/365, probs = c(0.05, 0.95)) ##95% CI ##Make graphs to assess how sample size impacts variance plot(outi, 100*outR, ylim=c(10, 100), cex=4, pch=16, col="cornflowerblue", cex.axis=3) da<-smooth.spline(outi, 100*outR, df=5) lines(da$x, da$y, col="cornflowerblue", lwd=2) abline(h=quantile(100*outR, probs = c(0.05, 0.95)), lty=2, col="cornflowerblue") quantile(100*outR, probs = c(0.05, 0.95)) ##95% CI ##Make graphs to assess how sample size impacts test statistics plot(outi, -log10(outP3_2), cex=4, col="cornflowerblue", cex.axis=3, pch=16, ylim=c(0, 10)) abline(h=-log10(0.05)) TOG<-cbind.data.frame(outi, outP3_2) subset<-subset(TOG, TOG$outP3_2<0.05) points(subset[ ,1], -log10(subset[ ,2]), col="cornflowerblue", cex=4, pch=16) ########################################################################## ###Apply dynamic time warping to the PFC growth curves in humans and chimpanzees ##Format data test<-cbind.data.frame(Chimp$Human_Age_days_post_conception, log10(Chimp$Human_PFC_white_matter)) test1<-cbind.data.frame(Human$Human_Age_days_post_conception, log10(Human$Human_PFC_white_matter)) ##Extract corresponding time points across humans and chimpanzees to yeild the same number of PFC measurements across the two species plot(test1[,1], test1[ ,2]) da<-smooth.spline(test1[,1], test1[ ,2]) lines(da$x, da$y) Translated<-predict(da, test[ ,1]) ##interpolate values from the smooth spline to compare PFC at equivalent ages acrosss the two species ##Apply dynamic time warping to find corresponding time poins from PFC growth curves Chimp_Human<-cbind.data.frame(Translated, test[ ,2]) alignment<-dtw(Chimp_Human[ ,2], Chimp_Human[ ,3], keep=TRUE); Human_alignment<-cbind.data.frame(1:nrow(Chimp_Human), test[ ,1]) Chimp_alignment<-cbind.data.frame(1:nrow(Chimp_Human), test[ ,1]) Cbind1<-cbind.data.frame(alignment$index1, alignment$index2) colnames(Chimp_alignment)[1]<-c("Chimp") colnames(Human_alignment)[1]<-c("Human") colnames(Cbind1)<-c("Human", "Chimp") merge_1<-merge(Cbind1, Chimp_alignment, by="Chimp") merge_2<-merge(merge_1, Human_alignment, by="Human") ##Fit a smooth spline through time points captured from dynamic time warping da<-smooth.spline(merge_2[ ,3]~merge_2[ ,4]) Translated<-predict(da, unique(merge_2[ ,3])) ##Plot the data, dynamic time warping plot(Translated$x, Translated$y, log="xy", col="plum4", cex=4, pch=16, xlim=c(300, 20000), ylim=c(300, 20000), cex.axis=3, cex.lab=3, xlab="Human", ylab="Chimpanzee") x<-c(100, 20000) y<-c(100, 20000) lines(x, y, col="grey")