You could make the amplitude of the hemodynamic response function for each stimulus cycle depend on the RT. Probably the thing to do is to have 2 regressors:
(1) the normal one with equal HRF amplitude in each stim cycle
(2) one where the amplitude of the HRF in each stim cycle is modulated by (RT-mean RT).
Subtracting the mean RT from the amplitude should make the 2 regressors independent, and then the amplitude of #1 is the average activation and the amplitude of #2 is the amount correlated with changes in RT.

If you are using 'simple regression' with a fixed hemodynamic response function, then the amplitude of #1 will give the average response to the stimulus, while the amplitude of #2 will give the response to changes in the RT.

Thus, if you had stimuli at time steps 3,9,11, you could
code the constant HR amplitude response from waver using something like

waver -GAM -inline 0 0 0 1 0 0 0 0 0 1 0 1 -TR 2.5 > a1.1D

and then if you had RT of 1.2, 1.0, and 0.5 for these 3 stimuli, you would subtract the mean (0.9) and then encode

waver -GAM -inline 0 0 0 0.3 0 0 0 0 0 0.1 0 -0.4 -TR 2.5 > a2.1D

You would probably want to do a GLT test with a1.1D and a2.1D grouped together, and then look at the locations where the GLT F statistic was large (indicating some activation) and also at the same time where the a2.1D t-statistic by itself was large (indicating activation level sensitive to RT). You can perform this conjunction using 3dcalc, as usual.

If you are doing deconvolution, then similar remarks apply, but you skip the waver step and just input the 0-1 timeseries for the first stimulus file and the (RT-mean RT) time series for the second stimulus file.