Engineering spatiotemporal patterns: information encoding, processing, and controllability in oscillator ensembles

Abstract The ability to finely manipulate spatiotemporal patterns displayed in neuronal populations is critical for understanding and influencing brain functions, sleep cycles, and neurological pathologies. However, such control tasks are challenged not only by the immense scale but also by the lack of real-time state measurements of neurons in the population, which deteriorates the control performance. In this paper, we formulate the control of dynamic structures in an ensemble of neuron oscillators as a tracking problem and propose a principled control technique for designing optimal stimuli that produce desired spatiotemporal patterns in a network of interacting neurons without requiring feedback information. We further reveal an interesting presentation of information encoding and processing in a neuron ensemble in terms of its controllability property. The performance of the presented technique in creating complex spatiotemporal spiking patterns is demonstrated on neural populations described by mathematically ideal and biophysical models, including the Kuramoto and Hodgkin-Huxley models, as well as real-time experiments on Wein bridge oscillators.


Introduction
Effective regulation of synchronization patterns formed by a neural population is of practical importance to the treatment of neurological pathologies that often result from disruptions in these synchronization patterns. For example, motor dysfunctions in Parkinson's disease are linked to excessive neural synchronization [1,2], sleep disorders are caused by disturbances in the synchronized circadian rhythms [3], and disruption of clustered synchronous states in the brain is associated with cognitive dysfunctions [4]. Controlling neuron populations, possibly with a relatively small number of electrodes, is also critical for retinal implants, where the activation of a large number of sensory neurons in parallel and, more importantly, in the correct sequence, is required to transmit the visual information to the visual cortex [5].
Processing or encoding information using a network of neurons requires controlling the timing of neuronal spikes, which is described by complex mathematical models that make analysis and control design insurmountable. The complexity of the description of neural spiking can be significantly reduced by focusing only on the phase of the spiking dynamics, described using phase models [6]. These phase reduction approaches provide a direct link between the biophysical and computational models and weakly coupled phase oscillator models, enabling several insights into the nature of synchronous activity at the neuronal level, and have become a staple tool in computational neuroscience [7]. As a result of the effectiveness and simplicity of phase models, we utilize these simplified models for control design and analysis.
Robust and optimal control of either a single neuron or a population of neurons has been the subject of active research. This led to the development of novel control techniques for asymptotic and phase-selective entrainment of neural oscillators with periodic waveforms [8][9][10][11][12]; modulation of the spiking rates of neurons using charge-balanced and minimum-power controls [13,14]; synchronization of neural populations [15,16]; and desynchronization of globally coupled neurons [17][18][19][20][21]. Instead of regulating the whole population of neurons, selective spiking of neurons is also achieved by designing stimuli that induce firing of a subset of the population while inhibiting the rest [22]. On the other hand, the problem of suppressing hypersynchronous oscillations of the whole population during seizures is considered in [23].
The primary contribution of our work is to transform the challenging task of controlling dynamic structures in an ensemble of nonlinear oscillators into a tracking problem and to present a computationally tractable iterative algorithm to solve the formulated tracking problem. This tracking formulation does not require real-time feedback and is able to produce any desired phase configuration while simultaneously offering flexibility to accommodate the trade-off between different objectives, such as tracking errors and input energy. Moreover, we uncover the relationship between the information encoding capacity of a neuronal network and its controllability property and establish the conditions under which a network loses its encoding capacity. Specifically, we show that arbitrary phase patterns can be created by using a common control input (a stimulus) without requiring any specific coupling structure in a controllable neuron oscillator network, while only limited types of patterns are possible for a partially controllable network, which in turn reduces its information encoding capacity.

Phase reduction
The dynamics of a nonlinear oscillator, such as the Hodgkin-Huxley (HH) neuron model [24] or the limit cycle model for circadian rhythms [25], are often described by a set of coupled ordinary differential equations (ODEs) exhibiting a stable limit cycle. Consider a generic time-invariant dynamical model of an oscillating system, given by where x t n Î ( )  is the state vector and u t Î ( )  is the stimulus input. A system exhibiting an attractive, be reduced to a one-dimensional system, described by a nonlinear phase equation [11,[26][27][28] given by where θ ä [0, 2π) is the phase variable, f and z are realvalued functions, and  u Î Ì  is the external control [29][30][31] taken from the set of admissible control functions  . The function f, also referred to as the instantaneous frequency, represents the baseline dynamics of the oscillator in the absence of the control u, and z describes the response of the phase to u applied at a given phase θ and is referred to as the phase response curve (PRC) [32,33].

Dynamic control of phase patterns
We formulate the phase pattern formation problem as an optimal tracking problem. Such a formulation offers the advantage of controlling the path to the desired phase pattern while simultaneously ensuring the optimality of the designed control input. Specifically, we consider a system involving a population of n-coupled neurons, as described by the phase model )  , and u t Î ( )  is the common control input broadcast to the ensemble. The drift f(Θ) = ( f 1 (Θ), L ,f n (Θ)) contains the self-and coupling-dynamics of each oscillator in the network and Z(Θ) = (z 1 (θ 1 ), L ,z n (θ n )) with z i (θ i ) being the PRC of oscillator i for i = 1,K,n. The objective is to find an input stimulus u(t) that uses minimum energy to drive the system to follow a desired reference trajectory  t , , ) . To this end, we first construct a desired reference Θ d (t) that results in the desired phase pattern (figure 1). After that, a control input that minimizes a quadratic cost functional, a combination of control energy and tracking error e(t) = Θ d (t) − Θ(t), is obtained by solving the dynamic optimization problem, where  :  The optimal tracking problem, equation (4), can be solved using the principle of optimal control theory by considering its associated value function, V(t, θ). We derive the following solution to the minimization problem (4) (refer to the Supplementary Materials for the detailed analysis): where the optimal trajectory Θ * (t) is the trajectory obtained by driving the system using the control input in equation (5) as follows and P t gt ,  satisfy a set of coupled differential equations given by ). The set of differential equations, in equations (5), (6), (7), and (8), characterizes the optimal control solution and is in general analytically and numerically intractable. To effectively solve this set of differential equations, a computational algorithm is also developed (algorithm 1). The developed iterative algorithm is computationally tractable and finds convergent tracking controls for the desired phase pattern if the ensemble system in (3) is controllable. Algorithm 1. Algorithmic description of the iterative method.

Simulation and experimental details
We apply the proposed method to control the phase pattern formation in a network of HH neurons and Kuramoto oscillators [24,29] (refer to the Supplementary Materials for the models and parameters used for simulations). The HH neuron model describes the generation and propagation of action potentials in a squid giant axon based on the dynamic interplay between ionic conductances and electrical activity [24], and serves as a canonical example for neural ensemble dynamics. Kuramoto oscillators are also frequently used to study the interaction of cortical brain areas [7,34].
To demonstrate the applicability of our method in real-world experimental settings and robustness against measurement noise, we test our method on a network of four Wien bridge oscillators. The experimental set-up consisted of four Wien bridge oscillators, buffer circuits connected to the input and output of each oscillator, and an Arduino Due microcontroller that measured the output voltages and was also used to apply the control input to the oscillators. The experimental setup along with the measured PRC of each oscillator is depicted in figure 2. The PRCs of the oscillators were experimentally measured by applying brief pulses and observing the phase shift of each oscillator.
Given that the Arduino Due only measures positive voltages, buffers and offset circuitries were used to enable the measurement of the full waveform from each oscillator. As for the control input u(t), we used both analog outputs of the Arduino Due with a difference amplifier circuit to reproduce both positive and negative voltages. The control input was designed in Matlab and applied to the hardware using Matlab/ Simulink, which was running the Arduino in real-time ( figure 3).

Results
The proposed iterative algorithm is a unified approach to controlling limit cycle oscillators, which offers flexibility in designing optimal controls for different objectives. For example, when tracking a specified reference trajectory is not required, our method can be used to design minimum energy controls to achieve a desired final state by setting the penalty matrix Q = 0. For such objectives, an arbitrary small terminal error, j, can be achieved at the expense of input energy by tuning the matrices R and F, respectively. The flexibility to regulate the terminal error, tracking error, and applied input energy, by simply tuning the respective penalty matrices, enables our method for phase assignment [12], and control of spike timing [35,36]. Furthermore, this method can be applied for designing real-time feedback controls by pre-computing the feedback gains and storing them in lookup tables as it is done in flight control [37].

Phase pattern formation in a network of HH oscillators
Control of the spatiotemporal patterns in neural networks using the minimum control energy is desired since the application of large inputs could be harmful, e.g. a strong electric field could harm neurons or some brain tissues. We illustrate the use of the proposed iterative tracking method by considering the problem of assigning a phase difference of Δθ ref = π/2 rad between successive neurons in an ensemble of four neurons. Figure 4(a) illustrates the final phase pattern at time t f with a system of four HH neurons with natural frequencies (ω 1 , ω 2 , ω 3 , ω 4 ) = (0.4207, 0.4264, 0.4322, 0.4379) rad/ms, and the neurons are ordered from the slowest (neuron 1, red circle in figure 4(a)) to the fastest (neuron 4, purple circle in figure 4(a)). We accomplish the displayed phase pattern, figure 4(a), by formulating the phase assignment problem as a tracking problem where each neuron tracks a predefined trajectory designed such that the desired phase pattern is achieved at the time t f = 40 × T 0 (T 0 = 14.75 ms). A naive yet instinctive way to design the reference trajectories is to define a linear phase reference trajectory for the slowest neuron as θ 1, When tracking or assigning relative phase differences between heterogeneous oscillators, it is best to consider the system of phase differences, given by where Δθ j,ref is the reference trajectory of the phase difference between the j th and j 1 th -( ) oscillators. This formulation results in a control with smaller energy and improved tracking accuracy. This is evident in   , and F = 0.5I, where I is a n × n identity matrix with n denoting the system dimension, for both scenarios, i.e. whether regulating the phase of each oscillator (n = 4) or the phase differences between the oscillators (n = 3). Figure 4(e) reveals that the frequency of the slowest oscillator (red dots) remains constant after applying the control u g , unlike the case when the control u f is applied ( figure 4(d)). The underlying reason is that we define a linear phase reference trajectory θ 1,ref (t), i.e. constant frequency, for the slowest neuron while designing u g . No such constraints on the slowest oscillator's phase trajectory are placed to design u f . In the first case, u g control, the frequency of the first oscillator is kept constant, which causes the remaining oscillators to speed up more than the second case, u f control, to get to the proper phase configuration in the allocated time, t f = 590 ms.
It is a well-established fact that increasing the spiking frequency of HH neurons demands more control energy than decreasing it. This can also be observed by examining the PRCs of the HH neurons, which are asymmetric with respect to the x-axis (see  In addition to the uniform phase assignment, we consider desynchronization and cluster-formation in a population of 100 HH neurons with uniformly distributed frequencies in [0.4261, 0.4326] rad/ms (see figure 5). Specifically, first, we design a control input to uniformly desynchronize the population from identical initial phases in time t f = 40 × T 0 , where T 0 = 14.75 ms. Once the population is desynchronized, we design another control input that drives the desynchronized ensemble into 3 uniformly distributed clusters in a predefined time 12 × T 0 . The control design is carried out with 12 HH neurons, and the resulting control input is then applied to 100 HH neurons, with the results shown in figure 5(g). Specifically, we select 12 HH neurons with 4 neurons taken from each frequency group, (ω 1l , ω 1h ) = (0.42605, 0.42825) rad/ms, (ω 2l , ω 2h ) = (0.42831, 0.43011) rad/ms and (ω 3l , ω 3h ) = (0.43017, 0.43258) rad/ms, such that the boundaries of frequency groups are included in the selected neurons of each group. This is to ensure that the designed control input works as expected when applied to 100 HH neurons. The penalty matrices are taken as R = 1, Q t I t t , and F = 2I for both uniform desynchronization and cluster formation. These numerical simulations illustrate that the proposed control strategy could generate complex patterns from arbitrary initial conditions in a large population of nonlinear oscillators.

Information processing and encoding in networks of oscillators
It is well-accepted that neuronal networks encode environmental information into firing patterns [38][39][40][41][42]. For instance, olfactory systems process scent cues by encoding them into complex spatiotemporal patterns [43], and brain functions, such as memory and information processing, are determined by sequential patterns [44,45]. Therefore, control algorithms for the effective synthesis of stimuli that can create arbitrary spatiotemporal firing patterns will be crucial for neural coding and information processing. The prospect of controlling the firing sequence of neurons also has applications in the entertainment and medical industries. For example, being able to stimulate the neurons in the olfactory system to simulate scents in the era of 3D videography will enhance the entertainment experience, and the ability to activate the neurons responsible for transferring visual images to the brain will aid in the treatment of vision impairment [5,46]. Consider the all-to-all coupled network of 12 sinusoidal PRC neurons shown in figure 6(b), we design a control input that forms 3 different clusters such that neurons (1,3,5,7,9,11) are in one cluster (c1), while the neurons (2,4,8) and (6,10,12) form the remaining 2 clusters (c2 and c3, respectively) and the neurons in c2 and c3 have a phase difference of −π/4 and π/2 with the neurons in c1. The designed control input and the phase difference trajectories, within successive neurons, are shown in figures 6(e) and (a). The penalty matrices are taken as R = 4, Q t I t t 0.08 f = ( ) , and F = 10I where t f = 61.8 s. Using this synchronization structure, we portray two different schemes that could be utilized to, hypothetically, encode information into the spiking patterns of a neuronal network using the relative phase differences (figure 6(c)) or the phase of firing ( figure 6(d)). Firstly, the encoding depicted in figure 6(c), which resembles a color barcode, can be utilized to encode information as various color patterns according to the relative phases of the oscillators. Secondly, the information can also be encoded into the phase of the firing of each neuron, as shown in figure 6(d). These schemes reveal that the encoding or storage capacity of oscillatory networks remains significant, even with a small number of neurons. Now, we show that our control framework can be directly used to encode information by producing spatial patterns in a globally coupled network of 25 heterogeneous SNIPER oscillators, z a 1 cos where a i is a model-dependent constant [32]. SNIPER oscillators describe neurons near a SNIPER bifurcation (i.e. a saddle-node bifurcation on a periodic orbit) [36]. Figure 7 displays the pattern formation from a uniform initial phase distribution tothe pattern 'W' (figure 7(c), middle panel) and then switch to the 'U' pattern (figure 7(c), right panel) using the derived control sequence depicted in figure 7(d) (the two letters W and U stand for Washington University). These letter patterns are obtained by assigning the same phase, e.g. 0rad, to all the oscillators forming the letter pattern, while assigning the remaining oscillators a phase equal to π rad. The trajectories of the phase differences in the top right panel of figure 7(b) show how the π phase difference between oscillators is assigned following a linear reference trajectory for 23 seconds. Once the oscillators are controlled to the right configuration, the last 8 seconds of the control signals continued to maintain the phase differences for both patterns, respectively. The oscillator phases on a unit circle are displayed in figure 7(a), where the first panel depicts the initial oscillator phases and the second, third, and fourth panels display the uncontrolled final phases, and the final phases corresponding to pattern 'W' and 'U.  Our developed broadcast control strategy is particularly important for applications involving manipulating ensemble systems with limited control sources. For example, the number of electrodes needed can be significantly reduced for retinal implants to transmit visual information to the visual cortex. In particular, using one electrode may be sufficient to control all the neurons in a given surface area. Such capability is quantified by the controllability property, which we will elaborate on in the next section.

Controllability and information encoding
The tasks depicted in figure 6 and figure 7 consisted of designing a control capable of steering the population of oscillators from an initial state (pattern) x 0 to a desired final state (pattern) x 1 . This was possible because x 1 was reachable from x 0 by the application of a control input, as evident by the numerical simulations in these figures. The collection of points that can be reached from x 0 forms the reachable set, denoted  x 0 . Consider a nonlinear system of the form, ) . Note that if the input gains, z 1 , z 2 , and z 3 , are equal, and the natural frequencies, ω 1 , ω 2 , and ω 3 , are identical, then this network is not controllable. This can be shown by evaluating the recursive Lie brackets Z ad f k , that are  As the uncontrollability of a network leads to the reduction in the capability of forming diverse phase patterns, this may be compensated by using networks of greater size, composed of more oscillators, to enhance encoding capacity. This is analogous to a digital system that encodes information as 0ʼs and 1ʼs, requiring a sufficient number of bits to represent data. Similarly, the loss of controllability in the brain circuitry could lead to the loss of brain functions due to the inability of the neural network to reproduce the precise spiking pattern that contains the necessary information.

Real time control of Wien bridge oscillators
We apply the developed iterative algorithm to four heterogeneous Wien bridge oscillators to form two clusters, each consisting of two oscillators. For convenience, we refer to the slowest and the fastest oscillator as oscillators 1 and 4, respectively. We perform two different phase assignment tasks: For the first task, we aim to form the clusters {1, 4} and {2, 3} with π rad apart ( figure 8(a)). For the second task, we form the clusters {1, 2} and {3, 4} with π rad apart ( figure 8(c)).
The reference trajectory in these two examples had three parts: (1) a constant section that maintained the oscillators synchronized for 10 sec, (2) a ramp section that placed the oscillators in the right phase configuration (clusters), and (3) another constant section to maintain the clusters. Additionally, we apply a 5-volt resetting pulse at the start of the experiment to ensure that each oscillator has identical phases since we assume similar initial phases while estimating the control input.
Note, however, that designing the control using the phase difference equation (9) only requires the knowledge of the initial phase differences between the oscillators. In contrast, if the control design is done using equation (3), one needs to know the initial phases of the oscillators to design the control input. In the case of unknown initial phases, applying any synchronizing control or a resetting pulse will suffice.

Discussion
In this paper, we develop a unified iterative control framework for controlling spatiotemporal patterns in an ensemble of limit-cycle oscillators. We illustrate that the optimal tracking control can be precisely computed by knowing the phase dynamics of the system. We show the versatility and efficacy of the control algorithm through various numerical simulations and experiments. Further, we characterize the controllability of neuronal ensembles using tools from geometric control theory and show its implications on the information encoding capacity in neural networks.
One of the main ideas of our work is to transform the phase assignment problem into a tracking problem. By recasting the phase assignment problem as a tracking problem, we gain control over the phase trajectories at each moment and the time required to attain the desired phase pattern, which is useful for rapid cardiac resynchronization and fast jet lag recovery. Furthermore, unlike the existing methods for phase assignment [12,45], which require sufficient non-linearity of the PRC to realize the desired phase pattern, our method merely requires reachability from the initial to the target phase configuration.
The control algorithm presented here has the potential to make an impact on a wide range of applications, ranging from the treatment of neurological pathologies such as Parkinson's disease and epilepsy [50,51] to the design of neurocomputers [52,53]. The applicability of our open-loop control methodology to various applications can be attributed to two factors: first, the robustness of our methods, as proved by experimental results, and second, the fact that feedback information of individual units is not required to synthesize a control rule. The presented broadcast control strategy will also have a significant impact on clinical applications, such as retinal implants, by reducing the number of required electrodes for transmitting visual information to the visual cortex via optic neurons. Although the core of this work considered neuron oscillator phase models, the methods and analysis presented here are applicable to any limit-cycle oscillators, e.g. pulsating cardiac cells generating heartbeats, which are often modeled as reaction-diffusion systems [33]. Furthermore, due to the algorithm's iterative nature, our suggested control method is computationally tractable, making it suitable for large networks. In each iteration, the computation cost is just  n M 2 ( )to solve three ordinary differential equations, where n is the system dimension and M is the number of time steps. All the numerical experiments are implemented in Matlab on a single workstation with Xeon Gold 6144 3.5 GHz processor and 192 GB memory, and the simulation time is within 2-3 H for all the examples.
It is also essential to acknowledge and discuss some limitations of our study. First, we design the control waveform based on phase models and then apply the designed control input to the original high-dimensional system. Suppose the designed control input is not weak. In that case, the phase models will not accurately reflect the actual phase dynamics of the original system, and the output obtained from applying the control input to the original system will not be as desired. However, by properly regulating the input and state penalty matrices, sufficiently weak input can be obtained to ensure the validity of phase models at the cost of small phase assignment errors. Second, while doing controllability analysis, we neglect the connection between oscillator units. Even though our results could be applied to weakly coupled oscillator ensembles, a detailed investigation of the influence of coupling on the network's information encoding capability, i.e. controllability property, is still required, which forms the foundation for future works.