Modeling the pupillometric signal#

The algorithms in this notebook, together with simulation studies showing their advantage over more traditional approaches, are described in

Mittner, M., & Groot, J. M. (2025). A novel method for modeling tonic and phasic pupil dynamics in humans. Behavior Research Methods, 57(9), 233. https://doi.org/10.3758/s13428-025-02755-7

Please cite this paper if you use pupil.estimate_baseline() or pupil.estimate_response().

[1]:

import pypillometry as pp
from pypillometry import u
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats

It can be useful to think about the pupillometric signal as being composed of different components. One common assumption, based on the finding that the pupil reflects activity in the norepinephrinergic system, is that it consists of slow, tonic (baseline) and faster, phasic (response) fluctuations.

To illustrate this, we simulate data according to such a model: event-induced responses are superimposed on a slow, baseline-like component. The response to each event is a scaled version of a “response kernel” (Hoeks & Levelt, 1993), which can be plotted with pp.signal.plot_prf():

[2]:

plt.figure(figsize=(6, 3))
pp.signal.plot_prf(npar=10.1, tmax=930);

We create 100 events with random inter-stimulus intervals between 0.6 and 1.4 seconds (a “fast” paradigm in which the responses to successive events overlap), a slowly fluctuating baseline (pp.signal.fake_pupil_baseline()), a response of random amplitude to each event (using the design matrix from pp.signal.pupil_build_design_matrix()) and some measurement noise:

[3]:

rng = np.random.default_rng(12)
fs = 50       # sampling rate (Hz)
ntrials = 100

event_onsets = 2000 + np.cumsum(rng.uniform(600, 1400, ntrials))  # in ms
duration = event_onsets[-1] + 5000                                   # in ms

# time (ms) and baseline
tx, real_baseline, _ = pp.signal.fake_pupil_baseline(duration/1000, fs, mean=1000, amplitude=100,
                                              freq=0.1, seed=12)
X = pp.signal.pupil_build_design_matrix(tx, event_onsets, fs, npar=10.1, tmax=930)
real_response_coef = rng.uniform(50, 150, ntrials)
real_response = real_response_coef @ X
noise = rng.normal(0, 5, len(tx))

faked = pp.EyeData(left_pupil=real_baseline + real_response + noise, time=tx,
                   event_onsets=event_onsets, event_labels=["event"]*ntrials,
                   name="simulated")

plt.figure(figsize=(15, 5))
faked.plot.plot_pupil(show_events=False)
plt.plot(tx/1000, real_baseline, label="baseline")
plt.vlines(event_onsets/1000, *plt.ylim(), color="grey", alpha=0.5)
plt.legend();

The orange line is the tonic, baseline component. At each event (grey lines), a scaled response kernel is added. The challenge is to extract both the size of the response and the baseline value at each event when only using the measured pupil data (the blue line).

Traditional approaches for disentangling tonic and phasic components#

One common way to analyse pupillometric data on the trial-by-trial level is to extract the average signal just before an event (as a measure of the baseline) and some time after the event (as a measure of the pupil’s response). In pypillometry, this is implemented by stat_per_event(), which summarizes the signal in a time window around each event.

For example, the following code extracts

the average signal from 200 ms before each event until the event (0) as baseline
the average signal from 800 ms until 1200 ms after each event as a measure of the response

The baseline is then subtracted from the response:

[4]:

baseline = faked.stat_per_event(event_select="event", window=(-200, 0) * u.ms)
response = faked.stat_per_event(event_select="event", window=(800, 1200) * u.ms)
baseline

[4]:

{'left_pupil': array([ 755.3445703 ,  910.5008777 ,  932.7298659 , 1032.40897831,
        1064.03141911, 1056.86661384, 1091.24004593, 1054.12046947,
        1070.96222088, 1139.40891598, 1113.25888123, 1195.99348093,
        1199.65930305, 1300.49149841, 1316.62737324, 1282.12409623,
        1246.72110993, 1136.01067965, 1059.99835976, 1044.56272194,
         956.21648324,  859.03635638,  903.51045264,  873.00710252,
         876.18155969,  932.69087178,  966.2978947 ,  967.53196025,
        1042.23556955, 1000.24557847,  958.35198931,  921.75006445,
         990.53000669,  951.96188817,  994.59993086,  997.73060282,
         938.00568988,  977.57391566, 1021.453641  , 1081.30504165,
        1094.97402141, 1124.03325259, 1108.26405136, 1074.26110607,
        1018.94383666,  949.14281559,  968.54194546,  869.84203316,
         938.66220767, 1070.52615893, 1098.13718577, 1177.17278248,
        1180.59678379, 1114.84777511, 1060.63578694, 1005.37472371,
         980.93643544, 1016.52835625, 1028.4051062 , 1051.84671225,
        1142.5630315 , 1166.94380007, 1234.42000908, 1202.93942607,
        1164.63545442, 1153.9282239 , 1118.07155382, 1140.35593043,
        1219.56578743, 1148.34374671, 1237.91395018, 1177.00180163,
        1044.6882333 ,  956.96898122,  935.12733253,  969.56222201,
        1020.253291  , 1013.06523706, 1091.89257587, 1113.22617499,
        1156.46164584, 1189.72445557, 1153.85375071, 1144.32535071,
        1078.47147651, 1110.182412  , 1058.94757395, 1006.90644063,
        1025.38538538, 1005.60095871,  992.93209897,  991.09570187,
         951.00503493, 1040.14332618,  977.77872015, 1004.7527711 ,
        1038.7322798 , 1015.5552674 , 1097.69076972, 1020.27321193])}

The result contains one array per eye and variable, with one value per event. We subtract the baseline from the response:

[5]:

baseline = baseline["left_pupil"]
response = response["left_pupil"] - baseline

stat_per_event() supports selecting specific events (event_select=) and any summary function (statfct=, default is numpy.mean). Masked data (e.g., blinks) is ignored when computing the statistic, so the function has to handle masked arrays (e.g., np.ma.median instead of np.median):

[6]:

faked.stat_per_event(event_select="event", window=(-200, 0) * u.ms, statfct=np.ma.median)

[6]:

{'left_pupil': array([ 755.92091005,  907.89529851,  930.73411829, 1031.65582806,
        1063.32439519, 1056.42293438, 1090.89687451, 1050.88269661,
        1069.58743809, 1139.02064069, 1111.96922201, 1195.85132123,
        1200.58884988, 1304.19762758, 1318.58902327, 1281.16290544,
        1247.59004759, 1133.00943923, 1059.81610238, 1045.0975211 ,
         957.29146249,  859.37766491,  903.91447365,  872.71620589,
         875.11773972,  935.05155369,  966.75632829,  968.71623455,
        1042.46240916, 1001.13380517,  957.60786568,  922.57282414,
         990.97282501,  951.78999637,  992.8116932 ,  994.68552031,
         936.65360014,  976.94021008, 1021.28482066, 1081.64026888,
        1097.19348532, 1123.37870692, 1109.40938278, 1073.47117931,
        1020.16660628,  949.9561078 ,  970.5497937 ,  870.9889643 ,
         938.98842965, 1070.78287978, 1098.14762076, 1179.25007842,
        1180.88396895, 1116.67619587, 1058.05325367, 1006.53086871,
         981.29214479, 1016.15942405, 1026.96828037, 1053.22477922,
        1142.1563809 , 1166.98151226, 1235.29960898, 1203.72369645,
        1168.51365154, 1154.16128504, 1116.87413336, 1138.41680091,
        1216.94732524, 1146.40525942, 1237.81622373, 1179.40856231,
        1047.62349375,  957.71144323,  934.87009215,  973.40229597,
        1019.21802694, 1009.47026062, 1087.92048071, 1114.74169578,
        1153.38104173, 1193.25887528, 1156.48598797, 1144.41650203,
        1075.88090747, 1111.50299268, 1061.46855166, 1004.52165797,
        1023.30285072, 1007.96107893,  992.75167008,  990.87522947,
         949.595824  , 1040.70773824,  977.13194747, 1006.10751207,
        1042.4995676 , 1013.38451585, 1097.11168414, 1022.60628766])}

Advanced methods for tonic/phasic component estimation#

The simple method above is appealing for its simplicity but has severe limitations. Most importantly, multiple overlapping pupillary responses can “look like” baseline fluctuations when added together, thereby artificially inflating baseline estimates, particularly when events are spaced closely in time (“fast paradigms”). For that reason, we developed specialized algorithms to disentangle tonic and phasic components of the pupillometric signal.

The algorithm uses an iterative procedure to remove an initial estimate of the responses from the signal to continue to estimate the underlying baseline. Details about how this algorithm works and which parameters it supports are available in Mittner & Groot (2025).

The functionality is implemented in pupil.estimate_baseline() and pupil.estimate_response(). The response estimation depends on the estimated baseline, hence pupil.estimate_baseline() should be called first. Both functions are slow for high sampling rates, so it is a good idea to low-pass filter and downsample the data first (here, the simulated data already has a sampling rate of 50 Hz):

[7]:

d = faked | pp.pupil.estimate_baseline() | pp.pupil.estimate_response()

pp: 14:34:37 | INFO     | pupil_response:598 | optimizing both npar and tmax, might take a while...
pp: 14:34:38 | INFO     | log_progress:601 | npar=9.25, tmax=992.01
pp: 14:34:39 | INFO     | log_progress:601 | npar=7.68, tmax=976.09
pp: 14:34:39 | INFO     | log_progress:601 | npar=6.66, tmax=959.13
pp: 14:34:39 | INFO     | log_progress:601 | npar=6.89, tmax=955.95
pp: 14:34:39 | INFO     | log_progress:601 | npar=6.91, tmax=953.34
pp: 14:34:40 | INFO     | log_progress:601 | npar=6.90, tmax=952.99
pp: 14:34:40 | INFO     | log_progress:601 | npar=6.90, tmax=952.96
pp: 14:34:40 | INFO     | log_progress:601 | npar=6.90, tmax=952.98
pp: 14:34:41 | INFO     | log_progress:601 | npar=6.90, tmax=952.98

Here, we allowed the two shape parameters of the response kernel, npar and tmax, to vary freely together with the amplitudes of the responses. This allows an individualized shape of the pupillary response for each subject (which appears reasonable given the results in Hoeks & Levelt, 1993) but may take a long time to optimize and can result in pathological solutions. In that case, one or both of the parameters can be fixed, for example to reasonable group-level values (e.g., pupil.estimate_response(npar=10.1, tmax=930)).

After running these methods, the baseline is stored as a new variable baseline and the predicted response as variable response:

[8]:

d.variables

[8]:

['pupil', 'baseline', 'response']

The estimated parameters of the response are stored in d.params["response"]:

[9]:

d.params["response"]["left_pupil"]

[9]:

{'npar': np.float64(6.898073863824489),
 'npar_free': True,
 'tmax': np.float64(952.9813212623851),
 'tmax_free': True,
 'coef': array([183.28911879,  84.04643122,  77.54227549,  53.06493822,
         80.82517569,  98.19065567,  66.771559  ,  60.47232813,
        150.46244084, 156.96295404,  42.39546353, 110.77393718,
         86.08681697,  61.37768494,  68.38481793, 108.44439133,
        127.15683717,  76.50062899, 117.40338136,  73.73638378,
         83.78826984, 161.86545469,  95.14261821, 106.27273525,
         88.61519282,  94.76137818, 112.60369117, 153.57435172,
         83.21178757,  81.47609056, 107.79893068, 183.3302553 ,
         66.75597556,  85.36239623,  73.62776758,  96.02169245,
        120.18382614, 130.98046227,  60.55223607,  84.67406329,
         76.09955041,  54.48380836,  72.6847929 , 155.25686911,
        141.13478619, 163.9118957 ,  48.92529894,  80.84413828,
        157.44370732, 157.71781099, 206.24919954, 126.34508741,
        121.80487169, 167.62832281, 187.02334815, 191.95811714,
        118.65565098,  99.18828573,  83.30767876,  83.4939494 ,
         73.41970949,  97.41950475,  75.12146785,  84.91750688,
         85.60071323,  48.83913909, 106.40510061, 143.88505571,
         32.72912238,  85.63759271,  97.54719046,  86.56839365,
        102.06446948,  58.10358953,  99.21977287, 126.88673397,
         90.02978459,  68.66468092,  39.91344642, 100.61478132,
         86.54381202, 118.91362733, 119.97832558, 153.76445372,
        177.40643497, 120.67078531,  84.15320822, 106.65904961,
         87.83733635,  82.1974781 , 144.08721123, 149.99290614,
        157.01548329, 100.83537451,   0.        , 137.45475764,
        157.59419856, 161.5473934 ,  84.83611519, 186.01864865]),
 'bounds': {'npar': (1, 20), 'tmax': (100, 2000)},
 'valid': array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]),
 'observed': array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]),
 'min_valid': 0.5}

We can plot the estimated baseline and the full model (baseline+response) together with the data and the real baseline:

[10]:

plt.figure(figsize=(15, 5))
plt.plot(d.tx/1000, d["left", "pupil"], label="data")
plt.plot(d.tx/1000, real_baseline, label="real baseline")
plt.plot(d.tx/1000, d["left", "baseline"], label="estimated baseline")
plt.plot(d.tx/1000, d["left", "baseline"] + d["left", "response"], label="model (baseline+response)")
plt.xlabel("seconds")
plt.legend();

The estimated baseline is more wiggly than the real one but gives a considerably better estimate than the traditional method. The overall fit of the model is excellent (which is to be expected for simulated data).

We can quantify how well the baseline estimation works relative to the traditional method by comparing both to the “ground truth” which is available for simulated data. We get the real and estimated baseline at each event onset with stat_per_event():

[11]:

# store the real baseline as a variable so that we can use stat_per_event()
d.data["left_realbaseline"] = real_baseline

at_onset = (0, 1000/fs) * u.ms   # one sample at the event onset
real_baseline_ev = d.stat_per_event(event_select="event", window=at_onset)["left_realbaseline"]
novel_baseline = d.stat_per_event(event_select="event", window=at_onset)["left_baseline"]
novel_response = d.params["response"]["left_pupil"]["coef"]

traditional_baseline = baseline
traditional_response = response

and compare them by means of the correlation between the estimated and the real values:

[12]:

print("Traditional method:")
print("Baseline: Corr(trad, real) =", scipy.stats.pearsonr(traditional_baseline, real_baseline_ev)[0])
print("Response: Corr(trad, real) =", scipy.stats.pearsonr(traditional_response, real_response_coef)[0])
print("\nNovel method:")
print("Baseline: Corr(nov, real) =", scipy.stats.pearsonr(novel_baseline, real_baseline_ev)[0])
print("Response: Corr(nov, real) =", scipy.stats.pearsonr(novel_response, real_response_coef)[0])

Traditional method:
Baseline: Corr(trad, real) = 0.9552897127692387
Response: Corr(trad, real) = 0.4355020295417933

Novel method:
Baseline: Corr(nov, real) = 0.9695993602751325
Response: Corr(nov, real) = 0.6044710168795537

The correlations are higher for the novel method than for the traditional method. More detailed simulation studies are reported in Mittner & Groot (2025).

The parameters of the baseline estimation are described in the API documentation of pp.signal.baseline_envelope_iter_bspline() and, in more detail, in Mittner & Groot (2025).

The parameters of the response estimation are described in the API documentation of pp.signal.pupil_response().