SOCR ≫ TCIU Website ≫ TCIU GitHub ≫

1 Introduction

This TCIU Appendix implements the Kime-Phase Tomography (KPT) algorithm to map repeated fMRI time-series measurements to 2D kime-surfaces. We test the algorithm on synthetic fMRI data with ON (stimulus) and OFF (rest) conditions.

2 Required Libraries

# Load required libraries
library(plotly)
library(dplyr)
library(tidyr)
library(circular)
library(MASS)
library(pracma)  # For numerical integration
library(ggplot2)
library(viridis)

# Set random seed for reproducibility
set.seed(42)

3 Simulate fMRI Data

# ---------------------------- 
# Simulate fMRI Data
# ---------------------------- 
simulate_fmri_data <- function(n_participants = 15, n_runs = 20, tr = 2, duration = 300) {
  # Hemodynamic Response Function (HRF)
  hrf <- function(t) {
    a1 <- 6; a2 <- 12; b1 <- 0.9; b2 <- 0.9; c <- 0.35
    t^(a1-1)*b1^a1*exp(-b1*t)/gamma(a1) - c*(t^(a2-1)*b2^a2*exp(-b2*t)/gamma(a2))
  }
  
  # Experimental design parameters
  n_timepoints <- duration/tr
  block_length <- 30  # 30s blocks
  n_blocks <- duration/block_length
  
  # Storage arrays
  data <- list(
    ON = array(0, dim = c(n_participants, n_runs, n_timepoints)),
    OFF = array(0, dim = c(n_participants, n_runs, n_timepoints))
  )
  
  # Generate data for each participant and run
  set.seed(123)
  for(p in 1:n_participants) {
    for(r in 1:n_runs) {
      # Generate block design
      design <- rep(c(0,1), each = block_length/(2*tr))
      if(length(design) < n_timepoints) {
        design <- c(design, rep(design[length(design)], n_timepoints - length(design)))
      } else {
        design <- design[1:n_timepoints]
      }
      
      # Generate HRF response
      hrf_values <- hrf(seq(0, block_length/tr - 1, length.out = block_length/tr))
      hrf_response <- numeric(n_timepoints)
      
      # Convolve design with HRF
      for(i in 1:n_timepoints) {
        if(design[i] == 1) {
          end_idx <- min(i + length(hrf_values) - 1, n_timepoints)
          hrf_response[i:end_idx] <- hrf_response[i:end_idx] + 
            hrf_values[1:(end_idx - i + 1)]
        }
      }
      
      # Add participant-specific noise and drift
      noise <- cumsum(rnorm(n_timepoints, sd = 0.1))
      drift <- 0.01 * seq(n_timepoints)
      
      # Store responses
      data$ON[p, r, ] <- hrf_response + noise + drift + rnorm(n_timepoints, sd = 0.2)
      data$OFF[p, r, ] <- noise + drift + rnorm(n_timepoints, sd = 0.2)
    }
  }
  return(data)
}

# Generate synthetic fMRI data
fmri_data <- simulate_fmri_data(n_participants = 10, n_runs = 20, tr = 2, duration = 300)

3.1 fMRI Simulated Data Visualization

# ---------------------------- 
# Visualize Original fMRI Time-Series
# ---------------------------- 

# Function to create interactive time-series plot
create_timeseries_plot <- function(fmri_data, n_display = 5) {
  # Extract dimensions
  n_participants <- dim(fmri_data$ON)[1]
  n_runs <- dim(fmri_data$ON)[2]
  n_timepoints <- dim(fmri_data$ON)[3]
  time_vec <- seq(0, (n_timepoints-1)*2, by = 2)  # TR = 2
  
  # Randomly select indices to display
  set.seed(42)
  selected_indices <- sample(1:(n_participants * n_runs), n_display)
  
  # Create plot
  fig <- plot_ly()
  
  # Add ON condition traces
  for(i in 1:n_display) {
    idx <- selected_indices[i]
    participant <- ((idx - 1) %/% n_runs) + 1
    run <- ((idx - 1) %% n_runs) + 1
    
    fig <- fig %>% add_trace(
      x = time_vec,
      y = fmri_data$ON[participant, run, ],
      type = 'scatter',
      mode = 'lines',
      name = paste0("ON - P", participant, "R", run),
      line = list(color = 'red', width = 2),
      opacity = 0.7,
      visible = (i <= 3)  # Show only first 3 by default
    )
  }
  
  # Add OFF condition traces
  for(i in 1:n_display) {
    idx <- selected_indices[i]
    participant <- ((idx - 1) %/% n_runs) + 1
    run <- ((idx - 1) %% n_runs) + 1
    
    fig <- fig %>% add_trace(
      x = time_vec,
      y = fmri_data$OFF[participant, run, ],
      type = 'scatter',
      mode = 'lines',
      name = paste0("OFF - P", participant, "R", run),
      line = list(color = 'blue', width = 2),
      opacity = 0.7,
      visible = (i <= 3)  # Show only first 3 by default
    )
  }
  
  # Create buttons for navigation
  button_list <- list()
  
  # Add "Show All" button
  visible_all <- rep(TRUE, 2 * n_display)
  button_list[[1]] <- list(
    method = "restyle",
    args = list("visible", as.list(visible_all)),
    label = "Show All"
  )
  
  # Add individual selection buttons
  for(i in 1:n_display) {
    visible_vec <- rep(FALSE, 2 * n_display)
    visible_vec[c(i, i + n_display)] <- TRUE  # Show ON and OFF for this index
    
    idx <- selected_indices[i]
    participant <- ((idx - 1) %/% n_runs) + 1
    run <- ((idx - 1) %% n_runs) + 1
    
    button_list[[i + 1]] <- list(
      method = "restyle",
      args = list("visible", as.list(visible_vec)),
      label = paste0("P", participant, "R", run)
    )
  }
  
  # Update layout with buttons
  fig <- fig %>% layout(
    title = "fMRI Time-Series: ON (Red) vs OFF (Blue) Conditions",
    xaxis = list(title = "Time (seconds)"),
    yaxis = list(title = "BOLD Signal"),
    updatemenus = list(
      list(
        type = "buttons",
        direction = "left",
        xanchor = "left",
        yanchor = "top",
        x = 0.01,
        y = 1.15,
        buttons = button_list
      )
    ),
    annotations = list(
      list(
        text = "Select Time-Series:",
        showarrow = FALSE,
        x = 0,
        y = 1.2,
        xref = "paper",
        yref = "paper",
        xanchor = "right",
        yanchor = "bottom"
      )
    )
  )
  
  return(fig)
}

# Create and display the time-series plot
ts_plot <- create_timeseries_plot(fmri_data, n_display = 10)
ts_plot

3.2 Simulated fMRI Data Summary Statistics

# ---------------------------- 
# Time-Series Statistics
# ---------------------------- 

# Compute statistics for each condition
compute_ts_stats <- function(data_array) {
  n_series <- dim(data_array)[1] * dim(data_array)[2]
  n_timepoints <- dim(data_array)[3]
  
  # Reshape to matrix
  reshaped <- matrix(data_array, nrow = n_series, ncol = n_timepoints)
  
  # Compute statistics
  mean_ts <- colMeans(reshaped)
  sd_ts <- apply(reshaped, 2, sd)
  
  # Compute temporal SNR
  tsnr <- mean_ts / sd_ts
  
  return(list(
    mean = mean_ts,
    sd = sd_ts,
    tsnr = tsnr,
    global_mean = mean(mean_ts),
    global_sd = mean(sd_ts)
  ))
}

stats_on <- compute_ts_stats(fmri_data$ON)
stats_off <- compute_ts_stats(fmri_data$OFF)

# Create comparison plot
time_vec <- seq(0, (length(stats_on$mean)-1)*2, by = 2)

stats_plot <- plot_ly() %>%
  # Mean traces
  add_trace(x = time_vec, y = stats_on$mean, type = 'scatter', mode = 'lines',
            name = 'ON Mean', line = list(color = 'red', width = 3)) %>%
  add_trace(x = time_vec, y = stats_off$mean, type = 'scatter', mode = 'lines',
            name = 'OFF Mean', line = list(color = 'blue', width = 3)) %>%
  # Confidence bands
  add_trace(x = c(time_vec, rev(time_vec)), 
            y = c(stats_on$mean + 2*stats_on$sd, rev(stats_on$mean - 2*stats_on$sd)),
            type = 'scatter', mode = 'none', fill = 'toself',
            fillcolor = 'rgba(255,0,0,0.2)', name = 'ON ±2SD', showlegend = FALSE) %>%
  add_trace(x = c(time_vec, rev(time_vec)), 
            y = c(stats_off$mean + 2*stats_off$sd, rev(stats_off$mean - 2*stats_off$sd)),
            type = 'scatter', mode = 'none', fill = 'toself',
            fillcolor = 'rgba(0,0,255,0.2)', name = 'OFF ±2SD', showlegend = FALSE) %>%
  layout(
    title = "Average Time-Series with Confidence Bands",
    xaxis = list(title = "Time (seconds)"),
    yaxis = list(title = "BOLD Signal"),
    hovermode = 'x unified'
  )

stats_plot
# Print summary statistics
cat("\n=== Time-Series Summary Statistics ===\n")
## 
## === Time-Series Summary Statistics ===
cat(sprintf("ON Condition - Global Mean: %.3f (SD: %.3f)\n", 
            stats_on$global_mean, stats_on$global_sd))
## ON Condition - Global Mean: 1.481 (SD: 0.792)
cat(sprintf("OFF Condition - Global Mean: %.3f (SD: %.3f)\n", 
            stats_off$global_mean, stats_off$global_sd))
## OFF Condition - Global Mean: 0.783 (SD: 0.792)
cat(sprintf("Mean Difference (ON-OFF): %.3f\n", 
            stats_on$global_mean - stats_off$global_mean))
## Mean Difference (ON-OFF): 0.699
cat(sprintf("Average tSNR - ON: %.3f, OFF: %.3f\n", 
            mean(stats_on$tsnr), mean(stats_off$tsnr)))
## Average tSNR - ON: 1.796, OFF: 0.896

4 Implement KPT Core Functions

4.1 Phase Distribution Functions

# Wrapped Laplace distribution density
wrapped_laplace_density <- function(theta, mu, beta) {
  # Normalize to [-pi, pi]
  theta <- ((theta - mu + pi) %% (2*pi)) - pi + mu
  
  # Wrapped Laplace density
  if(beta <= 0) return(rep(0, length(theta)))
  
  dens <- beta / (2 * sinh(beta)) * exp(beta * cos(theta - mu))
  return(dens)
}

# Sample from wrapped Laplace distribution
sample_wrapped_laplace <- function(n, mu, beta) {
  # Sample from linear Laplace
  u <- runif(n)
  linear_samples <- ifelse(u < 0.5,
                          mu + log(2*u)/beta,
                          mu - log(2*(1-u))/beta)
  
  # Wrap to [-pi, pi]
  wrapped_samples <- ((linear_samples + pi) %% (2*pi)) - pi
  return(wrapped_samples)
}

# Phase distribution parameters as functions of time
mu_t <- function(t, t_max = 150) {
  # Location parameter varies sinusoidally with time
  pi * sin(2 * pi * t / t_max)
}

beta_t <- function(t, t_max = 150) {
  # Scale parameter varies with time
  2 + sin(4 * pi * t / t_max)
}

4.2 Harmonic Expansion and Moment Computation

# Compute empirical moments
compute_empirical_moments <- function(s_values, theta_values, n_max = 10) {
  n_samples <- length(s_values)
  moments <- complex(length = 2*n_max + 1)
  
  for(n in -n_max:n_max) {
    moments[n + n_max + 1] <- mean(s_values * exp(-1i * n * theta_values))
  }
  
  return(moments)
}

# Estimate harmonic coefficients f_n(t)
estimate_harmonic_coefficients <- function(time_series_collection, n_max = 10) {
  n_time_points <- dim(time_series_collection)[3]
  n_total_series <- dim(time_series_collection)[1] * dim(time_series_collection)[2]
  
  # Reshape to (n_total_series, n_time_points)
  reshaped_data <- matrix(time_series_collection, 
                         nrow = n_total_series, 
                         ncol = n_time_points)
  
  # Initialize coefficient matrix
  f_coeffs <- matrix(0, nrow = n_time_points, ncol = 2*n_max + 1)
  
  # For each time point, estimate mean signal
  for(t_idx in 1:n_time_points) {
    # Simple estimate: f_0(t) = mean signal, others = 0 for now
    f_coeffs[t_idx, n_max + 1] <- mean(reshaped_data[, t_idx])
    
    # Add some harmonic structure (for demonstration)
    for(n in 1:n_max) {
      f_coeffs[t_idx, n_max + 1 + n] <- 0.1 * f_coeffs[t_idx, n_max + 1] * 
        exp(-n/3) * cos(2*pi*n*t_idx/n_time_points)
      f_coeffs[t_idx, n_max + 1 - n] <- 0.1 * f_coeffs[t_idx, n_max + 1] * 
        exp(-n/3) * cos(2*pi*n*t_idx/n_time_points)
    }
  }
  
  return(f_coeffs)
}

4.3 KPT Reconstruction Algorithm

# KPT reconstruction function
kpt_reconstruct <- function(time_series_collection, n_max = 10, n_phase_samples = 100) {
  n_participants <- dim(time_series_collection)[1]
  n_runs <- dim(time_series_collection)[2]
  n_time_points <- dim(time_series_collection)[3]
  
  # Step 1: Estimate harmonic coefficients
  f_coeffs <- estimate_harmonic_coefficients(time_series_collection, n_max)
  
  # Step 2: For each time point, sample phases and compute moments
  reconstructed_phases <- list()
  empirical_moments_all <- matrix(0, nrow = n_time_points, ncol = 2*n_max + 1)
  
  for(t_idx in 1:n_time_points) {
    t_val <- t_idx * 2  # TR = 2
    
    # Get all measurements at this time point
    s_values <- as.vector(time_series_collection[, , t_idx])
    
    # Sample phases from wrapped Laplace
    mu <- mu_t(t_val)
    beta <- beta_t(t_val)
    theta_samples <- sample_wrapped_laplace(length(s_values), mu, beta)
    
    # Compute empirical moments
    moments <- compute_empirical_moments(s_values, theta_samples, n_max)
    empirical_moments_all[t_idx, ] <- moments
    
    # Store phase distribution parameters
    reconstructed_phases[[t_idx]] <- list(
      t = t_val,
      mu = mu,
      beta = beta,
      theta_samples = theta_samples,
      s_values = s_values
    )
  }
  
  return(list(
    f_coeffs = f_coeffs,
    empirical_moments = empirical_moments_all,
    phase_data = reconstructed_phases,
    n_max = n_max
  ))
}

4.4 Kime-Surface Generation

# Generate kime-surface from KPT reconstruction
generate_kime_surface <- function(kpt_result, n_grid = 50) {
  # Create polar grid
  r_max <- max(sapply(kpt_result$phase_data, function(x) x$t))
  r_grid <- seq(0, r_max, length.out = n_grid)
  theta_grid <- seq(-pi, pi, length.out = n_grid)
  
  # Initialize surface
  surface_z <- matrix(0, nrow = n_grid, ncol = n_grid)
  
  # Interpolate intensity values
  for(i in 1:n_grid) {
    for(j in 1:n_grid) {
      r <- r_grid[i]
      theta <- theta_grid[j]
      
      if(r > 0) {
        # Find nearest time point
        t_idx <- which.min(abs(sapply(kpt_result$phase_data, function(x) x$t) - r))
        
        # Compute expected intensity using harmonic expansion
        z_val <- Re(kpt_result$f_coeffs[t_idx, kpt_result$n_max + 1])  # f_0(t)
        
        # Add harmonic contributions
        for(n in 1:kpt_result$n_max) {
          z_val <- z_val + 2 * Re(kpt_result$f_coeffs[t_idx, kpt_result$n_max + 1 + n] * 
                                  exp(1i * n * theta))
        }
        
        surface_z[i, j] <- z_val
      }
    }
  }
  
  # Convert to Cartesian coordinates
  x_grid <- outer(r_grid, cos(theta_grid))
  y_grid <- outer(r_grid, sin(theta_grid))
  
  return(list(
    x = x_grid,
    y = y_grid,
    z = surface_z,
    r_grid = r_grid,
    theta_grid = theta_grid
  ))
}

5 Apply KPT to fMRI Data

# Apply KPT to ON and OFF conditions
cat("Applying KPT to ON condition...\n")
## Applying KPT to ON condition...
kpt_on <- kpt_reconstruct(fmri_data$ON, n_max = 8)

cat("Applying KPT to OFF condition...\n")
## Applying KPT to OFF condition...
kpt_off <- kpt_reconstruct(fmri_data$OFF, n_max = 8)

# Generate kime-surfaces
surface_on <- generate_kime_surface(kpt_on, n_grid = 60)
surface_off <- generate_kime_surface(kpt_off, n_grid = 60)

6 Visualize Results

6.1 Interactive 3D Kime-Surfaces

# Create interactive 3D plot for ON condition
fig_on <- plot_ly(
  x = surface_on$x,
  y = surface_on$y,
  z = surface_on$z,
  type = "surface",
  colorscale = list(
    c(0, "rgb(20,20,80)"),
    c(0.25, "rgb(20,80,140)"),
    c(0.5, "rgb(20,140,200)"),
    c(0.75, "rgb(80,200,120)"),
    c(1, "rgb(255,220,80)")
  ),
  contours = list(
    z = list(
      show = TRUE,
      usecolormap = TRUE,
      highlightcolor = "#ff0000",
      project = list(z = TRUE)
    )
  ),
  lighting = list(
    ambient = 0.4,
    diffuse = 0.6,
    specular = 0.2,
    roughness = 0.85
  )
) %>%
  layout(
    title = "fMRI ON Condition - Kime-Surface",
    scene = list(
      xaxis = list(title = "X (t·cos θ)"),
      yaxis = list(title = "Y (t·sin θ)"),
      zaxis = list(title = "Z (Intensity)"),
      camera = list(
        eye = list(x = 1.5, y = 1.5, z = 1.5)
      )
    )
  )

# Create interactive 3D plot for OFF condition
fig_off <- plot_ly(
  x = surface_off$x,
  y = surface_off$y,
  z = surface_off$z,
  type = "surface",
  colorscale = list(
    c(0, "rgb(40,20,60)"),
    c(0.25, "rgb(60,40,100)"),
    c(0.5, "rgb(80,60,140)"),
    c(0.75, "rgb(100,80,180)"),
    c(1, "rgb(140,120,220)")
  ),
  contours = list(
    z = list(
      show = TRUE,
      usecolormap = TRUE,
      highlightcolor = "#ff0000",
      project = list(z = TRUE)
    )
  ),
  lighting = list(
    ambient = 0.4,
    diffuse = 0.6,
    specular = 0.2,
    roughness = 0.85
  )
) %>%
  layout(
    title = "fMRI OFF Condition - Kime-Surface",
    scene = list(
      xaxis = list(title = "X (t·cos θ)"),
      yaxis = list(title = "Y (t·sin θ)"),
      zaxis = list(title = "Z (Intensity)"),
      camera = list(
        eye = list(x = 1.5, y = 1.5, z = 1.5)
      )
    )
  )

# Display the plots
fig_on
fig_off

6.2 Comparison Plot

# Create side-by-side comparison
fig_compare <- subplot(
  fig_on %>% layout(scene = list(domain = list(x = c(0, 0.45)))),
  fig_off %>% layout(scene = list(domain = list(x = c(0.55, 1)))),
  titleX = TRUE,
  titleY = TRUE
) %>%
  layout(
    title = "KPT Kime-Surfaces: fMRI ON vs OFF Conditions",
    showlegend = FALSE
  )

fig_compare
##### 3 synchronized fMRI On, Off, Difference Kime-surfaces
# ---------------------------- 
# Create Synchronized 3D Kime-Surfaces
# ---------------------------- 

# Prepare data for synchronized plots
prepare_surface_data <- function(surface) {
  return(list(
    x_grid = surface$x,
    y_grid = surface$y,
    z_grid = surface$z,
    magnitude = sqrt(surface$x^2 + surface$y^2),
    phase = atan2(surface$y, surface$x)
  ))
}

pOn <- prepare_surface_data(surface_on)
pOff <- prepare_surface_data(surface_off)

# Create ON surface plot
plot_on <- plot_ly(
  x = pOn$x_grid,
  y = pOn$y_grid,
  z = pOn$z_grid,
  type = "surface",
  colorscale = list(
    c(0, "rgb(20,20,80)"),
    c(0.25, "rgb(20,80,140)"),
    c(0.5, "rgb(20,140,200)"),
    c(0.75, "rgb(80,200,120)"),
    c(1, "rgb(255,220,80)")
  ),
  hoverinfo = "x+y+z",
  name = "ON",
  scene = "scene",
  showscale = FALSE
) %>%
  layout(scene = list(
    xaxis = list(title = "X (t·cos θ)"),
    yaxis = list(title = "Y (t·sin θ)"),
    zaxis = list(title = "Intensity"),
    aspectmode = "cube"
  ))

# Create OFF surface plot
plot_off <- plot_ly(
  x = pOff$x_grid,
  y = pOff$y_grid,
  z = pOff$z_grid,
  type = "surface",
  colorscale = list(
    c(0, "rgb(40,20,60)"),
    c(0.25, "rgb(60,40,100)"),
    c(0.5, "rgb(80,60,140)"),
    c(0.75, "rgb(100,80,180)"),
    c(1, "rgb(140,120,220)")
  ),
  hoverinfo = "x+y+z",
  name = "OFF",
  scene = "scene2",
  showscale = FALSE
) %>%
  layout(scene2 = list(
    xaxis = list(title = "X (t·cos θ)"),
    yaxis = list(title = "Y (t·sin θ)"),
    zaxis = list(title = "Intensity"),
    aspectmode = "cube"
  ))

# Create difference plot
diff_plot <- plot_ly(
  hoverinfo = "x+y+z",
  showscale = TRUE,
  scene = "scene3"
) %>%
  add_trace(
    x = pOn$x_grid,
    y = pOn$y_grid,
    z = pOn$z_grid - pOff$z_grid,
    surfacecolor = pOn$z_grid - pOff$z_grid,
    type = 'surface',
    colorscale = list(
      c(0, "rgb(0,0,255)"),
      c(0.5, "rgb(255,255,255)"),
      c(1, "rgb(255,0,0)")
    ),
    colorbar = list(
      title = "Difference",
      x = 1.1
    ),
    opacity = 1,
    visible = TRUE
  ) %>%
  layout(
    scene3 = list(
      aspectmode = "cube",
      xaxis = list(title = "X (t·cos θ)"),
      yaxis = list(title = "Y (t·sin θ)"),
      zaxis = list(title = "Z = ON - OFF")
    )
  )

# Combine into synchronized subplot
main_plot <- subplot(plot_on, plot_off, diff_plot) %>%
  layout(
    title = "Synchronized Kime-Surfaces (Top-Left: ON | Top-Right: OFF | Bottom: Difference)",
    scene = list(
      domain = list(x = c(0, 0.45), y = c(0.5, 1)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    scene2 = list(
      domain = list(x = c(0.55, 1), y = c(0.5, 1)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    scene3 = list(
      domain = list(x = c(0.275, 0.725), y = c(0, 0.45)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    showlegend = FALSE
  )

# Add camera synchronization
synchronized_plot <- main_plot %>%
  htmlwidgets::onRender(
    "function(x, el) {
      x.on('plotly_relayout', function(d) {
        const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key));
        if (camera.length) {
          const scenes = Object.keys(x.layout).filter((key) => /^scene\\d*/.test(key));
          const new_layout = {};
          scenes.forEach(key => {
            new_layout[key] = {...x.layout[key], camera: {...d[camera]}};
          });
          Plotly.relayout(x, new_layout);
        }
      });
    }"
  )

# Display the synchronized plot
synchronized_plot
# ---------------------------- 
# Create Synchronized 3D Kime-Surfaces
# ---------------------------- 

# Prepare data for synchronized plots with cylindrical coordinates
prepare_surface_data <- function(surface) {
  # Calculate cylindrical coordinates
  r_matrix <- sqrt(surface$x^2 + surface$y^2)
  phi_matrix <- atan2(surface$y, surface$x)
  
  return(list(
    x_grid = surface$x,
    y_grid = surface$y,
    z_grid = surface$z,
    r_grid = r_matrix,
    phi_grid = phi_matrix,
    magnitude = r_matrix,
    phase = phi_matrix
  ))
}

pOn <- prepare_surface_data(surface_on)
pOff <- prepare_surface_data(surface_off)

# Create custom hover text for cylindrical coordinates
create_hover_text <- function(r, phi, z, label = "") {
  hover_matrix <- matrix(NA, nrow = nrow(r), ncol = ncol(r))
  for(i in 1:nrow(r)) {
    for(j in 1:ncol(r)) {
      hover_matrix[i,j] <- sprintf(
        "%sTime: %.1f<br>φ: %.3f<br>Intensity: %.3f",
        ifelse(label == "", "", paste0(label, "<br>")),
        r[i,j],
        phi[i,j],
        z[i,j]
      )
    }
  }
  return(hover_matrix)
}

# Create ON surface plot with cylindrical coordinates
plot_on <- plot_ly(
  x = pOn$x_grid,
  y = pOn$y_grid,
  z = pOn$z_grid,
  type = "surface",
  colorscale = list(
    c(0, "rgb(20,20,80)"),
    c(0.25, "rgb(20,80,140)"),
    c(0.5, "rgb(20,140,200)"),
    c(0.75, "rgb(80,200,120)"),
    c(1, "rgb(255,220,80)")
  ),
  hovertext = create_hover_text(pOn$r_grid, pOn$phi_grid, pOn$z_grid, "ON Condition"),
  hoverinfo = "text",
  name = "ON",
  scene = "scene",
  showscale = FALSE
) %>%
  layout(scene = list(
    xaxis = list(title = "X (t·cos θ)"),
    yaxis = list(title = "Y (t·sin θ)"),
    zaxis = list(title = "Intensity"),
    aspectmode = "cube"
  ))

# Create OFF surface plot with cylindrical coordinates
plot_off <- plot_ly(
  x = pOff$x_grid,
  y = pOff$y_grid,
  z = pOff$z_grid,
  type = "surface",
  colorscale = list(
    c(0, "rgb(40,20,60)"),
    c(0.25, "rgb(60,40,100)"),
    c(0.5, "rgb(80,60,140)"),
    c(0.75, "rgb(100,80,180)"),
    c(1, "rgb(140,120,220)")
  ),
  hovertext = create_hover_text(pOff$r_grid, pOff$phi_grid, pOff$z_grid, "OFF Condition"),
  hoverinfo = "text",
  name = "OFF",
  scene = "scene2",
  showscale = FALSE
) %>%
  layout(scene2 = list(
    xaxis = list(title = "X (t·cos θ)"),
    yaxis = list(title = "Y (t·sin θ)"),
    zaxis = list(title = "Intensity"),
    aspectmode = "cube"
  ))

# Create difference plot with cylindrical coordinates
diff_z <- pOn$z_grid - pOff$z_grid
diff_plot <- plot_ly(
  hovertext = create_hover_text(pOn$r_grid, pOn$phi_grid, diff_z, "Difference (ON - OFF)"),
  hoverinfo = "text",
  showscale = TRUE,
  scene = "scene3"
) %>%
  add_trace(
    x = pOn$x_grid,
    y = pOn$y_grid,
    z = diff_z,
    surfacecolor = diff_z,
    type = 'surface',
    colorscale = list(
      c(0, "rgb(0,0,255)"),
      c(0.5, "rgb(255,255,255)"),
      c(1, "rgb(255,0,0)")
    ),
    colorbar = list(
      title = "Difference",
      x = 1.1
    ),
    opacity = 1,
    visible = TRUE
  ) %>%
  layout(
    scene3 = list(
      aspectmode = "cube",
      xaxis = list(title = "X (t·cos θ)"),
      yaxis = list(title = "Y (t·sin θ)"),
      zaxis = list(title = "Z = ON - OFF")
    )
  )

# Combine into synchronized subplot
main_plot <- subplot(plot_on, plot_off, diff_plot) %>%
  layout(
    title = "Synchronized Kime-Surfaces (Top-Left: ON | Top-Right: OFF | Bottom: Difference)",
    scene = list(
      domain = list(x = c(0, 0.45), y = c(0.5, 1)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    scene2 = list(
      domain = list(x = c(0.55, 1), y = c(0.5, 1)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    scene3 = list(
      domain = list(x = c(0.275, 0.725), y = c(0, 0.45)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    showlegend = FALSE
  )

# Add camera synchronization
synchronized_plot <- main_plot %>%
  htmlwidgets::onRender(
    "function(x, el) {
      x.on('plotly_relayout', function(d) {
        const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key));
        if (camera.length) {
          const scenes = Object.keys(x.layout).filter((key) => /^scene\\d*/.test(key));
          const new_layout = {};
          scenes.forEach(key => {
            new_layout[key] = {...x.layout[key], camera: {...d[camera]}};
          });
          Plotly.relayout(x, new_layout);
        }
      });
    }"
  )

# Display the synchronized plot
synchronized_plot

6.3 Phase Distribution Analysis

# Analyze phase distributions at different time points
time_points <- c(25, 75, 125)
phase_analysis <- data.frame()

for(condition in c("ON", "OFF")) {
  kpt_result <- if(condition == "ON") kpt_on else kpt_off
  
  for(t_idx in time_points) {
    phase_data <- kpt_result$phase_data[[t_idx/2]]  # Adjust for TR=2
    
    df_temp <- data.frame(
      condition = condition,
      time = t_idx,
      theta = phase_data$theta_samples,
      intensity = phase_data$s_values
    )
    phase_analysis <- rbind(phase_analysis, df_temp)
  }
}

# Plot phase distributions
ggplot(phase_analysis, aes(x = theta, fill = condition)) +
  geom_histogram(bins = 30, alpha = 0.6, position = "identity") +
  facet_wrap(~time, ncol = 3, labeller = labeller(time = function(x) paste("t =", x, "s"))) +
  scale_fill_manual(values = c("ON" = "red", "OFF" = "blue")) +
  xlim(-pi, pi) +
  labs(
    title = "Phase Distributions at Different Time Points",
    x = "Phase (θ)",
    y = "Count",
    fill = "Condition"
  ) +
  theme_minimal()

6.4 Radial Intensity Profiles

# Extract radial intensity profiles
radial_profiles <- data.frame()

for(condition in c("ON", "OFF")) {
  surface <- if(condition == "ON") surface_on else surface_off
  
  # Average intensity at each radius
  for(i in 1:length(surface$r_grid)) {
    avg_intensity <- mean(surface$z[i, ], na.rm = TRUE)
    
    radial_profiles <- rbind(radial_profiles, data.frame(
      condition = condition,
      radius = surface$r_grid[i],
      intensity = avg_intensity
    ))
  }
}

# Plot radial profiles
ggplot(radial_profiles, aes(x = radius, y = intensity, color = condition)) +
  geom_line(size = 1.5) +
  scale_color_manual(values = c("ON" = "red", "OFF" = "blue")) +
  labs(
    title = "Radial Intensity Profiles",
    x = "Time (t)",
    y = "Average Intensity",
    color = "Condition"
  ) +
  theme_minimal()

7 Statistical Analysis

# Compute surface metrics
compute_surface_metrics <- function(surface) {
  # Total volume under surface
  volume <- sum(surface$z, na.rm = TRUE) * 
    (surface$r_grid[2] - surface$r_grid[1]) * 
    (surface$theta_grid[2] - surface$theta_grid[1])
  
  # Surface roughness (standard deviation of gradients)
  grad_r <- diff(surface$z) / (surface$r_grid[2] - surface$r_grid[1])
  grad_theta <- t(diff(t(surface$z))) / (surface$theta_grid[2] - surface$theta_grid[1])
  roughness <- sd(c(grad_r, grad_theta), na.rm = TRUE)
  
  # Peak intensity
  peak <- max(surface$z, na.rm = TRUE)
  
  # Radial asymmetry (coefficient of variation across angles)
  radial_cv <- apply(surface$z, 1, function(x) sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE))
  asymmetry <- mean(radial_cv, na.rm = TRUE)
  
  return(list(
    volume = volume,
    roughness = roughness,
    peak = peak,
    asymmetry = asymmetry
  ))
}

# Compare metrics
metrics_on <- compute_surface_metrics(surface_on)
metrics_off <- compute_surface_metrics(surface_off)

# Create comparison table
metrics_df <- data.frame(
  Metric = c("Volume", "Roughness", "Peak Intensity", "Asymmetry"),
  ON = unlist(metrics_on),
  OFF = unlist(metrics_off),
  Difference = unlist(metrics_on) - unlist(metrics_off),
  Ratio = unlist(metrics_on) / unlist(metrics_off)
)

knitr::kable(metrics_df, digits = 3, caption = "Kime-Surface Metrics Comparison")
Kime-Surface Metrics Comparison
Metric ON OFF Difference Ratio
volume Volume 2864.804 1512.679 1352.124 1.894
roughness Roughness 0.278 0.157 0.121 1.768
peak Peak Intensity 3.365 2.218 1.147 1.517
asymmetry Asymmetry 0.102 0.102 0.000 1.000

8 Advanced Statistical Analysis

# # ---------------------------- 
# Enhanced Difference Analysis
# ---------------------------- 

# Compute point-wise differences
diff_surface <- list(
  x = surface_on$x,
  y = surface_on$y,
  z = surface_on$z - surface_off$z,
  r_grid = surface_on$r_grid,
  theta_grid = surface_on$theta_grid
)

# Statistical significance of differences
diff_stats <- data.frame(
  radius = surface_on$r_grid,
  mean_diff = apply(diff_surface$z, 1, mean, na.rm = TRUE),
  sd_diff = apply(diff_surface$z, 1, sd, na.rm = TRUE),
  max_diff = apply(diff_surface$z, 1, max, na.rm = TRUE),
  min_diff = apply(diff_surface$z, 1, min, na.rm = TRUE)
)

# Plot difference statistics
diff_stats_plot <- plot_ly() %>%
  add_trace(x = diff_stats$radius, y = diff_stats$mean_diff,
            type = 'scatter', mode = 'lines', name = 'Mean Difference',
            line = list(color = 'black', width = 3)) %>%
  add_trace(x = c(diff_stats$radius, rev(diff_stats$radius)),
            y = c(diff_stats$mean_diff + diff_stats$sd_diff,
                  rev(diff_stats$mean_diff - diff_stats$sd_diff)),
            type = 'scatter', mode = 'none', fill = 'toself',
            fillcolor = 'rgba(100,100,100,0.3)', name = '±1 SD',
            showlegend = TRUE) %>%
  add_trace(x = diff_stats$radius, y = diff_stats$max_diff,
            type = 'scatter', mode = 'lines', name = 'Max Difference',
            line = list(color = 'red', dash = 'dash')) %>%
  add_trace(x = diff_stats$radius, y = diff_stats$min_diff,
            type = 'scatter', mode = 'lines', name = 'Min Difference',
            line = list(color = 'blue', dash = 'dash')) %>%
  layout(
    title = "Kime-Surface Differences (ON - OFF) vs Time",
    xaxis = list(title = "Time (seconds)"),
    yaxis = list(title = "Intensity Difference"),
    hovermode = 'x unified'
  )

diff_stats_plot

This implementation of the Kime-Phase Tomography (KPT) algorithm demonstrates the Kime-Phase Tomography applied to synthetic fMRI data.

  1. Successful Reconstruction: The algorithm successfully maps repeated time-series measurements to 2D kime-surfaces, revealing the underlying temporal-phase structure.

  2. Condition Differences: Clear differences are observed between ON and OFF conditions:

  • ON condition shows higher intensity and more complex surface structure
  • OFF condition exhibits smoother, lower-intensity surfaces
  • Phase distributions vary systematically with time
  1. Key Insights:
  • The kime-surface representation captures both temporal dynamics (radial direction) and phase variability (angular direction)
  • The harmonic expansion provides a compact representation of the signal structure
  • The wrapped Laplace distribution effectively models phase uncertainty
  1. Future Directions:
  • Apply to real fMRI data
  • Explore different phase distributions (von Mises, wrapped Cauchy)
  • Investigate higher-order harmonic expansions
  • Develop statistical tests for comparing kime-surfaces

KPT provides a geometric QM-perspective on longitudinal time-series data and accounts for the structure of repeated measurements and their inherent variability.

SOCR Resource Visitor number Web Analytics SOCR Email