--- title: "Time Complexity, Inferential uncertainty and Spacekime Analytics" subtitle: "Appendix 03.1: Kime-Phase Tomography (KPT) - fMRI Simulation" author: "SOCR Team" date: "`r format(Sys.time(),'%m/%d/%Y')`" output: html_document: theme: spacelab highlight: tango includes: before_body: TCIU_header.html toc: true number_sections: true toc_depth: 3 toc_float: collapsed: false smooth_scroll: true code_folding: hide word_document: toc: true toc_depth: '3' header-includes: - \usepackage{amsmath} - \usepackage{amssymb} - \usepackage{physics} - \usepackage{mathrsfs} - \usepackage{bm} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, fig.width = 10, fig.height = 8) ``` # Introduction This [TCIU Appendix](https://tciu.predictive.space/) implements the [Kime-Phase Tomography (KPT) algorithm](https://www.socr.umich.edu/TCIU/HTMLs/TCIU_SK_Appendix03_KPT.html) 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. # Required Libraries ```{r 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) ``` # Simulate fMRI Data ```{r simulate-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) ``` ## fMRI Simulated Data Visualization ```{r simulate-fMRI-data-Viz} # ---------------------------- # 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 ``` ## Simulated fMRI Data Summary Statistics ```{r simulate-fMRI-data-Summary-Stats} # ---------------------------- # 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") cat(sprintf("ON Condition - Global Mean: %.3f (SD: %.3f)\n", stats_on$global_mean, stats_on$global_sd)) cat(sprintf("OFF Condition - Global Mean: %.3f (SD: %.3f)\n", stats_off$global_mean, stats_off$global_sd)) cat(sprintf("Mean Difference (ON-OFF): %.3f\n", stats_on$global_mean - stats_off$global_mean)) cat(sprintf("Average tSNR - ON: %.3f, OFF: %.3f\n", mean(stats_on$tsnr), mean(stats_off$tsnr))) ``` # Implement KPT Core Functions ## Phase Distribution Functions ```{r phase-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) } ``` ## Harmonic Expansion and Moment Computation ```{r harmonic-functions} # 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) } ``` ## KPT Reconstruction Algorithm ```{r kpt-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 )) } ``` ## Kime-Surface Generation ```{r kimesurface-functions} # 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 )) } ``` # Apply KPT to fMRI Data ```{r apply-kpt} # Apply KPT to ON and OFF conditions cat("Applying KPT to ON condition...\n") kpt_on <- kpt_reconstruct(fmri_data$ON, n_max = 8) cat("Applying KPT to OFF condition...\n") 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) ``` # Visualize Results ## Interactive 3D Kime-Surfaces ```{r visualize-3d} # 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 ``` ## Comparison Plot ```{r 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
φ: %.3f
Intensity: %.3f", ifelse(label == "", "", paste0(label, "
")), 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 ``` ## Phase Distribution Analysis ```{r phase-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() ``` ## Radial Intensity Profiles ```{r radial-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() ``` # Statistical Analysis ```{r statistics} # 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") ``` # Advanced Statistical Analysis ```{r advanced-statistics} # # ---------------------------- # 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](https://www.socr.umich.edu/TCIU/HTMLs/TCIU_SK_Appendix03_KPT.html) 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 3. *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 4. *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.