Figure 2

NOTE: Run all the code included on the “Base Code” page prior to running the following.

Setup

# working directory
setwd("")  # enter appropriate working directory

# load eigenworms [Broekmans et al. 2016]
ew = read.csv("eigenworms.csv", header=F, sep=",")

# specify EDM parameters/ranges
Es = seq(1,8,1)
thetas = seq(0,3,0.5)
tau = 1
tp = 1
lib = c(1,1000)
pred = c(1001,2000)
maxE = 8

N2 Wildtype foraging worms [Broekmans et al 2016]

# working directory
setwd("")  # enter appropriate working directory

# filenames
filenames = c(
  '1.txt',
  '2.txt',
  '3.txt',
  '4.txt',
  '5.txt',
  '6.txt',
  '7.txt',
  '8.txt',
  '9.txt',
  '10.txt',
  '11.txt',
  '12.txt'
)

datasections = rbind(
  c(10001,NA),  #worm1
  c(10001,NA),  #worm2
  c(8501,NA),   #worm3
  c(10501,NA),  #worm4
  c(10501,NA),  #worm5
  c(10501,NA),  #worm6
  c(10501,NA),  #worm7
  c(10501,NA),  #worm8
  c(11001,NA),  #worm9
  c(11001,NA),  #worm10
  c(11001,NA),  #worm11
  c(11001,NA)   #worm12
)

N2out = sensitivity_main(filenames, datasections, Es, thetas, tau, tp, lib, pred, sep="", maxE)

LSJ1

# working directory
setwd("")  # enter appropriate working directory

filenames_LSJ1 = c('1 300 LS51 on food L_2010_11_26__16_25___3___13_features.csv',
                   '2 300 LS51 on food L_2011_02_17__11_36___3___3_features.csv',
                   '3 300 LS51 on food R_2010_11_26__11_47___3___5_features.csv',
                   '4 300 LS51 on food L_2011_02_17__12_48___3___7_features.csv',
                   '5 300 LS51 on food L_2010_11_26__10_38_21___7___1_features.csv',
                   '6 300 LS51 on food R_2011_02_25__10_31_11___7___1_features.csv',
                   '8 300 LSJ1 on food L_2011_03_10__10_52_16___6___1_features.csv',
                   '10 300 LSJ1 on food R_2011_03_07__11_42_35___6___3_features.csv',
                   '11 300 LSJ1 on food R_2011_04_12__12_30_20___6___5_features.csv',
                   '12 300 LSJ1on food R_2011_03_15__12_37_50___6___5_features.csv')

datasections_LSJ1 = rbind(
  c(1001,NA), #worm1
  c(1001,NA), #worm2
  c(1001,NA), #worm3
  c(1001,NA), #worm4
  c(1001,NA), #worm5
  c(1001,NA), #worm6
  c(1001,NA), #worm8
  c(1001,NA), #worm10
  c(1001,NA), #worm11
  c(1001,NA)  #worm12
)

# test data read
# data = read.csv(filenames_LSJ1[6], header=F, sep=",")
# data = data[datasections_LSJ1[6,1]:datasections_LSJ1[6,2],1:5]

LSJ1out = sensitivity_main(filenames_LSJ1, datasections_LSJ1, Es, thetas, tau, tp, lib, pred,
                           sep=",", maxE)

dpy-20

# working directory
setwd("")  # enter appropriate working directory

filenames_dpy20 = c('1 dpy-20 (e1282)IV on food L_2011_08_04__11_30_01___2___4_features.csv',
                    '2 dpy-20 (e1282)IV on food L_2011_08_19__12_00_59___6___5_features.csv',
                    '3 dpy-20 (e1282)IV on food R_2011_08_04__12_51_11___6___6_features.csv',
                    '4 dpy-20 (e1282)IV on food R_2011_08_09__12_58_43___6___8_features.csv',
                    '5 dpy-20 (e1282)IV on food R_2011_08_24__11_21_59___6___5_features.csv',
                    '6 dpy-20 (e1282)IV on food R_2011_08_25__15_40_16___6___9_features.csv',
                    '7 dpy-20 (e1282)IV on food L_2011_08_04__10_52_51___8___2_features.csv',
                    '8 dpy-20 (e1282)IV on food L_2011_08_24__12_07_53___8___7_features.csv',
                    '9 dpy-20 (e1282)IV on food L_2011_08_25__11_23_10___8___1_features.csv',
                    '10 dpy-20 (e1282)IV on food R_2011_08_09__10_54_50___8___2_features.csv',
                    '11 dpy-20 (e1282)IV on food R_2011_08_19__12_43_48___8___7_features.csv',
                    '12 dpy-20 (e1282)IV on food L_2011_08_09__10_33_28___7___1_features.csv')

datasections_dpy20 = rbind(
  c(1001,NA), #worm1
  c(1001,NA), #worm2
  c(1001,NA), #worm3
  c(1001,NA), #worm4
  c(1001,NA), #worm5
  c(1001,NA), #worm6
  c(1001,NA), #worm7
  c(1001,NA), #worm8
  c(1001,NA), #worm9
  c(1001,NA), #worm10
  c(1001,NA), #worm11
  c(1001,NA)  #worm12
)

# test data read
# data = read.csv(filenames_dpy20[1], header=F, sep=",")

dpy20out = sensitivity_main(filenames_dpy20, datasections_dpy20, Es, thetas, tau, tp, lib, pred,
                            sep=",", maxE)

octr-1

# working directory
setwd("")  # enter appropriate working directory

filenames_octr1 = c('1 tag-24 (ok371)X on food R_2010_01_22__15_03_36___2___12_features.csv',
                    '2 tag-24 (ok371)X on food L_2010_01_28__13_19_44___2___9_features.csv',
                    '3 tag-24 (ok371)X on food L_2010_01_26__11_26___3___3_features.csv',
                    '4 tag-24 (ok371)X on food R_2010_01_28__13_19___3___9_features.csv',
                    '5 tag-24 (ok371)X on food L_2010_01_22__15_02_47___1___12_features.csv',
                    '6 tag-24 (ok371)X on food R_2010_01_28__13_19_32___1___9_features.csv',
                    '7 tag-24 (ok371)X on food R_2010_01_22__15_04_07___4___12_features.csv',
                    '8 tag-24 (ok371)X on food R_2010_01_26__11_27_13___4___3_features.csv',
                    '9 tag-24 (ok371)X on food R_2010_01_28__13_20_16___4___10_features.csv',
                    '10 tag-24 (ok371)X on food L_2010_01_22__15_05_37__10_features.csv',
                    '11 tag-24 (ok371)X on food L_2010_01_28__13_21_42__1_features.csv',
                    '12 tag-24 (ok371)X on food R_2010_01_26__11_28_44__2_features.csv')

datasections_octr1 = rbind(
  c(1001,NA), #worm1
  c(1001,NA), #worm2
  c(1001,NA), #worm3
  c(1001,NA), #worm4
  c(1001,NA), #worm5
  c(1001,NA), #worm6
  c(1001,NA), #worm7
  c(1001,NA), #worm8
  c(1001,NA), #worm9
  c(1001,NA), #worm10
  c(1001,NA), #worm11
  c(1001,NA)  #worm12
)

# test data read
# data = read.csv(filenames_octr1[1], header=F, sep=",")

octr1out = sensitivity_main(filenames_octr1, datasections_octr1, Es, thetas, tau, tp, lib, pred,
                            sep=",", maxE)

unc-80

# working directory
setwd("")  # enter appropriate working directory

filenames_unc80 = c('1 unc-80 (e1069) on food L_2010_04_14__14_57_39___2___12_features.csv',
                    '2 unc-80 (e1069)V on food L_2009_12_09__16_36_22___2___15_features.csv',
                    '3 unc-80 (e1069)V on food R_2009_12_10__13_43_33___2___8_features.csv',
                    '4 unc-80 (e1069) on food L_2010_04_14__14_57___3___12_features.csv',
                    '5 unc-80 (e1069)V on food R_2009_12_09__16_36___3___15_features.csv',
                    '6 unc-80 (e1069)V on food R_2009_12_10__13_43___3___8_features.csv',
                    '7 unc-80 (1069)V on food R_2009_12_10__13_42_50___1___9_features.csv',
                    '8 unc-80 (e1069) on food L_2010_04_14__14_56_45___1___12_features.csv',
                    '9 unc-80 (e1069)V on food L_2009_12_09__16_35_38___1___15_features.csv')

datasections_unc80 = rbind(
  c(1001,NA), #worm1
  c(1001,NA), #worm2
  c(1001,NA), #worm3
  c(1001,NA), #worm4
  c(1001,NA), #worm5
  c(1001,NA), #worm6
  c(1001,NA), #worm7
  c(1001,NA), #worm8
  c(1001,NA)  #worm9
)

unc80out = sensitivity_main(filenames_unc80, datasections_unc80, Es, thetas, tau, tp, lib, pred,
                            sep=",", maxE)

Figure-Specific Functions

Main Function

sensitivity_main <- function(filenames, datasections, Es, thetas, tau, tp, lib, pred, sep, maxE){
  
  # generate placeholder for output
  main_output = array(NA, dim=c(length(Es), length(thetas), length(filenames)))
  
  # loop over files
  for (i in c(1:length(filenames))) {
    print(filenames[i])
    print(datasections[i,])
    data = read.csv(filenames[i], header=F, sep=sep)
    data = data[datasections[i,1]:dim(data)[1],1:5]
    main_output[,,i] = run_sensitivity_loop(data, lib, pred, Es, thetas, tau, tp, maxE)
  }
  return(main_output)
}

Function to Test Sensitivity Looping over Range of Es and thetas

run_sensitivity_loop <- function(data, lib, pred, Es, thetas, tau, tp, maxE){
  ###inputs###
  #data: 5-dimensional eigenworm data
  #lib: vector of length 2 denoting start and end rows for library within data
  #pred: vector of length 2 denoting start and end rows for prediction within data
  #Es: vector of E values to test
  #thetas: vector of theta values to test
  
  ###returns###
  #output: array of size length(Es) by length (thetas) containing mean of...
  #...rms prediction error over the interval specified by pred for each of...
  #...the values of E and theta specified by the user.
  ###
  
  # setup output
  output = array(NA, dim=c(length(Es), length(thetas)))
  
  # generate base embedding
  base_embedding = make_embed_skipna(data, lib, pred, maxE, tau, tp)
  
  # set initial outer counter
  i = 1
  
  # outer loop: E
  for(E in Es){
    
    # make embedding
    matrix = base_embedding[,1:((E+1)*5)]  # assumes 5 eigenworms
    
    # set initial inner counter
    j = 1
    
    # inner loop: theta
    for(theta in thetas){
      
      # make prediction
      new_pred <- make_pred(E, matrix, theta, lib, pred)
      observations_total = new_pred$obs
      predictions_total = new_pred$pred
      
      # calculate error
      errors <- edm_error_noangles(observations_total, predictions_total)
      
      # store mean error in output
      output[i,j] = mean(errors, na.rm=TRUE)
      
      # update inner counter
      j = j + 1
    }
    
    # update outer counter
    i = i + 1
  }
  
  return(output)
}

Function to Generate Valid Embeddings without NaN Values

# generate embedding that doesn't include any vectors with NA values
make_embed_skipna <- function(data, lib, pred, E, tau, tp){
  
  # generate a standard embedding
  embed = make_embed(data, E, tau, tp)
  
  # remove all rows of the embedding that contain NA values
  for(i in c(1:dim(embed)[2])){
    embed = embed[!is.na(embed[,i]),]
  }
  
  # generate library embedding
  liblength = lib[2] - lib[1] + 1
  lib_embed = embed[1:liblength,]
  
  # generate prediction interval embedding with no library overlap
  predlength = pred[2] - pred[1] + 1
  pred_embed = embed[(liblength+(E*tau)+tp+1):(liblength+(E*tau)+tp+predlength),]
  
  # concatenate embeddings
  output = rbind(lib_embed, pred_embed)
  
  return(output)
}

Function to Calculate RMS Error on Eigenworm Coeffs (not angles)

# calculate error from raw eigenworm coefficients
edm_error_noangles <- function(observations_total, predictions_total){
  errors = (rowMeans((observations_total - predictions_total)^2))^(0.5)
  return(errors)
}

Function to Generate Time Series of Worm Body Angle RMS Error Looping over E and theta Sequence (Note: This uses a standard continuous embedding.)

run_sensitivity_simple <- function(data, lib, pred, Eandtheta, tau, tp){
  
  # setup output
  output = array(NA, dim=c(dim(Eandtheta)[1],1000))
  
  for(i in c(1:dim(Eandtheta)[1])){
    
    E = Eandtheta[i,1]
    theta = Eandtheta[i,2]
    
    # make embedding
    matrix1 <- make_embed(data[lib[1]:lib[2],], E, tau, tp)
    matrix2 <- make_embed(data[pred[1]:pred[2],], E, tau, tp)
    matrix = rbind(matrix1, matrix2)
    
    # make prediction
    new_pred <- make_pred(E, matrix, theta, lib, pred)
    observations_total = new_pred$obs
    predictions_total = new_pred$pred
    
    # calculate error
    output[i,] <- edm_error(pred, observations_total, ew, predictions_total)
  }
  
  return(output)
}

Plotting

Fig. 2 Plot

# main figure 2
m = matrix(c(1:8), nrow=2, ncol=4, byrow=TRUE)
layout(m)
par(mai=c(0.35,0.35,0.35,0.1), omi=c(0,0,0,0), cex=1, mgp=c(1.5,0.5,0))

eplots(N2out)
eplots(octr1out)
eplots(unc80out)
eplots(dpy20out)

thetaplots(N2out)
thetaplots(octr1out)
thetaplots(unc80out)
thetaplots(dpy20out)

Fig. S1 Plot

# figure 2 supplement 1
m = matrix(c(1:8), nrow=2, ncol=4, byrow=TRUE)
layout(m)
par(mai=c(0.35,0.35,0.35,0.1), omi=c(0,0,0,0), cex=1, mgp=c(1.5,0.5,0))

eplots(LSJ1out)
thetaplots(LSJ1out)

Fig. S2 Plot

Eandtheta = rbind(
  c(1,0),
  c(3,1),
  c(5,2),
  c(7,3)
)

# working directory
setwd("")  # enter appropriate working directory for N2 WT foraging worms

# load eigenworms [Broekmans et al. 2016]
ew = read.csv("eigenworms.txt", header=F, sep="")  # adjust as appropriate

# filenames
filenames = c(
  '1.txt',
  '2.txt',
  '3.txt',
  '4.txt',
  '5.txt',
  '6.txt',
  '7.txt',
  '8.txt',
  '9.txt',
  '10.txt',
  '11.txt',
  '12.txt'
)

datasections = rbind(
  c(10001,12000), #worm1
  c(10001,12000), #worm2
  c(8501,10500),  #worm3
  c(10501,12500), #worm4
  c(10501,12500), #worm5
  c(10501,12500), #worm6
  c(10501,12500), #worm7
  c(10501,12500), #worm8
  c(11001,13000), #worm9
  c(11001,13000), #worm10
  c(11001,13000), #worm11
  c(11001,13000)  #worm12
)


# calculate the errors for first 3 worms
worm = 1
data = read.csv(filenames[worm], header=F, sep="")
data = data[datasections[worm,1]:datasections[worm,2],1:5]
w1_detailsens = run_sensitivity_simple(data, lib, pred, Eandtheta, tau, tp)

worm = 2
data = read.csv(filenames[worm], header=F, sep="")
data = data[datasections[worm,1]:datasections[worm,2],1:5]
w2_detailsens = run_sensitivity_simple(data, lib, pred, Eandtheta, tau, tp)

worm = 3
data = read.csv(filenames[worm], header=F, sep="")
data = data[datasections[worm,1]:datasections[worm,2],1:5]
w3_detailsens = run_sensitivity_simple(data, lib, pred, Eandtheta, tau, tp)


# plot supp fig
m = matrix(c(1,2,3), nrow=3, ncol=1, byrow=TRUE)
layout(m)
par(mai=c(0.3,0.3,0.35,0.1), omi=c(0,0,0,0), cex=1, mgp=c(1.5,0.5,0))
cols = c(rgb(1,0,0,0.7), 'orange', rgb(0,0,0.5,1))

# first plot
y = w1_detailsens[1,]
plot(y, type='o', col='gray', pch='', lty=1, lwd=1, xlab='', ylab='')

for (i in c(2:4)) {
  y = w1_detailsens[i,]
  lines(y, col=cols[i-1], lty=1, lwd=1)
}

# second plot
y = w2_detailsens[1,]
plot(y, type='o', col='gray', pch='', lty=1, lwd=1, xlab='', ylab='')

for (i in c(2:4)) {
  y = w2_detailsens[i,]
  lines(y, col=cols[i-1], lty=1, lwd=1)
}

# third plot
y = w3_detailsens[1,]
plot(y, type='o', col='gray', pch='', lty=1, lwd=1, xlab='', ylab='')

for (i in c(2:4)) {
  y = w3_detailsens[i,]
  lines(y, col=cols[i-1], lty=1, lwd=1)
}

Functions for Plotting

eplots <- function(indata) {
  xdata = Es  # same sequence of E values used to generate the data
  linecolor = rgb(0,0,0.7,0.5)
  linecolor2 = rgb(0,0,0.7,1)
  symbol = '.'
  smallwidth = 0.25
  bigwidth = 3
  thetaindex = 5
  
  y = indata[,thetaindex,1]
  plot(xdata, y, type='o', col=linecolor, pch=symbol, lty=1, lwd=smallwidth, xlab='',
       ylab='', ylim=c(min(indata[,thetaindex,]),max(indata[,thetaindex,])))
  
  for (i in c(2:dim(indata)[3])) {
    y = indata[,thetaindex,i]
    points(xdata, y, col=linecolor, pch=symbol)
    lines(xdata, y, col=linecolor, lty=1, lwd=smallwidth)
  }
  
  y = rowMeans(indata[,thetaindex,])
  lines(xdata, y, col=linecolor2, lty=1, lwd=bigwidth)
}


thetaplots <- function(indata){
  xdata = thetas  # same sequence of thetas used to generate the data
  linecolor = rgb(0.7,0,0,0.5)
  linecolor2 = rgb(0.7,0,0,1)
  symbol = '.'
  smallwidth = 0.25
  bigwidth = 3
  eindex = 5
  
  y = indata[eindex,,1]
  plot(xdata, y, type='o', col=linecolor, pch=symbol, lty=1, lwd=smallwidth, xlab='',
       ylim=c(min(indata[eindex,,]),max(indata[eindex,,])), ylab='')
  
  
  for (i in c(2:dim(indata)[3])) {
    y = indata[eindex,,i]
    points(xdata, y, col=linecolor, pch=symbol)
    lines(xdata, y, col=linecolor, lty=1, lwd=smallwidth)
  }
  
  y = rowMeans(indata[eindex,,])
  lines(xdata, y, col=linecolor2, lty=1, lwd=bigwidth)
}