Skip to content
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
84 changes: 53 additions & 31 deletions hrf_estimation.md
Original file line number Diff line number Diff line change
Expand Up @@ -8,10 +8,10 @@ touch: true

---

<!-- .slide: data-background-image="images/AdobeStock_477448716 [Converted].png" data-background-opacity="0.5" -->
<!-- slide bg="[[AdobeStock_477448716 [Converted].png]]" data-background-opacity="0.5"-->

::: block <!-- element class="title-border" -->
<h2>Estimating the Hemodynamic Response Function</h2>
## Estimating the Hemodynamic Response Function
Lynne J Williams<br>
March 28, 2022
:::
Expand All @@ -20,8 +20,8 @@ March 28, 2022

### Why estimate HRFs?

- gives better estimation of activation/deactivation [Boch et al., 2021]
- remember BOLD is a proxy for neural activation/deactivation
- gives better estimation of activation/deactivation <span style="font-size:small">[Boch et al., 2021]</span>
- remember BOLD is a proxy for neural activation/deactivation
+ allows for inter-voxel variability
- different HRF rise times for different brain regions/tasks/populations
- children's HRFs rise faster than adults'
Expand All @@ -33,35 +33,41 @@ March 28, 2022

<split even gap=2>
::: block
<h4>Model-based</h4>
Assume shape of HRF (e.g., double gamma, Poisson, radial-basis, inverse logit, etc.) <span class="small-font">
[Lei et al., 2012] </span><br>
<img src="https://i.stack.imgur.com/4jAuf.png">
<br>
<div class="small-font" style="text-align: center;">Double Gamma distribution [<a href="https://i.stack.imgur.com/4jAuf.png">Stack Exchange</a>]</div>:::
#### Model-based

Assume shape of HRF (e.g., double gamma, Poisson, radial-basis, inverse logit, etc.) <span style="font-size:small">
[Lei et al., 2012] </span>

![Image](https://i.stack.imgur.com/4jAuf.png)

Double Gamma distribution <!-- element align="center" style="font-size:smaller; margin:0"--> [[Stack Exchange](https://i.stack.imgur.com/4jAuf.png)]

:::

::: block
<h4>Model-free</h4>
Estimate voxel-specific HRF and use the estimated HRF to deconvolve BOLD signal <span class="small-font"> [Lei et al., 2012; Wink et al., 2008; Zhang et al., 2012] </span><br>
<img src="https://www.frontiersin.org/files/Articles/110699/fncom-09-00054-HTML/image_m/fncom-09-00054-g001.jpg">
<!--<img src="./images/fncom-09-00054-g0001.jpeg">--><br>
<div class="small-font" style="text-align: center;">Estimations of HRF in fMRI where black box is input [Rosa et al, 2015]</div>
#### Model-free

Estimate voxel-specific HRF and use the estimated HRF to deconvolve BOLD signal <span style="font-size:small"> [Lei et al., 2012; Wink et al., 2008; Zhang et al., 2012] </span>

![[fncom-09-00054-g0001.jpeg]]
Estimations of HRF in fMRI where black box is input<!-- element align="center" style="font-size:smaller; margin:0"--> [Rosa et al, 2015]
:::
</split>


--

### The BOLD response is assumed to be linear time-invariant
<div style="margin-top: 25%; text-align: center;">HRF can be represented as a linear combination of basis functions </div>
HRF can be represented as a linear combination of basis functions <!-- element align="center" drag="100 10" drop="0 25"-->

note: Here's the problem. We want a mathematical description of a curve or any other data distributed over space, time, and other types of continuum. Flexibility is a central issue since we usually cannot say in advance how complex the curve will be, or specify certain of its characteristics. We don't have the time or patience to search some handbook of known functions for one that looks like what we want to study, either. Moreover, however it is that we design our mathematical function, we will want to do any computation that is needed to fit the data quickly and with a minimum of programming.

--

### What are basis functions [Ramsay & Silverman, 1997]?
### What are basis functions ? <span style="font-size:small">[Ramsay & Silverman, 1997]</span>

A set of basic functional building blocks that can be stacked on top of one another to generate curves distributed over time or space (or other continuum).
$$f(t) = a_1\phi_{1}(t)+a_2\phi_{2}(t)+\ldots+a_k\phi_k(t)$$
$f(t) = a_1\phi_{1}(t)+a_2\phi_{2}(t)+\ldots+a_k\phi_k(t)$
This is a linear combination of building blocks, where $\phi_k(t)$ is the $k$th function in our toolbox and $a_k$ is a constant. $\phi_k(t)$ is called a basis function.

note: We need, therefore, a set of basic functional building blocks that can be stacked on top of one another so as to have the features that we need. Mathematical Lego, in other words. Since this is mathematics, we use the symbol $\phi_k(t)$ to stand for the $k$th function in our toy box, and we call this the $k$th _basis function_.
Expand All @@ -76,27 +82,41 @@ In math-speak, this is a "linear combination". The construction of an actual fun

### Why basis functions?

<div style="text-align: center;"><img src="https://ieeexplore.ieee.org/mediastore_new/IEEE/content/media/10/6238368/6210369/6210369-fig-1-source-large.gif" width=40%><!--<img src="images/Illustration-of-the-double-gamma-function-and-its-partial-derivatives.png" width=40%>--><br> Illustration of the double-gamma function and its partial derivatives. [Barré et al., 2012]</div>
::: block

![[Illustration-of-the-double-gamma-function-and-its-partial-derivatives.png|380]]

Illustration of the double-gamma function and its partial derivatives. <span style="font-size:small">[Barré et al., 2012]</span>

<div style="text-align: center;"><img src="https://ieeexplore.ieee.org/mediastore_new/IEEE/content/media/10/6238368/6210369/6210369-fig-4-source-large.gif" width=40%><!--<img src="images/Modeled-haemodynamic-response-block-paradigm-noisefree-black-FIR-blue-IIR.png" width="40%">--><br> Modeled haemodynamic response (block paradigm): noisefree (black), FIR (blue), IIR (Magenta), linear regression (green), FOP (red). [Barré et al., 2012]</div>
![[Modeled-haemodynamic-response-block-paradigm-noisefree-black-FIR-blue-IIR.png|380]]

Modeled haemodynamic response (block paradigm): noisefree (black), FIR (blue), IIR (Magenta), linear regression (green), FOP (red). <span style="font-size:small">[Barré et al., 2012]</span>

::: <!-- element align="center"-->

--

### Now to model the HRF

<div style="text-align: center;">Let's say we have an HRF that looks like this:<br><img src="https://theclevermachine.files.wordpress.com/2012/12/hrfestimation-fir-representation-b.png?w=2000&h"> <!--<img src="images/hrfestimation-fir-representation-b.jpg">--></div>
<!-- slide align="center"-->

Let's say we have an HRF that looks like this:
![[hrfestimation-fir-representation-b.jpg]]

note: We would like to be able to model the HRF as a weighted combination of simple basis functions. The simplest set of basis functions is the FIR basis, which is a series of ![H](https://s0.wp.com/latex.php?latex=H&bg=ffffff&fg=4e4e4e&s=-1&c=20201002) distinct unit-magnitude (i.e. equal to one) impulses, each of which is delayed in time by ![t = 1 \dots H](https://s0.wp.com/latex.php?latex=t+%3D+1+%5Cdots+H&bg=ffffff&fg=4e4e4e&s=-1&c=20201002) TRs.

--

### Now to model the HRF
<!-- slide align="center"-->

<split left="1" right="4" gap=2>
<div style="text-align: center;"><img src="https://theclevermachine.files.wordpress.com/2012/12/hrfestimation-fir-representation-b.png?w=2000&h"><!--<img src="images/hrfestimation-fir-representation-b.jpg">--></div>
![[hrfestimation-fir-representation-b.jpg]]

<div style="text-align: center;">Let's model that HRF using Finite Impulse Respose basis set:<br><img src="https://theclevermachine.files.wordpress.com/2012/12/hrfestimation-fir-representation-a.png?w=2000&h="><!--<img src="images/hrfestimation-fir-representation-a.jpg">--></div>
::: block
Let's model that HRF using Finite Impulse Respose basis set:
![[hrfestimation-fir-representation-a.jpg]]
:::
</split>

note: Each of the basis functions ![b_t](https://s0.wp.com/latex.php?latex=b_t&bg=ffffff&fg=4e4e4e&s=0&c=20201002) has an unit impulse that occurs at time ![t = 1 \dots 20](https://s0.wp.com/latex.php?latex=t+%3D+1+%5Cdots+20&bg=ffffff&fg=4e4e4e&s=0&c=20201002); otherwise it is equal to zero. Weighting each basis function ![b_t](https://s0.wp.com/latex.php?latex=b_t&bg=ffffff&fg=4e4e4e&s=-1&c=20201002) with the corresponding value of the HRF at each time point ![t](https://s0.wp.com/latex.php?latex=t&bg=ffffff&fg=4e4e4e&s=-1&c=20201002), followed by a sum across all the functions gives the target HRF in the first plot above. The FIR basis model makes no assumptions about the shape of the  HRF–the weight applied to each basis function can take any value–which allows the model to capture a wide range of HRF profiles.
Expand All @@ -108,15 +128,16 @@ Given an experiment where various stimuli are presented to a subject and BOLD re
### Now to model the HRF

4 voxels:
<div style="text-align: center;"><img src="https://theclevermachine.files.wordpress.com/2012/12/hrfestimation-fir-simulated-voxels1.png?w=2000&h="><!--<img src="images/hrfestimation-fir-simulated-voxels1.jpg">--></div>
![[hrfestimation-fir-simulated-voxels1.jpg]]

note: For this example we revisit a simulation of voxels with 4 different types of tuning. One voxel is strongly tuned for visual stimuli (such as a light), the second voxel is weakly tuned for auditory stimuli (such as a tone), the third is moderately tuned for somatosensory stimuli (such as warmth applied to the palm), and the final voxel is unselective (i.e. weakly and equally selective for all three types of stimuli). We simulate an experiment where the blood-oxygen-level dependent (BOLD) signals evoked  in each voxel by a series of stimuli consisting of nonoverlapping lights, tones, and applications of warmth to the palm, are measured over ![T=330](https://s0.wp.com/latex.php?latex=T%3D330&bg=ffffff&fg=4e4e4e&s=-1&c=20201002) fMRI measurments (TRs).

--

### Now to model the HRF

<div style="text-align: center;"><img src="https://theclevermachine.files.wordpress.com/2012/12/hrfestimation-fir-design-matrix.png?w=2000&h="><!--<img src="images/hrfestimation-fir-design-matrix.jpg">--></div>
<!-- slide align="center"-->
![[hrfestimation-fir-design-matrix.jpg]]

note: Now let’s estimate the HRF of each voxel to each of the ![C = 3](https://s0.wp.com/latex.php?latex=C+%3D+3&bg=ffffff&fg=4e4e4e&s=-1&c=20201002) stimulus conditions using an FIR basis function model. To do so, we create a design matrix composed of successive sets of delayed impulses, where each set of impulses begins at the onset of each stimulus condition. For the ![[T \times C]](https://s0.wp.com/latex.php?latex=%5BT+%5Ctimes+C%5D&bg=ffffff&fg=4e4e4e&s=0&c=20201002)-sized stimulus onset matrix ![D](https://s0.wp.com/latex.php?latex=D&bg=ffffff&fg=4e4e4e&s=-1&c=20201002), we calculate an ![[T \times HC]](https://s0.wp.com/latex.php?latex=%5BT+%5Ctimes+HC%5D&bg=ffffff&fg=4e4e4e&s=0&c=20201002) FIR design matrix ![X_{FIR}](https://s0.wp.com/latex.php?latex=X_%7BFIR%7D&bg=ffffff&fg=4e4e4e&s=-1&c=20201002), where ![H](https://s0.wp.com/latex.php?latex=H&bg=ffffff&fg=4e4e4e&s=-1&c=20201002) is the assumed length of the HRF we are trying to estimate. The code for creating and displaying the design matrix for an assumed HRF length ![H=16](https://s0.wp.com/latex.php?latex=H%3D16&bg=ffffff&fg=4e4e4e&s=-1&c=20201002)

Expand All @@ -135,24 +156,25 @@ then we can use the Ordinary Least Squares (OLS) to solve the for ![\beta_{HRF}
--

### Now to model the HRF
<div style="text-align: center;"><img src="https://theclevermachine.files.wordpress.com/2012/12/hrfestimation-fir-hrf-estimates.png?w=2000&h="><!--<img src="images/hrfestimation-fir-hrf-estimates.jpg">--></div>
<!-- slide align="center"-->
![[hrfestimation-fir-hrf-estimates.webp]]

note: Once determined, the resulting ![[CH \times V]](https://s0.wp.com/latex.php?latex=%5BCH+%5Ctimes+V%5D&bg=ffffff&fg=4e4e4e&s=-1&c=20201002) matrix of weights ![\hat \beta_{FIR}](https://s0.wp.com/latex.php?latex=%5Chat+%5Cbeta_%7BFIR%7D&bg=ffffff&fg=4e4e4e&s=0&c=20201002) has the HRF of each of the ![V=4](https://s0.wp.com/latex.php?latex=V%3D4&bg=ffffff&fg=4e4e4e&s=-1&c=20201002) different voxels to each stimulus condition along its columns. The first ![H](https://s0.wp.com/latex.php?latex=H&bg=ffffff&fg=4e4e4e&s=-2&c=20201002) (1-16) of the weights along a column define the HRF to the first stimulus (the light). The second ![H](https://s0.wp.com/latex.php?latex=H&bg=ffffff&fg=4e4e4e&s=-1&c=20201002) (17-32) weights along a column determine the HRF to the second stimulus (the tone), etc… Below we parse out these weights and display the resulting HRFs for each voxel

Here we see that estimated HRFs accurately capture both the shape of the HRF and the selectivity of each of the voxels. For instance, the HRFs estimated from the responses of first voxel indicate strong tuning for the light stimulus. The HRF estimated for the light stimulus has an amplitude that is approximately 4 times that of the true HRF. This corresponds with the actual tuning of the voxel (compare this to the value of ![\beta(1,1)](https://s0.wp.com/latex.php?latex=%5Cbeta%281%2C1%29&bg=ffffff&fg=4e4e4e&s=0&c=20201002)). Additionally, time delay till the maximum value (time-to-peak) of the HRF to the light is the same as the true HRF. The first voxel’s HRFs estimated for the other stimuli are essentially noise around baseline. This (correctly) indicates that the first voxel has no selectivity for those stimuli. Further inspection of the remaining estimated HRFs indicate accurate tuning and HRF shape is recovered for the other three voxels as well.

--

<h3 style="margin-top: 25%;">PROS and CONS</h3>
### PROS and CONS <!-- element align="center" drag="100 10" drop="0 25"-->

--

### Tools

- rsHRF [Wu et al., 2021]
- FLOBS (FIMRIB Linear Optimal Basis Sets) [Woolrich et al., 2004]
- SPM Optimal Basis Sets [Ashburner et al., 2021]
- hrf_estimation [Pedregosa & Eickenberg, 2015]
- rsHRF <span style="font-size:small">[Wu et al., 2021]</span>
- FLOBS (FIMRIB Linear Optimal Basis Sets) <span style="font-size:small">[Woolrich et al., 2004]</span>
- SPM Optimal Basis Sets <span style="font-size:small">[Ashburner et al., 2021]</span>
- hrf_estimation <span style="font-size:small">[Pedregosa & Eickenberg, 2015]</span>

--

Expand Down