Spiking neurons as predictive controllers of linear systems

The control of a generic $K$-dimensional linear dynamical system through an $N$-dimensional control signal, as illustrated in Fig. 1A, can be expressed as:

$$\mathbf{\dot{x}}=\mathbf{Ax}+\mathbf{B}\mathbf{u}\,,$$

(1)

where $\mathbf{x}\in\mathbb{R}^{K}$ is the system state, $\mathbf{A}\in\mathbb{R}^{K\times K}$ is the system matrix, describing the system’s uncontrolled dynamics, $\mathbf{B}\in\mathbb{R}^{K\times N}$ is the input matrix which maps the control input onto the system, and $\mathbf{u}\in\mathbb{R}^{N}$ is the control input.

In classical control, $\mathbf{u}$ is typically a continuous signal, as in LQR control, where the optimal control input is computed to minimize a quadratic cost function (Fig. 1B). For control systems that make use of neural networks, the control signal can instead be derived from neural activity. In spiking networks, a filtered-spike signal can be employed, where the firing rates of neurons determine the value of the control input for each timestep (Fig. 1C). Here, the spikes are aiding the network’s internal representation of the systems, but the control input still approximates a continuous control signal \citepboerlinPredictive2013, deneveEfficient2016, huangOptimizing2017, huangDynamical2018, guoNeural2021, slijkhuisClosedform2023, mooreneuron2024. In our work, we instead consider a sparse and discrete control signal, where the input is made of brief pulses or ‘kicks’ (Fig. 1D) \citepyangImpulsive1999, yangStability2007, raffObservers2007, wangImpulsive2014, ouyangImpulsive2020, petriAnalysis2024. The system then becomes:

$$\mathbf{\dot{x}}=\mathbf{Ax}+\mathbf{B}\mathbf{s}\,,$$

(2)

where $\mathbf{s}\in\mathbb{R}^{N}$ describes the on/off state of each of $N$ neurons in the network. A spike signal is emitted every time the voltage $V_{i}$ of neuron $i$ reaches its threshold $T_{i}$, resulting in the spike train $s_{i}(t)=\sum_{j}\delta(t-t_{j})$ with $t_{j}$ indexing the spike times. Importantly, here $\delta$ is the Dirac delta function, which instantaneously assumes the value one at spike times $t_{j}$ and is zero otherwise. This is different from other forms of impulsive control like bang-bang control or chattering control \citepfelgenhauerstability2003, where the instantaneous input is maintained for a certain period of time.

To formalize our control approach, we define a loss function that expresses the discrepancy between the current state of the linear dynamical system and a desired target state, by taking the squared $L^{2}$ norm of their difference. Specifically, we consider the loss:

$$L=\|\mathbf{x}(t+f)-\mathbf{z}(t+f)\|^{2}+\mu||\mathbf{s}(t)||^{2}\,,$$

(3)

where $\mathbf{x}(t+f)$ is the system state some time $f$ in the future, $\mathbf{z}(t+f)\in\mathbb{R}^{K}$ the target state, and $\mu$ is a regularization parameter that determines the spiking cost (Sec. 2.5). We consider both reactive control ($f=0$) and predictive control ($f>0$). Classically, such a loss might be minimized by taking the gradient with respect to certain system parameters, or by finding the optimal control input that minimizes this loss for some future time window, as in standard LQR control. However, the use of a spiking control signal makes both of these approaches challenging to solve mathematically. Here, we instead take inspiration from the neuroscientific theory of spike coding networks. SCN theory assumes spikes should only happen when they improve on an underlying coding loss and derives the required spiking dynamics. From this principle, spiking neural networks that represent their inputs \citepboerlinPredictive2013, deneveEfficient2016 and perform a variety of computations \citepbrendelLearning2020, guoNeural2021, slijkhuisClosedform2023 can be defined.

Following this same principle, but applied to the problem of controlling dynamical systems, we require that a spike should only happen if it improves the control loss (Eq. (3)). According to these SCNs-inspired approaches, whenever a neuron in the network spikes, it (directly or indirectly) influences the state of the plant, resulting in a different loss value. Therefore, we can compute two different losses $L_{\text{s}}$ and $L_{\text{ns}}$, representing the distance to the target in the presence or absence of a spike, respectively. In order to control the plant’s dynamics, the network compares the two losses to determine wether inputting a spike is advantageous to reach the control target. The inequality $L_{\text{s}}<L_{\text{ns}}$ is describing how our loss can be improved by the control input, decreasing the distance towards the target and helping our network solve the control task. These losses are describing a spiking condition, reminiscent of a voltage-threshold inequality $\mathbf{V}>\mathbf{T}$ (where $\mathbf{V}$ is the voltage of the neurons, and $\mathbf{T}$ is their firing threshold) that allows biological neural networks to produce action potentials. In Sec. 4.2, we show that one can indeed derive voltage and thresholds values from our loss comparison. As a consequence, in our framework, the values of both voltage and threshold are composed exclusively from quantities that describe the (current and future) states of the plant. In particular, the voltage of neuron $i$ is expressed by:

$$V_{i}=\mathbf{b}_{i}^{T}\mathbf{A}_{f}{{}^{T}}(\mathbf{\mathbf{z}}-\mathbf{A}_{f}\mathbf{x}_{0})\,,$$

(4)

where $\mathbf{b}_{i}$ is the $i$th column of the matrix $\mathbf{B}$, $\mathbf{A}_{f}\equiv e^{\mathbf{A}f}\in\mathbb{R}^{K\times K}$ is the matrix exponential of the state matrix $\mathbf{A}$ times the future time window $f$, and $\mathbf{x}_{0}=\mathbf{x}(t)$ is the current state of the system. The threshold of neuron $i$ is described as:

$$T_{i}=\frac{\mathbf{A}_{f}\mathbf{b}_{i}\mathbf{b}_{i}^{T}\mathbf{A}_{f}{{}^{T}}}{2}+\mu\,.$$

(5)

A detailed derivation of the values of $\mathbf{V}$ and $\mathbf{T}$ can be found in Sec. 4.2.4. Although a cost $\mathbf{C}$ on each dimensions of the system can be introduced to further modulate spiking activity (Sec. 4.2), the values of $V_{i}$ and $T_{i}$ in Eq. (4) and Eq. (5) do not consider this weighted case (assuming $\mathbf{C}$ is an identity matrix), for simplicity. Further developing this simple spiking rule leads to a recurrent network of integrate-and-fire neurons, illustrated in Fig. 2, and described by:

$$\mathbf{\dot{V}}=-\mathbf{V}+\mathbf{G}(\dot{\mathbf{z}}+\mathbf{z})-\mathbf{F}\mathbf{x}_{0}-\mathbf{\Omega s}\,,$$

(6)

where $\mathbf{V}\in\mathbb{R}^{N}$ represents the membrane potentials, $\mathbf{G}\in\mathbb{R}^{N\times K}$ are the forward connections that encode the target $\mathbf{z}$, $\mathbf{F}\in\mathbb{R}^{N\times K}$ contains the forward connections that encode the current system’s state $\mathbf{x_{0}}$, $\mathbf{\Omega}\in\mathbb{R}^{N\times N}$ represents recurrent connections, and $\mathbf{s}$ are the spikes. Interestingly, the obtained recurrent network has low-rank connectivity with respect to the dimensions of the downstream system (Sec. 4.2.4). This network structure is not assumed, but rather derived from each neuron as a greedy spiking condition given our loss comparison (Sec. 3.1.1). Consequently, although the network obtained is recurrent, each neuron computes the losses without access to network-level information (Sec. 3.1.1). For a full derivation of the network, see Sec. 4.2.

For this proof-of-concept, we enforce an asynchronous firing scheme: at each timestep, all neurons independently compute their local loss comparisons, but only one is allowed to spike, even if several exceed threshold. The spiking neuron can be chosen randomly from the suprathreshold set, or, as we do here, selected as the one with the highest voltage above threshold. While for clarity we keep this firing policy in all the shown examples, it is not strictly required for achieving control. This is because networks that operate on local spiking rules can distribute their activity very effectively among neurons \citepcalaimRobustness2021, which means no specific neuron is more relevant for the computation of the network. Such a simple asynchronous policy is only feasible in zero-delay systems, where the new state of the plant is inputted in the network without any time delay. In a system with delays, we could not impose this asynchronous firing. However, we can employ different methods that allow multiple neurons to spike while still being compatible with our local spiking rule, for instance, by scaling each neuron’s threshold according to its past firing activity.

We first consider a physically unconstrained system with reactive control (Fig. 4A). In this case, the system can be successfully made to track a target (Fig. 4A top right). This is achieved by continuously spiking directly towards the target, which is always the better choice since the neurons do not consider future states of the plant (as shown in Fig. 3A). In this example, the control is exerted by progressively kicking the mass in the same position as the evolving target, essentially ‘teleporting’ it without affecting its velocity. Of course, this is only possible because we are considering a physically unconstrained system, where the mass is allowed to be moved instantly in any direction of the state space. Also note that because of this property, the system ends up in a physically impossible state with a non-zero velocity (green curve in upper right plot), despite remaining in place according to the position variable. We therefore ask if a reactive control method could work on physically constrained systems. We employ reactive control to our physically constrained example, running the simulation of an SMD. In this context, applying a spiking control rule means directly influencing the velocity component of the state $\mathbf{x}$, while the position component evolves as a consequence. This reflects a key constraint of physical systems: the position cannot be instantaneously changed, as the mass cannot move from one position to another in no time. Instead, we constrain spiking control to act on the simulated physical system by applying a force to only kick the velocity, guaranteeing realistic state evolution in accordance with Newtonian mechanics. As shown in Fig. 4B, a reactive approach fails to control a physically constrained system. This is because influencing the velocity of the mass has no instant effects on its position. To appreciate these effects, the network should be able to observe future states of the plant. In short, spiking up or down in state space does not immediately affect the loss, and no spike is produced.

We introduce our predictive spiking algorithm to a non-constrained system. In this scenario, the controller is allowed to spike any dimension of the system, but it is also able to predict the plant’s future states, and it produces spikes that improve the loss $f$ seconds from the current state. From Fig. 4C, we can see how when exerting the control, the position is sometimes directly moved towards the target (as in the reactive case), but optimizing the local losses often equates to kicking the system in its velocity, exploiting the intrinsic dynamics of the plant. Our predictive algorithm exploits the intrinsic dynamics of the system even in the absence of physical constraints, although in some instances, the mass is still instantaneously moved towards the target. Finally, we apply the same predictive algorithm to the SMD simulation. In this scenario, the spiking network accounts for future states of the system, and is constrained to only spike the velocity either up or down, adjusting to the changing target position. The trajectory approaches the target and the loss $f$ seconds after a spike is minimized. If the position then departs from the target, the losses of the neurons will change their value, eventually causing the neuron to reach its threshold. At this point, the velocity will be either spiked up or down, causing the network to trace the target correctly. As shown in Fig. 4D, our predictive spiking rule introduces an element of anticipation, making the control strategy more robust in physically constrained dynamical systems, and allowing for the control of the SMD where a reactive control would fail.

Different combinations of these meta-parameters can give different results in the control of the plant. Quantifying the optimality of the control strategy in comparison to other paradigms (standard LQR, filtered-spikes) is very difficult. As described in Sec. 3.1.1, each neuron in our network follows a greedy spiking rule given the control set, but like in other approaches, this does not necessarily produce the best overall controller. It is of interest to vary some control parameters, considering the cases where the changes in target value are traced by the system’s state (and therefore the SMD is successfully controlled). Within these cases, we identified different regimes of control for our network. In these explorations we decided to maintain fixed the number of neurons ($N=2$) and the kick strength $\mathbf{B}$. We therefore vary the future time window $f$, influencing the predicted state based on which the network computes the spiking decision, and $\mu$, the global scaling factor for threshold. As shown in Fig. 6A, networks with a lower $\mu$ tend to spike more, using more energy, as defined by the effected acceleration onto the system (for this system, this equates to the sum of the total spike count of the simulation, with each spike scaled by the corresponding value in the matrix $\mathbf{B}$). Given the same $f$, these networks also tend to have a lower error compared to high $\mu$ networks (Fig. 6B), since a lower threshold allows to adjust more often to the target position. In general, it emerges from Fig. 6B how a relatively wide range of parameters allows for the control of the plant, giving comparable error measures. Only more extreme parameters combinations (a high threshold scaling, especially with small $f$) result in the network not tracing the target, yielding significantly higher cumulative error.

In our work, we derive a spiking neural network (SNNs) of leaky integrate-and-fire (LIF) neurons from control principles. A network of $N$ such neurons is described by their voltage dynamics

$$\mathbf{\dot{V}}(t)=-\lambda\mathbf{V}(t)+\mathbf{Fx}(t)+\mathbf{\Omega s}(t),$$

(7)

where $\mathbf{V}(t)\in\mathbb{R}^{N}$are the neurons’ voltages, $\lambda$ determines the membrane leak time-constant, $\mathbf{x}(t)\in\mathbb{R}^{K}$ are the inputs, $\mathbf{F}\in\mathbb{R}^{N\times K}$ are the forward weights, $\mathbf{\Omega}\in\mathbb{R}^{N\times N}$ are the recurrent weights, and $\mathbf{s}(t)\in\mathbb{R}^{N}$ are the neural spike trains. A spike is generated whenever a neuron’s voltage crosses its threshold, $T_{i}$, and a spike train is described as a sum of Dirac delta functions, $s_{i}(t)=\sum_{t_{j}}\delta(t-t_{j})$. Whenever a neuron spikes, its voltage is reset to a resting potential, $R_{i}$. We now will generally be leaving out the time notation $(t)$ for notional clarity, unless it’s important for a specific derivation step.

We want to control a linear dynamical system. We consider a $K$-dimensional linear dynamical system, over which we want to directly exert control using spike events. The system is described by the following equation:

$$\mathbf{\dot{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{s},$$

(8)

where $\mathbf{x}\in\mathbb{R}^{K}$ represents the state of the system (one value for each dimension of the system), $\mathbf{A}\in\mathbb{R}^{K\times K}$ is a matrix which defines the transformations in $\mathbf{x}$ that give rise to the dynamics (for instance, an oscillatory behavior), $\mathbf{B}\in\mathbb{R}^{K\times N}$ is a matrix that maps the spike events onto the system, and $\mathbf{s}\in\mathbb{R}^{N}$ are the spike trains, which represent our control signal.

Since it’s a linear system, we can solve it analytically (Sec. 4.1) to derive the state of the system in a future time step $f$, starting from the current state $\mathbf{x}_{0}$.

To control the system, we set a target state $\mathbf{\mathbf{z}}$ which holds one value to be reached for each dimension of the system. To quantify the distance from this target, we define a loss as the $L^{2}$ norm of the difference between the target and the predicted state:

$$L=||\mathbf{\mathbf{z}}-\mathbf{x}_{f}||^{2}+\mu||\mathbf{s}||^{2},$$

(12)

where $\mathbf{z}\in\mathbb{R}^{K}$ is the target state, and $\mu$ is a parameter that scales the cost for each spike. In the most general case, we can scale the part of this loss that describes the distance to the target. We introduce the matrix $\mathbf{C}$, which represents a cost on the dimensions of the system. This allows us to introduce a penalty on the values of the system’s dimensions, which will weigh on the input to the controller.

The loss then becomes:

$$L=(\mathbf{z}-\mathbf{x}_{f})^{T}\mathbf{C}(\mathbf{z}-\mathbf{x}_{f})+\mu||\mathbf{s}||^{2},$$

(13)

where $\mathbf{C}\in\mathbb{R}^{K\times K}$ is a cost matrix on each dimensions of the system. We can express this loss in the presence and absence of spikes.

As mentioned in Sec. 3.2, previous works conceptualized spiking neurons as capable of directly encoding dynamical variables, and derived a "greedy" spiking rule where a neuron $i$ fires whenever this results in a decrease of a cost function, which usually measures the distance between the dynamical variable $x$ and its estimate $\hat{x}$, or a target $z$. This is a fully local rule: a neuron’s decision to fire depends solely on its current membrane potential and its threshold, which in turn depend on information locally accessible to the neuron at that moment. The neuron does not need global information about the activity of the entire network. This rule guarantees that each spike emitted by a neuron contributes to improving the representation (or the control) of the dynamical variable, and has therefore a clear interpretation within the computing objective of the network. We start deriving our network based on this simple firing condition:

$$L_{\text{s}}\leq L_{\text{ns}}$$

(18)

from which:

$$\begin{split}&(\mathbf{\mathbf{z}}-\mathbf{A}_{f}\mathbf{x}_{0}-\mathbf{A}_{f}\mathbf{b}_{i})^{T}\mathbf{C}(\mathbf{\mathbf{z}}-\mathbf{A}_{f}\mathbf{x}_{0}-\mathbf{A}_{f}\mathbf{b}_{i})+\mu\\ &-(\mathbf{\mathbf{z}}-\mathbf{A}_{f}\mathbf{x}_{0})^{T}\mathbf{C}(\mathbf{\mathbf{z}}-\mathbf{A}_{f}\mathbf{x}_{0})\leq 0.\end{split}$$

(19)

Expanding the squares:

$$-2(\mathbf{A}_{f}\mathbf{b}_{i})^{T}\mathbf{C}(\mathbf{\mathbf{z}}-\mathbf{A}_{f}\mathbf{x}_{0})+\mu+(\mathbf{A}_{f}\mathbf{b}_{i})^{T}\mathbf{C}(\mathbf{A}_{f}\mathbf{b}_{i})\leq 0$$

(20)

and rearranging these terms, we obtain the inequality:

$$\frac{\mathbf{b}_{i}^{T}\mathbf{A}_{f}{{}^{T}}\mathbf{C}(\mathbf{A}_{f}\mathbf{b}_{i})}{2}+\mu\leq\mathbf{b}_{i}^{T}\mathbf{A}_{f}{{}^{T}}\mathbf{C}(\mathbf{\mathbf{z}}-\mathbf{A}_{f}\mathbf{x}_{0}).$$

(21)

On the left side of this inequality, all the factors are time independent. We can attribute these to the threshold of the neuron $i$,

$$T_{i}=\frac{\mathbf{b}_{i}^{T}\mathbf{A}_{f}{{}^{T}}\mathbf{C}(\mathbf{A}_{f}\mathbf{b}_{i})}{2}+\mu$$

(22)

On the right side, all the time dependent factors can be attributed to the voltage of that neuron:

$$V_{i}=\mathbf{b}_{i}^{T}\mathbf{A}_{f}{{}^{T}}\mathbf{C}(\mathbf{\mathbf{z}}-\mathbf{A}_{f}\mathbf{x}_{0}).$$

(23)

This gives us a spiking condition for neuron $i$: ${V_{i}}\geq T_{i}.$

If we do not pose any constraint on the dimensions of the system, we would make the matrix $\mathbf{C}$ an identity matrix $\mathbf{I}$, obtaining the voltage and threshold values described in Eq. (4) and Eq. (5) from Sec. 2.1.2.

Stepping back from the single neuron perspective, and considering instead the voltages of all the neurons in the network, we obtain:

$$\mathbf{V}=\mathbf{B}^{T}\mathbf{A}_{f}{{}^{T}}\mathbf{C}(\mathbf{\mathbf{z}}-\mathbf{A}_{f}\mathbf{x}_{0}).$$

(24)

From there, we pose:

$$\mathbf{G}=\mathbf{B}^{T}\mathbf{A}_{f}{{}^{T}}\mathbf{C},$$

(25)

obtaining:

$$\mathbf{V}=\mathbf{G}(\mathbf{\mathbf{z}}-\mathbf{A}_{f}\mathbf{x}_{0}).$$

(26)

We now consider the rate of change of $\mathbf{V}$, as well as the rate of change of all time-dependent variables:

$$\dot{\mathbf{V}}=\mathbf{G}(\dot{\mathbf{\mathbf{z}}}-\mathbf{A}_{f}\dot{\mathbf{x}}_{0}),$$

(27)

the change of the current state of the system over time $\dot{\mathbf{x}}_{0}$ can be expressed by Eq. (8), obtaining:

$$\dot{\mathbf{x}}_{0}=\mathbf{A}\mathbf{x}_{0}+\mathbf{Bs}.$$

(28)

Substituting this into Eq. (27) yields:

	$\displaystyle\dot{\mathbf{V}}$	$\displaystyle=\mathbf{G}(\dot{\mathbf{z}}-\mathbf{A}_{f}(\mathbf{A}\mathbf{x}_{0}+\mathbf{Bs}))$		(29)
		$\displaystyle=\mathbf{G}(\dot{\mathbf{z}}-\mathbf{A}_{f}\mathbf{A}\mathbf{x}_{0}+\mathbf{A}_{f}\mathbf{Bs})$		(30)
		$\displaystyle=\mathbf{G}\dot{\mathbf{z}}-\mathbf{G}\mathbf{A}_{f}\mathbf{A}\mathbf{x}_{0}-\mathbf{G}\mathbf{A}_{f}\mathbf{Bs}.$		(31)

We then add and subtract a $\mathbf{V}$ on the right side of the equation. We leave the $-\mathbf{V}$ as is, while we substitute the $+\mathbf{V}$ using the definition of $\mathbf{V}$ from Eq. (26), obtaining:

$$\dot{\mathbf{V}}=-\mathbf{V}+\mathbf{G}\dot{\mathbf{z}}-\mathbf{G}\mathbf{A}_{f}\mathbf{A}\mathbf{x}_{0}+\mathbf{G}\mathbf{z}-\mathbf{G}\mathbf{A}_{f}\mathbf{x}_{0}-\mathbf{G}\mathbf{A}_{f}\mathbf{Bs};$$

(32)

we can rewrite this as

$$\dot{\mathbf{V}}=-\mathbf{V}+\mathbf{G}(\dot{\mathbf{z}}+\mathbf{z})-\mathbf{G}\mathbf{A}_{f}\mathbf{(A+I)}\mathbf{x}_{0}-\mathbf{G}\mathbf{A}_{f}\mathbf{Bs}.$$

(33)

We can trace this back to the form of a recurrent network of integrate and fire neurons:

$$\mathbf{\dot{V}}=-\mathbf{V}+\mathbf{G}(\dot{\mathbf{z}}+\mathbf{z})-\mathbf{F}\mathbf{x}_{0}-\mathbf{\Omega s},\,$$

(34)

where $\mathbf{F}=\mathbf{G}\mathbf{A}_{f}\mathbf{(A+I)}$, and $\mathbf{\Omega}=\mathbf{G}\mathbf{A}_{f}\mathbf{B}=\mathbf{B}^{T}\mathbf{A}_{f}{{}^{T}}\mathbf{C}\mathbf{A}_{f}\mathbf{B}.$

This full network equation describes the changes in $\mathbf{V}\in\mathbb{R}^{N}$ , which represents the membrane potentials. Here, $\mathbf{G}\in\mathbb{R}^{N\times K}$ are the target input weights: forward connections that encode for the target $\mathbf{z}$. Considering the changes in target value over time $\dot{\mathbf{z}}$ as well as the actual value of the target $\mathbf{z}$ allows for the network to trace rapid changes of the target. Only encoding for the current value of $\mathbf{z}$ would still allow to control the system, but with a slower tracing of the target. The matrix $\mathbf{F}\in\mathbb{R}^{N\times K}$ represents the state input weights, that contain the forward connections to encode the current system’s state $\mathbf{x_{0}}$, $\mathbf{\Omega}\in\mathbb{R}^{N\times N}$ are the recurrent weights, and $\mathbf{s}$ are the spikes. Starting from a greedy spiking condition, we obtained a network of leaky integrate-and-fire (LIF) neurons in the form of Eq.  (7). Such network has a voltage decay term, forward weights that map the target onto the network, forward weights that encode the system state in input, recurrent connectivity, and spikes, and is able to control our linear dynamical system. For a $K$-dimensional system to control and a network of $N$ neurons, $\mathbf{B}$ is a $N\times K$ matrix, while $\mathbf{A}_{f}$ and $\mathbf{C}$ are $K\times K$ matrices. Thus, $\mathbf{\Omega}=\mathbf{B}^{T}\mathbf{A}_{f}{{}^{T}}\mathbf{C}\mathbf{A}_{f}\mathbf{B}$ is of rank $K$. Hence, our spiking control solution uses a network with low-rank connectivity (low-rank with respect to the system that is controlling).