The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Hemodynamic Response Functions

Bradley R. Buchsbaum

2025-09-11

Introduction to Hemodynamic Response Functions (HRFs)

A hemodynamic response function (HRF) models the temporal evolution of the fMRI BOLD (Blood-Oxygen-Level-Dependent) signal in response to a brief neural event. Typically, the BOLD signal peaks 4-6 seconds after the event onset and then returns to baseline, often with a slight undershoot.

fmrihrf provides tools to define, manipulate, and visualize various HRFs commonly used in fMRI analysis.

Pre-defined HRF Objects

fmrihrf includes several pre-defined HRF objects, which are essentially functions with specific attributes defining their type, number of basis functions (nbasis), and effective duration (span).

Let’s look at two common examples: the SPM canonical HRF (HRF_SPMG1) and a Gaussian HRF (HRF_GAUSSIAN).

# SPM canonical HRF (based on difference of two gamma functions)
print(HRF_SPMG1)
#> function (t, P1 = 5, P2 = 15, A1 = 0.0833) 
#> {
#>     ifelse(t < 0, 0, exp(-t) * (A1 * t^P1 - 1.274527e-13 * t^P2))
#> }
#> <bytecode: 0x12149de20>
#> <environment: namespace:fmrihrf>
#> attr(,"class")
#> [1] "HRF"      "function"
#> attr(,"name")
#> [1] "SPMG1"
#> attr(,"nbasis")
#> [1] 1
#> attr(,"span")
#> [1] 24
#> attr(,"param_names")
#> [1] "P1" "P2" "A1"
#> attr(,"params")
#> attr(,"params")$P1
#> [1] 5
#> 
#> attr(,"params")$P2
#> [1] 15
#> 
#> attr(,"params")$A1
#> [1] 0.0833

# Gaussian HRF
print(HRF_GAUSSIAN)
#> function (t, mean = 6, sd = 2) 
#> {
#>     stats::dnorm(t, mean = mean, sd = sd)
#> }
#> <bytecode: 0x1215240a8>
#> <environment: namespace:fmrihrf>
#> attr(,"class")
#> [1] "HRF"      "function"
#> attr(,"name")
#> [1] "gaussian"
#> attr(,"nbasis")
#> [1] 1
#> attr(,"span")
#> [1] 24
#> attr(,"param_names")
#> [1] "mean" "sd"  
#> attr(,"params")
#> attr(,"params")$mean
#> [1] 6
#> 
#> attr(,"params")$sd
#> [1] 2

These objects are functions themselves, so you can evaluate them at specific time points:

time_points <- seq(0, 25, by = 0.1)

# Evaluate the HRFs
y_spm <- HRF_SPMG1(time_points)
y_gauss <- HRF_GAUSSIAN(time_points)

# Manually scale each to peak at 1.0 for easier shape comparison
y_spm_scaled <- y_spm / max(y_spm)
y_gauss_scaled <- y_gauss / max(y_gauss)

# Combine into a data frame for plotting
plot_df <- data.frame(
  Time = time_points,
  SPM_Canonical = y_spm_scaled,
  Gaussian = y_gauss_scaled
) %>% 
  pivot_longer(-Time, names_to = "HRF_Type", values_to = "Response")

# Plot using ggplot2
ggplot(plot_df, aes(x = Time, y = Response, color = HRF_Type)) +
  geom_line(linewidth = 1) +
  labs(title = "Comparison of SPM Canonical and Gaussian HRFs",
       subtitle = "HRFs manually scaled to peak at 1.0 for comparison",
       x = "Time (seconds)",
       y = "BOLD Response (normalized)",
       color = "HRF Type") +
  theme_minimal()

Note that the span attribute (e.g., 24 seconds) indicates the approximate time window over which the HRF is non-zero.

Modifying HRF Parameters with gen_hrf

The gen_hrf function is a flexible way to create new HRF functions, often by modifying the parameters of existing ones.

For example, the hrf_gaussian function takes mean and sd arguments. We can use gen_hrf to create Gaussian HRFs with different peak times (mean) and widths (sd).

# Create Gaussian HRFs with different parameters using gen_hrf
# Note: hrf_gaussian is the underlying function, not the HRF object HRF_GAUSSIAN
hrf_gauss_7_3 <- gen_hrf(hrf_gaussian, mean = 7, sd = 3, name = "Gaussian (Mean=7, SD=3)")
hrf_gauss_5_2 <- gen_hrf(hrf_gaussian, mean = 5, sd = 2, name = "Gaussian (Mean=5, SD=2)")
hrf_gauss_4_1 <- gen_hrf(hrf_gaussian, mean = 4, sd = 1, name = "Gaussian (Mean=4, SD=1)")

# Evaluate the new HRFs
vals1 <- hrf_gauss_7_3(time_points)
vals2 <- hrf_gauss_5_2(time_points)
vals3 <- hrf_gauss_4_1(time_points)

# Combine for plotting
plot_df_params <- data.frame(
  Time = time_points,
  `Mean=7, SD=3` = vals1,
  `Mean=5, SD=2` = vals2,
  `Mean=4, SD=1` = vals3
) %>%
  pivot_longer(-Time, names_to = "Parameters", values_to = "Response")

# Plot
ggplot(plot_df_params, aes(x = Time, y = Response, color = Parameters)) +
  geom_line(linewidth = 1) +
  labs(title = "Gaussian HRFs with Different Parameters",
       x = "Time (seconds)",
       y = "BOLD Response",
       color = "Parameters") +
  theme_minimal()

gen_hrf can also directly incorporate lags and durations (see later sections).

Modeling Event Duration with block_hrf

fMRI events often have a duration (e.g., a stimulus presented for several seconds). The block_hrf function (or gen_hrf with a width argument) modifies an HRF to model the response to a sustained event of a specific width (duration). Internally, it convolves the original HRF with a boxcar function of the specified width.

The precision argument controls the sampling resolution used for this convolution.

# Create blocked HRFs using the SPM canonical HRF with different durations
hrf_spm_w1 <- block_hrf(HRF_SPMG1, width = 1)
hrf_spm_w2 <- block_hrf(HRF_SPMG1, width = 2)
hrf_spm_w4 <- block_hrf(HRF_SPMG1, width = 4)

# Evaluate
resp_w1 <- hrf_spm_w1(time_points)
resp_w2 <- hrf_spm_w2(time_points)
resp_w4 <- hrf_spm_w4(time_points)

# Combine for plotting
plot_df_blocked <- data.frame(
  Time = time_points,
  `Width=1s` = resp_w1,
  `Width=2s` = resp_w2,
  `Width=4s` = resp_w4
) %>%
  pivot_longer(-Time, names_to = "Duration", values_to = "Response")

# Plot
ggplot(plot_df_blocked, aes(x = Time, y = Response, color = Duration)) +
  geom_line(linewidth = 1) +
  labs(title = "SPM Canonical HRF for Different Event Durations",
       subtitle = "Using block_hrf()",
       x = "Time (seconds)",
       y = "BOLD Response",
       color = "Duration") +
  theme_minimal()

Normalization

By default, longer durations lead to higher peak responses (assuming summation, see next section). Setting normalize=TRUE in block_hrf (or gen_hrf) rescales the response so the peak amplitude is approximately 1, regardless of duration.

# Create normalized blocked HRFs
hrf_spm_w1_norm <- block_hrf(HRF_SPMG1, width = 1, normalize = TRUE)
hrf_spm_w2_norm <- block_hrf(HRF_SPMG1, width = 2, normalize = TRUE)
hrf_spm_w4_norm <- block_hrf(HRF_SPMG1, width = 4, normalize = TRUE)

# Evaluate
resp_w1_norm <- hrf_spm_w1_norm(time_points)
resp_w2_norm <- hrf_spm_w2_norm(time_points)
resp_w4_norm <- hrf_spm_w4_norm(time_points)

# Combine for plotting
plot_df_blocked_norm <- data.frame(
  Time = time_points,
  `Width=1s` = resp_w1_norm,
  `Width=2s` = resp_w2_norm,
  `Width=4s` = resp_w4_norm
) %>%
  pivot_longer(-Time, names_to = "Duration", values_to = "Response")

# Plot
ggplot(plot_df_blocked_norm, aes(x = Time, y = Response, color = Duration)) +
  geom_line(linewidth = 1) +
  labs(title = "Normalized SPM Canonical HRF for Different Durations",
       subtitle = "Using block_hrf(normalize = TRUE)",
       x = "Time (seconds)",
       y = "BOLD Response",
       color = "Duration") +
  theme_minimal() +
  ylim(0, NA) # Ensure y-axis starts at 0

Modeling Saturation with summate

The summate argument in block_hrf controls whether the response accumulates over the duration (summate=TRUE, default) or saturates (summate=FALSE). Saturation implies that the response reaches a plateau and doesn’t increase further with longer stimulation.

# Create non-summating blocked HRFs
hrf_spm_w2_nosum <- block_hrf(HRF_SPMG1, width = 2, summate = FALSE)
hrf_spm_w4_nosum <- block_hrf(HRF_SPMG1, width = 4, summate = FALSE)
hrf_spm_w8_nosum <- block_hrf(HRF_SPMG1, width = 8, summate = FALSE)

# Evaluate
resp_w2_nosum <- hrf_spm_w2_nosum(time_points)
resp_w4_nosum <- hrf_spm_w4_nosum(time_points)
resp_w8_nosum <- hrf_spm_w8_nosum(time_points)

# Combine for plotting
plot_df_blocked_nosum <- data.frame(
  Time = time_points,
  `Width=2s` = resp_w2_nosum,
  `Width=4s` = resp_w4_nosum,
  `Width=8s` = resp_w8_nosum
) %>%
  pivot_longer(-Time, names_to = "Duration", values_to = "Response")

# Plot
ggplot(plot_df_blocked_nosum, aes(x = Time, y = Response, color = Duration)) +
  geom_line(linewidth = 1) +
  labs(title = "Non-Summating (Saturating) SPM HRF for Different Durations",
       subtitle = "Using block_hrf(summate = FALSE)",
       x = "Time (seconds)",
       y = "BOLD Response",
       color = "Duration") +
  theme_minimal()

We can combine summate=FALSE and normalize=TRUE:

# Create normalized, non-summating blocked HRFs
hrf_spm_w2_nosum_norm <- block_hrf(HRF_SPMG1, width = 2, summate = FALSE, normalize = TRUE)
hrf_spm_w4_nosum_norm <- block_hrf(HRF_SPMG1, width = 4, summate = FALSE, normalize = TRUE)
hrf_spm_w8_nosum_norm <- block_hrf(HRF_SPMG1, width = 8, summate = FALSE, normalize = TRUE)

# Evaluate
resp_w2_nosum_norm <- hrf_spm_w2_nosum_norm(time_points)
resp_w4_nosum_norm <- hrf_spm_w4_nosum_norm(time_points)
resp_w8_nosum_norm <- hrf_spm_w8_nosum_norm(time_points)

# Combine for plotting
plot_df_blocked_nosum_norm <- data.frame(
  Time = time_points,
  `Width=2s` = resp_w2_nosum_norm,
  `Width=4s` = resp_w4_nosum_norm,
  `Width=8s` = resp_w8_nosum_norm
) %>%
  pivot_longer(-Time, names_to = "Duration", values_to = "Response")

# Plot
ggplot(plot_df_blocked_nosum_norm, aes(x = Time, y = Response, color = Duration)) +
  geom_line(linewidth = 1) +
  labs(title = "Normalized, Non-Summating SPM HRF for Different Durations",
       subtitle = "Using block_hrf(summate = FALSE, normalize = TRUE)",
       x = "Time (seconds)",
       y = "BOLD Response",
       color = "Duration") +
  theme_minimal() +
  ylim(0, NA)

Modeling Temporal Shifts with lag_hrf

Sometimes, the hemodynamic response might be delayed or advanced relative to the event onset. The lag_hrf function (or gen_hrf_lagged) shifts an existing HRF in time by a specified lag (in seconds). A positive lag delays the response, while a negative lag advances it.

# Create lagged versions of the Gaussian HRF
hrf_gauss_lag_neg2 <- lag_hrf(HRF_GAUSSIAN, lag = -2)
hrf_gauss_lag_0 <- HRF_GAUSSIAN # Original (lag=0)
hrf_gauss_lag_pos3 <- lag_hrf(HRF_GAUSSIAN, lag = 3)

# Evaluate
resp_lag_neg2 <- hrf_gauss_lag_neg2(time_points)
resp_lag_0 <- hrf_gauss_lag_0(time_points)
resp_lag_pos3 <- hrf_gauss_lag_pos3(time_points)

# Combine for plotting
plot_df_lagged <- data.frame(
  Time = time_points,
  `Lag=-2s` = resp_lag_neg2,
  `Lag= 0s` = resp_lag_0,
  `Lag=+3s` = resp_lag_pos3
) %>%
  pivot_longer(-Time, names_to = "Lag", values_to = "Response")

# Plot
ggplot(plot_df_lagged, aes(x = Time, y = Response, color = Lag)) +
  geom_line(linewidth = 1) +
  labs(title = "Gaussian HRF with Different Temporal Lags",
       subtitle = "Using lag_hrf()",
       x = "Time (seconds)",
       y = "BOLD Response",
       color = "Lag") +
  theme_minimal()

Combining Lag and Duration

We can combine lag_hrf and block_hrf using the pipe operator (%>%) from dplyr (or magrittr).

# Create HRFs that are both lagged and blocked
hrf_lb_1 <- HRF_GAUSSIAN %>% lag_hrf(1) %>% block_hrf(width = 1, normalize = TRUE)
hrf_lb_3 <- HRF_GAUSSIAN %>% lag_hrf(3) %>% block_hrf(width = 3, normalize = TRUE)
hrf_lb_5 <- HRF_GAUSSIAN %>% lag_hrf(5) %>% block_hrf(width = 5, normalize = TRUE)

# Evaluate
resp_lb_1 <- hrf_lb_1(time_points)
resp_lb_3 <- hrf_lb_3(time_points)
resp_lb_5 <- hrf_lb_5(time_points)

# Combine for plotting
plot_df_lagged_blocked <- data.frame(
  Time = time_points,
  `Lag=1, Width=1` = resp_lb_1,
  `Lag=3, Width=3` = resp_lb_3,
  `Lag=5, Width=5` = resp_lb_5
) %>%
  pivot_longer(-Time, names_to = "Settings", values_to = "Response")

# Plot
ggplot(plot_df_lagged_blocked, aes(x = Time, y = Response, color = Settings)) +
  geom_line(linewidth = 1) +
  labs(title = "Gaussian HRFs with Combined Lag and Duration",
       subtitle = "Using lag_hrf() %>% block_hrf()",
       x = "Time (seconds)",
       y = "BOLD Response",
       color = "Settings") +
  theme_minimal() +
  ylim(0, NA)

Alternatively, gen_hrf can apply lag and width directly:

# Using gen_hrf directly
hrf_lb_gen_3 <- gen_hrf(hrf_gaussian, lag = 3, width = 3, normalize = TRUE)
resp_lb_gen_3 <- hrf_lb_gen_3(time_points)

# Compare (should be very similar to hrf_lb_3 from piped version)
# plot(time_points, resp_lb_3, type = 'l', col = 2, lwd = 2, main = "Piped vs gen_hrf")
# lines(time_points, resp_lb_gen_3, col = 1, lty = 2, lwd = 2)
# legend("topright", legend = c("Piped", "gen_hrf"), col = c(2, 1), lty = c(1, 2), lwd = 2)

Multivariate HRFs: Basis Sets

Instead of assuming a fixed HRF shape, we can model the response using a linear combination of multiple basis functions. This allows for more flexibility in capturing variations in HRF shape across brain regions or individuals. The resulting HRF function returns a matrix where each column corresponds to a basis function.

SPM Basis Sets

fmrihrf provides pre-defined HRF objects for the SPM canonical HRF plus its temporal derivative (HRF_SPMG2), and additionally its dispersion derivative (HRF_SPMG3).

# SPM + Temporal Derivative (2 basis functions)
print(HRF_SPMG2)
#> function (t) 
#> {
#>     do.call(cbind, lapply(xs, function(f) f(t)))
#> }
#> <bytecode: 0x13018a908>
#> <environment: 0x130180ef0>
#> attr(,"class")
#> [1] "SPMG2_HRF" "HRF"       "function" 
#> attr(,"name")
#> [1] "SPMG2"
#> attr(,"nbasis")
#> [1] 2
#> attr(,"span")
#> [1] 24
#> attr(,"params")
#> list()
resp_spmg2 <- HRF_SPMG2(time_points)

# SPM + Temporal + Dispersion Derivatives (3 basis functions)
print(HRF_SPMG3)
#> function (t) 
#> {
#>     do.call(cbind, lapply(xs, function(f) f(t)))
#> }
#> <bytecode: 0x130228b38>
#> <environment: 0x13021bcf0>
#> attr(,"class")
#> [1] "SPMG3_HRF" "HRF"       "function" 
#> attr(,"name")
#> [1] "SPMG3"
#> attr(,"nbasis")
#> [1] 3
#> attr(,"span")
#> [1] 24
#> attr(,"params")
#> list()
resp_spmg3 <- HRF_SPMG3(time_points)

# Plot SPMG2
matplot(time_points, resp_spmg2, type = 'l', lty = 1, lwd = 1.5,
        xlab = "Time (seconds)", ylab = "BOLD Response",
        main = "SPM + Temporal Derivative Basis Set (HRF_SPMG2)")
legend("topright", legend = c("Canonical", "Temporal Deriv."), col = 1:2, lty = 1, lwd = 1.5)


# Plot SPMG3
matplot(time_points, resp_spmg3, type = 'l', lty = 1, lwd = 1.5,
        xlab = "Time (seconds)", ylab = "BOLD Response",
        main = "SPM + Temporal + Dispersion Derivative Basis Set (HRF_SPMG3)")
legend("topright", legend = c("Canonical", "Temporal Deriv.", "Dispersion Deriv."), col = 1:3, lty = 1, lwd = 1.5)

B-Spline Basis Set

The hrf_bspline function generates a B-spline basis set. We typically use it within gen_hrf to create an HRF object. Key parameters are N (number of basis functions) and degree.

# B-spline basis with N=5 basis functions, degree=3 (cubic)
hrf_bs_5_3 <- gen_hrf(hrf_bspline, N = 5, degree = 3, name = "B-spline (N=5, deg=3)")
print(hrf_bs_5_3)
#> function (t) 
#> hrf(t, ...)
#> <bytecode: 0x127e322e8>
#> <environment: 0x130a5c400>
#> attr(,"class")
#> [1] "HRF"      "function"
#> attr(,"name")
#> [1] "B-spline (N=5, deg=3)"
#> attr(,"nbasis")
#> [1] 5
#> attr(,"span")
#> [1] 24
#> attr(,"params")
#> list()
resp_bs_5_3 <- hrf_bs_5_3(time_points)
matplot(time_points, resp_bs_5_3, type = 'l', lty = 1, lwd = 1.5,
        xlab = "Time (seconds)", ylab = "BOLD Response",
        main = "B-spline Basis Set (N=5, degree=3)")


# B-spline basis with N=10 basis functions, degree=1 (linear -> tent functions)
hrf_bs_10_1 <- gen_hrf(hrf_bspline, N = 10, degree = 1, name = "Tent Set (N=10)")
print(hrf_bs_10_1)
#> function (t) 
#> hrf(t, ...)
#> <bytecode: 0x127e322e8>
#> <environment: 0x1305dff70>
#> attr(,"class")
#> [1] "HRF"      "function"
#> attr(,"name")
#> [1] "Tent Set (N=10)"
#> attr(,"nbasis")
#> [1] 10
#> attr(,"span")
#> [1] 24
#> attr(,"params")
#> list()
resp_bs_10_1 <- hrf_bs_10_1(time_points)
matplot(time_points, resp_bs_10_1, type = 'l', lty = 1, lwd = 1.5,
        xlab = "Time (seconds)", ylab = "BOLD Response",
        main = "Tent Function Basis Set (B-spline, N=10, degree=1)")

Sine Basis Set

The hrf_sine function creates a basis set using sine waves of different frequencies.

hrf_sin_5 <- gen_hrf(hrf_sine, N = 5, name = "Sine Basis (N=5)")
print(hrf_sin_5)
#> function (t) 
#> hrf(t, ...)
#> <bytecode: 0x127e322e8>
#> <environment: 0x120199ae8>
#> attr(,"class")
#> [1] "HRF"      "function"
#> attr(,"name")
#> [1] "Sine Basis (N=5)"
#> attr(,"nbasis")
#> [1] 5
#> attr(,"span")
#> [1] 24
#> attr(,"params")
#> list()
resp_sin_5 <- hrf_sin_5(time_points)
matplot(time_points, resp_sin_5, type = 'l', lty = 1, lwd = 1.5,
        xlab = "Time (seconds)", ylab = "BOLD Response",
        main = "Sine Basis Set (N=5)")

Half-Cosine Basis Set (FLOBS-like)

The hrf_half_cosine function implements the basis set described by Woolrich et al. (2004), often used in FSL’s FLOBS (FMRIB’s Linear Optimal Basis Sets). It uses four half-cosine functions to model initial dip, rise, fall/undershoot, and recovery.

# Use default parameters from Woolrich et al. (2004)
# Note: hrf_half_cosine itself returns a single HRF shape, not a basis set directly.
# To use as a basis in modeling, you'd typically include it alongside derivatives
# or other basis functions. Here we just visualize the shape.
resp_half_cos <- hrf_half_cosine(time_points)
plot(time_points, resp_half_cos, type = 'l', lwd = 1.5,
     xlab = "Time (seconds)", ylab = "BOLD Response",
     main = "Half-Cosine HRF Shape (Woolrich et al., 2004)")

Other HRF Shapes

Gamma HRF

The hrf_gamma function uses the gamma probability density function.

````

hrf_gam <- gen_hrf(hrf_gamma, shape = 6, rate = 1, name = "Gamma (shape=6, rate=1)")
print(hrf_gam)
#> function (t) 
#> hrf(t, ...)
#> <bytecode: 0x127e322e8>
#> <environment: 0x12159b6a8>
#> attr(,"class")
#> [1] "HRF"      "function"
#> attr(,"name")
#> [1] "Gamma (shape=6, rate=1)"
#> attr(,"nbasis")
#> [1] 1
#> attr(,"span")
#> [1] 24
#> attr(,"params")
#> list()
resp_gam <- hrf_gam(time_points)
plot(time_points, resp_gam, type = 'l', lwd = 1.5,
     xlab = "Time (seconds)", ylab = "BOLD Response",
     main = "Gamma HRF")

Mexican Hat Wavelet HRF

The hrf_mexhat function uses the Mexican hat wavelet (second derivative of a Gaussian).

hrf_mh <- gen_hrf(hrf_mexhat, mean = 6, sd = 1.5, name = "Mexican Hat (mean=6, sd=1.5)")
print(hrf_mh)
#> function (t) 
#> hrf(t, ...)
#> <bytecode: 0x127e322e8>
#> <environment: 0x121d25e40>
#> attr(,"class")
#> [1] "HRF"      "function"
#> attr(,"name")
#> [1] "Mexican Hat (mean=6, sd=1.5)"
#> attr(,"nbasis")
#> [1] 1
#> attr(,"span")
#> [1] 24
#> attr(,"params")
#> list()
resp_mh <- hrf_mh(time_points)
plot(time_points, resp_mh, type = 'l', lwd = 1.5,
     xlab = "Time (seconds)", ylab = "BOLD Response",
     main = "Mexican Hat Wavelet HRF")

Inverse Logit Difference HRF

The hrf_inv_logit function creates an HRF shape by subtracting two inverse logit (sigmoid) functions, allowing control over rise and fall times.

hrf_il <- gen_hrf(hrf_inv_logit, mu1 = 5, s1 = 1, mu2 = 15, s2 = 1.5, name = "Inv. Logit Diff.")
print(hrf_il)
#> function (t) 
#> hrf(t, ...)
#> <bytecode: 0x127e322e8>
#> <environment: 0x1231f5730>
#> attr(,"class")
#> [1] "HRF"      "function"
#> attr(,"name")
#> [1] "Inv. Logit Diff."
#> attr(,"nbasis")
#> [1] 1
#> attr(,"span")
#> [1] 24
#> attr(,"params")
#> list()
resp_il <- hrf_il(time_points)
plot(time_points, resp_il, type = 'l', lwd = 1.5,
     xlab = "Time (seconds)", ylab = "BOLD Response",
     main = "Inverse Logit Difference HRF")

Creating Custom Basis Sets with gen_hrf_set

The gen_hrf_set function allows you to combine any set of HRF functions into a single multivariate HRF object (a basis set).

For example, we can create a basis set from a series of lagged Gaussian HRFs:

# Create a list of lagged Gaussian HRFs
lag_times <- seq(0, 10, by = 2)
list_of_hrfs <- lapply(lag_times, function(lag) {
  lag_hrf(HRF_GAUSSIAN, lag = lag)
})

# Combine them into a single HRF basis set object
hrf_custom_set <- do.call(gen_hrf_set, list_of_hrfs)
print(hrf_custom_set) # Note: name is default 'hrf_set', nbasis is 6
#> function (t) 
#> {
#>     do.call(cbind, lapply(xs, function(f) f(t)))
#> }
#> <bytecode: 0x1242a3a70>
#> <environment: 0x1242a2ce0>
#> attr(,"class")
#> [1] "HRF"      "function"
#> attr(,"name")
#> [1] "hrf_set"
#> attr(,"nbasis")
#> [1] 6
#> attr(,"span")
#> [1] 34
#> attr(,"params")
#> list()

# Evaluate and plot
resp_custom_set <- hrf_custom_set(time_points)
matplot(time_points, resp_custom_set, type = 'l', lty = 1, lwd = 1.5,
        xlab = "Time (seconds)", ylab = "BOLD Response",
        main = "Custom Basis Set (Lagged Gaussians)")

Creating Empirical HRFs

From a Single Measured Response (gen_empirical_hrf)

If you have a measured or estimated hemodynamic response profile (e.g., from deconvolution), you can turn it into an HRF function using gen_empirical_hrf. It uses linear interpolation between the provided points.

# Simulate an average measured response profile
sim_times <- 0:24
set.seed(42) # For reproducibility
sim_profile <- rowMeans(replicate(20, {
  h <- HRF_SPMG1 %>% lag_hrf(lag = runif(n = 1, min = -1, max = 1)) %>%
                    block_hrf(width = runif(n = 1, min = 0, max = 2))
  h(sim_times)
}))

# Normalize profile to max = 1 for better visualization
sim_profile_norm <- sim_profile / max(sim_profile)

# Create the empirical HRF function from the normalized profile
emp_hrf <- gen_empirical_hrf(sim_times, sim_profile_norm)
print(emp_hrf)
#> function (v) 
#> .approxfun(x, y, v, method, yleft, yright, f, na.rm)
#> <bytecode: 0x122f517a8>
#> <environment: 0x11735f318>
#> attr(,"class")
#> [1] "HRF"      "function"
#> attr(,"name")
#> [1] "empirical_hrf"
#> attr(,"nbasis")
#> [1] 1
#> attr(,"span")
#> [1] 24
#> attr(,"params")
#> list()

# Evaluate and plot (using a finer time grid for interpolation)
fine_times <- seq(0, 24, by = 0.1)
resp_emp <- emp_hrf(fine_times)

# Plot the interpolated curve with the original points
plot(fine_times, resp_emp, type = 'l', lwd = 1.5,
     xlab = "Time (seconds)", ylab = "BOLD Response",
     main = "Empirical HRF from Simulated Average Profile")
points(sim_times, sim_profile_norm, pch = 16, col = "red", cex = 1) # Show original points

Empirical Basis Set via PCA

You can create an empirical basis set by applying dimensionality reduction (like PCA) to a collection of observed or simulated HRFs.

# 1. Simulate a matrix of diverse HRFs
set.seed(123) # for reproducibility
n_sim <- 50
sim_mat <- replicate(n_sim, {
  hrf_func <- HRF_SPMG1 %>%
              lag_hrf(lag = runif(1, -2, 2)) %>%
              block_hrf(width = runif(1, 0, 3))
  hrf_func(sim_times)
})

# Show a sample of simulated HRFs to illustrate variability 
matplot(sim_times, sim_mat[, 1:10], type = 'l', col = scales::alpha("gray", 0.7), lty = 1,
        xlab = "Time (seconds)", ylab = "Response", 
        main = "Sample of Simulated HRF Profiles")


# 2. Perform PCA on the transpose (each column = one HRF, each row = one time point)
pca_res <- prcomp(t(sim_mat), center = TRUE, scale. = FALSE)
n_components <- 3

# Print variance explained by top components
variance_explained <- summary(pca_res)$importance[2, 1:n_components]
cat("Variance explained by top", n_components, "components:", 
    paste0(round(variance_explained * 100, 1), "%"), "\n")
#> Variance explained by top 3 components: 67.1% 29.9% 2.7%

# Extract the top principal components
pc_vectors <- pca_res$rotation[, 1:n_components]

# 3. Convert principal components into HRF functions
# Create list to store our empirical HRF functions
list_pc_hrfs <- list()

for (i in 1:n_components) {
  # Get the PC vector
  pc_vec <- pc_vectors[, i]
  
  # Start at 0 (shift so first value is 0)
  pc_vec_zeroed <- pc_vec - pc_vec[1]
  
  # Normalize peak to 1 (or -1 if the peak is negative)
  max_abs <- max(abs(pc_vec_zeroed))
  pc_vec_norm <- pc_vec_zeroed / max_abs
  
  # Create empirical HRF function
  list_pc_hrfs[[i]] <- gen_empirical_hrf(sim_times, pc_vec_norm)
}

# 4. Combine PC HRFs into a basis set using gen_hrf_set
emp_pca_basis <- do.call(gen_hrf_set, list_pc_hrfs)
print(emp_pca_basis)
#> function (t) 
#> {
#>     do.call(cbind, lapply(xs, function(f) f(t)))
#> }
#> <bytecode: 0x1242a3a70>
#> <environment: 0x130b6d890>
#> attr(,"class")
#> [1] "HRF"      "function"
#> attr(,"name")
#> [1] "hrf_set"
#> attr(,"nbasis")
#> [1] 3
#> attr(,"span")
#> [1] 24
#> attr(,"params")
#> list()

# 5. Evaluate and plot the basis functions
resp_pca_basis <- emp_pca_basis(sim_times)

# Create a prettier plot of the PCA basis functions
pc_df <- as.data.frame(resp_pca_basis)
names(pc_df) <- paste("PC", 1:n_components)
pc_df$Time <- sim_times

# Use ggplot for a nicer visualization
pc_df_long <- pivot_longer(pc_df, -Time, names_to = "Component", values_to = "Value")

ggplot(pc_df_long, aes(x = Time, y = Value, color = Component)) +
  geom_line(linewidth = 1.2) +
  scale_color_brewer(palette = "Set1") +
  labs(title = "Empirical Basis Set from PCA",
       subtitle = paste0("First ", n_components, " Principal Components"),
       x = "Time (seconds)",
       y = "Component Value") +
  theme_minimal() +
  theme(legend.position = "right")

This empirical basis set can then be used in regression models just like any other pre-defined or custom basis set.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.