To enable a neural mass model to mimic the effect of habituation, we introduce the schematic control structure in . This structure stipulates the existence of an auxiliary, slow control variable, denoted s
, which emulates a saturation (hence the notational choice) and relaxation dynamics in the neural mass. Recalling the intuition highlighted in Sect. 2.2, the purpose of the control variable s
is thus to capture both the cumulative contribution to habituation of previous response cycles, and the dynamic forgetting
of that contribution while no new stimulus is applied. The variable s
is therefore defined as the output of a stable dynamic system that integrates the amplitude of the response to past stimuli. In turn, as a new stimulus is applied, at a time ts
, parameters of the neural mass model are adjusted as functions of s
). The class of parameter scheduling methods to which this mechanism belongs is common in nonlinear system identification (Nelles 2001
, p. 195). It builds on classical ideas of time constant separation and inertial manifold embedding that are common in models of myriad physical and chemical phenomena (Gorban et al. 2007
); in particular, these ideas are fundamental in the closely related, low order models of firing patterns in single neurons (Ghigliazza and Holmes 2004
; Shilnikov and Rulkov 2004
), and in low order models in fluid dynamics (Stuart 1971
; Tadmor et al. 2011
), the second author’s other area of interest.
Fig. 5 A generic block diagram of habituation control. The predicted habituation in the response of a neural mass model to a sensory stimulus is determined by adjustments of model parameters, as functions of the habituation control variable s. That parameter (more ...)
The following discussion is used to motivate our choice of the structure of the ordinary differential equation (ODE) governing the auxiliary state, s. The convention we use, in what follows, is that s ≥ 0, that s(ts) = 0 is associated with no habituation, and that habituation increases as s(ts) grows. To be able to serve in the role we envision, the spread of values of s(ts) needs to match the spread of levels of evoked response, observed in our data. Specific examples include:
- Following an isolated stimulus4 at the time ts, the value of s(ts + 0.5 s) should be sufficiently small to be associated with no habituation. This will match the absence of noticeable habituation in our EEG response data, when the ISI is 1 and 0.5 s (1 and 2 Hz stimulus trains).
- The value of s(ts + τ), following an isolated stimulus, should be sufficiently high to cause, at least a moderate habituation, when the ISI is τ = 0.33 s or shorter. This is dictated by the observed habituation in the response to the second stimulus, at stimulus frequencies of 3 Hz and above.
- The value of s(ts + τ), following an isolated stimulus, should be sufficiently high to prevent any discernible response, when the ISI is τ = 0.166 s or shorter. This is dictated by the lack of noticeable evoked response to the next stimulus, for stimulus frequencies of 6 Hz and above.
- Beyond these, easily stated requirements, s should be able to integrate the effect of frequent, small responses, to generate the alternating response amplitude observed in the data for the shorter range of ISI lengths.
These requirements translate into constraints on the decay time constant of a hypothesized first-order ODE, governing, the evolution of s. A detailed examination of our data revealed that they could not be matched by the response to a first order, linear time invariant ODE. Rather, the data suggest that the desired dynamics should provide for a far shorter decline time constant, at large values of s, than at small values of s. Our choice of a nonlinear system structure that meets these characteristics is
where the input, u
), is the rectified (i.e., the absolute value of the) “P1–N1–P2” portion of the EEG signal. Indeed, the effective time constant, τ
), in this formulation, decays as a function of s
(for details about how this control structure meets the qualitative behavior of habituation refer Sect. C.2). Model parameters were selected as gs
= 20 and τs
= 0.5 s. We shall motivate the guidelines for this selection shortly.
Before we turn to our next task, the determination of model parameter dependencies on values of s, we use and to illustrate the required properties of s and our considerations in making the choice of its governing dynamics.
Fig. 6 The rectified EEG signal, as predicted neural mass model (top panels) is the input to the ODE (2) that governs the evolution of the auxiliary state s, shown in the bottom panels. Results are shown for the first 2 s of stimulus trains with ISI of 0.25 (more ...)
Fig. 7 The solid curve is a 1-s trajectory of s, from Eq. 2, driven by the un-habituated EEG response to an isolated stimulus. The + marker indicates samples of s(t) at the time of arrival of subsequent stimuli. Eight data records were used, each characterized (more ...)
shows trajectories of the rectified P1–N1–P2 waveforms (top panels), and of the auxiliary state s
, (bottom panels), when the former are used as inputs in (2
), initiated at s
(0) = 0. The P1–N1–P2 waveforms are extracted from response data for stimulus trains with ISI of 0.25 s (4 Hz, left panels) and 0.125 s (8 Hz, right panels), using a best-fit procedure we elaborate in Sect. 4. In essence, during each response cycle, the parameters of the neural mass model are optimized for best fit of the predicted P1–N1–P2 waveform to the data signal.5
Having no preceding stimulus, the response to the first stimulus is large, driving the value of s
, sampled at the end of that cycle, to be the largest of the values sampled along the respective train. By the same token, the habituation effects of the first response cycle is the most pronounced, and the response to the subsequent stimulus, in each of the two cases, is the smallest. Indeed, no response is detected during the second stimulus cycle, when an ISI of 0.125 s is used (8 Hz case). Subsequent variations in the response amplitudes continue to correlate with variations of s
), the evaluations of s
) at the time of a new stimulus, at the end of each response cycle.
The data used to produce were obtained by the same procedure as in . The figure shows the continuous 1-s long trajectory of s
, following the application of the first, isolated stimulus in a 4-s stimulus train (solid/red curve). The + signs mark the sampled s
) for eight 4-s response data records, corresponding to each of the stimulus frequencies, from 1 to 8 Hz. (The dotted rectilinear lines are included as visualization aids, and are not
continuous time trajectories of s
.) This figure illustrates the requirements from a successful parameter choice in (2
): The sampled value of s
are intended to guide the adjustment of the neural mass model parameters, and thus, to determine the level of habituation in the response waveform, generated by the model, following each stimulus. The vertical separation of the sampled values, in , should therefore be ample to robustly distinguish between the corresponding response intensities, in subsequent cycles. The faster-than-exponential decay of the response to the initial stimulus (the continuous solid curve) is needed to achieve that separation. While, the selection of the gain, gs
, would have been immaterial to the desired stretch in a linear equation, the implied nonlinearity makes a successful choice dependent on both gs
and of τs
, which were found by a parameter scan.
Our next task is to determine the precise functional dependence of the neural mass model parameters on the auxiliary, saturation state, s.