4. Applications¶

In NPDS Toolbox we can apply different control strategies to change neural population oscillators synchrony and investigate the dynamics of some of the common neural models according to the values of their parameters. In this section, we describe different options of the toolbox for these tasks.

As we mentioned before, this toolbox has two main pages that can be accessed by running the following command in MATLAB.

run NPDSToolbox.m

Running the above command opens the main page of the package, which is related to the control of synchronous neuronal populations (See Fig. 4.1).

_images/Main1.png — Fig. 4.1 The main window of NPDS Toolbox¶

On this page, using $\texttt{Tools}$, you can access another main page of the package called NeuronDynamics, which is related to the study of neuron dynamics (see Fig. 4.2).

_images/Neurondynamics.png — Fig. 4.2 The graphical user interface showing the dynamics of the included neural models.¶

We will now describe each of them in detail.

4.1. Neural populations (de)synchronization¶

The main innovation of the NPDS toolbox is that not only are there some famous and predefined control strategies as proportional and bang–bang control laws, but researchers can also define their own control laws and examine their effectiveness in controlling the phase distribution.

In this part, users can set the desired control input (the area marked with number 2 in red in Figure Fig. 4.1), then, choose the type of neural oscillators with their period and determine if they have noise (the area marked with number 3 in red in Figure Fig. 4.1). Afterwards, they should specify the initial and final distributions (number 4 in red in Figure Fig. 4.1).These distributions can be chosen from predefined ones (Von-Mises or uniform) or defined manually in $\texttt{user\_defined\_initial\_dist}$ and $\texttt{user\_defined\_final\_dist}$ files. Moreover, the number of oscillators and simulation time can be defined by users (number 5 in red in Figure Fig. 4.1). Finally, the numerical approach for the simulation should be selected (number 6 in red in Figure Fig. 4.1). The details of the proposed numerical methods are described in previous sections. At the top of the toolbox home screen, you have access to two $\texttt{Tools}$ and $\texttt{Help}$ options, which you can use to check the dynamics of the neurons and the information about the toolbox, respectively. Also below these two options is an advanced toolbar that can be used to configure simulated results. After setting all the above, you can click the $\texttt{Start}$ option at the bottom right of the page to start the simulation and see the results of the phase distribution control of the specified neuronal population in the middle panel of the page as an animation in different time units and whenever you want youcan hold down the simulation using the $\texttt{Pause}$ option at the bottom right of the window. In the following, the properties of each panel in the window is describe in details.

4.1.1. The display panels¶

When the simulation starts, four figures appear in the main window like Fig Fig. 4.3 which show the simulation results over time. In the top left panel the black and red lines show the probability distribution $\rho(\theta,t)$ and desired probability distribution $\rho_f(\theta,t)$, respectively. For example, here in this figure, from the right panel, the initial and final distributions are selected $\texttt{Von-Mises}$ and $\texttt{uniform}$, respectively (we will discuss this panel in detail later). The top right panel shows the control input $u(t)$, for example $\texttt{Proportional}$ is selected here. The bottom left and right panels demonstrate the error ($\int_0^{2\pi}\big(\rho(\theta,t)-\rho_f(\theta,t)\big)^2$) and the phase oscillators, respectively. Using the Pause button, we can stop the simulation whenever we want and resume it again.

_images/1.png — Fig. 4.3 Plots of simulation results over the time.¶

At the end of the simulation (see Figure Fig. 4.4), in the bottom panel some information are reported such as the type of phenomenon (desynchronization, synchronization, clustering, etc.), the final error value, the variance of the initial and final distributions, and the initial and final population variances. The obtained results can be saved as figures and matrix forms using the toolbar buttons. In addition, there are other features in the toolbar such as zoom in, zoom out, add legends, etc.

_images/2.png — Fig. 4.4 An example of the final results.¶

4.1.2. Control inputs¶

In the Control type panel (top right in Figure Fig. 4.5), the control strategy is designed. There is three types of control input: $\texttt{Proportional}$ [MM19b][ Moayeri2021], $\texttt{Bang-Bang}$ [MFM18] and $\texttt{User-defined}$.

In the $\texttt{Proportional}$ control law we have:

(4.2)¶\[U_p(t)=max(min(U(max),-KI(t)),U(min)),\]

where the oscillators should be noise-free. On the other hand we have in noisy oscillators

(4.2)¶\[U_p(t)=max(min(U(max),-KI(t)-\frac{G(t)}{I(t)}),U(min)),\]

in which

(4.3)¶\[I(t)=2\int_0^{2\pi}\Big(\frac{\partial\rho(\theta,t)}{\partial \theta}-\frac{\partial \rho_f(\theta,t)}{\partial \theta}\Big)\mathcal{Z}(\theta)\rho(\theta,t)d\theta,\]

(4.4)¶\[G(t)=-\mathcal{B}\int_0^{2\pi}\Big(\frac{\partial\rho(\theta,t)}{\partial \theta}-\frac{\partial \rho_f(\theta,t)}{\partial \theta}\Big)\frac{\partial \rho(\theta,t)}{\partial \theta}d\theta.\]

where $K$, $U(min)$ (lower bound) and $U(max)$ (upper bound) are free parameters that are determined by the user (see top right in Figure Fig. 4.4).

_images/Main1_1.png — Fig. 4.5 Take a closer look at the $\texttt{control type}$ panel options.¶

In Fourier decomposition method, $I(t)$ and $G(t)$ are computed as follows:

(4.5)¶\[I(t)=\sum_{k=1}^{N-1}\Big(\big(A_k(t)-\hat{A}_k(t)\big)^2+\big(B_k(t)-\hat{B}_k(t)\big)^2\Big),\]

(4.6)¶\[G(t)=-\frac{\mathcal{B}^2}{2}\sum_{k=1}^{N-1}k^2\Big(A_k(t)\big(A_k(t)-\hat{A}_k(t)\big)+B_k(t)\big(B_k(t)-\hat{B}_k(t)\big)\Big).\]

Actually, we can not choose large values for upper and lower bounds of control input because large electrical pulses may damage tissues in practice. Additionally, large values for $K$ make proportional controllers like bang-bang one. On the other hand, if we choose small values for $K$ it just decreases $L^2-$ cannot be negative and $U(min)$ should be less than $U(max)$.

In $\texttt{Bang-Bang}$ and $U(max)$ according to the sign of $I(t)$:

(4.7)¶\[U_b(t)=\begin{cases} U(min), & I(t)>0, \\ U(max), & I(t)<0.\end{cases}\]

Finally, by choosing $\texttt{User-defined}$ option, a Matlab file is opened. In this file, the user should define his control law in the body of $\texttt{user\_defined\_control(varagin)}$ function according to the appropriate arguments. The arguments of $\texttt{user\_defined\_control}$ is depends on the numerical method and the existence of noise. If the numerical approach is a method other than Fourier decomposition, the arguments are the spatial nodes, $\rho(\theta,t)$, $\rho_f(\theta,t)$, $\mathcal{Z}(\theta)$, $\mathcal{Z}'(\theta)$, $\omega$, the number of time steps, and $\Delta t$. If there is noise in the system, the intensity of noise, $\mathcal{B}$, $\frac{\partial\rho(\theta,t)}{\partial\theta}$, and $\frac{\partial\rho_f(\theta,t)}{\partial\theta}$ are added to input arguments. But when the numerical approach is Fourier decomposition method, the input arguments are spatial nodes, $\mathcal{Z}(\theta)$, $\mathcal{Z}'(\theta)$, $\omega$, number of time steps, $\Delta t$, $A_k$, $B_k$, $\Tilde{A}_k$, $\Tilde{B}_k$, $\mathcal{I}_k^A(t)$ and $\mathcal{I}_k^B(t)$. Moreover, for the noisy systems, $\mathcal{B}$, intensity, and $k$ are also considered as inputs. Now, the control law should be defined using any of these input arguments.

For instance, assume we choose the spectral method as our numerical approach for simulation and a continuous control law as $U(t)=SI(t)$ where $S$ is constant. So, we have:

function u=user_defined_control(varagin)
if nargin ==9
    domain=varargin{1};phi=varargin{2};phi_f=varargin{3};
    prc=varargin{4};dprc=varargin{5};error=varargin{6};
    omega=varargin{7};iteration_number=varargin{8};
    dt=varargin{9};
elseif nargin == 12
    domain=varargin{1};prc=varargin{2};dprc=varargin{3};
    omega=varargin{4};iteration_number=varargin{5};
    dt=varargin{6};
    bk=varargin{7};bfk=varargin{8};ak=varargin{9};
    afk=varargin{10};Ika=varargin{11};Ikb=varargin{12};
elseif nargin == 13
    domain=varargin{1};phi=varargin{2};phi_f=varargin{3};
    prc=varargin{4};dprc=varargin{5};error=varargin{6};
    omega=varargin{7};iteration_number=varargin{8};
    dt=varargin{9};
    B=varargin{10};dphi=varargin{11};dphi_f=varargin{12};
    intensity=varargin{13};
elseif nargin == 15
    domain=varargin{1};prc=varargin{2};dprc=varargin{3};
    omega=varargin{4};iteration_number=varargin{5};
    dt=varargin{6};
    B=varargin{7};intensity=varargin{8};bk=varargin{9};
    bfk=varargin{10};ak=varargin{11};afk=varargin{12};
    Ika=varargin{13};Ikb=varargin{14};k=varargin{15};
end
%User-defined control law
S=20;
I=(trapz(domain,(phi-phif).*Zp'.*phi));
u=S*I;

4.1.3. Phase response curve (PRC)¶

On the toolbox homepage, in the second right panel from the top, we have the $\texttt{Neurons property}$ settings (see Figure Fig. 4.6).

In fact, in the neurons property panel, the neural model is chosen from the available models i.e. Hodgkin-Huxley (HH), Fitzhugh-Nagumo (FHN), Rose-Hindmarsh (RH), or Thalamic model. Each model has its own PRC function described in the following table. PRC is experimentally measurable by perturbing an oscillatory neuron at different phases and determining the change in spike timing [WM14]. Different neural models and populations may lead to different PRC which affects the control model.

Model	PRC
HH	$\texttt{PRC\_interpol(domain,'Hodgkin-Huxley')}$
FHN	$-43.41215711 \sin(\theta)$
RH	$\frac{1-\cos(\theta)}{2\pi}$
Thalamic	$\texttt{PRC\_interpol(domain,'Thalamic')}$

_images/Main1_2.png — Fig. 4.6 Take a closer look at the $\texttt{Neurons property}$ panel options.¶

Since HH and Thalamic PRCs have not any specific formula, we calculate PRC by interpolation from data. Their data are stored in $\texttt{matrices/Thalamic.mat}$ and $\texttt{matrices/Thalamic.mat}$ directories.

In addition, in this panel you can specify the period of the oscillators ( $\tau$), the presence or absence of noise. If there exists noise in the system (as Figure Fig. 4.7), the intensity option is available to the user.

_images/4.png — Fig. 4.7 An example of noisy neural oscillators.¶

In fact, if we want to add a Gaussian white noise $\sqrt{2D}\eta(t)$ with with zero mean and variance $2D$ variance to the model, it is enough to select the $\texttt{noise}$ option.

4.1.4. Initial and final distributions¶

In distribution control, the oscillators move from an initial probability distribution to the desired distribution over time.

In this toolbox (see Figure Fig. 4.8), users can utilize predefined $\texttt{Von-Mises}$ or $\texttt{uniform}$ distribution as initial and final distributions. Also, they can define their desired distributions in $\texttt{user\_defined\_initial\_dist}$ and :math:`texttt{user_defined_final_dist`$ functions for initial and final distributions, respectively.

_images/Main1_3.png — Fig. 4.8 Take a closer look at the $\texttt{Distribution}$ panel options.¶

In particular, $texttt{Von-Mises}$ distribution is defined as [BF79]:

(4.9)¶\[U\rho(\theta,t)=\frac{\exp(\kappa\cos(\theta+\pi))}{2\pi\mathcal{I}_0(\kappa)},\]

where $\mathcal{I}_0(\kappa)$ is the modified Bessel function of first kind of order 0. For such a distribution, the mean is $\theta_0$, and the variance is $1 -\mathcal{I}_1(\kappa)/\mathcal{I}_0(\kappa)$, in which $\mathcal{I}_1(\kappa)$ is the modified Bessel function of first kind of order 1. $\theta_0$ is interpreted as location in toolbox. Moreover the variance decreases as $\kappa$ (concentration in toolbox) increases. If we consider $\kappa=0$, the $\texttt{uniform}$ distribution is obtained. Figure Fig. 4.9 is an example of synchronization where the final distribution is uniform instead of the initial one.

If one choose $\texttt{User-defined$\_$initial}$ or $\texttt{User-defined$\_$final}$ options, the corresponding Matlab file is opened; then, a new distribution can be defined in it. For example, assume we tend to cluster synchronized neurons. Thus, we choose the initial distribution as $\texttt{uniform}$. For final distribution, we should choose $\texttt{User-defined$\_$final}$. Then, $\texttt{user\_defined\_final\_dist.m}$ is created. Our desired final distribution is a bi-modal distribution, which can be realized as a sum of two uni-modal $\texttt{von-Mises}$ distributions; so we have [MM19b]:

(4.9)¶\[\rho_f(\theta,t)=\frac{\exp(\kappa\cos(\theta+\frac{\pi}{2}))+\exp(\kappa\cos(\theta+\frac{3\pi}{2}))}{4\pi\mathcal{I}_0(\kappa)}.\]

It is enough to do the following in the user-defined function:

function [dist,diff_dist]=user_defined_final_dist(domain,omega,i,dt)
kappa=52;
dist=(exp(kappa*cos(domain-pi/2-omega*i*dt))...
+exp(kappa*cos(domain-3*pi/2-omega*i*dt)))...
/(4*pi*besseli(0,kappa));
diff_dist=(-kappa*exp(kappa*sin(domain-pi/2-omega*i*dt))...
-kappa*exp(kappa*sin(domain-3*pi/2-omega*i*dt)))...
/(4*pi*besseli(0,kappa));

$\texttt{omega*i*dt}$ expression in equations is due to the traveling of the wave over time.

_images/3.png — Fig. 4.9 An example of synchronization.¶

4.2. Neural dynamics¶

There is an option in the toolbox for those who are interested in examining the dynamics of neural models in the toolbox. For this purpose, choose the $\texttt{Neuron Dynamics}$ from $\texttt{Tools}$ as Figure Fig. 4.10 or use short key $\texttt{Crtl+D}$. Moreover, it possible to run this part directly on the MATLAB command lineby $\texttt{run NeuronDynamics.m}$ command. By selecting this item a window like Figure Fig. 4.11 is loaded. As mentioned earlier, by loading this part, the solution of the Hodgkin-Huxley dynamical model is computed automatically with default values. By adjusting any parameter or type of model, the solution is recomputed.

_images/toolss.png — Fig. 4.10 Select neuron dynamics tool to load that part where the dynamics of the model can be investigated.¶

_images/dynamic_main.png — Fig. 4.11 The main window of Neuron Dynamic part.¶

In the following, we describe the panels of this part and the details of the existing models.

4.2.1. The display panels¶

There are two main figures in this part. the left panel shows the time plots of the selected functions. The time plot of membrane potential is demonstrated as a default and you can add the other functions using the checkboxes below the figure. Moreover, the scale mode option scales the functions to the interval $[0,1]$ in such a way their details can be seen better.

The right panel displays the phase portrait of the two or three selected variables. For the three variables phase portrait, you should select the 3D checkbox first; then select the three desired variables. One can hide this figure by clicking on the phase portrait option in the toolbar. Moreover, in the toolbar, there are vector field and stream options that can be added to this figure. Note that one of these two options can be added each time.

At the bottom of these figures, there is a panel that represents some information of the simulation such as the neuron model, indicating whether the system is periodic or not and the CPU time for solving the system.

The right side panel contains the values of parameters of the selected model and slide bars to determine initial conditions and the amount of external current. this panel option is changed depending on the selected model.

The toolbox at the top of the window has other options as well. We can save the results and figures from there easily. Moreover, there is an option called “Set domain”. By clicking on it, the panel on the right side of the window is changed (see Figure Fig. 4.12). In this part, you can change the upper and lower bounds of initial conditions for variables, external current, and the duration of the simulation. In the end, we should click on OK to apply changes.

_images/11.png — Fig. 4.12 The set domain option where determine the upper and lower bounds of variables, external currents, and the time duration.¶

The remaining sections describe the existing dynamical models in the toolbox and show an example of each.

4.2.2. Hodgkin-Huxley model¶

The Hogkin-Huxley (HH) [HH52] is the pioneer model in computational neuroscience. This model describes the action potential in the squid giant axon by the kinetics of voltage-dependent ion channels i.e. sodium and potassium in the cell membrane.

The properties of an excitable cell are described by the following system of ordinary differential equations:

(4.11)¶\[C\dot{v}=I- \underbrace{g_K n^4(v-v_K)}_{I_{K}}-\underbrace{g_{Na}m^3h(v-v_{Na})}_{I_{Na}}-\underbrace{g_L(v-v_L)}_{I_{L}},\\ \dot{n}=\frac{n_{\infty}(v)-n}{\tau_n(v)},\\ \dot{m}=\frac{m_{\infty}(v)-m}{\tau_m(v)},\\ \dot{h}=\frac{h_{\infty}(v)-h}{\tau_h(v)},\]

for $p=(n,m,h)$, we have:

(4.11)¶\[p_{\infty}=\frac{\alpha_p}{\alpha_p+\beta_p},~~\tau_p=\frac{1}{\alpha_p+\beta_p},\]

where $\alpha_p$ and $\beta_p$ are called rate coefficient which depend on membrane potential, and we have:

(4.12)¶\[p_{\infty}=\frac{\alpha_p}{\alpha_p+\beta_p},~~\tau_p=\frac{1}{\alpha_p+\beta_p}, \alpha_n(v)=\frac{0.01(10-v)}{\exp\big(\frac{10-v}{10}\big)-1},~~\beta_n(v)=0.125\exp\big(\frac{-v}{80}\big),\\ \alpha_m(v)=\frac{0.1(25-v)}{\exp\big(\frac{25-v}{10}\big)-1},~~\beta_n(v)=4\exp\big(\frac{-v}{18}\big),\\ \alpha_h(v)=\exp\big(\frac{-v}{20}\big),~~\beta_h(v)=\frac{1}{\exp\big(\frac{30-v}{10}\big)+1}.\]

where $v(t)$ is the electrical potential across the membrane, $m(t)$ and $h(t)$ are the activation and inactivation variable of sodium channel and $n(t)$ is activation variable of potassium channel. Also, $I$ is an external current applying to the membrane, $C$ is the membrane capacity, $g_K$, $g_{Na}$ and $g_L$ are maximal conductance of the ion channels. In fact, $g_L$ corresponds to the leak channel which is carried mostly $Cl^-$ ions. Also, $v_k$, $v_{Na}$ and $v_L$ are the Nernst equilibrium potentials of the ion channels [Izh07].

The results of system HH_eq are obtained by using $\texttt{ode23tb}$ Matlab function. Figure Fig. 4.13 displays an example for HH model with respected parameters and initial condition values.

_images/6.png — Fig. 4.13 An Example of an HH dynamical system.¶

4.2.3. Fitzhugh-Nagumo model¶

Fitzhugh-Nagumo (FHN) model is a simplified 2D version of the HH model. This model is described as [Izh07][IF06]:

(4.13)¶\[\dot{v}=v-\frac{v^3}{3}-w+I,\\ \dot{w}=a(v+b-cw).\]

In this model $v(t)$ mimics the membrane potential and $w(t)$ is called the recovery variable showing activation of an outward current. Also, $I$ is applied current ,and $a$, $b$ and $c$ are constants. This system is solved by $\texttt{ode45}$ function in Matlab. An example of this dynamical system is represented in Figure Fig. 4.14. It should be noted that in this example the obtained results are in scale mode.

_images/8.png — Fig. 4.14 An example of an FHN dynamical system.¶

4.2.4. Rose-Hindmarsh model¶

The Rose-Hindmarsh (RH) model [RH89] is proposed to study the spike bursting phenomena in the membrane potential of a single neuron. This model is as follows:

(4.14)¶\[\dot{x}=y-ax^3+bx^2-z+I,\\ \dot{y}=c-dx^2-y,\\ \dot{z}=r\big(s(x-X_r)-z\big),\]

where $x(t)$ is dimensionless variable related to the membrane potential. Also, $y(t)$ is called the spiking variable and measures the rate of transportation of sodium and potassium ions. On the other hand, $z(t)$ corresponds to adaptation current. $a,~b,~c,~d,~r,~s,~X_r$ are model parameters. Also, $I$ is the applied current to the neuron [HR84]. This nonlinear system is solved by $\texttt{ode15s}$ Matlab function. Figure Fig. 4.15 demonstrates an example of RH dynamical model.

_images/5.png — Fig. 4.15 an example of an RH dynamical system.¶

4.2.5. Thalamic neuron model¶

The neurons in the thalamus can be modeled as [RT04]:

(4.15)¶\[\dot{v}=\frac{1}{C}\bigg(-\underbrace{g_L(v-e_L)}_{I_{L}}-\underbrace{g_{Na}m_{\infty}^3h(v-e_{Na})}_{I_{Na}}-\underbrace{g_K((0.75(1-h))^4)(v-e_K)}_{I_{K}}-\underbrace{g_T p_{\infty}^2(v-e_T)}_{I_{T}}+I\bigg),\\ \dot{h}=\frac{h_{\infty}-h}{\tau_h},\\ \dot{r}=\frac{r_{\infty}-r}{\tau_r},\]

where

(4.16)¶\[m_{\infty}=\frac{1}{1+\exp(\frac{-v-37}{7})},~~p_{\infty}=\frac{1}{1+\exp(\frac{-v-60}{6.2})},\\ h_{\infty}=\frac{1}{1+\exp(\frac{v+41}{4})},\\ \tau_h=\frac{1}{\alpha_h+\beta_h},\\ \alpha_h=0.128\exp(\frac{-v-46}{18}),~~\beta_h=\frac{4}{1+\exp(\frac{-v-23}{5})},\\ r_{\infty}=\frac{1}{1+\exp(\frac{v+84}{4})},\\ \tau_r=28+\exp(\frac{-v-25}{10.5}).\]

Here $v(t)$ is membrane voltage, and $h(t)$, $r(t)$ are the gating variables of the neuron describing the modulation of the flow of ions across the neural membrane. Also $g_L$, $e_L$, $g_{Na}$, $e_{Na}$, $g_K$, $e_K$, $g_T$, $e_T$ are the parameters of the model. On the other hand, $I_L$, $I_{Na}$, $I_K$, and $I_T$ are leak, sodium, potassium and calcium ions channels spiking currents. Moreover, $C$ and $I$ are the membrane capacitance and injected current, respectively. In order to simulate this dynamical system, we use $\texttt{ode23tb}$ Matlab function. An example for this model is provided in Figure Fig. 4.16).

_images/7.png — Fig. 4.16 An example of an Thalamic neuron dynamical system.¶

Model	PRC
HH	\(\texttt{PRC\_interpol(domain,'Hodgkin-Huxley')}\)
FHN	\(-43.41215711 \sin(\theta)\)
RH	\(\frac{1-\cos(\theta)}{2\pi}\)
Thalamic	\(\texttt{PRC\_interpol(domain,'Thalamic')}\)