Thalamus Specific c-Fos Response

Introduction

The purpose of this markdown is to walk though our process for calculating a c-Fos score at the regional level, and provide the code supporting the visualizations we provide in our manuscript https://docs.google.com/document/d/166X4o_6HegeS0uUHRa4ahKedkaVEMeiRagCm7T9nOVU/edit?usp=sharing. Specifically, here we demonstrate how to generate Figures 7k-r. First, we load relevant packages.

library(nat)
library(plyr)
library(cocoframer)
library(purrr)
library(rgl)
library(RNifti)
library(tidyverse)
library(ggplot2)
library(ggrepel)
library(RColorBrewer)
library(tidyr)
library(plotly)
library(dplyr)
library(reshape2)
library(data.table)
'%notin%' = Negate('%in%')

Anatomy of the Thalamus

Our first step is to understand where the thalamus resides within the brain, and where the substructures of interest to us are located. First, we load in the Allen Brain Atlas mesh, and translate the whole brain mesh slightly, so that the Slideseq data and the ABA atlas are within the same reference space.

#Load the meshes
structures <- c("root", "TH", "PVT", "RE", "Xi")
mesh_list <- map(structures, ccf_2017_mesh)
names(mesh_list) <- structures
#Translate the ABA to fit the Slideseq space
for(x in structures){
  mesh_list[[x]]$vb[2,] = mesh_list[[x]]$vb[2,] - 500
  mesh_list[[x]]$vb[3,] = mesh_list[[x]]$vb[3,] - 500

}

Plot the whole brain, as well as the structure of the thalamus within the brain.

fig = plot_ly(type = 'mesh3d',  x=mesh_list$root$vb[1,], y=mesh_list$root$vb[2,],z=mesh_list$root$vb[3,],i=mesh_list$root$it[1,]-1,j=mesh_list$root$it[2,]-1,k =mesh_list$root$it[3,]-1, facecolor = rep("lightgray", length(mesh_list$root$vb[2,])*2), opacity = 0.1) %>% add_trace(type = 'mesh3d', x=mesh_list$TH$vb[1,], 
 y=mesh_list$TH$vb[2,],z=mesh_list$TH$vb[3,], i=mesh_list$TH$it[1,]-1,j=mesh_list$TH$it[2,]-1, k =mesh_list$TH$it[3,]-1, facecolor = rep("darkslateblue",length(mesh_list$TH$vb[2,])*2 ),opacity = 1)

fig <- fig %>%  layout(title = "Thalamus",scene = list(xaxis = list(title = '', showgrid = F , zerolinecolor = '#ffff', showticklabels=FALSE),yaxis = list(title = '', showgrid = F,zerolinecolor = '#ffff', showticklabels=FALSE),zaxis = list(title = '', showgrid = F, zerolinecolor = '#ffff', showticklabels=FALSE)))
fig

Now, let’s look at some of the structures within the thalamus

fig = plot_ly(type = 'mesh3d', x=mesh_list$TH$vb[1,], y=mesh_list$TH$vb[2,],z=mesh_list$TH$vb[3,], i=mesh_list$TH$it[1,]-1, j=mesh_list$TH$it[2,]-1, k =mesh_list$TH$it[3,]-1, facecolor = rep("darkslateblue",length(mesh_list$TH$vb[2,])*2 ), opacity = 0.1, showlegend = T) %>% add_trace(type = 'mesh3d', x=mesh_list$PVT$vb[1,], y=mesh_list$PVT$vb[2,],z=mesh_list$PVT$vb[3,], i=mesh_list$PVT$it[1,]-1,j=mesh_list$PVT$it[2,]-1, k =mesh_list$PVT$it[3,]-1, facecolor = rep("lightpink",length(mesh_list$PVT$vb[2,])*2 ),opacity = 1, showlegend = T) %>% add_trace(type = 'mesh3d', x=mesh_list$RE$vb[1,],  y=mesh_list$RE$vb[2,],z=mesh_list$RE$vb[3,],  i=mesh_list$RE$it[1,]-1,j=mesh_list$RE$it[2,]-1, k =mesh_list$RE$it[3,]-1, facecolor = rep("darkseagreen",length(mesh_list$RE$vb[2,])*2 ),opacity = 1, showlegend = T)  %>% add_trace(type = 'mesh3d', x=mesh_list$Xi$vb[1,], y=mesh_list$Xi$vb[2,], z=mesh_list$Xi$vb[3,], i=mesh_list$Xi$it[1,]-1, j=mesh_list$Xi$it[2,]-1, k =mesh_list$Xi$it[3,]-1, facecolor = rep("orange",length(mesh_list$Xi$vb[2,])*2 ),opacity = 1, showlegend = T) 
fig <- fig %>%  layout(title = "Thalamus",
                       scene = list(xaxis = list(title = '', showgrid = F , zerolinecolor = '#ffff', showticklabels=FALSE),
                                    yaxis = list(title = '', showgrid = F,zerolinecolor = '#ffff', showticklabels=FALSE),
                                    zaxis = list(title = '', showgrid = F, zerolinecolor = '#ffff', showticklabels=FALSE)))
fig

Visualization of the Thalamus c-Fos and Slideseq Data

Now that we know where the thalamus and its sub anatomical structures lie within the brain, we now want to visualize the Slideseq beads in conjunction with the thalamus structures.

#Read in both the Slideseq and c-Fos data
slideseq = readRDS("Data/Vascular/SupportFiles/100um_Celltypes_in_each_Voxel_and_totalsums.rds")
cfos_data = readRDS("Data/cFos/DataFiles/CFOS_data_registeredto_100UmABA_OrganizedbyTimepoint.rds")
#Filter for Slide-seq voxels that confidentally contain atleast one voxel
slideseq = filter(slideseq, slideseq$total_sums >= 1)
joint = left_join(slideseq, cfos_data)
#Reduce to only voxels from the Thalamus
voxelsTH = readRDS("Data/cFos/SupportFiles/Voxels_TH.rds")
joint = filter(joint, joint$Voxel_ID %in% voxelsTH)

Now we can plot the Slide-seq thalamic beads on top of the ABA thalamus mesh to insure correspondence.

ssFig = plot_ly(type = "scatter3d", mode = "markers",data = joint,x = ~Voxel_X,y = ~Voxel_Y,z = ~Voxel_Z, marker = list(size = 3, color = 'black', opacity = 1))  %>% add_trace(type = 'mesh3d', x=mesh_list$TH$vb[1,],y=mesh_list$TH$vb[2,],z=mesh_list$TH$vb[3,], i=mesh_list$TH$it[1,]-1,j=mesh_list$TH$it[2,]-1,k =mesh_list$TH$it[3,]-1,facecolor = rep("darkslateblue",length(mesh_list$TH$vb[2,])*2 ),opacity = 0.05) %>% add_trace(type = 'mesh3d', x=mesh_list$PVT$vb[1,],
 y=mesh_list$PVT$vb[2,],z=mesh_list$PVT$vb[3,],i=mesh_list$PVT$it[1,]-1, j=mesh_list$PVT$it[2,]-1, k =mesh_list$PVT$it[3,]-1,facecolor = rep("lightpink",length(mesh_list$PVT$vb[2,])*2 ),opacity = 0.3) %>% add_trace(type = 'mesh3d', x=mesh_list$RE$vb[1,],y=mesh_list$RE$vb[2,],z=mesh_list$RE$vb[3,],i=mesh_list$RE$it[1,]-1,j=mesh_list$RE$it[2,]-1,k=mesh_list$RE$it[3,]-1,facecolor =rep("darkseagreen",length(mesh_list$RE$vb[2,])*2 ), opacity = 0.3)  %>% add_trace(type = 'mesh3d', x=mesh_list$Xi$vb[1,],y=mesh_list$Xi$vb[2,],z=mesh_list$Xi$vb[3,],  i=mesh_list$Xi$it[1,]-1,j=mesh_list$Xi$it[2,]-1,k =mesh_list$Xi$it[3,]-1, facecolor = rep("orange",length(mesh_list$Xi$vb[2,])*2 ), opacity = 0.3)
ssFig <- ssFig %>%  layout(title = "Thalamus Voxels", scene = list(xaxis = list(title = '', showgrid = F , zerolinecolor = '#ffff', showticklabels=FALSE),yaxis = list(title = '', showgrid = F,zerolinecolor = '#ffff', showticklabels=FALSE),zaxis = list(title = '', showgrid = F, zerolinecolor = '#ffff', showticklabels=FALSE)))
ssFig

Examining the Number of c-Fos Expressing Cells within the Thalamus

Before calculating a c-Fos score for the molecular types present within the thalamus, we want to examine the voxelized c-Fos data that is present within the thalamus. Beginning with Time point 30 minutes after foot shock, we plot the voxels of the Thalamus that overlap with the Slide-seq data. Each dot is a 100 um voxel, and the color of each dot reflects the number of c-Fos expressing cells within that voxel. For T30:

#Filter out voxels with 0 expressing c-Fos cells for visability.
cfosT30 = filter(joint, joint$Change_T30 > 0)
cfosT30$Colors = as.numeric(cut(cfosT30$Change_T30,breaks=50))
voxelcolors = grDevices::palette(viridis::viridis(option="A",n=50,direction = -1))
t30 = plot3d(cfosT30$Voxel_X,cfosT30$Voxel_Y,cfosT30$Voxel_Z,col = cfosT30$Colors,xlab="Anterior-Posterior",ylab="Inferior-Superior",zlab="Left-Right",size=5, add = TRUE)
t30
view3d(theta = -2, phi = 8)

For T90:

close3d()
cfosT90 = filter(joint, joint$Change_T90 > 0)
cfosT90$Colors = as.numeric(cut(cfosT90$Change_T90,breaks=50))
voxelcolors = grDevices::palette(viridis::viridis(option="A",n=50,direction = -1))
t90 = plot3d(cfosT90$Voxel_X,cfosT90$Voxel_Y,cfosT90$Voxel_Z,col = cfosT90$Colors,xlab="Anterior-Posterior",ylab="Inferior-Superior",zlab="Left-Right",size=5, add = TRUE)
t90
view3d(theta = -2, phi = 8)

For T180:

close3d()
#Run for each timepoint
cfosT180 = filter(joint, joint$Change_T180 > 0)
cfosT180$Colors = as.numeric(cut(cfosT180$Change_T180,breaks=50))
voxelcolors = grDevices::palette(viridis::viridis(option="A",n=50,direction = -1))
t180 = plot3d(cfosT180$Voxel_X,cfosT180$Voxel_Y,cfosT180$Voxel_Z,col = cfosT180$Colors,xlab="Anterior-Posterior",ylab="Inferior-Superior",zlab="Left-Right",size=5, add = TRUE)
t180
view3d(theta = -2, phi = 8)

Calculating a Thalamus Specific Score

Next we calculate a thalamus specifc c-Fos score for the molecular cell types within the thalamus.

close3d()
#For all instances where we have a Slideseq voxel, but do not have a c-Fos score, we set these values to zero
joint[is.na(joint)] <- 0
time = c("Change_T180", "Change_T90", "Change_T30", "Timepoint_0")
for(point in time){
  print(point)
  running_list = c()
  celltypes = colnames(joint)[1:1937]
  for (type in celltypes){
    ratio = joint[,type]/joint[, "total_sums"]
    cfos_count = joint[,point]
    products = ratio*cfos_count
    cfos_score = sum(products)/sum(joint[,type] != 0)
    newAdd =  data.frame(cfos_score, type, sum(joint[,type] != 0))
    colnames(newAdd) = c("CFOSScore", "CellType", "VoxelsOccupied")
    if(length(running_list) == 0){
      running_list = newAdd
    } else{
      running_list = rbind(running_list, newAdd)
    }
  }
  saveRDS(running_list, paste0("/nfs/turbo/umms-welchjd/BRAIN_initiative/Final_integration_workflow/Markdowns/Data/cFos/Results/Thalamus_specific_CFOS_expression_score_", point, ".RDS"))
}

Onset of Experimental Condition within the Thalamus

After calculating a thalamus-specific c-Fos score, we want to ensure that the score shows molecular cell type specificity, and reflects the onset of the experimental condition (footshock).

#Load in C-Fos activation scores for each timepoint
cfos_t30 = readRDS("Data/cFos/Results/Thalamus_specific_CFOS_expression_score_Change_T30.RDS")
colnames(cfos_t30)[1] = "CFOSScore_T30"

cfos_t90 = readRDS("Data/cFos/Results/Thalamus_specific_CFOS_expression_score_Change_T90.RDS")
colnames(cfos_t90)[1] = "CFOSScore_T90"

cfos_t180 = readRDS("Data/cFos/Results/Thalamus_specific_CFOS_expression_score_Change_T180.RDS")
colnames(cfos_t180)[1] = "CFOSScore_T180"

cfos_t0 = readRDS("Data/cFos/Results/Thalamus_specific_CFOS_expression_score_Timepoint_0.RDS")
colnames(cfos_t0)[1] = "CFOSScore_T0"

allFOS = left_join(cfos_t0, cfos_t30)
allFOS = left_join(allFOS, cfos_t90)
allFOS = left_join(allFOS, cfos_t180) 

#Filter for molecular cell types occupying greater than 10 voxels within the thalamus
allFOS = filter(allFOS, allFOS$VoxelsOccupied > 10)
#Remove non-neuronal cell typrs
allFOS = separate(allFOS, CellType, into = c("Broad"), remove = FALSE, sep = "_")
nonneuronal = c("Pit" , "Astro" , "Ependymal" ,"Tanycyte" , "Oligo", "OPC", "Endo" ,"Fibro" ,"Micro","Macro") 
allFOS = filter(allFOS, allFOS$Broad %notin% nonneuronal) 

#We want to examine the specificity of the c-Fos response at the regional level. In order to do this, we label each molecular cell type with the thalamuc region that most of its beads fall into. If the cell type has a significant presence in multiple thalamic regions, we label it "TH", for thalamus. 
annies = readRDS("/nfs/turbo/umms-welchjd/BRAIN_initiative/Final_integration_workflow/Markdowns/Data/cFos/SupportFiles/TH_celltype_Location_annotations.rds")
allFOS = left_join(allFOS, annies)

#Prepare the data for plotting
allFOS = select(allFOS, -c("VoxelsOccupied", "Broad"))
allFOS <- melt(setDT(allFOS), id.vars = c("CellType","DeepCCF"), variable.name = "TimePoint")
allFOS$TimePoint = factor(sub("CFOSScore_", "", allFOS$TimePoint), levels=c("T0", "T30", "T90", "T180"))
allFOS = allFOS %>% mutate(label = if_else(TimePoint == "T180" & value > 4.5, as.character(CellType), NA_character_))


p = allFOS %>%
  ggplot(aes(x=TimePoint, y=value, group=CellType, colour=DeepCCF)) +
  geom_line(linewidth = 2) + theme_bw() + xlab("Time") + ylab("cFOS Score") + ggtitle("Thalamic c-Fos Response Over Time")  +
  geom_label_repel(aes(label = label), nudge_x = 1, na.rm = TRUE,max.overlaps = 25, show.legend = FALSE) + ylim(-2, 7)
p

Spatial Distribution of High Scoring Cell Types

After examining the c-Fos score of cell types across time within the thalamus (see above plot), we identify that there does seem to be some regional specificity in the molecular cell types with the highest c-Fos scores. We will examine the spatial distribution of the top five highest scoring molecular cell types at each experimental timepoint. Note: The code in the following code chunk consumes a considerable amount of memory. Therefore, while we provide the code, we recommend skipping the following segment of code, and resuming with the pre-saved object in the next code segment when running on a personal laptop.

beads = readRDS("/Data/cFos/SupportFiles/All_Beads_CellTypeProportions_Filtered.rds")
beads = as.matrix(beads) %>% data.frame
meta = readRDS("Data/cFos/SupportFiles/Slideseq_beads_CCF_Registered_TopStruct.rds")
#Reduce to beads in the thalamus
meta = filter(meta, meta$TopStruct == "TH")
beads = beads[meta$BeadBarcode,]
#Set only for cell types confidentally assigned
beads[beads <= .3] <- 0
beads[beads > .3] <- 1
saveRDS(beads, "Data/cFos/SupportFiles/SlideseqBeads_Thalamus_confidentCelltypes.rds")

Resume here with the pre-saved object

meta = readRDS("Data/cFos/SupportFiles/Slideseq_beads_CCF_Registered_TopStruct.rds")
#Read in beads
beads = readRDS("Data/cFos/SupportFiles/SlideseqBeads_Thalamus_confidentCelltypes.rds")

cfos_t30 = readRDS("Data/cFos/Results/Thalamus_specific_CFOS_expression_score_Change_T30.RDS")
cfos_t30 = filter(cfos_t30, cfos_t30$VoxelsOccupied > 10)
high_t30 = cfos_t30[order(cfos_t30$CFOSScore, decreasing = TRUE),][1:5,]
cellTypes = unique(high_t30$CellType)
#Select only our top scoring cell types
beads30 = beads[, cellTypes]
beads30$BeadBarcode = rownames(beads30)
beads30 <- melt(setDT(beads30), id.vars = c("BeadBarcode"), variable.name = "CellType")
#Filter out all beads that are not confidentally assigned one of our cell types of interest
beads30 = filter(beads30, beads30$value != 0)
beads30 = left_join(beads30, meta)

## Joining with `by = join_by(BeadBarcode)`

palette = c("deepskyblue",  "darkgreen", "springgreen", "plum", "darkorange")
ssFig = plot_ly(type = "scatter3d",mode = "markers",data =beads30,x = ~CCF_X, y = ~CCF_Y,z = ~CCF_Z, color = ~CellType,colors = palette, size =3) %>% add_trace(type = 'mesh3d', x=mesh_list$TH$vb[1,], y=mesh_list$TH$vb[2,],z=mesh_list$TH$vb[3,],  i=mesh_list$TH$it[1,]-1, j=mesh_list$TH$it[2,]-1,  k =mesh_list$TH$it[3,]-1, facecolor = rep("darkslateblue",length(mesh_list$TH$vb[2,])*2 ),opacity = 0.01) %>% add_trace(type = 'mesh3d', x=mesh_list$PVT$vb[1,],  y=mesh_list$PVT$vb[2,], z=mesh_list$PVT$vb[3,],i=mesh_list$PVT$it[1,]-1, j=mesh_list$PVT$it[2,]-1,k =mesh_list$PVT$it[3,]-1, facecolor = rep("lightpink",length(mesh_list$PVT$vb[2,])*2 ),opacity = 0.02) %>% add_trace(type = 'mesh3d', x=mesh_list$RE$vb[1,], y=mesh_list$RE$vb[2,],z=mesh_list$RE$vb[3,],  i=mesh_list$RE$it[1,]-1, j=mesh_list$RE$it[2,]-1, k =mesh_list$RE$it[3,]-1, facecolor = rep("darkseagreen",length(mesh_list$RE$vb[2,])*2 ),opacity = 0.02)  %>% add_trace(type = 'mesh3d', x=mesh_list$Xi$vb[1,], y=mesh_list$Xi$vb[2,],z=mesh_list$Xi$vb[3,], i=mesh_list$Xi$it[1,]-1,j=mesh_list$Xi$it[2,]-1,                                                       k=mesh_list$Xi$it[3,]-1,  facecolor = rep("orange",length(mesh_list$Xi$vb[2,])*2 ), opacity = 0.02) 
ssFig <- ssFig %>%  layout(title = "Timepoint30", scene = list(xaxis = list(title = '', showgrid = F , zerolinecolor = '#ffff', showticklabels=FALSE),yaxis = list(title = '', showgrid = F,zerolinecolor = '#ffff', showticklabels=FALSE), zaxis = list(title = '', showgrid = F, zerolinecolor = '#ffff', showticklabels=FALSE)))
ssFig

cfos_t90 = readRDS("Data/cFos/Results/Thalamus_specific_CFOS_expression_score_Change_T90.RDS")
cfos_t90 = filter(cfos_t90, cfos_t90$VoxelsOccupied > 10)
high_t90 = cfos_t90[order(cfos_t90$CFOSScore, decreasing = TRUE),][1:5,]
cellTypes = unique(high_t90$CellType)
#Select only our top scoring cell types
beads90 = beads[, cellTypes]
beads90$BeadBarcode = rownames(beads90)
beads90 <- melt(setDT(beads90), id.vars = c("BeadBarcode"), variable.name = "CellType")
#Filter out all beads that are not confidentally assigned one of our cell types of interest
beads90 = filter(beads90, beads90$value != 0)
beads90 = left_join(beads90, meta)

## Joining with `by = join_by(BeadBarcode)`

palette = c("darkorchid", "gold", "darkseagreen", "deeppink", "firebrick")
ssFig = plot_ly(type = "scatter3d",mode = "markers",data = beads90,x = ~CCF_X, y = ~CCF_Y,z = ~CCF_Z, color = ~CellType,colors = palette, size =3)  %>% add_trace(type = 'mesh3d', x=mesh_list$TH$vb[1,], y=mesh_list$TH$vb[2,],z=mesh_list$TH$vb[3,],  i=mesh_list$TH$it[1,]-1, j=mesh_list$TH$it[2,]-1,  k =mesh_list$TH$it[3,]-1, facecolor = rep("darkslateblue",length(mesh_list$TH$vb[2,])*2 ),opacity = 0.01) %>% add_trace(type = 'mesh3d', x=mesh_list$PVT$vb[1,],  y=mesh_list$PVT$vb[2,], z=mesh_list$PVT$vb[3,],i=mesh_list$PVT$it[1,]-1, j=mesh_list$PVT$it[2,]-1,k =mesh_list$PVT$it[3,]-1, facecolor = rep("lightpink",length(mesh_list$PVT$vb[2,])*2 ),opacity = 0.02) %>% add_trace(type = 'mesh3d', x=mesh_list$RE$vb[1,], y=mesh_list$RE$vb[2,],z=mesh_list$RE$vb[3,],  i=mesh_list$RE$it[1,]-1, j=mesh_list$RE$it[2,]-1, k =mesh_list$RE$it[3,]-1, facecolor = rep("darkseagreen",length(mesh_list$RE$vb[2,])*2 ),opacity = 0.02)  %>% add_trace(type = 'mesh3d', x=mesh_list$Xi$vb[1,], y=mesh_list$Xi$vb[2,],z=mesh_list$Xi$vb[3,], i=mesh_list$Xi$it[1,]-1,j=mesh_list$Xi$it[2,]-1,                                                       k=mesh_list$Xi$it[3,]-1,  facecolor = rep("orange",length(mesh_list$Xi$vb[2,])*2 ), opacity = 0.02) 
ssFig <- ssFig %>%  layout(title = "Timepoint90", scene = list(xaxis = list(title = '', showgrid = F , zerolinecolor = '#ffff', showticklabels=FALSE),yaxis = list(title = '', showgrid = F,zerolinecolor = '#ffff', showticklabels=FALSE), zaxis = list(title = '', showgrid = F, zerolinecolor = '#ffff', showticklabels=FALSE)))
ssFig

cfos_t180 = readRDS("Data/cFos/Results/Thalamus_specific_CFOS_expression_score_Change_T180.RDS")
cfos_t180 = filter(cfos_t180, cfos_t180$VoxelsOccupied > 10)
high_t180 = cfos_t180[order(cfos_t180$CFOSScore, decreasing = TRUE),][1:5,]
cellTypes = unique(high_t180$CellType)
#Select only our top scoring cell types
beads180 = beads[, cellTypes]
beads180$BeadBarcode = rownames(beads180)
beads180 <- melt(setDT(beads180), id.vars = c("BeadBarcode"), variable.name = "CellType")
#Filter out all beads that are not confidentally assigned one of our cell types of interest
beads180 = filter(beads180, beads180$value != 0)
beads180 = left_join(beads180, meta)

## Joining with `by = join_by(BeadBarcode)`

palette = c("darkorchid", "deepskyblue", "gold", "firebrick", "mediumblue")
ssFig = plot_ly(type = "scatter3d", mode = "markers",data = beads180, x = ~CCF_X, y = ~CCF_Y, z = ~CCF_Z, color = ~CellType,colors = palette, size =3) %>% add_trace(type = 'mesh3d', x=mesh_list$TH$vb[1,], y=mesh_list$TH$vb[2,],z=mesh_list$TH$vb[3,],  i=mesh_list$TH$it[1,]-1, j=mesh_list$TH$it[2,]-1,  k =mesh_list$TH$it[3,]-1, facecolor = rep("darkslateblue",length(mesh_list$TH$vb[2,])*2 ),opacity = 0.01) %>% add_trace(type = 'mesh3d', x=mesh_list$PVT$vb[1,],  y=mesh_list$PVT$vb[2,], z=mesh_list$PVT$vb[3,],i=mesh_list$PVT$it[1,]-1, j=mesh_list$PVT$it[2,]-1,k =mesh_list$PVT$it[3,]-1, facecolor = rep("lightpink",length(mesh_list$PVT$vb[2,])*2 ),opacity = 0.02) %>% add_trace(type = 'mesh3d', x=mesh_list$RE$vb[1,], y=mesh_list$RE$vb[2,],z=mesh_list$RE$vb[3,],  i=mesh_list$RE$it[1,]-1, j=mesh_list$RE$it[2,]-1, k =mesh_list$RE$it[3,]-1, facecolor = rep("darkseagreen",length(mesh_list$RE$vb[2,])*2 ),opacity = 0.02)  %>% add_trace(type = 'mesh3d', x=mesh_list$Xi$vb[1,], y=mesh_list$Xi$vb[2,],z=mesh_list$Xi$vb[3,], i=mesh_list$Xi$it[1,]-1,j=mesh_list$Xi$it[2,]-1,                                                       k=mesh_list$Xi$it[3,]-1,  facecolor = rep("orange",length(mesh_list$Xi$vb[2,])*2 ), opacity = 0.02) 
ssFig <- ssFig %>%  layout(title = "Timepoint180", scene = list(xaxis = list(title = '', showgrid = F , zerolinecolor = '#ffff', showticklabels=FALSE),yaxis = list(title = '', showgrid = F,zerolinecolor = '#ffff', showticklabels=FALSE), zaxis = list(title = '', showgrid = F, zerolinecolor = '#ffff', showticklabels=FALSE)))
ssFig