SOCR ≫

TCIU Website ≫

TCIU GitHub ≫

Time Complexity, Inferential uncertainty and Spacekime Analytics

Appendix 03.1: Kime-Phase Tomography (KPT) - fMRI Simulation

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