fMRI BOLD Signal Response

import matplotlib.pyplot as plt
import numpy as np
from math import factorial

np.random.seed(1111)

# HRF data
time = np.linspace(0, 30, 50)
# Double gamma HRF
a1, a2, b1, b2, c = 6, 16, 1, 1, 1/6
hrf = (time**a1 * np.exp(-time/b1)) / (b1**(a1+1) * factorial(int(a1))) -       c * (time**a2 * np.exp(-time/b2)) / (b2**(a2+1) * factorial(int(a2)))
hrf = hrf / np.max(np.abs(hrf)) + np.random.normal(0, 0.08, 50)

# Smooth curve
t_smooth = np.linspace(0, 30, 200)
hrf_smooth = (t_smooth**a1 * np.exp(-t_smooth/b1)) / (b1**(a1+1) * factorial(int(a1))) -              c * (t_smooth**a2 * np.exp(-t_smooth/b2)) / (b2**(a2+1) * factorial(int(a2)))
hrf_smooth = hrf_smooth / np.max(np.abs(hrf_smooth))

fig, ax = plt.subplots(figsize=(10, 7), facecolor='#ffffff')
ax.set_facecolor('#ffffff')

ax.yaxis.grid(True, color='#f0f0f0', linewidth=1, zorder=1)
ax.xaxis.grid(True, color='#f0f0f0', linewidth=1, zorder=1)

ax.fill_between(t_smooth, hrf_smooth - 0.1, hrf_smooth + 0.1, color='#F527B0', alpha=0.12, linewidth=0, zorder=2)
ax.plot(t_smooth, hrf_smooth, color='#F527B0', linewidth=2.5, zorder=3)
ax.scatter(time, hrf, c='#276CF5', s=50, alpha=0.85, edgecolors='white', linewidths=0.5, zorder=4)

# Zero line
ax.axhline(y=0, color='#888888', linestyle='-', linewidth=0.8, alpha=0.5)

# Peak annotation
peak_idx = np.argmax(hrf_smooth)
ax.scatter([t_smooth[peak_idx]], [hrf_smooth[peak_idx]], c='#F5276C', s=100, marker='v', zorder=5)
ax.text(t_smooth[peak_idx] + 0.5, hrf_smooth[peak_idx] + 0.08, f'peak: {t_smooth[peak_idx]:.1f}s', fontsize=9, color='#F5276C')

# Undershoot
ax.annotate('undershoot', xy=(18, -0.15), fontsize=9, color='#666666')

for spine in ['top', 'right']:
    ax.spines[spine].set_visible(False)
for spine in ['bottom', 'left']:
    ax.spines[spine].set_color('#cccccc')

ax.set_xlabel('Time (seconds)', fontsize=12, color='#333333', fontweight='500', labelpad=10)
ax.set_ylabel('BOLD Signal (a.u.)', fontsize=12, color='#333333', fontweight='500', labelpad=10)
ax.set_title('Hemodynamic Response Function', fontsize=15, color='#1a1a1a', fontweight='bold', pad=20, loc='left')
ax.tick_params(colors='#666666', labelsize=10, length=0)

plt.tight_layout()
plt.show()