The nervous system encodes physical quantities as neural activity. In insects, the firing rate of sensory neurons encode the stretch of chordotonal organs (Field and Matheson, 1998) and the strain of campaniform sensilla (Zill et al., 2004), among other physical quantities. Typical robot controllers perform operations on these signals to provide meaningful information for control actions. These operations may include the subtraction of measured and reference values, differentiation or integration of error values, or gain adjustments. Neural systems perform these same operations, but in a transformed space. The exact transformation that nervous systems use is not known, but for reliable behavior, it is necessary that sensory information is mapped to neural activity in apredictable way. Thus, we map any sensory input, θ, to an applied current (7)Iapp=R⋅θ−θminθmax−θmin, where R is the “operating range” specified in the previous section. Figure 1B illustrates such a transformation within a diagram of a neural feedback loop that controls the position of a motor. The purpose of this paper is not to analyze how this particular network functions; for a detailed analysis of this network and its function, see Szczecinski et al. (2017b). Instead, the purpose is to show how R and other functional values (time constants, gains, etc.) may be used to constrain neural and synaptic parameter values.

Figure 1B (i,ii) graphically illustrate the mapping in equation (7), and Figure 1B (iii) graphically illustrates the inverse relationship (i.e., neural value to mechanical value). If a sensory neuron has only this applied current and leak current, equation (6) shows that (8)CmdUdt+U=R⋅θ−θminθmax−θmin.

This means that the sensory neuron acts as a low-pass filter with time constant τ = C_m. It is trivial to show that, when the neuron is at equilibrium (i.e., dU/dt = 0), (9)U∗=R⋅θ−θminθmax−θmin, where the superscript “∗” specifies the equilibrium value. (Throughout this manuscript, the equilibrium activation of neuron U will be referred to as U*, and the neuron itself will be referred to as U.) Equation (9) means that the neuron’s activation above its rest potential encodes the sensory signal. In addition to perceiving sensory information, commands must be issued in the same transformation. Thus, we map the commanded sensory quantity, θ_comm, to the commanded neural activation, U_comm, with the inverse function of equation (9), (10)θcomm=θmin+UcommR⋅(θmax−θmin).

In this way, the nervous system may specify an intended motion, such as the rotation of a joint, encoded in neural activity. In our synthetic nervous systems, R specifies how mechanical quantities and neural activation are related. Thus, the tuning of every functional subnetwork described in this work relies on R, which the designer specifies before tuning the rest of the network. Two other parameters are critical for tuning these subnetworks: the amplification of synaptic transmission, k_syn (discussed in Sec. 3.1), and the synaptic reversal potential, ΔE_s. From these values, biological parameters such as synaptic conductance and neural tonic drive can be directly calculated. This makes network design intuitive, enabling the designer to select biological parameter values based on functional ones.

The goal of a signal transmission pathway is to cause the postsynaptic neuron’s voltage to be some ratio of the presynaptic neuron’s voltage. We call this ratio k_syn. The U_pre,i terms in the denominator of the right hand side of equation (13) mean that k_syn changes as U_pre changes, so we approximate k_syn as U_post/U_pre when the presynaptic neuron is fully activated (i.e., U_pre = R). The steady-state response of a neuron with a single synaptic input and no applied current can be written based on equation (13), as below: (14)Upost∗=gsR⋅Upre⋅ΔEs1+gsR⋅Upre.

To find k_syn for this synapse, we first divide both sides of equation (14) by U_pre, (15)Upost∗Upre=gsR⋅Upre⋅ΔEsUpre⋅1+gsR⋅Upre.

Next, we want to find k_syn, which can be calculated for any value of U_pre. To simplify analysis and improve the clarity of this derivation, we set find k_syn when U_pre = R. Then, we show how to set parameter values to keep k_syn nearly constant, even as U_pre changes. Making this substitution, (16)Upost∗R=ksyn=gsR⋅R⋅ΔEsR⋅1+gsR⋅R.

Finally, reducing R/R terms reveals (17)ksyn=gs⋅ΔEsR⋅(1+gs).

Rearranging to solve for g_s, (18)gs=ksyn⋅RΔEs−ksyn⋅R.

Because g_s must be positive, and the numerator of equation (18) is always positive, equation (18) is also subject to the constraint (19)ΔEs>ksyn⋅R.

Equation (18) will be used to tune addition, subtraction, multiplication, and division networks (Secs. 3.3 through 3.6).

We may also use synapses to modulate a neuron’s sensitivity to other inputs. Based on equation (13), the steady-state response of a neuron with only an applied current I_app is simply (20)Upost∗=Iapp, if we set G_m = 1. For example, this is the case for a sensory neuron that receives applied current proportional to a sensor’s state, such as a joint angle (Figure 1B), muscle stretch, or touch sensor. However, the nervous system may need to actively increase or reduce the sensitivity of the sensory neuron depending on context. Hyperpolarizing or depolarizing the neuron, however, would cause sensory information to be truncated (i.e., V_pre < E_lo). We can change the sensitivity of this neuron without losing sensory information by adding a synaptic input to the response from equation (20): (21)Upost∗=gsR⋅Upre⋅ΔEs+Iapp1+gsR⋅Upre.

To quantify how U_pre modulates Upost∗ for a given I_app, we introduce the parameter c_syn, which quantifies this degree of modulation. We define c_syn as Upost∗∕Upre, the same as k_syn, but with the understanding that U_pre will decrease Upost∗ in this case. Dividing both sides of equation (21) by U_pre and using the definition of c_syn, (22)Upost∗Upre=csyn=gsR⋅Upre⋅ΔEs+IappUpre⋅1+gsR⋅Upre.

As in the previous section, we will solve for c_syn when U_pre = R to simplify analysis. Making this substitution and reducing R/R terms, (23)csyn⋅R=gs⋅ΔEs+R1+gs.

Multiplying both sides by the denominator of the right hand side and expanding, (24)csyn⋅R+csyn⋅R⋅gs=gs⋅ΔEs+R.

Collecting g_s terms on the left hand side, (25)csyn⋅R⋅gs−gs⋅ΔEs=R−csyn⋅R.

Solving equation (25) for g_s, (26)gs=csyn⋅R−RΔEs−csyn⋅R.

Just as in Sec. 3.1, g_s > 0 depends only on R, which the designer specifies beforehand, ΔE_s, which is limited by biological constraints, and c_syn, which the designer picks based on network function. ΔE_s should be negative, and as close to 0 as possible to minimize hyperpolarization of the postsynaptic neuron. Equation (26) will be used to tune division and multiplication networks (Secs. 3.5 and 3.6).

A subnetwork that approximates linear addition of the form Upost∗=ksyn⋅(Upre,1+Upre,2) may underlie positive feedback mechanisms, which increase motor neuron activation proportional to sensory inputs such as force sensing organs (Zill et al., 2004), or used to sum sensory signals from different body segments (Mittelstaedt, 1957). We construct such a network by using two Signal Transmission pathways as presented in Sec. 3.1.

Let us rewrite equation (13) here, for clarity: (27)U∗=∑i=1ngs,iR⋅Upre,i⋅ΔEs,i+Iapp1+∑i=1ngs,iR⋅Upre,i.

This equation shows U_pre,i in both the numerator and denominator. To capture addition, we wish to minimize the impact of U_pre,i on the denominator. This is accomplished by minimizing g_s. However, if g_s = 0, then the network will not function at all. Therefore, we instead maximize ΔE_s, which yields a small g_s (equation (18)). Mathematically, there is no limit on ΔE_s, but synaptic potentials are limited in biological systems. In our work, we choose the reversal potential of calcium (E_s = 134 mV), which yields ΔE_s = E_s − E_r = 134 − (−60) = 194 mV, and specify R = 20 mV. To design a pathway where k_syn = 1, for example, we plug these values into equation (18), which gives g_s = 115 nS. The contour plots in Figure 2A show that the network matches the ideal behavior very closely over the operating range U_sum ∈ [0, R]. These design constraints are summarized in Table 1.

A subnetwork that approximates linear subtraction of the form Upost∗=ksyn⋅(Upre,1−Upre,2) may underlie negative feedback mechanisms, which are important for controlling many parameters in locomotion (Pearson, 1993; Peterka, 2003; Buschmann et al., 2015). Just as in the previous section, equation (18) is used to find g_s for each pathway.

Designing a subtraction network requires that we pay attention to how the two synapses affect one another. Since the reversal potentials of hyperpolarizing ion channels are not much more negative than typical resting potentials, larger g_s,2 values are required to transmit information than for depolarizing ion channels. This makes it harder to minimize g_s like we did in the previous section. Equation (13) enables us to constrain g_s,2 such that when U_pre,1 = R and U_pre,2 = R, Upost∗=0. Starting with the neuron response in equation (13) for two synaptic currents and no applied current, (28)Upost∗=gs,1∕R⋅Upre,1⋅ΔEs,1+gs,2∕R⋅Upre,2⋅ΔEs,21+gs,1∕R⋅Upre,1+gs,2∕R⋅Upre,2.

Substituting in U_pre,1 = R, U_pre,2 = R, and Upost∗=0, (29)0=gs,1∕R⋅R⋅ΔEs,1+gs,2∕R⋅R⋅ΔEs,21+gs,1∕R⋅R+gs,2∕R⋅R (30)0=gs,1⋅ΔEs,1+gs,2⋅ΔEs,21+gs,1+gs,2 (31)0=gs,1⋅ΔEs,1+gs,2⋅ΔEs,2 (32)gs,2=ΔEs,1ΔEs,2⋅−gs,1.

Substituting equation (18) for g_s,1, (33)gs,2=ΔEs,1ΔEs,2⋅−ksyn⋅RΔEs,1−ksyn⋅R.

To be physically realizable, g_s,2 > 0. Because g_s,1 > 0 and ΔE_s,1 > 0, g_s,2 > 0 if and only if ΔE_s,2 < 0. Thus, it is critical that ΔE_s,2 < 0.

Just as for the addition network, we minimize g_s,1 by maximizing ΔE_s,1. If R = 20 mV and k_syn = 1, then g_s,1 = 115 nS and ΔE_s,1 = 194 mV. To tune g_s,2, we first select ΔE_s,2 = −40 mV, then we solve equation (33) to find g_s,2 = 558 nS. These design constraints are summarized in Table 1, and Figure 2B graphically shows the accuracy of the subtraction network.

A subnetwork that approximates division of the form (34)Upost∗=Upre,11+1−csyncsyn⋅R⋅Upre,2 replicates the function of GABA synapses that regulate activity in the brain. A key reason for this behavior is that the reversal potential of GABA-ergic synapses is about equal to the resting potential of the postsynaptic neuron (Trappenberg, 2009). Equation (26) is used to find g_s for the division pathway.

The synapse from U_pre,1 to U_post is tuned as an excitatory Signal Transmission pathway with k = 1, as in Sec. 3.1. In our work, R = 20 mV, ΔE_s,1 = 194 mV, and equation (18) tells us that g_s,1 = 115 nS. Such a small g_s ensures that the signal from U_pre,1 to U_post is transmitted without greatly affecting the sensitivity of U_post to inputs. That is, the effect of U_pre,1 on the denominator of Upost∗ is very nearly 0.

The synapse from U_pre,2 to U_post is tuned as a Signal Modulation pathway, as analyzed in Sec. 3.2. Setting ΔE_s,2 = 0 will eliminate U_pre,2’s influence on the numerator of Upost∗. Substituting this case into equation (26) and reducing, (35)gs,2=1−csyncsyn, where Upost∗=csyn⋅R when U_pre,1 = U_pre,2 = R, their maximal value. Equation (35) also reveals that since g_s,2 > 0, 0 < c_syn < 1.

The steady-state response of the network is the result of these two synaptic inputs, as written in equation (28). Substituting equation (35), and specifying that k_syn,1 = 1, Upost∗ simplifies to (36)Upost∗=gs,1∕R⋅ΔEs,1⋅1Upre,1+gs,2∕R⋅ΔEs,2⋅0Upre,21+gs,1⋅0Upre,1∕R+1−csyncsyn⋅Upre,2∕R≈Upre,11+1−csyncsyn⋅R⋅Upre,2

In our network, we wished Upost∗=1 when U_pre,2 = R, so we set c_syn = 1/R = 0.05, which makes g_s,2 = 19 μS. When c_syn is close to 0, U_pre,2 can strongly reduce U_post’s sensitivity to inputs. When c_syn is close to 1, U_pre,2 can only weakly reduce U_post’s sensitivity to inputs. Figure 2C shows that this network performs the intended division of the signals. Table 1 summarizes these design constraints.

A subnetwork that approximates multiplication of the form Upost∗=Upre,1⋅Upre,2∕R can be used to control the gain of a sensory feedback loop, a frequently observed characteristic of neural systems that control locomotion (Cruse, 1981; Gabriel and Büschges, 2007) and posture (Peterka and Loughlin, 2004).

A multiplication network can be assembled by replacing the Modulatory Pathway in the division network with two identical Modulatory Pathways in series, connected into a disinhibitory network (see Figure 2D). This works because the product of two numbers, a⋅b = a/(1/b). However, tuning the Modulatory Pathway for the multiplication network differs from tuning the division network. This is because the right-side pathway of the network in Figure 2D must make Upost∗=0, no matter how active U_pre,1 becomes (because a⋅0 = 0, no matter the value of a). Thus, according to equation (22), c_syn = 0, unlike the division network, for which 0 < c_syn < 1. Solving equation (26) when c_syn = 0 reveals that (37)gs,2=−R∕ΔEs,2.

To solve for g_s,2, we must first select ΔE_s,2. If ΔE_s,2 = 0 like for the division network, then equation (37) divides by 0. If ΔE_s,2 > 0, then g_s,2 < 0, which is physically not realizable. Therefore, we must choose a value ΔE_s,2 < 0. The more negative ΔE_s,2 is, the more small-amplitude signals are clipped; however, the less negative it is, the larger g_s,2 must be. Therefore, g_s,2 is the limiting factor to maintain biological realism. We have chosen g_s,2 = 20 μS and R = 20 mV, making ΔE_s,2 = −1.

Now that we have designed one of the Modulatory synapses, we can calculate the response of the complete multiplication network seen in Figure 2D, which includes two identical Modulatory Pathways in series. When U_pre,2 is inactive, then it does not inhibit U_inter, which is tonically active. In this case, U_inter’s activity completely desensitizes U_post to inputs. When U_pre,2 is active, then it inhibits U_inter. In this case, U_inter is hyperpolarized, and cannot desensitize U_post to inputs. To show that this is the case, let us find the full response of the system. We first calculate Uinter∗, which has one Modulatory Pathway input and a tonic applied current I_app = R. Its response is the same as in equation (21), with the constraint from equation (37), which causes terms to cancel: (38)Uinter∗=gs,2R⋅Upre,2⋅−Rgs,2+R1−Upre,2ΔEs,2=R−Upre,21−Upre,2ΔEs,2.

U_post has two presynaptic neurons, U_pre,1 and U_inter. The synapse from U_pre,1 is a Signal Transmission synapse, and the synapse from U_inter is a Signal Modulation synapse. Its response is found via equation (13), (39)Upost∗=gs,3R⋅Uinter⋅ΔEs,3+gs,1R⋅Upre,1⋅ΔEs,11+gs,3R⋅Uinter+gs,1R⋅Upre,1.

We showed in Sec. 3.3 that equation (18) can be used to design a synapse that transmits the presynaptic neuron’s activity to the postsynaptic neuron, while minimizing its impact on the denominator of the postsynaptic neuron’s steady-state response, Upost∗. This enables us to approximate U_pre,1’s effect on Upost∗ as an applied current I_app ≈ U_pre,1. Making this substitution in equation (39), (40)Upost∗≈gs,3R⋅Uinter⋅ΔEs,3+Upre,11+gs,3R⋅Uinter.

Because we previously specified that the Modulatory Pathways are identical, we can apply the constraint from equation (37), (41)Upost∗=Upre,1−Uinter1−UinterΔEs,3.

We can now substitute equation (38) for U_inter, (42)Upost∗=Upre,1−R−Upre,21−Upre,2ΔEs,21−1ΔEs,3⋅R−Upre,21−Upre,2ΔEs,2.

This expression can be simplified. First, as noted previously, synapses 2 and 3 are identical, so ΔE_s,2 = ΔE_s,3 = ΔE_s. Second, we can multiply the first term in both the numerator and denominator by the factor (1 − U_pre,2/ΔE_s), which enables us to combine terms. Performing these simplifications, (43)Upost∗=Upre,1−Upre,1⋅Upre,2∕ΔEs−R+Upre,21−Upre,2∕ΔEs−R∕ΔEs+Upre,2∕ΔEs, (44)Upost∗=−Upre,1⋅Upre,2∕ΔEs+Upre,1+Upre,2−R1−R∕ΔEs.

Equation (44) contains a lot of information about how the multiplication network functions. First, U_post’s response indeed contains a term that multiplies U_pre,1 and U_pre,2. When ΔE_s = −1, then Upost∗scales with U_pre,1⋅U_pre,2 in a 1:1 fashion. Second, the numerator will be ≤ 0 if either U_pre,1 = 0 or U_pre,2 = 0, Upost∗≤0. This is because U_pre,1 and U_pre,2 must each be less than or equal to R. If either input is greater than R, then their synaptic inputs to U_post will saturate (see equation (4)), preventing this condition from being violated. Third, the denominator does not depend on the input values. Technically, because of the approximation made in equation (40), the denominator does change slightly with U_pre,1. However, with our chosen values of R (20), ΔE_s (−1), and g_s,1 (0.115), this change is less than 1%, justifying this approximation. Figure 2D demonstrates that this network multiplies the two inputs.

Table 1 summarizes the function, component pathways, constraint equations, and free parameters of each network from this section. This analysis enables direct construction and parameter selection for functional subnetworks that can be assembled into more complex networks capable of performing real-time robotic control (e.g., Szczecinski and Quinn (2017) and Szczecinski et al. (2017a)). Additionally, one of the key advantages to using dynamic neural systems for motor control is the handling of time varying signals. The next section examines how the dynamics of these neurons can be exploited to perform calculus on signals.

One dynamic response neural systems are known to utilize is differentiation of signals. Early examination of neural networks led to the discovery of the Reichardt detector network (Reichardt, 1961), an autocorrelation network with delays that approximates the differential of an incoming signal. Other examples include human balance, which relies on feedback proportional to the position, velocity, and acceleration of the center of mass (Peterka, 2003; Safavynia and Ting, 2012). Also, positive velocity feedback plays an important role in insect muscle control (Cruse, 1981).

We have developed differentiation networks based on the Reichardt detector network, shown in Figure 3A. We can understand its function by examining a neuron’s response to a ramp input, I_app = A⋅t, where A is an arbitrary slope of the ramp. The response of the network should be a step with a magnitude proportional to A, as shown in Figure 3B. Inserting this applied current into equation (6), a single neuron’s response is (45)Cm⋅dUdt=−U+A⋅t (46)Cm⋅dUdt+U=A⋅t.

The response of the neuron, U(t), is the sum of the particular and homogeneous solutions to equation (46), U_p(t) and U_h(t), respectively. Simulating the dynamics of equation (46) suggests that the particular solution is a ramp of slope A, which lags behind the input with a time constant C_m. To confirm this, we can substitute a candidate solution and its derivative into equation (46), and check for equality. The result is the particular (i.e., steady-state) response, (47)Up(t)=A⋅(t−Cm)

This means that if the same I_app were injected into neurons with different C_m values, and then their outputs were subtracted from one another with a network from Sec. 3.4, the network would perform a finite-difference approximation of the derivative of I_app, once the transient response decays (illustrated in Figures 3A,B).

Calculating the homogeneous solution, U_h(t), informs us how quickly the transient response decays. The homogeneous solution to first-order linear equation like equation (46) is well-known, U_h(t) = b⋅exp(− t/C_m). The constant b is found by plugging the initial condition into the full response, U(t) = U_p(t) + U_h(t), (48)b=A⋅Cm.

To tune this network, the response of U_post is written as the difference between neuron U_pre,1 with C_m,1 and neuron U_pre,2 with C_m,2 > C_m,1, (49)Upost(t)=Upre,1(t)−Upre,2(t)=A⋅t−A⋅Cm,1⋅1−exp(−t∕Cm,1)−A⋅t−A⋅Cm,2⋅1−exp(−t∕Cm,2).

Canceling the terms that are linear in t and expanding, (50)Upost(t)=A⋅(Cm,2−Cm,1)+A⋅Cm,1⋅exp(−t∕Cm,1)−Cm,2⋅exp(−t∕Cm,2).

Properly tuning a differentiator network requires tuning C_m,1 and C_m,2 to obtain the intended gain of the network, k_d, and an appropriately high cutoff frequency, ω_c. Equation (50) reveals how these may be tuned. First, the steady-state response of this network to a ramp input defines k_d = (C_m,2 − C_m,1). Second, the cutoff frequency ω_c = 1/τ_d quantifies the frequency of incoming signals (i.e., I_app = A⋅sin(ω⋅t)) above which the network’s response has less than half the energy of a lower-frequency signal. This is especially useful because although differential calculations amplify high-frequency noise, this network filters out noise with a frequency ω > ω_c. Because C_m,2 > C_m,1, the time constant τ_d = C_m,2.

Figure 3C shows contours of k_d and τ_d as C_m,1 and C_m,2 change. The plots show that increasing C_m,2 relative to C_m,1 increases k_d, which may be valuable for amplifying signals. However, this also increases τ_d, making ω_c impractically low, which will cause the network’s output to lag behind the input substantially. The contour for k_d = 1 is drawn on the contour of τ_d, showing that the smallest τ_d achievable for this gain value is 1,000 ms, which would filter out all incoming signals for which ω > ω_c = 1/(1 s) = 1 rad/s (0.159 Hz).

We can gain further insight into tuning τ_d using our FeedbackDesign tool (Szczecinski et al., 2017b). Figure 3D shows Bode plots for this network’s response, given two different values for C_m,2. When C_m,2 = 1,000 nF, like in Figure 3B, the network functions properly for inputs with ω < 1 rad/s, as predicted in the previous paragraph. Lowering C_m,2 to 50 nF increases ω_c to 20 rad/s (3.18 Hz). Lowering C_m,2 also lowers the magnitude response as a function of ω, that is, it decreases k_d. To regain this lost gain, we may increase k_syn in the subtraction network. Figure 4 shows simulation data that explores this tradeoff. Table 2 lists how to use τ_d and k_d to tune the entire differentiation network.

Our neuron model is a leaky integrator, which means that the membrane voltage will integrate an applied current, but “leak” current to return to its resting potential. As a result, data cannot be stored in individual neurons, because neurons only have one stable equilibrium point. A network that is constructed to have a marginally stable equilibrium curve (or subspace) will not leak. A network will have this property if the determinant of the Jacobian matrix is 0, or in other words, if it is not full rank (Khalil, 2002). Instead of leaking, it will maintain its activation when no external currents are applied; when currents are applied, the state of the system will change continuously. This is analogous to the position of a box on a table with friction; it will remain wherever it is placed indefinitely, unless an external force is applied. In this section, we expand on previous work (Szczecinski et al., 2017a) to show how to construct a network that is marginally stable by applying constraints to reduce the rank of its Jacobian matrix; demonstrate that such a network can be used to integrate signals over time; and relate the integration rate, k_i, to the parameter values of the network, such that U1˙=ki⋅Iapp.

Marginally stable networks are hypothesized to play an important role in navigation (Haferlach et al., 2007) and the regulation of muscle forces in posture (Lévy and Cruse, 2008). Some memory models use carefully tuned self-excitation to cancel the leak current with excitatory synaptic current (Seung et al., 2000). In a similar vein, our network uses self-disinhibition (Figure 5A) to produce a line-attractor network in which a continuum of marginally stable equilibrium states exist. Simulation data in Figure 5B shows that stimulating U₁ with an applied constant current u causes U₁ to increase at an apparently constant rate, and when u is removed, neither U₁ nor U₂ leak to their rest potentials. This is the behavior of an integrator, as described in the previous paragraph.

Let us write the response of the integrator network as shown in Figure 5A to find its equilibrium states. Each neuron has leak current, synaptic current, and a constant applied current. Let all parameter values be symmetrical between the two neurons. We make the same substitutions as before; U = V − E_r, E_r = E_lo, ΔE_s = E_s − E_r, and R = E_hi − E_lo. If I_app = R, (51)Cm⋅dU1dt=−U1+gs⋅U2R⋅(ΔEs−U1)+R (52)Cm⋅dU2dt=−U2+gs⋅U1R⋅(ΔEs−U2)+R.

Moving dynamical terms to the left hand side, and applied current to the right hand side, (53)dU1dt+1CmU1−gs⋅U2R⋅(ΔEs−U1)=RCm (54)dU2dt+1CmU2−gs⋅U1R⋅(ΔEs−U2)=RCm

Solving equation (53) when dU₁/dt = 0 reveals the equilibrium curve (55)U2=R⋅(U1−R)gs⋅(ΔEs−U1).

Solving equation (54) when dU₂/dt = 0 reveals the equilibrium curve (56)U1=R⋅(U2−R)gs⋅(ΔEs−U2), which can be algebraically rearranged to be the same as equation (55) as long as g_s and ΔE_s are constrained such that (57)gs⋅ΔEs=−R.

Multiplying both sides of equation (56) by the denominator of the right hand side, and expanding, (58)gs⋅ΔEs⋅U1−g⋅U1⋅U2=R⋅U2−R2.

Collecting multiples of U₂ and applying equation (57), (59)U2=R⋅(U1−R)gs⋅(ΔEs−U1).

Thus, equations (55) and (56) are the same equilibrium curve if g_s and ΔE_s satisfy equation (57). This curve, drawn on the phase-space diagram in Figure 5C, describes every equilibrium state that this network can have. In other words, a [U₁, U₂] pair is an equilibrium state of the system if and only if it satisfies equation (55). In the coming paragraph, we will use eigenvalue analysis to show that this network always functions as an integrator, as long as equation (57) is satisfied.

To find the system’s eigenvalues, let us write equations (53) and (54) together in matrix form, (60)U1˙U2˙+1Cm⋅1+U2⋅gsR−gsR⋅(ΔEs−U1)−gsR⋅(ΔEs−U2)1+U1⋅gsR⋅U1U2=1Cm⋅RR, in which the square matrix is J, the system Jacobian. Because J contains U₁ and U₂ terms, it is not constant, but still describes the stability of the system, given specific values of U₁ and U₂. To construct a marginally stable equilibrium subspace for the network, we must show that J has insufficient rank (i.e., the rows are identical) when U₁ and U₂ are at equilibrium (i.e., equation (55) is satisfied). However, the rows are identical, no matter the values of U₁ and U₂, if we apply the constraint from equations (57) to (60), (61)U1˙U2˙+1Cm⋅1+gsR⋅U21+gsR⋅U11+gsR⋅U21+gsR⋅U1⋅U1U2=1Cm⋅RR.

Thus, the system will always have one null direction, and we do not need to calculate J for specific equilibrium conditions to determine the system’s stability. To make notation more compact, let us define (62)a=1+gs∕R⋅U1 (63)b=1+gs∕R⋅U2

These expressions let us write equation (61) as simply (64)U1˙U2˙+b∕Cma∕Cmb∕Cma∕Cm⋅U1U2=R∕CmR∕Cm.

Plotting the simulation data of the network’s forced response from Figure 5B on a phase-space diagram (Figure 5C) suggests that u causes U₁ and U₂ to change in such a way that the state of the system (x→(t), blue) moves tangent to the equilibrium curve (dashed violet), with some constant distance away from it. These curves do not overlap because the forced response is not the same as the equilibrium condition while the external current u is applied. Motion in the X₂ direction is resisted by the neural dynamics, much how a spring resists the translation of an object with an applied force.

Nonetheless, these direction-dependent responses suggest that the state can be generalized into two decoupled degrees of freedom in the phase-space: unresisted, marginally stable motion parallel to the equilibrium curve (X₁, green in Figure 5C); and resisted, stable motion away from the equilibrium curve (X₂, red). The natural coordinates, x→=[U1,U2]T, are transformed into generalized coordinates, q→=[q1,q2]T, by a matrix X comprised of the eigenvectors of J. This same transformation matrix is used to transform J into the generalized coordinate system, yielding J_q. J_q is diagonal, decoupling the dynamics of the generalized coordinates and enabling us to quantify how quickly x→ moves parallel to the equilibrium curve.

Appendix shows the calculation of X, with q₁ representing the marginally stable mode and q₂ representing the stable mode. Using X, we can transform the system into generalized coordinates. First, we write the dynamics from equation (64) in a compact format. (65)x→˙+Jx→=F→, where J is the square matrix in equation (64) and (66)F→=R∕Cm+u∕CmR∕Cm.

The generalized coordinates, q→, are defined as (67)x→=Xq→.

To transform equation (65) into generalized coordinates, premultiply both sides of equation (65) by X⁻¹, (68)q→˙+Jqq→=Q→, where J_q = X⁻¹JX and Q→=X−1F→. The top and bottom rows of equation (68) are decoupled because J_q is a diagonal matrix. Furthermore, Jqi,i=λi, meaning that Jq1,1=0, so the system simplifies even further.

To find the particular solution of this system, we can guess the form of q_p,1 and q_p,2, and substitute those in to equation (68). We observe that q˙1(t)=B⋅u in steady state, where B is a constant that relates q˙1(t) and u. q₁(t) would be the integral of q˙1(t), but because the top row of J_q is zeros, it will not appear in the particular solution, and thus need not be explicitly included. We also observe that q˙2(t)=0 in steady state, so q₂(t) = D, a constant. We can calculate Q→=X−1F using X⁻¹, which is calculated in Appendix (equation (A9)). Solving for the particular solution of this system, q→˙p(t), (69)q→˙p(t)+Jqqp→=B⋅u0+000a+bCm⋅qp,1D=a⋅dCm⋅(a+b)⋅uR2Cm+b2Cm⋅(a+b)⋅u. (70)B=a⋅dCm⋅(a+b), where d is defined in equation (A5). B describes how quickly q_p,1 varies with u, but we want to know how quickly U₁ varies with u. Therefore, we use equation (67) to transform q˙→p=[B⋅u,0]T into natural coordinates to find x→˙, (71)xp→˙=Xqp→˙ (72)U˙1,p(t)U˙2,p(t)=1∕d1∕2−b∕(ad)1∕2⋅a⋅dCm⋅(a+b)⋅u0 (73)U˙1,p(t)=aCm⋅(a+b)⋅u (74)ki=aCm⋅(a+b).

Recall that a and b are functions of U₁ and U₂, respectively. This means that k_i, the integral gain of the network, is not a constant. To place bounds on k_i, let us substitute equations (62) and (63) into equation (74), (75)ki=1+gs∕R⋅U1Cm⋅(2+gs∕R⋅(U1+U2)).

We can now plug in different values of U₁ and U₂ to see how k_i varies. Using equations (55) and (56), we find that the most extreme cases are when [U₁, U₂] = [0, R] and [U₁, U₂] = [R, 0]. We can plug these cases into equation (75) to find the minimum and maximum values for k_i, (76)ki,min=1+gs∕R⋅0Cm⋅(2+gs∕R⋅(0+R))=1Cm⋅(2+gs) and (77)ki,max=1+gs∕R⋅RCm⋅(2+gs∕R⋅(R+0))=1+gsCm⋅(2+gs).

The difference between k_i,min and k_i,max: (78)ki,range=1+gsCm⋅(2+gs)−1Cm⋅(2+gs)=gsCm⋅(2+gs).

To find the mean rate of integration, we can calculate k_i,mean = (k_i,min + k_i,max)/2, (79)ki,mean=12⋅Cm.

This is the same value of k_i obtained from computing k_i when U₁ = U₂. This simple expression is a useful relationship for tuning the integrator network. One may select C_m to obtain the intended mean integration rate, and then minimize the variation of the integration rate by minimizing g_s, as long as equation (57) is satisfied.

Figure 5D graphically demonstrates how k_i,mean and k_i,range determine C_m and g_s. Just as in equation (79), k_i,mean is a function only of C_m. Therefore, the contour only shows vertical lines. The value of k_i,range is minimized by decreasing either g_s or Cm−1 (i.e., increasing C_m). Figure 6 shows simulation data of the integrator’s response to a step input with eight different parameter value combinations. In every case, the change in U₁ is bounded by the values of k_i,mean and k_i,range. As shown in Figure 5D, increasing C_m decreases the integration rate, and increasing g_s increases the variation in the integration rate.

Table 2 summarizes the design approach for this integrator network. The mean and range of the integration rate are free parameters that are determined by the intended network performance. Using these values and the constraint in equation (57), the neurons’ C_m value and the synapses’ g_s and ΔE_s values can be fully specified.

Some of the calculations in this paper are based on approximations, which lead to inaccuracies in the calculations of the subnetworks. One example is that the subtraction network does not produce linear output. This non-linearity occurs because the reversal potentials of synapses are rarely much lower than the resting potentials of neurons, requiring large values of g_s,2 to build a subtractor where k_syn = 1. A large g_s,2 value increases U_pre,2’s effect on the denominator of Upost∗’s response, causing the synaptic input to reduce U_post’s sensitivity to inputs. This is particularly noticeable in the differentiator’s response (Figure 4), especially as k_syn increases.

Another example of a simplification we made is that our calculation of k_i only used the particular solution of the system. This means that a transient response also exists, which we did not compute. In addition, k_i is a function of U₁ and U₂. This means that k_i is not a constant for this network. However, the impact of U₁ and U₂ on k_i can be minimized by minimizing g_s and maximizing R, as we showed in Sec. 4.2.

However, the developed networks are not intended to act as perfect analogs to their mathematical counterparts. These networks are intended to act as representations of real neural circuits, which likely do not act as perfect adders, multipliers, differentiators, etc. Dynamic and transient effects are a real part of biological control systems, and effective neural controllers have developed around these idiosyncrasies and have likely evolved to even exploit many of these aspects. In spite of these issues, the methods in this paper are valuable. Our recent robotics work (Szczecinski and Quinn, 2017; Szczecinski et al., 2017a), as well as related work in progress, is proof of the effectiveness of this approach.

The methods in this paper enable the direct construction of networks that perform arithmetic and dynamic calculations. Why bother building neural networks just to recreate mathematical operators? We believe there are several reasons to take this approach. From a neurobiology perspective, the constraints that we have identified may help explain why certain structures are common in the nervous system (David Friel, personal correspondence). For instance, mutually inhibitory parallel pathways are common in the thoracic control of insect locomotion (Büschges and Wolf, 1995), which may function as subtraction networks in negative feedback loops. As another example, networks in the retina of the rabbit are selectively sensitive to motion in one direction or the other (Barlow and Levick, 1965). Such a network could be constructed by using adjacent cells in the retina as inputs to differentiator networks. This would be consistent with both the function of direction-sensitivity, as well as the laterally inhibitive structure. Even though such consistency does not guarantee that the animal’s nervous system functions precisely this way, the design methods in this paper may aid in understanding the function of neural networks found in animals.

Additionally, the constraints that we identified may be used to constrain parameter values in large brain models. Rather than using global search techniques to understand the dynamics of a large pool of neurons, we believe it may be faster to begin with a number of functional subnetworks, and then use local search techniques to tune the connections between them. In this way, the designer is certain that parts of the network perform specific, useful computations, rather than naively optimizing a large network (Haferlach et al., 2007; Agmon and Beer, 2013; Izquierdo and Beer, 2013). The end result is something like a genetic program, but in a neuroscience context.