--- title: "Scaling and development of elastic mechanisms: the tiny strikes of larval mantis shrimp" author: "J.S.Harrison" date: "5/26/2020" output: html_document --- ```{r} #Load in all libraries necessary to run code library(plyr) library(readr) library(reshape2) library(scales) library(dplyr) library(pracma) library(ggplot2) ``` ```{r} #This is the method by which we converted the .mdf output from MTrackJ into '.csv', and converted the raw positional data into angular displacement #Within the Dryad submission there are the raw '.mdf' files and the processed '.csv' files read.one.mdf <- function(file) { ## mdf â output file of MTrackJ ## MTrackJ â a plugin of ImageJ useful for digitizing ## return 'track' â format for digitized points track.df <- read.table(file, skip = 4, fill = T) d <- data.frame() for (i in 1:dim(track.df)[1]) { line <- track.df[i,] if (line[1] == "Track") { track.id <- as.numeric(as.character(line[1, 2])) next } if (line[1] == "Point") { point.id <- as.numeric(as.character(line[1, 2])) point.x <- as.numeric(as.character(line[1, 3])) point.y <- as.numeric(as.character(line[1, 4])) frame.id <- as.numeric(as.character(line[1, 6])) d <- rbind( d, data.frame( track.id = track.id, point.id = point.id, point.x = point.x, point.y = point.y, frame.id = frame.id ) ) } } return(d) } get.rotation <- function(track, target.v = c(1,2), ref.v = c(3,4)) { ## current method to obtain angle made by two lines ## The the information for the two lines must be specified by four points. ## The four points are consisted of two points for one target line ("target.v"), and ## the other two points for one reference line ("ref.v"). focus.id <- c(target.v, ref.v) trackt <- data.frame() for (i in focus.id) { trackt <- rbind(trackt, subset(track, track.id == i)) } dmin <- ddply(trackt, .(track.id), subset, frame.id == min(frame.id)) dmax <- ddply(trackt, .(track.id), subset, frame.id == max(frame.id)) start.frame <- max(dmin$frame.id) end.frame <- min(dmax$frame.id) d <- subset(track, frame.id >= start.frame & frame.id <= end.frame) p1 <- subset(d, track.id == target.v[1], select = c("point.x", "point.y")) p2 <- subset(d, track.id == target.v[2], select = c("point.x", "point.y")) p3 <- subset(d, track.id == ref.v[1], select = c("point.x", "point.y")) p4 <- subset(d, track.id == ref.v[2], select = c("point.x", "point.y")) v1 <- p2 - p1 v0 <- p4 - p3 theta <- atan2(v1$point.y, v1$point.x) * 180 / pi + 180 theta0 <- atan2(v0$point.y, v0$point.x) * 180 / pi + 180 for(i in 1:length(theta)){ if(theta[i] > 270 & theta[i] < 360){ theta[i] <- theta[i] - 360 } } for(i in 1:length(theta0)){ if(theta0[i] > 270 & theta0[i] < 360){ theta0[i] <- theta0[i] - 360 } } rot <- theta - theta0 rot <- rot - rot[1] return(rot) } ``` ```{r} #Here is the code for calculating the kinematics of individual strikes from rotation data saved as .csv files setwd("~/Desktop/G.falcatus_DevelopmentandKinematicsProject/DryadSubmission/StrikeKinematicsData/csv files/") rotationfolder<-"~/Desktop/G.falcatus_DevelopmentandKinematicsProject/DryadSubmission/StrikeKinematicsData/csv files/" #Replace the above code with where you have the Dryad files saved on your directory rotationfiles <- list.files(rotationfolder, pattern="*.csv") filestoanalyze <- lapply(rotationfiles, function(x) {read.table(file =x, header = T, sep =",")}) setwd("~/Desktop/G.falcatus_DevelopmentandKinematicsProject/DryadSubmission/") #Replace this code with where you have the files saved on your directory results=NULL for(j in 1:length(rotationfiles)){ rotationdata <- filestoanalyze[[j]] IDdata <- read.csv("Larval.ID_forkinematics.csv") scale_m = IDdata$pixelconv_pixels[j]*1000 st.move= IDdata$startmove_frame[j] end.move= IDdata$endmove_frame[j] fps=IDdata$fps[j] p = 1022 fps = 20000 mu = 0.00087 #Ns/m2, at 30 degrees celsius polyn=6 #This sets up the polynomial fitting model to the kinematic data smooth <- function(y){lm(y ~ poly(x=1:length(y), polyn, raw=TRUE))} rot <- rotationdata$Angle[st.move:end.move] model <- smooth(rot) predicted <- predict(model,data.frame(x=1:length(rot)), interval='confidence', level=0.95) time_s= seq(from=0, length.out= length(rot), by = 1/fps) plot(time_s, rot) #plots the rotation and the polynomial model to check fit points(time_s, predicted[,1], pch=".", col="red", size=2, lty=1) #Pulls derivatives of the polynomial model Y <- predicted[,1] X <- time_s dY <- diff(Y)/diff(X) dX <- X[-length(X)] ddY<- diff(dY)/diff(dX) ddX <- dX[-length(dX)] #Extracts values for individual strikes displacement.angle=rotationdata$Angle[st.move:end.move] displacement.rad<- max(displacement.angle)*(pi/180) #Converts degree values into radians duration.s=length(rot)/fps #Extracts peak velocities/accelerations from models and converts to radians peak.rot.vel_rad<- max(dY)*(pi/180) peak.rot.acc_rad <- max(ddY)*(pi/180) #This converts the rotational kinematics into linear kinematics peak.lin.vel <- peak.rot.vel_rad*IDdata$strikingbody[j]/1000 peak.lin.acc <- peak.rot.acc_rad*IDdata$strikingbody[j]/1000 #Reynold's number calculations Re <- (p*peak.lin.vel*(IDdata$strikingbody[j]/1000))/mu results <- rbind(results, data.frame(rotationfiles[j], IDdata$LarvalNumber_formanuscript[j], IDdata$Strike[j], IDdata$fps[j], IDdata$strikingbody_mm[j], IDdata$rotatemass_kg[j], duration.s, displacement.rad, peak.rot.vel_rad, peak.rot.acc_rad, peak.lin.vel, peak.lin.acc, Re)) } results ``` ```{r} #Extracts means, min, max, sd for each individual Ind <- factor(results$IDdata.LarvalNumber_formanuscript.j.) Ind.means = aggregate(results[,5:13],list(Ind), mean, na.rm = TRUE) #means for each individual Ind.sd = aggregate(results[,5:13],list(Ind), sd, na.rm = TRUE) #SD for each individual Ind.min = aggregate(results[,5:13],list(Ind), FUN = "min", na.rm = TRUE) Ind.max = aggregate(results[,5:13],list(Ind), FUN = "max", na.rm = TRUE) counts = data.frame(table(results$IDdata.ID.j.)) Ind.means #Pulls out overall means, min, max, sd for all larvae from individual means overallmean = sapply(Ind.means[,2:10], mean, na.rm = TRUE) overallsd = sapply(Ind.means[,2:10], sd, na.rm = TRUE) overallmin = sapply(Ind.min[,2:10], FUN = "min", na.rm = TRUE) overallmax = sapply(Ind.max[,2:10], FUN = "max", na.rm = TRUE) overallmean overallsd ``` ```{r} #This section looks at the larval kinematic data, starting with histograms of angular velocity and acceleration and #follows with OLS regressions. ggplot(Ind.means, aes(x=peak.rot.vel_rad)) + geom_histogram(binwidth = 100) #Relatively large bins, feel free to make them smaller ggplot(Ind.means, aes(x=peak.rot.acc_rad)) + geom_histogram(binwidth = 200000) #OLS regressions comparing striking body length, duration, and displacement with the angular kinematics (velocity and acceleration) ggplot(Ind.means, aes(x=IDdata.strikingbody_mm.j.,y=peak.rot.vel_rad)) + geom_point()+ geom_smooth(method="lm") m1 <- lm(Ind.means$peak.rot.vel_rad~Ind.means$IDdata.strikingbody_mm.j, data = Ind.means) anova(m1) summary(m1) ggplot(Ind.means, aes(x=IDdata.strikingbody_mm.j.,y=peak.rot.acc_rad)) + geom_point()+ geom_smooth(method="lm") m2 <- lm(Ind.means$peak.rot.acc_rad~Ind.means$IDdata.strikingbody_mm.j, data = Ind.means) anova(m2) summary(m2) ggplot(Ind.means, aes(x=duration.s,y=peak.rot.vel_rad)) + geom_point()+ geom_smooth(method="lm") m3 <- lm(Ind.means$peak.rot.vel_rad~Ind.means$duration.s, data = Ind.means) anova(m3) #Significant: p = 0.015 summary(m3) ggplot(Ind.means, aes(x=duration.s,y=peak.rot.acc_rad)) + geom_point()+ geom_smooth(method="lm") m4 <- lm(Ind.means$peak.rot.acc_rad~Ind.means$duration.s, data = Ind.means) anova(m4) summary(m4) ggplot(Ind.means, aes(x=displacement.rad,y=peak.rot.vel_rad)) + geom_point()+ geom_smooth(method="lm") m5 <- lm(Ind.means$peak.rot.vel_rad~Ind.means$displacement.rad, data = Ind.means) anova(m5) summary(m5) ggplot(Ind.means, aes(x=displacement.rad,y=peak.rot.acc_rad)) + geom_point()+ geom_smooth(method="lm") m6 <- lm(Ind.means$peak.rot.acc_rad~Ind.means$displacement.rad, data = Ind.means) anova(m6) summary(m6) ``` ```{r} #Measurement error from manual tracking of videos setwd("~/Desktop/G.falcatus_DevelopmentandKinematicsProject/DryadSubmission/KinematicErrorData") error_rotationfolder<-"~/Desktop/G.falcatus_DevelopmentandKinematicsProject/DryadSubmission/KinematicErrorData" #Replace the above code with where you have the Dryad files saved on your directory error_rotationfiles <- list.files(error_rotationfolder, pattern="*.csv") filestoanalyze <- lapply(error_rotationfiles, read_csv) rotation_error <- lapply(filestoanalyze, get.rotation) error_data<-as.data.frame(t(ldply(rotation_error, rbind()))) error_data<- mutate(error_data, Frame=1:102) kinematics= function(x){ rotation=x[1:95] time= seq(from=0, length.out= length(rotation), by = 1/fps) polyn=6 smooth <- function(y){lm(y ~ poly(x=1:length(y), polyn, raw=TRUE))} model <- smooth(rotation) predicted <- predict(model,data.frame(x=1:length(rotation)), interval='confidence', level=0.95) time_s= seq(from=0, length.out= length(rotation), by = 1/20000) Y <- predicted[,1] X <- time_s dY <- diff(Y)/diff(X) dX <- X[-length(X)] ddY<- diff(dY)/diff(dX) ddX <- dX[-length(dX)] return(c(Y=max(Y)*pi/180, dY=max(dY)*pi/180, ddY=max(ddY)*pi/180)) } Error <- as.data.frame(t(sapply(error_data, kinematics))) original.measurement <- data.frame(Y=1.603505,dY= 762.81722, ddY=846554.54) Error <- rbind(Error, original.measurement) Summary.Error<-summarise(Error, meanY=mean(Y), seY=sd(Y)/sqrt(length(Y)), meandY=mean(dY), sedY=sd(dY)/sqrt(length(dY)), meanddY=mean(ddY), seddY=sd(ddY)/sqrt(length(ddY))) Summary.Error ``` ```{r} #Measure uncertainty from instruments, used to set significant digits for reported numbers #These are limits of instrumentation used to collect kinematics #DEGREE #Minimum resolution of angle in FIJI depends on resolution of image. #Minimum resolution of angle here is based on average distance in pixels between points tracked on merus #SigFigs - report to 0.69 degrees or 0.01 radians uncertDegree = 0.69/2 uncertRad = (0.69/2)*(pi/180) #DISTANCE #Ruler recalibrated to 0.1 mm using a 0.02 mm scale #SigFigs - report to 0.1 mm #Uncertainty: ± 0.05 mm or 0.00005 m uncertDISTm = 0.00005 #MASS #Masses were measured with an analytical balance (resolution: 0.1 µg; XPE56, Mettler Toledo, Pleasant Prairie, WI, USA).] #Sig figs: report the mass to 0.001 mg or 1 µg #Uncertainty: ±0.0005 mg uncertMASSmg = 0.0005 #TIME #Time uncertainty is +/- 1/fps This is ACCURACY of frame rate provided by the calibration of the video camera, calculated as 1/framerate #Sig figs: Report strike duration in increments of 0.05 ms or 0.00005 s #Uncertainty: ±0.05 ms uncertTIMEs = 1/(20000) #Here we are calculating the uncertainty of our measurements based on average kinematics from all individuals #Rotation= Uncertainty ratio is 0.5%, report to 0.001 rad avg_reported_rotation = 9.364337e-01 ratio_rotation = uncertRad/avg_reported_rotation mag_rot = avg_reported_rotation* ratio_rotation ratio_rotation mag_rot #DURATION = avg_reported_duration = 8.747917e-03 ratio_duration = uncertTIMEs/avg_reported_duration #T/t #this is a uncertainty ratio mag_dur = avg_reported_duration* ratio_duration ratio_duration mag_dur #ROTATIONAL VELOCITY = Uncertainty ratio is 0.8%, report to 0.1 rad/s avg_reported_rot = 9.364337e-01 avg_reported_rot_vel = 2.793349e+02 ratio_rot_v = sqrt((uncertRad/avg_reported_rot)^2 + (uncertTIMEs/avg_reported_duration)^2) #this is uncertainty ratio mag_rot_v = avg_reported_rot_vel* ratio_rot_v #this is magnitude of uncertainty ratio_rot_v mag_rot_v #ROTATIONAL ACCELERATION = Uncertainty ratio is 1.3%, report to 100 rad/s^2 avg_reported_rot_acc = 4.690133e+05 uncerttimeSQ = 2*(uncertTIMEs/avg_reported_duration) ratio_rot_a = sqrt((uncertRad/avg_reported_rot)^2 + (uncerttimeSQ)^2) #this is uncertainty ratio mag_rot_a = avg_reported_rot_acc* ratio_rot_a ratio_rot_a mag_rot_a #SPEED = Uncertainty ratio is 4%, report to 0.001 m/s avg_distance = 9.364337e-01*0.001329250 #avg rotation multiplied by avg striking body length avg_reported_speed = 3.686288e-01 ratio_speed = sqrt((uncertDISTm/avg_distance)^2 + (uncertTIMEs/avg_reported_duration)^2) #this is uncertainty ratio mag_speed = avg_reported_speed* ratio_speed #this is the magnitude of uncertainty ratio_speed mag_speed #LINEAR ACCELERATION = Uncertainty ratio is 4%, report to 1 m/s^2 avg_reported_lin_acc = 6.192120e+02 ratio_lin_a = sqrt((uncertDISTm/avg_distance)^2 + (uncerttimeSQ)^2) #this is uncertainty ratio mag_lin_a = avg_reported_lin_acc* ratio_lin_a ratio_lin_a mag_lin_a ``` ```{r} #All code below this point were used to make plotsfor manuscript #Figures comparing larval mantis shrimp kinematics to adult mantis shrimp setwd("~/Desktop/G.falcatus_DevelopmentandKinematicsProject/DryadSubmission/") #Replace this code with where you have the Dryad files saved on your directory comparativedata <- read.csv("MantisShrimpComparativeKinematics.csv") base_breaks <- function(n = 10){ function(x) { axisTicks(log10(range(x, na.rm = TRUE)), log = TRUE, n = n) } } AngVelFig <- ggplot(data=comparativedata, aes(x=StrkBody,y=MaxAngVel, colour=Taxon, shape=Type, Size=2)) + geom_point(aes(size=1.5)) + scale_shape_manual(values=c(15, 17)) + scale_y_continuous(trans = log_trans(), breaks = base_breaks()) + scale_x_continuous(trans = log_trans(), breaks = base_breaks())+ theme_test()+ ylab("Angular Velocity (rad/s)") + xlab("Striking Body Length (m)") AngVelFig AngAccelFig <- ggplot(data=comparativedata, aes(x=StrkBody,y=MaxAngAcc, colour=Taxon, shape=Type)) + geom_point(aes(size=1.5)) + scale_shape_manual(values=c(15, 17)) + scale_y_continuous(trans = log_trans(), breaks = base_breaks()) + scale_x_continuous(trans = log_trans(), breaks = base_breaks())+ theme_test() + ylab("Angular Acceleration (rad/s^2)") + xlab("Striking Body Length (m)") AngAccelFig LinVelFig <- ggplot(data=comparativedata, aes(x=StrkBody,y=MaxLinVel, colour=Taxon, shape=Type)) + scale_shape_manual(values=c(15, 17)) + geom_point(aes(size=1.5))+ scale_y_continuous(trans = log_trans(), breaks = base_breaks()) + scale_x_continuous(trans = log_trans(), breaks = base_breaks())+ theme_test() LinVelFig + ylab("Speed (m/s)") + xlab("Striking Body Length (m)") LinAccelFig <- ggplot(data=comparativedata, aes(x=StrkBody,y=MaxLinAcc, colour=Taxon, shape=Type)) + geom_point(aes(size=1.5)) + scale_y_continuous(trans = log_trans(), breaks = base_breaks()) + scale_shape_manual(values=c(15, 17)) + scale_x_continuous(trans = log_trans(), breaks = base_breaks()) + theme_test() + ylab("Linear Acceleration (m/s^2)") + xlab("Striking Body Length (m)") LinAccelFig ``` ```{r} #Figure comparing larval mantis shrimp strike speed to other small pelagic swimmers setwd("~/Desktop/G.falcatus_DevelopmentandKinematicsProject/DryadSubmission/") #Replace this code with where you have the Dryad files saved on your directory Larval_Speeds <- read.csv("LarvalSpeedComparisons.csv") speed.data <- Larval_Speeds p1 <- ggplot(data=speed.data, aes(x=speed.data$Body.Length..mm., y=speed.data$Speed..mm.s., shape = Behavior, fill = Family)) + scale_shape_manual(values=c(15, 17)) + geom_point(size=4) + xlab("Total Length (mm)") + ylab("Swimming Speed (mm/s)") + theme_test() fColors <- speed.data$Color names(fColors) <- speed.data$Family head(fColors) p1 + aes(color = Family) + scale_color_manual(values = fColors) ```