System identification of the nociceptive function

Pain is essential for survival but persistent pain reduces quality of life. One clinical symptom is hyperalgesia which results in either a decrease in threshold or an increase in suprathreshold response\cite{Sandkühler}. Hyperalgesia can be a consequence of peripheral and/or central sensitization. This perturbs the nociceptive system that is then malfunctioning. Our goal is to identify the nociceptive pathway and to characterize which parts are (mal)functioning. To this end, we adopt a system identification approach by integrating input-output measurements for estimating the parameters in a mathematical model of the ascending nociceptive pathway. We consider data from a (binary) detection task with suitable variation in stimulus parameters. This is achieved by application of electrocutaneous stimulation using an intra-epidermal needle electrode. In each trial, current is delivered as a pulse train characterized by the amplitude and the temporal combination of the number of pulses, the inter-pulse interval and the pulse width. Then, the subject reports a binary response, perceived or not. From one trial to the next, the amplitude of the current is be adjusted using a staircase procedure \cite{Doll}. The set of stimulus-response pairs (SRP) is the input-output data for our identification method. In a previous study, we constructed a mathematical model of the ascending nociceptive pathway considering peripheral and central mechanisms. Analogous to the experimental paradigm, one can define the perception threshold (PT) for the model: the amplitude such that given the experimental temporal parameters and model parameters, half of the stimuli are detected \cite{Yang}. The model reproduces PT variations similar as in experiments and depends on the parameters of the model. Based on this we were able to estimate the model parameters given the PTs. We found, however, that the addition of (measurement) noise in the PTs led to large estimation errors. Indeed, a relative error of only 0.1\% in the PTs resulted in a 10\% relative error in the estimation. We suspect this is due to losing information when transforming the SRP into a set of PTs. Therefore, we aim to improve the estimation by using SRP directly. The staircase procedures can be formulated as Markov models \cite{Wetherill}. This describes the eventual distribution of the amplitude of applied currents. For the staircase procedures we consider, the Markov models are ergodic, hence, a single realization or measurement is informative if the number of trials is large enough. Based on the final asymptotic distribution and the SRP, we formulate a maximum likelihood estimator for identification of the model parameters. We show that on simulated data this direct approach with SRP outperforms the indirect method with PTs. We study the effect of the number of trials and the choice of stimulus parameters, i.e. the pulse width and the inter-pulse interval, on the estimation performance