A neural machine code and programming framework for the reservoir computer

From logical reasoning to mental simulation, biological and artificial neural systems possess an incredible capacity for computation. Such neural computers offer a fundamentally novel computing paradigm by representing data continuously and processing information in a natively parallel and distributed manner. To harness this computation, prior work has developed extensive training techniques to understand existing neural networks. However, the lack of a concrete and low-level machine code for neural networks precludes us from taking full advantage of a neural computing framework. Here we provide such a machine code along with a programming framework by using a recurrent neural network—a reservoir computer—to decompile, code and compile analogue computations. By decompiling the reservoir’s internal representation and dynamics into an analytic basis of its inputs, we define a low-level neural machine code that we use to program the reservoir to solve complex equations and store chaotic dynamical systems as random-access memory. We further provide a fully distributed neural implementation of software virtualization and logical circuits

Article https://doi.org/10.1038/s42256-023-00668-8 In this Article, we provide two such programming frameworksstate neural programming (SNP) and dynamic neural programming (DNP)-by constructing an analytic representational basis of the RC neurons, alongside two architectures: open loop and closed loop. Through SNP, we program RCs to perform operations and solve analytic equations. Through DNP, we program RCs to store chaotic dynamical systems as random-access memories (dRAM), virtualize the dynamics of a guest RC, implement neural logic AND, NAND, OR, NOR, XOR and XNOR, and construct neural logic circuits such as a binary adder, flip-flop latch and multivibrator circuit. Using these circuits, we define a simple scheme for game development on RC architectures by programming an RC to simulate a variant of the game 'pong'. Finally, we decompile computations from conventionally trained RCs as analytic functions.

Open-loop architecture with SNP for output functions
We begin with the simplest architecture, which is an open-loop neural computer architecture that treats the RNN as a function approximator. We conceptualize an RNN as a computing machine comprising N neurons r, which receive k inputs x and produce m outputs o (Fig. 1a). This machine typically supports the basic instructions of multiplication of the neural states and inputs by weights A, B and W, the addition of these weighted states with each other and some bias term d, the transformation of the resultant quantities through an activation function g, and evolution in time (equations (6) and (7)). The output is given by o = Wr (Fig. 1b). Hence, the weights B, A, d and W specify the set of instructions for the RNN to run, which we conceptualize as the low-level neural machine code (Fig. 1c). Unlike conventional computers, these To bring this analogy to reality, we seek a neural network with a simple set of governing equations that demonstrates many computer-like capabilities 20 . One such network is a reservoir computer (RC), which is a recurrent neural network (RNN) that receives inputs, evolves a set of internal states forward in time and outputs a weighted sum of its states 21,22 . True to its namesake, RCs can be trained to perform fundamental computing functions such as memory storage 23,24 and manipulation 25,26 , prediction of chaotic systems 22 and model-free control 27 . Owing to the simplicity of the governing equations, the theoretical mechanism of training is understood, and recent advances have dramatically shortened training requirements by using a more efficient and expanded set of input training data 28 . But can we skip the training altogether and program RCs without any sampling or simulation just as we do for silicon computers?
These ideas have a rich history of exploration, notably including the Neural Engineering Framework 29 and system hierarchies for brain-inspired computing 30 . The former defines the guiding principles-representation, transformation and dynamics-to implement complex computations and dynamics in neuron models 31 . The latter defines an extension of Turing completeness to neuromorphic completeness, and builds an interface between neuromorphic software and hardware for program-level portability. Our work sits at the intersection of these two areas by defining an extension of the former for RCs to program computations in existing RCs with full-rank connectivity, decompile the computations performed by conventionally trained RCs and omit any sampling of state space in the optimization procedure (Supplementary Section X). Combined with the substantial advances in experimental RC platforms 32 , it is now timely to formalize a programming framework to develop neural software atop RC hardware.
Source code matrix: O c (rotation by 90°) Programming matrix: R d (discrete time) which are all treated as variables. b, The low-level instructions supported by this RNN. c, We randomly instantiate B, A and d, and compile only the output matrix W. d, To convert the machine code into an algebraic form (that is, as a function of x), we decompile the machine code into the SNP matrix R c by first approximating the RNN state as a function of the powers and the time derivatives of the inputs x (equations (8)-(10)), and then by taking the Taylor series expansion of h. Hence, the (i, j) term of R c corresponds to the Taylor series coefficient of neuron i and term j. e, This SNP matrix is used to program a desired output (for example, a rotation about the x 3 axis) by filling in the entries of the source code matrix O c , which correspond to the coefficients in front of the desired output functions. f, This program is compiled by training an output matrix W c that maps R c to O c (equation (1)). The RNN output, W c r(t), successfully rotates the three-dimensional Thomas attractor input x(t). g, The discrete-time decompiler is also composed of a Taylor series expansion, but the function h is now from equations (11) and (12), and the time derivative terms are replaced by time lags. h,i, We program highpass filters of various cut-off frequencies by weighting the appropriate time-lagged terms (h), which yields high-passed versions of the inputs (i).
Article https://doi.org/10.1038/s42256-023-00668-8 instructions are simultaneously evaluated in parallel by the global set of neurons. We randomly instantiate B, A and d, and program W.
To define a programming framework, we take the approach of translating the machine code into a representation that is meaningful to the user. We choose that representation to be the input variables x, as the inputs are usually meaningful and interpretable in most applications. This translation from the low-level to the high-level programming matrix is referred to as decompiling, and involves writing the neural states r as a function h of the inputs x given the machine code B, A and d (Methods, equations (8)-(10)). To write h in a more understandable form, we perform a Taylor series expansion of h with respect to all of the input variables x x x,ẋ x x, ⋯, thereby writing the state of every neuron as a weighted sum of linear, quadratic and higher-order polynomial terms of x x x,ẋ x x, ⋯. The weights that multiply these terms are precisely the coefficients of the Taylor series expansion, and form an N × K matrix of coefficients R c , where K is the number of terms in the expansion (Fig. 1d). This matrix of coefficients R is our SNP matrix in the open-loop architecture.
Next, we define the source code, which is the set of programmable output functions given by the rowspace of R. This is because the output of our RNN is determined by a linear combination of neural states, W, which are weighted sums of the K expansion terms from the decompiler. Here, programming refers to specifying m outputs as a weighted sum of the K terms from the decompiler, which forms an m × K matrix O c : the source code matrix. In this example, the RNN receives three inputs: x 1 , x 2 and x 3 . To program a 90° rotation about x 3 , we specify the coefficients in matrix O c for three outputs o 1 , o 2 and o 3 (Fig. 1e). Finally, to compile the source code O c into machine code W, we solve (1) When we drive the RNN with the complex, chaotic time series x(t) from the Thomas attractor ( Fig. 1f, blue), the RNN output is a rotated attractor Wr(t) (Fig. 1f, orange). This process of decompiling, programming and compiling also holds in discrete-time RNNs ( Fig. 1g; see Methods, equations (11) and (12)). Our programming matrix now consists of polynomial powers of time-lagged inputs, which allows us to program operations such as highpass filters (Fig. 1h). After compiling the program O d into W, the RNN outputs (Fig. 1i, orange, yellow and purple) filter away the lower-frequency components of an input signal (Fig. 1i,

Closed-loop architecture with SNP to solve algorithms
To increase the computational power of these RNNs, we use the same SNP, but introduce a closed-loop neural architecture. The idea is to use SNP to program some output function of the inputs as W r r r = f f f(x x x, x x x), and then feed that output back into the inputs to define an equivalence relation between the outputs and the inputs to solve x x . This feedback modifies the internal recurrent weights in the RNN, allowing it to store and modify a history of its states and inputs to store and run algorithms. Hence, the closed-loop architecture is no longer solely a function approximator. We consider the same RNN as in Fig. 1a, except now with two sets of inputs, x x x ∈ ℝ n and x x x ∈ ℝ k , and two sets of outputs,  (13)). We feed one set of outputs back into the inputs such that ō o o =x x x , which will determine the internal connectivity of the RNN. Using SNP, we program an output ō o o =Wr r r =f f f(x x x, x x x) and perform feedback as A =Ā +BW (Fig. 2b,c and equation (14)).
As a demonstration, we program an RNN to solve the continuous Lyapunov equation, where X and X are 6 × 6 matrices such that the pre-programmed RNN receives n = 36 inputs as x x x and k = 36 inputs as x. Given X, the solution X to equation (2) is important for control theory 33 and neuroscience 34 , and is often referred to as the controllability Gramian. To program equation (2) into our RNN, we first decompile the neural states r into the SNP matrix R c with respect to our input variables x x x and x (Fig. 2d). Next, we fill in matrix O c with the coefficients of all constant, linear and quadratic terms from equation (2) (Fig. 2e). Then, we compile the program O c by solving for argminW ∥WR c − O c ∥, and define the recurrent connections of a new feedback RNN, A =Ā +BW , which evolves as equation (14). By driving this new RNN with a matrix X, the output converges to the solution X of equation (2) (Fig. 2f).
As a demonstration for discrete-time systems, we program an RNN to store a substantial time history of a stochastic, non-differentiable time series x t , and perform a short-time Fourier transform. Starting with our decompiled RNN states (Fig. 2g), we write a program, Ō d , to store time history across n = 49 inputs for x x x t , by defining a sample lag operator that shifts the state of all inputs down by one index (Fig. 2h and equation (15)). Then, using this lagged history, we write another program, O d , that outputs a short-time Fourier transform (equation (16)) with a sliding window of length n. We compile these two programs by minimizing ∥WR d −Ō d ∥ and ∥WR d − O d ∥, and define the recurrent connectivity of a new feedback RNN-where A =Ā +BW -according to equation (14). By driving this RNN with a stochastic sawtooth wave of amplitude 0.2, period of eight samples and noise that is uniformly distributed about 0 with a width of 0.1, the RNN matches the true short-time Fourier transform (for a performance comparison with conventional RCs and FORCE, see Supplementary Section VII).

Closed-loop RNN with DNP to simulate and virtualize
Here we define a second, dynamical programming framework, DNP, which allows explicit programming of time history for continuous-time RNNs (for an extended discussion, see Supplementary Section VIII). Building on SNP where we decompiled the neural states r, we now decompile the activation function g, which encodes both state and dynamic information through a rearrangement of equation (6). We substitute r r r ≈ h h h(x x x, x x x) into equation (13) as Hence, our DNP decompiler takes the Taylor series coefficients of g instead of h, allowing us to program not only output functions of the input states, but also the input time history. As a demonstration, we consider an RNN with 15 states, r r r ∘ ∈ ℝ 15 , which receives three inputs x x x ∈ ℝ 3 . We will use the closed-loop architecture where Ā ∘ , B ∘ and d ∘ are randomly initialized, which is decompiled according to equation (3) (Fig. 3a). Because our decompiled code consists of analytic state and time-derivative variables, we can compile W ∘ to map RNN states to input states, and RNN time derivatives to input time derivatives. Prior work has trained RNNs to simulate time-evolving systems by copying exemplars 22 or sampling the dynamical state space 35 . Here we achieve the same simulation without any samples in a chaotic Lorenz attractor that evolves according to ẋ x Here, 'programming' refers to the construction of a matrix Ō ∘ c comprising the coefficients of x x x + 1 γ f(x x x) preceding the variables 1,x 1 ,x 2 ,x 3 ,x 1x 2 , ⋯ of the program matrix G ∘ c (Fig. 3b). Once we compile this code and perform feedback by defining new connectivity A =Ā +BW , the evolution of the RNN simulates the Lorenz attractor (Fig. 3c).
More generally, DNP allows us to program systems of the form ẋ x x = f(x x x), which raises an interesting phenomenon. We program another Dynamical decompiler (matrix G°c) ) of the guest RNN with 14 neurons and 3 inputs is decompiled into DNP matrix G ∘ c by taking the Taylor series coefficients of g ∘ as equation (3). b, Therefore, the programs we write are not output functions of x x x, but are rather given by where ẋ x x = f(x x x) for a dynamical system. c, After compiling into W ∘ and performing feedback, the evolution of the guest RNN projects onto the programmed Lorenz  RNN r r r ∈ ℝ 2000 -the host-to emulate the dynamics of the feedback RNN r r r ∘ ∈ ℝ 15 -the guest-that itself was programmed to evolve about the Lorenz attractor. We decompile the host RNN using DNP (Fig. 3d), write the code of the guest RNN in the format r r r ∘ + 1 γṙ r r ∘ (Fig. 3e), and compile the code into matrix W (Fig. 3f)

Op-codes, composition and dynamic RAM
Here we extend more of the functionality of general purpose computers to RNNs. The first functionality is support for op-codes, which is typically a string of bits specifying which instruction to run. We add control inputs c as a string of 0s and 1s such that which pushes the pre-programmed RNN to different fixed points, thereby generating a unique SNP at each point (Fig. 4a). Then, we program different matrix operations (Fig. 4b), and simultaneously compile each source code at a different SNP (Fig. 4c) into matrix W . When we drive our RNN with matrices P and Q at different c, the RNN outputs each operation (Fig. 4d). The second functionality is the ability to compose more complicated programs from simpler programs. We note that, in SNP, the output is programmed and compiled to perform an operation on the inputs, such as a matrix multiplication and vector addition for a neural processing unit (NPU, Fig. 4e). By feeding these outputs into another NPU, we can perform a successive series of feedback operations to define and solve more complex equations, such as least-squares regression (NPU, Fig. 4e).
The third functionality is the random access of chaotic dynamical memories. The control inputs c drive the RNN to different fixed points, thereby generating unique DNPs G c 1 ,c 2 (Fig. 4f). By compiling a single matrix W that maps each DNP to a unique attractor (Fig. 4g,h), the feedback RNN with internal connectivity A =Ā +BW autonomously evolves about each of the four chaotic attractors at different values of c (Fig. 4i).

A logical calculus using recurrent neural circuits
This dynamical programming framework allows us to greatly expand the computational capability of our RNNs by programming neural implementations of logic gates. While prior work has established the ability of biological and artificial networks to perform computations, here we provide an implementation that makes full use of existing computing frameworks. We program logic gates into distributed RNNs by using a simple dynamical systeṁ where a, b and z are parameters. This particular system has the nice property of hysteresis, where when z = 0.1, the value of x converges to x = 0.1, but when z = −0.1, the value of x jumps discontinuously to converge at x = −0.1 (Fig. 5a). This property enables us to program logic gates (Fig. 5b). Specifically, by defining the variable z as a product of two input variables p and q, we can program in the dynamics in equation (5) to evolve to −0.1 or 0.1 for different patterns of p and q. These logic gates can now take full advantage of existing computing frameworks. For example, we can construct a full adder using neural circuits that take Boolean values p and q as the two numbers to be added, and a 'carry' value from the previous addition operation. The adder outputs the sum s and the output carry v. We show the inputs and  outputs of a fully neural adder in Fig. 5c, forming the basis of our ability to program neural logic units, which are neural analogues of existing arithmetic logic units. The emulation of these neural logic gates to circuit design extends even to recurrent circuit architectures. For example, the set-reset (SR) latch-commonly referred to as a flip-flop-is a circuit that holds persistent memory, and is used extensively in computer random-access memory (RAM). We construct a neural SR latch using two NOR gates with two inputs, p and q (Fig. 5d). When p = 0.1 is pulsed high, the output o = − 0.1 changes to low. When q is pulsed high, the output changes to high. When both p and q are held low, then the output is fixed at its most recent value (Fig. 5d). As another example, we can chain an odd number of inverting gates (that is, NAND, NOR and XOR) to construct a multivibrator circuit that generates oscillations (Fig. 5e). Because the output of each gate will be the inverse of its input, if p is high, then q is low and o is high. However, if we use o as the input to the first gate, then p must switch to low. This discrepancy produces constant fluctuations in the states of p, q and o, which generate oscillations that are offset by the same phase (Fig. 5e).

Game development and decompiling trained RNNs
To demonstrate the flexibility and capacity of our framework, we program a variant of the game 'pong' into our RNN as a dynamical system. We begin with the game design by defining the relevant variables and behaviours (Fig. 6a). The variables are the coordinates of the ball, x, y, the velocity of the ball, ẋ ,ẏ , and the position of the paddle, x p . Additionally, we have the variables that determine contact with the left, right and upper walls as c l , c r and c u , respectively, and the variable that determines contact with the paddle, c p . The behaviour that we want is for the ball to travel with a constant velocity until it hits either a wall or the paddle, at which point the ball should reverse direction.
Here we run into our first problem: how do we represent contact detection-a fundamentally discontinuous behaviour-using continuous dynamics? Recall that we have already done so to program logic gates in Fig. 5a by using the bifurcation of the cubic dynamical system in equation (5). Here we will use the exact same equation, except rather than changing the parameter z to shift the dynamics up and down (Fig. 5a), we will set the parameter b to skew the shape. As an example, for the right-wall contact c r , we will let b = x − x r (Fig. 6b). When the ball is to the left such that x < x r , then c r approaches 0. When the ball is to the right such that x > x r , then c r becomes non-zero. To set the velocity of the ball, we use the SR latch developed in Fig. 5d. When neither wall is in contact, then c r and c l are both low, and the latch's output does not change. When either the right or the left wall is in contact, then either c r or c l pulses the latch, producing a shift in the velocity (Fig. 6c). Combining these dynamical equations together produces the code for our pong program (Fig. 6d), and the time evolution of our programmed RNN simulates a game of pong (Fig. 6e).
To demonstrate the capacity of our programming framework beyond compiling programs, we decompile the internal representation of a reservoir that has been trained to perform an operation. We instantiate a reservoir with random input matrix B and random recurrent matrix A, drive the reservoir with a sum of sinusoids-thereby generating a time series r(t)-and train the output weights W to reconstruct highpass-filtered versions of the inputs (Fig. 6f). To understand what the reservoir has learned, we decompile the reservoir state given A and B into the SNP matrix R d according to equation (12), and find that the output WR d as an analytic function of the inputs and time derivatives closely matches the true filter coefficients (Fig. 6g).

Discussion
Neural computation exists simultaneously at the core of and the intersection between many fields of study. From the differential binding of neurexins in molecular biology 37 and neural circuits in neuroscience [38][39][40][41][42] , to the RNNs in dynamical systems 43 and neural replicas of computer architectures in machine learning 44 , the analogy between neural and silicon computers has generated growing and interdisciplinary interest. Our work provides one possible realization of this analogy by defining a dynamical programming framework for RNNs that takes full advantage of their natively continuous and analytic representation, their parallel and distributed processing, and their dynamical evolution.
This work also makes an important contribution to the increasing interest in alternative computation. A clear example is the vast array of systems-such as molecular 45 , DNA 46 and single photon 47 -that implement Boolean logic gates. Other examples include the design of materials that compute 48,49 and store memories 50,51 . Perhaps of most relevance are physical instantiations of reservoir computing in electronic, photonic, mechanical and biological systems 32 . Our work demonstrates the potential of alternative computing frameworks to be fully programmable, thereby shifting paradigms away from imitating silicon computer hardware 6 , and towards defining native programming frameworks that bring out the full computational capability of each system. to different Boolean logic gates, and we program these logic gate dynamics into our RNNs. c-e, By connecting these neural logic gates, we can form neural circuits that add Boolean numbers (c), store persistent Boolean states according to a SR latch with output o and hidden variable o ′ (d), and oscillate at a fixed phase difference due to the propagation delay of inversion operations (e).
Article https://doi.org/10.1038/s42256-023-00668-8 One of the main current limitations is the linear approximation of the RC dynamics. While prior work demonstrates substantial computational ability for RCs with largely fluctuating dynamics (that is, computation at the edge of chaos 52 ), the approximation used in this work requires that the RC states stay reasonably close to the operating points. While we are able to program a single RC at multiple operating points that are far apart, the linearization is a prominent limitation. Future extensions would use more advanced dynamical approximations into the bilinear regime using Volterra kernels 53 or Koopman composition operators 54 to better capture non-linear behaviours.
Finally, we report in the Supplementary Section XI an analysis of the gender and the racial makeup of the authors we cited in a Citation Diversity Statement.

Open-loop architecture with SNP
In our framework, we conceptualize an RNN comprising N neurons r r r ∈ ℝ N , which receive k inputs x x x ∈ ℝ k and produce m outputs o o o ∈ ℝ m . This machine has weights A ∈ ℝ N×N , B ∈ ℝ N×k , and W ∈ ℝ m×N , and some bias term d d d ∈ ℝ N×1 . If the RNN evolves in continuous time, the instructions are where 1/γ is a time constant. If the RNN evolves in discrete time, these instructions are r r r t+1 = g g g(Ar r r t + Bx x x t + d d d).
We decompile the neural states r as a function h of the inputs x given the machine code B, A and d in three steps. First, we linearize the dynamics in equation (6) about a stable fixed point r * and an operating point x * to yield 1 γṙ r r(t) ≈ A * r r r(t) + g g g(Ar r r * + Bx x x(t) where A * = (dg(Ar * + Bx * + d)∘A − I). Second, because our system is now linear, we can write the neural states as the convolution of the impulse response and the inputs as Third, to obtain r(t) as an algebraic function without an integral, we perform a Taylor series expansion of this convolution with respect to t to yield x r x p x 1 x + x/γ x + x/γ x′ + x′/γ y + y/γ y + y/γ y′ + y′/γ c r + c r /γ c l + c l /γ c u + c u /γ The code matrix Ō c (scaled for visualization) is shown. b, Contact detection with the right wall is implemented using a supercritical pitchfork bifurcation by scaling the b term in equation (5) by x − x r . When x < x r , the contact variable c r goes to 0. When x > x r , a bifurcation occurs and c r becomes non-zero. c, These contact variables are used to drive an SR latch whose output is the ball's velocity. d, An RNN simulating a playable game of pong in its head. e, The colour from blue to yellow represents the evolution of time. The bottom square is the movement of the paddle, and the circle is the movement of the marker. f, Conventional training of a reservoir by first driving it with an input time series to generate the reservoir time series r(t), and then training an output matrix W to reconstruct highpass-filtered versions of the input. g, The decompiled analytic outputs of the trained reservoir closely match the true highpass filter coefficients.
We provide a detailed analytical derivation of h in the Supplementary Sections I-III, and demonstrate the goodness of the approximations in Supplementary Sections IV-VI.
To decompile discrete-time RNNs, first we linearize equation (7) about a stable fixed point r * and operating point x * : r r r t+1 = A * r r r t + g g g(Ar r r * + Bx x x t + d d d) − dg g g(Ar r r * + Bx x x t + d d d) ∘ Ar r r * ⏟⎵⎵⎵⎵⎵⎵⎵⎵⎵⎵⎵⎵⎵⏟⎵⎵⎵⎵⎵⎵⎵⎵⎵⎵⎵⎵⎵⏟ u u u(x x x t ) , (11) where A * = dg(Ar * + Bx * + d)∘A. Second, we write r t+1 as the convolved sum of inputs which we Taylor series expand to yield the N × K coefficient matrix for K expansion terms.

Closed-loop architecture with SNP
For the closed-loop architecture with SNP, we begin with the pre-programmed RNNs, 1 γṙ r r = −r r r + g g g(Ā r r r +Bx x x + Bx x x + d d d), r r r t+1 = g g g(Ā r r r t +Bx x x t + Bx x x t + d d d), (13) for continuous-time and discrete-time systems, respectively (Fig. 2b). Using SNP, we program an output ō o o =Wr r r =f f f(x x x, x x x) and perform feedback as A =Ā +BW to yield 1 γṙ r r = −r r r + g g g((Ā +BW)r r r + Bx x x + d d d), r r r t+1 = g g g((Ā +BW)r r r t + Bx x x t + d d d), (14) for continuous-time and discrete-time systems, respectively. The sample lag operator is defined as which shifts the state of all inputs down by one index. The short-time Fourier transform with a sliding window of length n is defined as

Data availability
There are no data with mandated deposition used in the manuscript or supplement. All data in the main text and Supplementary Information are generated by the code that is publicly available online.

Code availability
All figures were directly generated in MATLAB from the code available on Code Ocean, available upon publication at https://codeocean.com/ capsule/7077387/tree (ref. 55).