Hi Marina,
The regression matrix should be first put together by 3dDeconvolve, and then passed to 3dREMLfit.
When running an analysis with afni_proc.py, this can be done using -regress_reml_exec (and probably -regress_3dD_stop).
To do it directly with 3dDeconvolve and 3dREMLfit, just run a normal 3dDeconvolve command, but add -x1D_stop.
It will output a 3dREMLfit script that can then be run.
- rick