# Closed-Loop Current Stimulation Feedback Control of a Neural Mass Model Using Reservoir Computing
## 2 Methods
### Training Set Generation  
For the rest of the paper, we use the notation \(\mathbf{X}(\mathbf{n})\) to denote a vector time-series \(\mathbf{X}\) that ends at time \(n\). Each signal \(\mathbf{I}(n)\), \(\mathbf{I}(n-\delta)\), \(\mathbf{p}(n)\), \(\mathbf{p}(n+\delta)\) has dimensions \(\mathbb{R}^{2\times 60}\), where 2 corresponds to the number of columns (i.e., \(\mathbf{p}\) is a concatenation of \(\mathbf{p}_{1}\) and \(\mathbf{p}_{2}\)) and 60 is the number of samples. Signal \(\mathbf{I}(n-\delta)\) is a delayed copy of \(\mathbf{I}(n)\), delayed by the time step \(\delta\), which was set to 0.01 s. The time sample size 60 (taken with a sampling period of 0.01 s) was chosen to capture one full cycle of frequencies greater than 1.67 Hz, which we found to be adequate to capture the model dynamics since the peak model frequency was around 5 Hz. The resulting input to the ESN \(\mathbf{u}(n)\) was a concatenation of \(\mathbf{I}(n-\delta)\), \(\mathbf{p}(n)\), \(\mathbf{p}(n+\delta)\) to yield a signal with size \(\mathbb{R}^{6\times 60}\). The output signal predicted by the ESN \(\mathbf{y}(n)\) was set to be \(\mathbf{I}(n)\), the current stimulation input to be applied that led to the future dynamics \(\mathbf{p}(n+\delta)\). A schematic describing the prediction procedure can be seen in Figure 3.  
To invert the dynamics of the JRNMM, the ESN must observe a sufficient number of \(\mathbf{I}(n)\) current stimulations and corresponding \(\mathbf{p}(n)\) PC dynamics. To do so, we developed the following paradigm to generate various current stimulation inputs to the model. For each of the two PC layers in each of the columns, the current input \(\mathbf{I}(n)\) was white noise with a noise power of 0.1 (height of the power spectral density). The noise was then bandpass filtered between 0.1 and 30 Hz. To allow for the ESN to observe both onset and offsets of current stimulation, the resulting bandpass filtered noise was multiplied by a unit amplitude pulse function. The pulse function had a period uniformly drawn from \([0.1,0.2,\ldots,1]\) seconds. The pulse had a width of \([10,20,\ldots,90]\) percent of the period, also uniformly drawn. Each pulse started with a time delay of 0.5 s. The pulse was then multiplied by a gain uniformly drawn from \([1,2,\ldots,10]\). Figure 4 shows a sample current stimulation input. The simulation was run 1000 times at a sampling period of 0.001 s with a simulation duration of 4 s, then downsampled to a sampling period of 0.01. Each iteration had randomized current inputs \(\mathbf{I}(n)\) across all iterations and PC layers, as well as randomized \(\mathbf{u}(n)\) inputs to the excitatory interneurons.  
Figure 3: The current stimulation \(\mathbf{I}(n)\) is injected into the JRNMM to yield the PC dynamics \(\mathbf{p}(n)\). The previous inputs \(\mathbf{I}(n-\delta)\), along with the current and future dynamics \(\mathbf{p}(n)\) and \(\mathbf{p}(n+\delta)\) were used to predict the input \(\mathbf{I}(n)\) which gave rise to the future dynamics. Random noise \(\mathbf{g}(n)\) is also input into the excitatory interneurons. Only one cortical column is shown; however, this process is carried out on the two-column model.