Modeling and Parameter Analysis of Basic Single Channel Neuron Mass Model for SSVEP

Abstract

While steady-state visual evoked potentials (SSVEPs) are widely used in brain–computer interfaces (BCIs) due to their robustness to rhythmic visual stimuli, their generation mechanisms remain poorly understood. Challenges such as experimental complexity, inter-subject variability, and limited physiological interpretability hinder the development of efficient BCI systems. This study employed a single-channel neural mass model (NMM) of V1 cortical dynamics to investigate the biophysical underpinnings of SSVEP generation. By systematically varying synaptic gain, time constants, and external input parameters, we simulated δ/α/γ band oscillations and analyzed their generation principles. The model demonstrates that synaptic gain controls oscillation amplitude and harmonic content, and time constants determine signal decay kinetics and frequency precision, while input variance modulates harmonic stability. Our results reveal how V1 circuitry generates frequency-locked SSVEP responses through excitatory–inhibitory interactions and dynamic filtering mechanisms. This computational framework successfully reproduces fundamental SSVEP characteristics without requiring multi-subject experimental data, offering new insights into the physiological basis of SSVEP-based brain–computer interfaces.

1. Introduction

The brain is an extremely complex organ containing billions of interconnected neurons. To understand and analyze brain activity, one can abstract and simplify the brain through mathematical models by accurately expressing the interactions and information transfer between neurons. Theoretically, different brain electrical activities can be simulated by adjusting the parameters of the model to mimic the brain.

The neuronal models widely used to model electroencephalograph (EEG) activity can be divided into two main categories. The first one focuses on the micro level, including Hodgkin–Huxley [1], Chay [2], Traub [3], and so on, which can be used to understand the physiological mechanisms of the brain from the perspective of neurons. However, the brain is a highly complex system comprising hundreds of billions of neurons and hundreds of trillions of synapses. These neurons connect to form a vast network; therefore, it is very difficult to capture such a complex network structure and dynamic characteristics when modeling at the neuron level.

The second model type treats a mass of neurons in the brain as a single unit whose activities are described by a set of aggregated parameters, i.e., the lumped parameter model. The neural mass model (NMM) is based on neural populations and uses the lumped state variable as the average approximate behavior of a specific neuronal mass [4,5]. Rather than modeling individual microscopic neurons, it captures the overall characteristics of the neural mass composed of specific cell types. The NMM also accounts for interactions between different neural masses, enabling the simulation of large-scale neural network dynamics at a macro level. Wilson and Cowan explored the excitatory and inhibitory interactions within local neuronal populations, deriving coupled nonlinear differential equations that describe these interactions through excitatory and inhibitory neuronal models [6]. Silva et al. introduced a population-based neuronal structure for simulating α rhythmic EEG signals consisting of interconnected excitatory and inhibitory neuronal populations, where each subpopulation is modeled by a first-order differential equation containing a static nonlinear function [7,8]. This model is physiologically grounded and is widely used as a basic single-channel neuron mass model. Jansen et al. further refined the NMM, demonstrating that a mathematical model of neuronal mass could effectively simulate EEG activity and generate visual evoked potential (VEP) waveforms. Moreover, they also found that the ability to simulate VEP activity was primarily related to excitatory connections within the cortex [9].

Steady-state visual evoked potentials (SSVEPs) refer to the EEG electrical signals evoked by the visual stimuli of specific frequencies. When the retina receives visual stimuli ranging from 3.5 Hz to 75 Hz, the brain generates electrical activity at the same frequency or harmonics. This neural response arises from the synchronization of cortical neurons in the primary visual cortex (V1), where periodic visual inputs drive γ-band oscillations (30–80 Hz) through thalamocortical pathways [10]. The resulting SSVEP amplitude is modulated by attentional mechanisms as subjects maintain a focused gaze on the stimuli [11].

Brain–computer interfaces (BCIs) based on SSVEP achieve high information transfer rates by monitoring subjects’ responses to light stimuli at multiple frequencies. The high signal-to-noise ratio (SNR) of SSVEP signals in the frequency domain enables robust detection of user intent, significantly improving BCI accuracy. These systems typically do not require complex cognitive tasks by the subjects, only the maintenance of attention to visual stimuli at specific frequencies. This attentional focus enhances neural entrainment in the visual cortex, further amplifying the SSVEP response through top-down modulation from the parietal and frontal regions [12]. Hence, the fatigue and cognitive burden of the subjects are minimized. Compared to other EEG signals, SSVEPs are relatively stable and reliable, and their system response speed is fast, making them suitable for scenarios requiring rapid feedback and real-time applications. For instance, the ability to encode multiple commands via distinct stimulus frequencies (e.g., 10 Hz, 12 Hz, 15 Hz) allows SSVEP-BCIs to achieve high-dimensional control without prolonged user training. This enhances the stability and availability of BCI systems. Consequently, SSVEPs have been widely used in the development of various brain-controlled intelligent devices [13], such as brain-controlled robotic arms [14], brain-controlled quadcopters [15], brain-controlled wall-climbing robots [16], brain-controlled prosthetics [17], brain-controlled wheelchairs [18], brain-controlled smart homes [19], and other systems [20,21].

However, despite the great potential of SSVEP-BCIs, several challenges remain, including experimental difficulties and the limited availability of subjects. Firstly, the experimental difficulties refer to the need to carry out complex experimental setups and extensive signal processing. In the SSVEP-BCI, subjects are often required to look at visual stimuli with different frequencies. Accurately presenting the stimuli and acquiring the signals during this process requires designing and executing the designated experiments within a highly controllable laboratory environment. This complexity may lead to technical and operational challenges in the experimental process, further affecting the advancement of research. Secondly, recruiting sufficient and diverse subjects is essential to validating and promoting the SSVEP-BCI technology in practical applications. However, the recruitment of subjects is limited by many economic factors such as time, money, and manpower. The types of subjects are also limited by age, sex, and health. This restricts the breadth and reliability of SSVEP-BCI research.

The study of a neural mass model by Jansen et al. [9] found that the NMM can be used to describe the interaction between neuron populations and is suitable for explaining neural activity in the cerebral cortex. In addition, by simulating the collective behaviors of neural populations, the NMM can pick up oscillations with specific frequencies produced in response to visual stimuli. This understanding of the neural mass provides an effective tool for more accurately understanding and modeling the neural activity of the brain when perceiving steady-state visual evoked potentials. The main contributions of this paper are as follows:

• (1). We explore the mechanism for generating SSVEPs based on the NMM, which helps to construct the SSVEP by artificially generating virtual signals similar to actual EEGs. As a result, it does not need an experimental setup, subject recruitment, and equipment configuration, and the simulation modeling method provides a more cost-effective avenue for SSVEP-BCI research, which will improve repeatability.

• (2). The principle of the neural mass model for single-channel SSVEP is analyzed in detail, where the three basic dynamic waveforms are generated by the model, and the effects of key parameters are analyzed on the simulated signals.

• (3). We provide current study limitations and further research directions, which help to clarify the understanding of the physiological significance and functionality of the current NMM.

The remainder of the paper is organized as follows: Section 2 provides the basic neural mass model employed in the construction of SSVEP and a detailed illustration of the structure of the model. The specific simulation results are provided in Section 3, complete with a thorough analysis of each key parameter’s influence. Section 4 provides the conclusion to the study, including a discussion of the limitations of the current study and future research directions.

2. Neural Mass Model

Human functions such as thinking, perception, emotion, memory, and movement are mediated by electrical signal transmission and interaction between neurons. A neuron is the basic constituent unit of the neural system, typically comprising structures such as the cell body, dendrites, and axon. These components work in concert to receive, integrate, and transmit neural signals. Neurons communicate with each other efficiently and complexly via specialized connection structures known as synapses. Within the neural system, a neural mass is a functional grouping of neurons connected through synaptic connections. The neural masses usually work together to perform a specific neural activity or process specific information. There are two main types of cells in a neural mass: pyramidal cells and non-pyramidal cells. Pyramidal cells’ shapes are usually triangular or conical. Their characteristic synaptic morphology allows them to play an important role in transmission and long-range connection. Pyramidal cells predominantly release excitatory neurotransmitters, which activate downstream neurons. Non-pyramidal cells focus more on local neural modulation, which includes many types of neurons. Some of them are excitatory, releasing excitatory neurotransmitters that activate neighboring neurons, thereby enhancing signal transmission, while others are inhibitory, releasing inhibitory neurotransmitters that suppress the activity of neighboring neurons. These inhibitory cells help regulate and maintain the balance and stability of the neural network. The interactions between excitatory and inhibitory processes form the basis of information processing within neuron masses. Through the collaborative work of both cell types, the neural system can achieve complex and precise signal transmission, thereby enabling the execution of different brain functions [6].

2.1. Single-Channel Basic Neuron Mass Model

In the single-channel basic neuron mass model described by Jansen et al. [9], the pyramidal cell subpopulation receives excitatory and inhibitory inputs from interneuron subpopulations within the same region. Simultaneously, the outputs of this pyramidal cell subpopulation provide feedback to the excitatory and inhibitory interneuron subpopulations. The excitatory interneuron subpopulations also receive input from other regions or external stimuli, as shown in Figure 1.

In the single-channel basic neuron mass model, Gaussian white noise n(t)∼N(μ, σ²) is used as the excitatory input to represent all external stimuli and signals in the uncertain region, while the inhibitory connections between these signals are ignored when modeling. As shown in Figure 2, C₁ and C₂ represent the connectivity coefficients between the pyramidal cell subpopulation and excitatory interneuron subpopulation, reflecting the proportional relationship between the average numbers of synapses in the two subpopulations. The connection between the pyramidal cell subpopulation and the inhibitory interneuron subpopulation is represented by two additional connectivity coefficients, C₃ and C₄, which correspondingly reflect the proportional relationship between the pyramidal cell subpopulation and the excitatory interneuron subpopulation [22].

As shown in Figure 2, single-channel basic neuron mass modeling includes a feedback mechanism, physiological parameters, and two information transformation processes. The concrete modeling method is as follows:

(1) Static nonlinear function. The average membrane voltage of the neuron mass is converted into the average pulse density of the action potential. This information transformation is represented by the following instantaneous S-function:

(1)$S(v)={\displaystyle \frac{2e_0}{1+e^{r\left(v_0-v\right)}}}$

where S(v) is a static nonlinear function representing the average pulse density of the action potential, v represents the pre-synaptic average membrane voltage, e₀ represents the maximum firing rate of the neuron mass, e is the base value of the natural logarithm function, r represents the curvature of the S-function, and v₀ is the postsynaptic membrane voltage corresponding to the firing rate e₀.

Figure 3 illustrates the trend of the mean pulse density as a function of the pre-synaptic average membrane voltage. When the pre-synaptic mean membrane voltage v < −5 mV, the mean pulse density S(v) = 0. When v > −5 mV, the mean pulse density S(v) increases gradually with the increase in the pre-synaptic mean membrane voltage. The growth rate increases and then decreases, reaching its maximum between 5 mV and 10 mV. Beyond 10 mV, the growth rate gradually slows down until the maximum average pulse density is reached, i.e., 2e₀ = 5, and then the static nonlinearity reaches the saturation state, and the average pulse density no longer changes with the increase in the pre-synaptic average membrane voltage. There is a nonlinear mapping relationship between the mean pulse density and the mean pre-synaptic membrane voltage; the S-function presents a nonlinear feature, which can verify the nonlinearity of the EEG signal.

(2) Dynamic linear function. The average pulse density of the pre-synaptic membrane is converted into excitatory or inhibitory postsynaptic membrane voltage, which occurs at the synapses of neurons. Equations (2) and (3) show the impulse responses of the linear function to excitability and inhibition, respectively.

(2)$h_e(t)=u(t){\eta }_{e}t{e}^{-t{\omega }_e}$

(3)$h_i(t)=u(t){\eta }_{i}t{e}^{-t{\omega }_i}$

In Equation (2), h_e(t) represents the excitatory postsynaptic membrane impulse response; u(t) represents the Heaviside function, which is used to describe the mutation or step change in the system behavior; ω_e represents the reciprocal of the mean time constant, which is used to describe the influence of the mean time constant of the excitatory membrane and the mean distribution delay of dendrites; and η_e is the product of the excitatory mean synaptic gain G_e and ω_e.

In Equation (3), h_i(t) represents the inhibitory postsynaptic membrane impulse response; ω_i represents the reciprocal of the mean time constant, which is used to describe the influence of the mean time constant of the inhibitory membrane and the mean distribution delay of dendrites; and η_i is the product of the inhibitory mean synaptic gain G_i and ω_i.

In the following section, we will investigate the effects of the average synaptic gain G_e and the reciprocal of the average time constant ω_e on the dynamic linear function.

The variation curve of h (i.e., the impulse response of the linear function) with average synaptic gain G_e is shown in Figure 4, where the average time constant ω_e = 0.01 s⁻¹, and the values of G_e are set to 3.25 mV, 4 mV, 7 mV, and 10 mV. For the same G_e, the value of h increases with time. h reaches its maximum value when the time is ω_e and then gradually decreases to zero. As the average synaptic gain G_e increases, the peak value of the impulse response function h also increases during the same time t; that is, the maximum value of the postsynaptic membrane voltage increases.

The variation curve of h with the average time constant ω_e is shown in Figure 5, where average synaptic gain G_e = 3.25 mV, and the values of ω_e are set to 20.00, 92.60, and 218.00 s⁻¹. For the same ω_e, the amplitude increases with time until ω_e, at which point h reaches its maximum value. After that, h gradually decreases with time until it reaches zero. If the linear function has different ω_e⁻¹ values, the time to reach the maximum value will change accordingly. With the increase in ω_e, the times will shift backward, and the linear function will reach the maximum value at ω_e. Thus, adjusting ω_e can regulate the sensitivity of the synaptic gain, and the delay time can be adjusted by this linear function.

The dynamic linear transformation function h(t) can be further represented by a first-order ordinary differential equation, as shown in Equation (4).

(4)$\left\{\begin{array}{l}\overline{z}(t)=z_1(t)\\ {\overline{z}}_1(t)=\eta x(t)-2\omega z_1(t)-{\omega }^{2}z(t)\end{array}\right.$

where x(t) and z(t) are the input and output signals of the dynamic linear function, respectively, z₁(t) is the first order derivative of z(t), ω represents the reciprocal of the mean time constant, η represents the product of the mean synaptic gain G_e and ω, and the excitability or inhibition of η and ω is determined by the properties of the interneuron. Equation (4) describes the potential conversion completed on the synapses, through which the input and output of the linear transformation process can be more intuitively understood, and the information transformation process of the dynamic linear transformation can be more clearly represented.

The input $x(t)$ in Equation (4) varies across subpopulations, as shown in Figure 2. For excitatory interneurons, the input $x_{excitatory}(t)$ integrates two critical components that drive the model’s dynamics. First, the external stimulus $n(t)$ is modeled as Gaussian white noise. Second, feedback signals from pyramidal cells are incorporated through the term $C_{1}S(y_0(t))$, where $y_0(t)$ represents the output potential of pyramidal cells.

(5)$x_{excitatory}(t)=n(t)+C_{2}S(C_1{y}_0(t))$

The resultant postsynaptic potential $y_1(t)$ reflects the excitatory subpopulation’s contribution to the network’s oscillatory behavior, directly influencing the generation of SSVEP’s fundamental frequency.

(6)$\left\{\begin{array}{l}{\overline{y}}_1(t)=y_4(t)\\ {\overline{y}}_4(t)={\eta }_e{x}_{excitatory}(t)-2{\omega }_e{y}_4(t)-{{\omega }_e}^2{y}_1(t)\end{array}\right.$

The input $x_{inhibitory}(t)$ in the inhibitory interneuron subpopulation is exclusively derived from feedback signals originating in pyramidal cells.

(7)$x_{inhibitory}(t)=C_{4}S\left(C_3{y}_0(t)\right)$

This input is processed through the dynamic linear response (Equation (4)), generating the inhibitory postsynaptic potential $y_2(t)$, governed by the following system of equations:

(8)$\left\{\begin{array}{l}{\overline{y}}_2(t)=y_5(t)\\ {\overline{y}}_5(t)={\eta }_i{x}_{inhibitory}(t)-2{\omega }_i{y}_5(t)-{{\omega }_i}^2{y}_2(t)\end{array}\right.$

Pyramidal cells receive input $x_{pyramidal}(t)$ from the difference between the postsynaptic excitatory and inhibitory subpopulations.

(9)$x_{pyramidal}(t)=S(y_1(t)-y_2(t))$

Here, $y_1(t)-y_2(t)$ reflects the dynamic balance between excitation and inhibition, which is the output of the NMM. The dynamics of pyramidal cells are described by the following:

(10)$\left\{\begin{array}{l}{\overline{y}}_0(t)=y_3(t)\\ {\overline{y}}_3(t)={\eta }_e{x}_{pyramidal}(t)-2{\omega }_e{y}_3(t)-{{\omega }_e}^2{y}_0(t)\end{array}\right.$

In Equations (5)–(10), y₃(t), y₄(t), and y₅(t) are first derivatives of y₀(t), y₁(t), and y₂(t), respectively; ω_e is the reciprocal of the average time constant, which is used to describe the effect of the mean time constant of the excitatory cell membrane and the mean distribution delay of dendrites; η_e = G_e * ω_e is the product of the excitatory mean synaptic gain and the reciprocal of the mean time constant; ω_i is the reciprocal of the average time constant, which is used to describe the effect of the mean time constant of the inhibitory cell membrane and the mean distribution delay of the dendrites; η_i = G_i * ω_i is the product of the inhibitory mean synaptic gain and the reciprocal of the mean time constant; n(t) represents signals and external stimuli from all variable regions; and C₁, C₂, C₃, and C₄ represent the connectivity coefficients between the pyramidal cell subpopulation and the interneuron subpopulation.

2.2. Parameters of a Single-Channel Basic Neuron Mass Model

To explore the effect of variable model parameters on outputs, it is necessary to identify which parameters in the single-channel basic neuron mass model are related and which parameters have standard values.

The relationships among C₁–C₄ have been investigated in several studies. First, C₁ and C₃ represent the number of synapses on the dendrites from the pyramidal cell mass to the excitatory and inhibitory interneuron masses, respectively. C₂ and C₄ are positively correlated with the number of synapses on the dendrites from the excitatory and inhibitory interneuron masses to the pyramidal cell mass, respectively. Typically, C₂ includes the synapses of a thalamic origin, where the researcher generally indicates the count as C₂’.

Braitenberg and Schüz studied the pyramidal cells in the visual cortex of mice and found that C₁ + C₃ = C₂ + C₂’ + C₄ [23]. White observed that synapses produced in the spine are located in excitatory cells; however, synapses produced in the axis may also be located in inhibitory cells [24]. Elhanany and White found that in the somatic motor cortex of mice, 87% of the synapses are in the spine, and the remaining 13% are in the axis [25], indicating that about 6.5% of the synapses formed by pyramidal cells are inhibitory; thus, we have C₃/(C₁ + C₃) = 6.5/100. Later, Liu et al. found that 80% of synapses on the dendrites of pyramidal cells in cats’ motor cortex are excitatory [26]; hence, (C₂ + C₂’)/[(C₂ + C₂’) + C₄] = 0.8.

Larkman’s study showed that most excitatory cells in the visual cortex are pyramidal [27]; that is, most cells in the excitatory feedback circuit are pyramidal. If the cells in the circuit are consistent with synapses, the pyramidal cells have the same number of synapses as other cells (such as non-pyramidal cells, astrocytes, or panariocytes) in the circuit. The number of synapses from the pyramidal cell to the excitatory interneuron mass should be the same as the number of synapses from the excitatory interneuron mass to the pyramidal cell, i.e., C₁ = C₂ + C₂’.

According to [28,29], about 20% of the asymmetric synapses in the fourth layer of the cortex are generated by thalamocortical terminals; thus, C₂’/(C₂ + C₂’) = 0.2, i.e., C₂’ = C₂/4.

The transformation of the above formulas gives us C₂/C₁ = 0.8 and C₃ = C₄. However, heterogeneity in experimental subjects (e.g., age, neurological status), inconsistent criteria for pyramidal cell selection (morphological vs. electrophysiological markers), and diversity in imaging techniques and analytical pipelines collectively contribute to synaptic count variability, complicating the establishment of universal synaptic scaling relationships in excitatory and inhibitory feedback loops. For example, we can obtain C₁ = 14.4C₃ because C₃/(C₁ + C₃) = 6.5/100; however, taking C₁ = C₂ + C₂’ and C₃ = C₄ into (C₂ + C₂’)/[(C₂ + C₂’) + C₄] = 0.8 will give C₁ = 4C₃. Given these differences, C₁ = 4C₃ is selected here, where we assume that the relationship is more suitable for the human user when producing the interested signals; hence, we get C₁ = C₂/0.8 = 4C₃ = 4C₄.

Indeed, the functions of the synapses are affected by many factors, such as number, structure, and characteristics. Changing variable C can affect these factors, further affecting the occurrence and evolution of synaptic phenomena. The flexibility of variable C enables it to adapt to diverse physiological environments and produce significant effects under different conditions, especially in regulating synaptic function. Here, when C = C₁, then C₂ = 0.8C, C₃ = 0.25C, and C₄ = 0.25C.

From the dynamic linear transformation function, we can see that the amplitude of the impulse response of the linear function is proportional to the excitatory average synaptic gain G_e and the inhibitory average synaptic gain G_i. Rotterdam et al. proposed G_e = 3.25 mV and G_i = 22 mV [30]; however, Dodt et al. found that some neuropeptides may change the amplitude of the postsynaptic membrane voltage [31], and thus G_e and G_i cannot be strictly limited. Since ω_e and ω_i represent the reciprocals of the excitability mean time constant and the inhibition mean time constant, respectively, ω_e and ω_i are inversely proportional to the time required for the postsynaptic membrane voltage to reach its maximum value. Because the average synaptic gain and the average time constant do not change easily over short periods, they can be set to fixed values [32]; therefore, ω_e = 100 s⁻¹ and ω_i = 50 s⁻¹ in this study.

Many substances can affect the excitability of neurons [33,34]; for example, the change in v₀ in the nonlinear function will affect the ignition threshold. Thus, it is necessary to determine the parameters involved in the nonlinear function. In most cases, v₀ = 6 mV, e₀ = 2.5 s⁻¹, and r = 0.56 mV⁻¹, which we adopted in our current study.

3. Results

According to the single-channel basic neuron mass model, the EEG signals with δ, α, and γ rhythms can be output by adjusting different physiological parameters in the model. To simulate the single-channel basic neuron mass model in this study, the Runge–Kutta solver in MATLAB/Simulink is used, and the basic sampling time is set to 1/1024 s. The external input is represented by Gaussian white noise n(t)~N(μ,σ²), where μ = 220 and σ² = 100 [22]. The standard physiological parameters of the EEG are taken as the model parameters and are set to the widely used values in the basic neuron mass model. The model parameters and their values are shown in Table 1.

By adjusting η_e, η_i, ω_e, and ω_i of the dynamic linear transformation function, the basic single-channel neuron mass model can output various types of simulated EEG signals with δ, α, or γ rhythms. When η_e, η_i, ω_e, ω_i, and external input n(t) are set to the typical values provided in Table 1, the neuron mass model outputs three basic dynamic waveforms that are significantly different in the time domain and frequency domain (Figure 6). In the following section, we will analyze the effect of each parameter on the model in detail.

3.1. The Effect of η_e on the Model

For the simulated EEG signals, μ = 220, σ² = 100, η_i = 1100 mV/s, ω_e = 100 s⁻¹, ω_i = 50 s⁻¹, the basic sampling time is adjusted to 1/1024 s, and other parameters are selected as shown in Table 1. To explore the influence of η_e on the dynamic behavior of the model, the value of the excitatory synaptic gain G_e is adjusted by changing the value of η_e. The simulation results are shown in Figure 7.

Based on Figure 7, different EEG signals can be generated by changing the value of η_e. When η_e = 92.59 mV/s, the model outputs a low-frequency signal with sporadic spikes only at the beginning. When η_e = 277.78 mV/s, the model output signal has an obvious α rhythm component, and the frequency is mainly concentrated at 10 Hz. At this time, both the values of G_e and η_e are close to the typical physiological values. When η_e = 462.96 mV/s, the simulated signal has a regular spike wave that is slightly shifted compared to the α rhythm. Specifically, both the amplitude peak and the normalized spectrum peak of the model’s output signal increase, and the frequency is still concentrated at 10 Hz. When η_e = 648.15 mV/s, the outputted EEG signals of the model are re-transformed into an α waveform, the amplitude peak is further increased, and the frequency is stable at 10 Hz. Therefore, we can conclude from Figure 7 that the amplitude peak of the model’s output signal increases with η_e, and there is no frequency shift.

3.2. The Effect of η_i on the Model

Similar to the discussion on the effects of η_e, this section explores the effects of changing the value of inhibitory synaptic gain η_i on the simulation results. Setting η_e = 325 mV/s and keeping other parameters unchanged gives the EEG simulation results shown in Figure 8.

As shown in Figure 8, when η_i = 500 mV/s, the output signals do not fluctuate significantly, there are sporadic spikes in the initial stage, and the normalized spectrum does not show a peak. When η_i is increased from 500 mV/s to 1000 mV/s, the output EEG signals are transformed into an α waveform (see the left column of (b)), and the frequency is concentrated at 10 Hz, as shown in the right column of Figure 8b. However, α rhythms disappear, and regular spikes appear when η_i is increased from 1000 mV/s to 1500 mV/s. The frequencies in the normalized spectrum are mainly concentrated in the range of 3–16 Hz, showing frequency harmonics; that is, the frequency peak appears at 5 Hz, 10 Hz, and 15 Hz and decreases in turn. When η_i is increased to 2000 mV/s, the α feature of the simulated EEG signals disappears, and regular spikes appear. In the normalized spectrum, the frequencies are mainly concentrated between 3 and 20 Hz, and frequency harmonics remain. The frequency peaks appear at 5 Hz, 10 Hz, 15 Hz, and close to 20 Hz. Moreover, the spectrum peaks also decrease, but their values are more obvious. Changing the value of η_i shows that, in addition to the main frequency, the simulated signal spectrum can generate frequency harmonics, which are integral multiples of the main frequency. This is similar to the frequency features of EEG signals generated by steady-state visual induction.

3.3. The Effects of External Inputs on the Model

3.3.1. The Effect of External Inputs’ Mean Value

In the single-channel basic neuron mass model, all the signals in the variable regions and external stimuli are represented by Gaussian white noise n(t)~N(μ,σ²) as excitatory input. μ is the mean value of n(t), indicating the average strength of external input. Therefore, the effect of external input on the output signals of a simulated single-cell mass can be explored by changing the value of μ. Setting σ² = 100, basic sampling time = 1/1024 s, η_e = 325 mV/s, η_i = 1100 mV/s, ω_e = 100 s⁻¹, and ω_i = 50 s⁻¹ in the single-channel basic neuron mass model gives the output signals shown in Figure 9.

As shown in Figure 9, when μ = 50, the output signals have more low-frequency components, and sporadic spikes appear only at the beginning. When μ is increased from 50 to 100, regular spikes appear in the signals, and their frequencies are mainly centered in the 1–11 Hz range. Furthermore, when μ is increased from 50 to 150, a large number of regular spikes with a relatively obvious α waveform sign appear in the simulated signals, and there is a spectrum peak at 7–9 Hz. The rhythm characteristics are strengthened when μ is increased to 200, and a spectral peak appears at 10 Hz.

3.3.2. The Effect of External Inputs’ Variance

In n(t), when the average strength of the external inputs is determined, the variance σ² can be used to represent their degree of fluctuation. The effects of σ² on EEG signals when μ = 220 and the other parameters are kept constant are shown in Figure 10.

As illustrated in Figure 10, when σ² is increased from 50 to 100, all the simulated signals of the model show α waveform characteristics, the frequencies are concentrated at 10 Hz, and the amplitude and spectrum do not change. Since the changes in the simulated signal are too weak, the value of σ² is increased from 100 to 3000. In this situation, the values of the amplitude of the simulated signals change slightly, but the normalized spectrum remains unchanged. The α waveform characteristic of the simulated signals is still obvious, and the frequencies are still concentrated at 10 Hz when σ² is increased from 3000 to 6000; however, the peak amplitude and frequency of the simulated signals change slightly. Further increasing σ² to 20,000 results in the α rhythm changing significantly (although the output signals still show an α rhythm) and the frequency peak shifting to the left.

There is no obvious characteristic change after adjusting the external inputs’ fluctuation degree when the average strength of external inputs is 220, considering that when μ = 220, the model outputs the standard α rhythm. To further evaluate the effect of σ² on the generated signals, μ is adjusted to 90, and the other parameters are kept consistent with those for μ = 220. Figure 11 shows the results.

The simulation results show that when μ = 90 and σ² = 100, the simulation signals of the model exhibit weak oscillations, and there is no frequency peak. When σ² is increased from 100 to 1000 and then from 1000 to 3000, the amplitude range of the model simulation signals’ oscillations increases significantly without frequency peaks. Sporadic spikes appear in the model simulation signals, and the spectrum changes when σ² is increased from 3000 to 6000. Furthermore, when σ² is increased to 20,000, a large number of spikes appear in the simulation signals, and frequency peaks appear in the spectrum, which are mainly concentrated at 5–10 Hz. In summary, when the average strength of external inputs is adjusted from 220 to 90, the simulation signals of the model change significantly with different σ² values.

3.4. The Effect of Mean Time Constant on the Model

The above analyses of the effects of η, μ, and σ² on the waveform and frequency of single-cell mass simulation signals are analyzed. Changing the average time constant ω⁻¹ affects the value of η and changes other parts of the dynamic linear transformation function in Equation (4).

When μ = 220, σ² = 100, and the basic sampling time is 1/1024 s, η_e represents the product of the average value of excitability synaptic gain G_e and the reciprocal of the average time constant ω_e. Keeping G_e = 3.25 mV, η_i = 1100 mV/s, and ω_i = 50 s⁻¹ to analyze the effect of excitability time constant ω_e⁻¹ on the output EEG signals of the model gives the results shown in Figure 12.

Figure 12 shows that as ω_e⁻¹ increases, the number of spikes of simulated EEG signals gradually decreases, with the peak values gradually increasing and moving to the left, and the spectrum peak value also gradually increasing. The excitability time constant affects the spike amplitude, spike number, spectrum peak value, and spectrum amplitude of the simulated EEG signals the model outputs.

In addition to the excitability time constant, the inhibitory time constant ω_i⁻¹ can also affect the simulated EEG signals. To explore this, the following parameters were set: μ = 220, σ² = 100, basic sampling time = 1/1024 s, G_i = 22 mV, η_e = 325 mV/s, and ω_e = 100 s⁻¹. η_i represents the product of the inhibitory mean synaptic gain G_i and the reciprocal average time constant ω_i. The simulation results are shown in Figure 13.

As shown in Figure 13, when ω_i⁻¹ = 0.010 s, the EEG signals outputted by the model show sporadic spikes at the beginning and then oscillate weakly. Moreover, there is no peak value in the normalized spectrum. When ω_i⁻¹ is increased to 0.015 s, the amplitude of EEG signals increases, and they show an α rhythm. The spike peaks gradually decrease, and the spectrum peak appears at 10–11 Hz. The amplitude of the EEG signals further increases when ω_i⁻¹ is increased to 0.020 s. At the same time, the number of signal spikes decreases, and the spectrum peak shifts to the left, but its value increases. When ω_i⁻¹ is further increased to 0.025 s, regular spikes appear in the output signals, α rhythm characteristics gradually disappear, signal amplitude increases, the number of spikes decreases, and spectrum harmonics are present. The spectrum peaks appear at 5 Hz, 10 Hz, and 15 Hz, and their values decrease successively.

4. Discussion and Conclusions

Based on the feedback mechanism and principle of the basic neuron mass model, this study introduced the static nonlinear information transformation process and dynamic linear information transformation process of the model, with a detailed explanation of the physiological parameters and their physical significance. Through simulation analysis of a single-channel basic NMM, we thoroughly investigated the effects of key parameters on EEG signal output characteristics and their normalized spectra.

The simulation results of the model demonstrated that the average synaptic gain G, the average time constant ω⁻¹, and the external input n(t) are the primary factors affecting EEG signal output. Specifically, the mean intensity μ and fluctuation degree σ² of the external input n(t) jointly influence the model’s output, with particularly pronounced effects when σ² undergoes significant changes. Furthermore, by adjusting the inhibitory product η_i and inhibitory time constant ω_i⁻¹, we observed the emergence of frequency harmonics in the spectral analysis of the simulated signals, occurring at integer multiples of the fundamental frequency. This phenomenon closely aligns with the characteristics of fundamental and harmonic responses in SSVEP signals, providing theoretical support for understanding the physiological mechanisms of SSVEP.

While this research employed a basic NMM to analyze SSVEP response mechanisms, several critical limitations must be acknowledged.

First, the current model relies on white Gaussian noise as the sole input perturbation, which oversimplifies the stochastic nature of neuronal activity. Real cortical dynamics are influenced by a mixture of white and pink noise (1/f noise), with pink noise dominating low-frequency components. This omission may restrict our ability to capture low-frequency phase coherence and stochastic resonance effects critical for SSVEP generation [35]. For instance, pink noise’s energy distribution could enhance thalamocortical entrainment of γ-band oscillations (30–80 Hz), amplifying SSVEP amplitude through nonlinear interactions [36]. Therefore, the absence of pink noise in our simulations may underestimate the stability and frequency specificity of SSVEP responses in realistic scenarios.

Second, although this study concentrated on narrowband behavior, its limited capability to capture multiple dynamic patterns constitutes a critical methodological constraint. While our model effectively captures SSVEP generation within local V1 circuits, it fails to account for the rich poly-frequency oscillatory dynamics observed in real brains. David et al. [37] proposed multi-dynamic neuron population models (multi-dynamic NMMs) that capture complex oscillatory dynamics in specific brain regions, such as the spatiotemporal coupling of γ-band rhythms and θ-α cross-frequency interactions, thereby establishing a more comprehensive dynamical framework for SSVEP analysis.

Third, this research focused on a single-channel basic NMM, which lacks the capacity to model interregional neural coupling mechanisms, such as functional interactions between the occipital visual cortex and frontal executive control regions. These cross-brain synchronization processes—exemplified by theta–alpha cross-frequency coupling during attentional modulation—not only directly influence SSVEP amplitude modulation but also contribute to the spatiotemporal organization of phase locking [6]. Neglecting these coupled dynamics risks fundamental misinterpretations of the hierarchical structure in SSVEP-related neural networks [38].

Currently, the development of SSVEP-based brain–computer interface (BCI) devices is constrained by experimental equipment, conditions, and individual subject variations. Acquiring real SSVEP signals involves high experimental costs, and insufficient data volume limits further research. Recent studies have attempted to generate SSVEP data consistent with real signal distributions through deep learning methods for data augmentation. However, these generation methods typically maintain mathematical distribution consistency with original data while lacking physiological interpretability and credibility, making it challenging to fully meet practical SSVEP signal modeling requirements.

Looking ahead, to address the inherent limitations of current NMMs and ensure the reliability and physiological significance of generated signals, future research should explore the development of hybrid generation frameworks that synergize systematically optimized neuron population models (NMMs) with deep learning techniques. The NMM could provide a signal generation framework based on neurophysiological mechanisms, ensuring the biological plausibility of generated signals. Meanwhile, deep learning methods could optimize generated signals using actual data distribution characteristics, enhancing the model’s ability to fit real signals. This integration would not only address SSVEP data insufficiency but also enable the real-time generation of biologically meaningful signals for practical BCI device control and optimization. This exploration direction presents new opportunities for SSVEP-based BCI technology development while promoting the convergence of neuroscience and artificial intelligence.

Author Contributions

Conceptualization, H.L.; methodology, Y.W.; software, Y.W.; writing—original draft preparation, D.G.; writing—review and editing, D.G.; validation, P.F.; funding acquisition, J.Q.; project administration, H.L. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All data generated or analyzed to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Single-channel basic neuron mass feedback model.

Figure 2. Block diagram of a single-channel basic neuron mass.

Figure 3. Static nonlinear function.

Figure 4. Changes in the impulse response of the linear function with average synaptic gain Ge.

Figure 5. Changes in the impulse response of the linear function with the average time constant ωe.

Figure 6. Simulated sequences and the normalized spectrum of δ, α, and γ.

Figure 7. Simulated EEG signal curves and the normalized spectrum with different ηe.

Figure 7. Simulated EEG signal curves and the normalized spectrum with different ηe.

Figure 8. Simulated EEG signal curves and the normalized spectrum with different ηi.

Figure 9. Simulated EEG signal curves and the normalized spectrum with different μ.

Figure 10. Simulated EEG signal curves and the normalized spectrum with different σ2 when μ = 220.

Figure 11. Simulated EEG signal curves and the normalized spectrum with different σ2 when μ = 90.

Figure 11. Simulated EEG signal curves and the normalized spectrum with different σ2 when μ = 90.

Figure 12. Simulated EEG signal curves and the normalized spectrum with different ωe−1.

Figure 12. Simulated EEG signal curves and the normalized spectrum with different ωe−1.

Figure 13. Simulated EEG signal curves and the normalized spectrum with different ωe−1.

Figure 13. Simulated EEG signal curves and the normalized spectrum with different ωe−1.

Table 1

Physiological significance and typical values of the model parameters.

Parameters	Physiological Significance	Typical Values
		δ Waveform	α Waveform	γ Waveform
η e	Ge * ωe	40 mV/s	325 mV/s	1630.43 mV/s
η i	Gi * ωi	300 mV/s	1100 mV/s	51,724.14 mV/s
ωe −1	Excitability time constant	0.05 s	0.0108 s	0.0046 s
ωi −1	Inhibitory time constant	0.05 s	0.02 s	0.0029 s
C1, C2	Excitatory mean synaptic connections	C1 = C, C2 = 0.8C(C = 135)
C3, C4	Inhibitory mean synaptic connections	C3 = 0.25CC4 = 0.25C
v0, e0, r	Parameters of the nonlinearfunction	v0 = 6 mV, 0 = 2.5 s−1, r = 0.56 mV−1

References

1. Hodgkin, A.L.; Huxley, A.F. A quantitatived escription of membrane current and its application to conduction and excitation in nerve. J. Physiol.; 1952; 117, pp. 500-544. [DOI: https://dx.doi.org/10.1113/jphysiol.1952.sp004764] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12991237]

2. Chay, T.R. Chaos in a three-variable excitable cell model. Phys. D Nonlinear Phenom.; 1985; 16, pp. 233-242. [DOI: https://dx.doi.org/10.1016/0167-2789(85)90060-0]

3. Traub, R.D.; Contreras, D.; Cunningham, M.O.; Murray, H.; LeBeau, F.E.N.; Roopun, A.; Bibbig, A.; Wilent, W.B.; Higley, M.J.; Whittington, M.A. Single-column thalamocortical network model exhibiting gamma oscillations, sleep spindles, and epileptogenic bursts. J. Neurophysiol.; 2005; 93, pp. 2194-2232. [DOI: https://dx.doi.org/10.1152/jn.00983.2004] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15525801]

4. Schelter, B.; Mader, M.; Mader, W.; Sommerlade, L.; Platt, B.; Lai, Y.C.; Grebogi, C.; Thiel, M. Overarching framework for data-based modelling. Europhys. Lett.; 2014; 105, 30004. [DOI: https://dx.doi.org/10.1209/0295-5075/105/30004]

5. Dong, E.; Liang, Z. The multi-frequency EEG rhythms modeling based on two-parameter bifurcation of neural mass model. Proceedings of the IEEE International Conference on Mechatronics and Automation; Tianjin, China, 3–6 August 2014; pp. 1564-1569.

6. Wilson, H.R.; Cowan, J.D. Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J.; 1972; 12, pp. 1-24. [DOI: https://dx.doi.org/10.1016/S0006-3495(72)86068-5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/4332108]

7. da Silva, F.H.L.; Hoeks, A.; Smits, H.; Zetterberg, L.H. Model of brain rhythmic activity. Kybernetik; 1974; 15, pp. 27-37. [DOI: https://dx.doi.org/10.1007/BF00270757]

8. da Silva, F.H.L.; Rotterdam, A.V.; Barts, P.; van Heusden, E.; Burr, W. Models of neuronal populations: The basic mechanisms of rhythmicity. Prog. Brain Res.; 1976; 45, pp. 281-308.

9. Jansen, B.H.; Zouridakis, G.; Brandt, M.E. A neurophysiologically-based mathematical model of flash visual evoked potentials. Biol. Cybern.; 1993; 68, pp. 275-283. [DOI: https://dx.doi.org/10.1007/BF00224863]

10. Tsoneva, T.; Garcia-Molina, G.; Desain, P. SSVEP phase synchronies and propagation during repetitive visual stimulation at high frequencies. Sci. Rep.; 2021; 11, 4975. [DOI: https://dx.doi.org/10.1038/s41598-021-83795-9]

11. Ding, J.; Sperling, G.; Srinivasan, R. Attentional modulation of SSVEP power depends on the network tagged by the flicker frequency. Cereb. Cortex; 2006; 16, pp. 1016-1029. [DOI: https://dx.doi.org/10.1093/cercor/bhj044]

12. Li, F.; Tian, Y.; Zhang, Y.; Qiu, K.; Tian, C.; Jing, W.; Liu, T.; Xia, Y.; Guo, D.; Yao, D. et al. The enhanced information flow from visual cortex to frontal area facilitates SSVEP response: Evidence from model-driven and data-driven causality analysis. Sci. Rep.; 2015; 5, 14765. [DOI: https://dx.doi.org/10.1038/srep14765] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26434769]

13. Fu, P.R.; Li, H.; Zhang, S. Review of latest noninvasive EEG-based robotic device. Proceedings of the IEEE/ASME (AIM) International Conference on Advanced Intelligent Mechatronics; Sapporo, Japan, 11–15 July 2022; pp. 11-15.

14. Lin, C.; Deng, X.; Yu, Z.L.; Gu, Z. A SSVEP-based BCI for controlling a 4-DOF robotic manipulator. Proceedings of the IEEE International Conference on Systems, Man and Cybernetics; Bari, Italy, 6–9 October 2019; pp. 2174-2179.

15. Ishizuka, K.; Kobayashi, N.; Saito, K. High accuracy and short delay 1ch-ssvep quadcopter-BMI using deep learning. J. Robot. Mechatron.; 2020; 32, pp. 738-744. [DOI: https://dx.doi.org/10.20965/jrm.2020.p0738]

16. Shao, L.; Zhang, L.; Belkacem, A.N.; Zhang, Y.; Chen, X.; Li, J.; Liu, H. EEG-controlled wall-crawling cleaning robot using SSVEP-based brain-computer interface. J. Healthc. Eng.; 2020; 2020, 6968713. [DOI: https://dx.doi.org/10.1155/2020/6968713] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32399166]

17. Jebri, A.; Madani, T.; Djouani, K. Neural adaptive integral-sliding-mode controller with a SSVEP-based BCI for exoskeletons. Proceedings of the International Conference on Advanced Robotics; Belo Horizonte, Brazil, 2–6 December 2019; pp. 87-92.

18. Na, R.; Hu, C.; Sun, Y.; Wang, S.; Zhang, S.; Han, M.; Yin, W.; Zhang, J.; Chen, X.; Zheng, D. An embedded lightweight SSVEP-BCI electric wheelchair with hybrid stimulator. Digit. Signal Process.; 2021; 16, 103101. [DOI: https://dx.doi.org/10.1016/j.dsp.2021.103101]

19. Gao, Q.; Zhao, X.; Yu, X.; Song, Y.; Wang, Z. Controlling of smart home system based on brain-computer interface. Technol. Health Care; 2018; 26, pp. 769-783. [DOI: https://dx.doi.org/10.3233/THC-181292]

20. Li, H.; Li, X.; Millán, J.d.R. Noninvasive EEG-based intelligent mobile robots: A systematic review. IEEE Trans. Autom. Sci. Eng.; 2024; pp. 1-25. [DOI: https://dx.doi.org/10.1109/TASE.2024.3441055]

21. Li, H.; Bi, L.; Li, X. Robust predictive control for EEG-based brain-robot teleoperation. IEEE Trans. Intell. Transp. Syst.; 2024; 25, pp. 9130-9140. [DOI: https://dx.doi.org/10.1109/TITS.2024.3359216]

22. Jansen, B.H.; Rit, V.G. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol. Cybern.; 1995; 73, pp. 357-366. [DOI: https://dx.doi.org/10.1007/BF00199471]

23. Braitenberg, V.; Schüz, A. Anatomy of the Cortex: Statistics and Geometry; Springer-Verlag Publishing: Berlin/Heidelberg, Germany, 1991.

24. Keller, A.; White, E.L. Triads: A synaptic network component in the cerebral cortex. Brain Res.; 1989; 496, pp. 105-112. [DOI: https://dx.doi.org/10.1016/0006-8993(89)91056-1]

25. Elhanany, E.; White, E.L. Intrinsic circuitry: Synapses involving the local axon collaterals of corticocortical projection neurons in the mouse primary somatosensory cortex. J. Comp. Neurol.; 1990; 291, pp. 43-54. [DOI: https://dx.doi.org/10.1002/cne.902910105]

26. Liu, X.B.; Zheng, Z.H.; Xi, M.C.; Wu, C.P. Distribution of synapses on an intracellularly labeled small pyramidal neuron in the cat motor cortex. Brain Struct. Funct.; 1991; 184, pp. 313-318. [DOI: https://dx.doi.org/10.1007/BF01673266] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/1793170]

27. Larkman, A.U. Dendritic morphology of pyramidal neurones of the visual cortex of the rat. J. Comp. Neurol.; 1991; 306, pp. 332-343. [DOI: https://dx.doi.org/10.1002/cne.903060209] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/1711059]

28. White, E.L. Termination of Thalamic Afferents in the Cerebral Cortex. Sensory-Motor Areas and Aspects of Cortical Connectivity; Springer: Berlin/Heidelberg, Germany, 1986; Volume 5, pp. 271-289.

29. Schnepel, P.; Kumar, A.; Zohar, M.; Aertsen, A.; Boucsein, C. Physiology and impact of horizontal connections in rat neocortex. Cereb. Cortex; 2015; 25, pp. 3818-3835. [DOI: https://dx.doi.org/10.1093/cercor/bhu265]

30. van Rotterdam, A.; da Silva, F.H.L.; van den Ende, J.; Viergever, M.A.; Hermans, A.J. A model of the spatial-temporal characteristics of the alpha rhythm. Bull. Math. Biol.; 1982; 44, pp. 283-305. [DOI: https://dx.doi.org/10.1007/BF02463252]

31. Dodt, H.U.; Pawelzik, H.; Zieglga, W. Actions of noradrenaline on neocortical neurons in vitro. Brain Res.; 1991; 545, pp. 307-311. [DOI: https://dx.doi.org/10.1016/0006-8993(91)91303-I]

32. Mina, F.; Benquet, P.; Pasnicu, A.; Biraben, A.; Wendling, F. Modulation of epileptic activity by deep brain stimulation: A model-based study of frequency-dependent effects. Front. Comput. Neurosci.; 2013; 7, 94. [DOI: https://dx.doi.org/10.3389/fncom.2013.00094]

33. McCormick, D.A.; Pape, H.C.; Williamson, A. Actions of norepinephrine in the cerebral cortex and thalamus: Implications for function of the central noradrenergic system. Prog. Brain Res.; 1991; 88, pp. 293-305.

34. Terakado, M. Adrenergic regulation of GABA release from presynaptic terminals in rat cerebral cortex. J. Oral Sci.; 2014; 56, pp. 49-57. [DOI: https://dx.doi.org/10.2334/josnusd.56.49] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24739708]

35. Morikawa, M.; Nakamichi, A. A simple model for pink noise from amplitude modulations. Sci. Rep.; 2023; 13, 8364. [DOI: https://dx.doi.org/10.1038/s41598-023-34816-2]

36. Hatayama, A.; Matsubara, A.; Nakashima, S.; Nishifuji, S. Effect of pink noise on EEG and memory performance in memory task. Proceedings of the 2021 IEEE 10th Global Conference on Consumer Electronics (GCCE); Kyoto, Japan, 12–15 October 2021; pp. 238-241.

37. David, O.; Friston, K.J. A Neural Mass Model for MEG/EEG: Coupling and Neuronal Dynamics. NeuroImage; 2003; 20, pp. 1743-1755. [DOI: https://dx.doi.org/10.1016/j.neuroimage.2003.07.015]

38. Zanto, T.P.; Rubens, M.T.; Thangavel, A.; Gazzaley, A. Causal role of the prefrontal cortex in top-down modulation of visual processing and working memory. Nat. Neurosci.; 2011; 14, pp. 656-661. [DOI: https://dx.doi.org/10.1038/nn.2773] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21441920]