Neuronal networks with NMDARs and lateral inhibition implement winner-takes-all

A neural circuit that relies on the electrical properties of NMDA synaptic receptors is shown by numerical and theoretical analysis to be capable of realizing the winner-takes-all function, a powerful computational primitive that is often attributed to biological nervous systems. This biophysically-plausible model employs global lateral inhibition in a simple feedback arrangement. As its inputs increase, high-gain and then bi- or multi-stable equilibrium states may be assumed in which there is significant depolarization of a single neuron and hyperpolarization or very weak depolarization of other neurons in the network. The state of the winning neuron conveys analog information about its input. The winner-takes-all characteristic depends on the nonmonotonic current-voltage relation of NMDA receptor ion channels, as well as neural thresholding, and the gain and nature of the inhibitory feedback. Dynamical regimes vary with input strength. Fixed points may become unstable as the network enters a winner-takes-all regime, which can lead to entrained oscillations. Under some conditions, oscillatory behavior can be interpreted as winner-takes-all in nature. Stable winner-takes-all behavior is typically recovered as inputs increase further, but with still larger inputs, the winner-takes-all characteristic is ultimately lost. Network stability may be enhanced by biologically plausible mechanisms.


INTRODUCTION
The winner-takes-all (WTA) function is an operation that is often assumed to take place in biological nervous systems, and it has been demonstrated to be a powerful computational primitive (Maass, 2000). Over the years, winner-takes-all networks have been modeled widely in the fields of computational brain science and artificial neural networks (Amari, 1972;Grossberg, 1973;Koch and Ullman, 1985;Rumelhart and Zipser, 1987;Yuille and Grzywacz, 1989;Coultrip et al., 1992;Ermentrout, 1992;Winder, 1999;Yuille and Geiger, 2003;Mao and Massaquoi, 2007;Handrich et al., 2009;Chen et al., 2013), and have also been implemented in various analog electronic circuits (Lazzaro et al., 1989;Andreou et al., 1991;Deweerth and Morris, 1995;Lau and Lee, 1998;Fish and Yadid-Pecht, 2001;Indiveri, 2001;Baishnab et al., 2010). These mathematical and silicon models for the most part achieve the WTA characteristic using some sort of common inhibitory feedback in conjunction with high gain and a strong nonlinearity, often in the form of a neural thresholding operation. In its strictest forms, the winner-takes-all amounts simply to selection of the largest from among a set of input signals; in "soft" implementations, it involves strong amplification of the difference between the largest and other such signals, with significant suppression of weaker signals within the set. Depending on the model, it may retain some analog information about the winning signal.
In this paper, I demonstrate a mechanism by which a neuronal network with global lateral inhibition, arranged in perhaps the simplest possible topology, is capable of realizing winner-takes-all. Although nonlinearities such as neural thresholding contribute to this capacity, it depends most critically on the unique electrical properties of NMDA receptor (NMDAR) ion channels, which are assumed to mediate excitatory input to the network. The NMDAR is a class of glutamatergic receptor for which N-methyl-D-aspartate (NMDA) is an agonist, and is found in many phyla and frequently associated with synapses. The current-voltage relationship of NMDAR ion channels is nonmonotonic under physiological conditions (Nowak et al., 1984;Jahr and Stevens, 1990), with a negative slope conductance regime due to (kinetically fast) magnesium blockade. This characteristic renders the NMDAR capable of supporting neural amplification (Shoemaker, 2011) and bistability (Lazarewicz et al., 2006;Shoemaker, 2011;Sanders et al., 2013) in conjunction with other membrane conductances. The primary finding here is that a WTA characteristic can be induced by high-gain regimes that result from such interactions, rather than requiring a high intrinsic or parametric gain in the feedback loops. In this respect, the model contrasts with other biologically-inspired WTA network models that in some way incorporate NMDARs (Winder, 1999;Handrich et al., 2009;Chen et al., 2013). It is significant because it represents a mechanism for WTA that is both simple and at the same time entirely plausible biophysically, relying on known characteristics of ubiquitous classes of synaptic receptors. It is also of interest due to the widespread distribution of neurons with glutamatergic synapses and lateral inhibition in many areas of the brain. Although not depicted, the outputs of neurons N 1 , N 2 , N 3 , . . . are assumed to also project to other parts of the nervous system.

THE WTA NETWORK MODEL
The WTA microcircuit described herein assumes a classical lateral inhibitory topology, with a set of competitive neurons that receive excitatory inputs via NMDA synapses, and global feedback inhibition via a common interneuron, as illustrated schematically in Figure 1.
The WTA characteristic in this model is conceived as resulting from strong amplification and nonlinear effects that can emerge from interactions between NMDARs and the other ion channels in the membrane (in this case, the parallel resting and inhibitory synaptic conductances). These can result in an electrically bistable membrane, in which the dependence of the transmembrane current on the membrane potential shows the classical N-shape that is associated with destabilizing first order and stabilizing higher odd-order terms in a canonical bistable system (illustrated in Figure 8 in the Methods Section). Such behavior is enhanced when the total non-NMDAR conductance has a reversal potential below typical neural resting potentials, and also when its current-voltage dependence is sub-linear or inward-rectifying (Shoemaker, 2011;Sanders et al., 2013). With this in mind, consideration is given to inhibitory feedback with several different characteristics: either mildly hyperpolarizing (e.g., as might be the case if chloride channels were associated with the inhibitory synapses) or more strongly hyperpolarizing (which might be expected with potassium channels), and having either ohmic or inward-rectifying channels. Gamma-aminobutyric acid (GABA)mediated synapses, which are prime candidates to implement the inhibition, occur in two major classes that reflect different combinations of these characteristics: the GABA B class involves Kir channels, which conduct potassium and are inward-rectifying (Mott and Lewis, 1994;Sodickson and Bean, 1996;Kaupmann et al., 1998;Fowler et al., 2007) (and thus might be expected to promote high-gain behavior), whereas GABA A receptor channels conduct chloride and are non-rectifying (Johnston, 1996;Olsen and DeLorey, 1999).
Development of a mathematical model (described in the Methods section) and numerical analysis are useful tools for analyzing the range of behaviors than can be expected from such a network. The model used herein is based on simple isopotential neurons. Neural signals are represented as continuous in time, with any spiking behavior assumed representable by mean spike rates that are in turn related to neural depolarization. Several different inhibitory characteristics are considered, including ohmic and inward rectifying channels with either weakly or strongly hyperpolarizing reversal potentials, in order to illustrate a range of behaviors associated with plausible inhibitory mechanisms. Magnitudes of synaptic conductances are specified in the text relative to resting membrane conductance. In addition to the inhibitory reversal potential, free parameters in the model include the magnitudes of the inputs to the network, and a loop gain parameter (defined in Methods) that applies to the inhibitory feedback. The instantaneous loop gain (i.e., the incremental gain around a feedback loop) of any particular circuit is state-dependent and can greatly exceed this loop gain constant in magnitude, but as a parameter it accounts for the strengths of the input and output synapses of the inhibitory interneuron and thus is useful to quantify the effectiveness of the feedback.

ANALYSIS OF STATIONARY EQUATIONS
Although in a biological system such a network would be expected to operate under dynamic conditions, a stationary analysisi.e., determination of the fixed points of the model's governing equations-can give significant insight into its functional characteristics. I undertake such an analysis in this section to characterize the range of stationary behaviors that can be expected, and in particular to determine the existence of multiple fixed point solutions for given levels of synaptic input. In the absence of time-dependence in the feedback loops, such solutions are indicative of bi-or multi-stable regimes, and I use these terms to refer to them throughout this section. However, whether such fixed points are in fact stable depends on the dynamical characteristics of the feedback pathway, which will be considered in the following section. Figure 2 depicts some of these dc characteristics, and illustrates the emergence of WTA behavior. For reference purposes, Figure 2A shows the dc input (expressed as relative NMDAR conductance Γ 1 ) vs. output (membrane potential) relationship for a single competitive neuron N 1 when that neuron is the only one in the network receiving excitatory input. The loop gain constant is the parameter. Membrane potential varies smoothly with input activation, and as might be expected, the slope of this relationship decreases as the strength of inhibitory feedback is increased. Figures 2B-F show the behavior of different configurations of the model when two of the competitive neurons receive excitatory input, and the network transitions from one "winner" to another as the first input exceeds the other. (It should be noted that this transitional behavior applies not just to the two-input case, but to any situation in which the two neurons with the largest inputs are the only ones that are able to reach a FIGURE 2 | (A) A dc input-output characteristic for the network when only a single competitive neuron N 1 is activated. The inhibitory feedback synapses are ohmic with reversal potential −90 mV. The abscissa is the normalized conductance Γ 1 of the input (NMDA) synapse onto the neuron, the ordinate is its membrane potential, and the parameter is the loop gain constant A L , which is varied from 0 to -5 in decrements of −1. (B-F) dc transition characteristics for several configurations of the network when two neurons receive excitatory input. The loop gain constant A L = −4 in each case. Input Γ 1 is swept in individual trials, and the parameter is the input Γ 2 to the second neuron N 2 , which is fixed during each sweep of Γ 1 but varied from trial to trial. Membrane potential V m 1 of N 1 is shown in black and V m 2 of N 2 in gray. In (B-D), the inhibitory channels are ohmic and Γ 2 is stepped from 0 to 35 in increments of 5. In (B) the inhibitory reversal potential is −70 mV; in (C,D) it is −90 mV. (C) shows results for upward sweeps of Γ 1 while (D) shows V m 1 for sweeps in both directions, which trigger transitions at juxtaposed limit points and illustrate the presence of hysteresis. Directions of transitions are indicated by thin arrows. The heavy arrow at lower right points to a small secondary hysteresis loop (not fully visible at this resolution). In (E,F) the inhibitory channels are inward-rectifying, and Γ 1 is swept upward in both cases. In (E) the reversal potential is −70 mV and Γ 2 is stepped from 0 to 35 in increments of 5; in (F) the reversal potential is −90 mV and Γ 2 is stepped from 0 to 17.5 in increments of 2.5. Resting potential in all cases is −60 mV. state of depolarization). When the inhibition is ohmic and mildly hyperpolarizing (Figure 2B), the transitions are gentle and the network cannot reasonably be characterized as winner-takes-all. In the remaining cases, however, there are input ranges for which the transitions are not only sharp but discontinuous. These discontinuities and the hysteretic effects that accompany them are associated with bistable regimes. Such a regime prevails over the range of Γ 1 values bounded by each hysteresis loop in Figure 2D. In these regions, fixed-point solutions exist in which either neuron may be significantly depolarized (i.e., in a "high" state), while the other is either hyperpolarized or depolarized to a lesser extent (in a "low" state). Which of the two solutions might be assumed under quasistatic conditions depends on the history of excitation of the system. The discontinuous jumps depicted in Figures 2C-F represent limit point transitions, which occur at fold bifurcations that correspond to the boundaries of bistable regimes. Outside of transitions between winners, the state of the winning neuron depends on its input level, and thus carries analog information about that input.
The fixed-point solutions in the model as formulated are shaped by a rich array of nonlinearities, including the NMDAR current-voltage relationship and the neural thresholding function, but also the rectifying characteristics of Kir channels (if present) and saturation of the membrane potential in a "winning" neuron. Numerical analysis shows their structure can be quite complex. Bistability is not observed at all when the inhibition is ohmic and mildly hyperpolarizing, as in Figure 2B. When bistability is possible, it requires a minimum value of the loop gain constant, and when the constant exceeds this minimum, it also requires some minimum level of input activation. This is because bistability is not supported by NMDAR conductance acting in parallel with the resting membrane conductance alone, and there must be sufficient activation of inhibition with characteristics that do allow bistability in order to Frontiers in Computational Neuroscience www.frontiersin.org February 2015 | Volume 9 | Article 12 | 3 overcome the effects of the resting membrane in a competitive neuron. Bistable behavior is not observed when only a single neuron received excitatory input. This is due to the dependence of the inhibitory current in such a neuron on its own activation, which is approximately second-order due to feedback. Such a supra-linear characteristic suppresses bistable behavior. However, in the other neurons in the network, this dependence is more nearly linear, because they are inhibited below their activation thresholds and do not participate in the feedback. For this reason, when the excitatory input to one of these "losing" neurons is sufficient to bring it to the point of depolarization, bistability may be induced and a limit point crossed, as seen in Figures 2C-F.
Another feature of significance seen in Figure 2 is the fact that the bistable regime with the WTA characteristic-i.e., in which the neurons are in opposing high and low states-ultimately collapses and vanishes when the inputs become sufficiently large. This is associated with the fact that at high input levels, saturation of the membrane potential in the winning neuron begins to limit the recruitment of inhibition, ultimately allowing the losing neuron to depolarize and start participating in the feedback process even when its input is no larger than that of the winner. The result is seen in Figures 2C,D, wherein the curves for the largest value of Γ 1 reflect monostability, and in addition are so gentle in transition that the winner-takes-all characteristic must be regarded as having vanished there as well.
The monostable and WTA bistable regimes described above do not exhaust the dynamical regimes possible in this model. Detailed examination shows the presence in the parameter space of secondary bistable and multistable regimes, which tend to appear as the loop gain and input amplitude are increased. By way of example, Figure 3 shows a portrait of the dynamical structure of the network with instantaneous feedback, for the case of ohmic inhibitory channels and inhibitory reversal potential −90 mV, in the parameter space spanned by the mean or common-mode NMDAR activation Γ C = (Γ 1 + Γ 2 )/2, the differential activation Γ D = (Γ 1 − Γ 2 )/2, and A L .
The lenticular regions in Figure 3 show that increasing the loop gain promotes the range of inputs for which bistability is present, as well as the magnitude of hysteresis. The presence of regimes in which two distinct stable "low" states are present (the falciform regions in Figure 3) is related to the neural thresholding function; it may be of academic interest but is of little relevance for the WTA function. (Although not visible at the resolution of the graph, a small hysteresis loop associated with such a regime is present in one of the traces in Figure 2D, and is indicated by a heavy arrow).
When inhibitory conductance is modeled as inward-rectifying, the dc behavior of the two-neuron network is qualitatively similar to the ohmic case in that there is a critical value of the loop gain parameter necessary to support bi-and multi-stable behavior, and above that value a bistable WTA regime appears and then disappears with increasing common-mode NMDAR activation. This bistable behavior is more pronounced: it is supported at higher inhibitory reversal potentials (c.f. Figures 2B,F) and at lower values of the loop gain parameter than when the inhibition is ohmic, and hysteresis is more extensive as well. The co-dimension three cusp bifurcation that marks its beginning occurs at [A L = 1.09, Γ C = 3.81, Γ D = 0], which may be compared to the values cited in Figure 3. In addition, the secondary bistable regimes described for the ohmic case-the green falciform regions in Figure 3-are also present. FIGURE 3 | Characterization of the stationary dynamical regimes of the two-neuron system, for ohmic inhibition with reversal potential −90 mV, in the space spanned by A L , Γ C = (Γ 1 + Γ 2 )/2, and Γ D = (Γ 1 − Γ 2 )/2. These regimes are symmetric about Γ D = 0 due to the physical symmetry of the network. Cross-sections of multistable regimes in Γ C − Γ D planes are depicted, from an oblique view at left, and in a projection onto the Γ C − Γ D plane at right. A cusp bifurcation of co-dimension three at [A L = 1.43, Γ C = 9.55, Γ D = 0] marks the beginning of multistability and is indicated by a star. Associated with the region that has blue lenticular cross-sections are two sets of stable equilibria, in which one neuron resides in a "low" state (either hyperpolarized or slightly depolarized) and the second in a "high" state (significantly depolarized), or vice-versa. This defines a WTA bistable regime. As Γ D increases, this region begins and then terminates at co-dimension two cusp bifurcations. Associated with the narrow regions that have green falciform cross-sections are two sets of stable equilibria, in which one neuron resides in one of two low states (one slightly hyperpolarized and the other slightly depolarized) and the second in one of two respective high states. Where a lenticular region intersects a single falciform region, the network is tristable with one neuron residing in one of two low states or a single high state, while the other resides respectively in one of two high states or a single low state. Where all three intersect, the network has four stable states, two high and two low for each neuron. Outside the depicted regions, the network is monostable.

Frontiers in Computational Neuroscience
www.frontiersin.org February 2015 | Volume 9 | Article 12 | 4 However, further complexity is introduced by the saturation of inhibitory synaptic current at depolarized membrane potentials. Significantly, as common-mode NMDAR activation increases, mono-and multi-stable regimes emerge in which both neurons can simultaneously assume significantly depolarized equilibrium states (a characteristic that precludes a single "winner"). Figure 4 illustrates the various possible dynamical regimes, for one particular choice of the loop gain parameter.
Characterization of the dc behavior of this network has been given for the case in which only a pair of competitive neurons receive non-zero input signals. Simulations show that as the number of active inputs increases, dynamical regimes of increasing complexity become possible. In particular, with inward-rectifying inhibition, multistable states are possible in which three or more competitive neurons are depolarized if their inputs are large enough and similar in magnitude. Nevertheless, with the appropriate loop gain, there are input ranges for which WTA behavior-significant depolarization of a single neuron with complete or substantial suppression of activation in all other neuronsprevails in networks with more than two inputs, and these ranges subsume multistable regimes in which any one neuron among a set with similar input levels may be a winner. As in the twoneuron case, hysteretic behavior is associated with these regimes and changes in winning states occur at limit points corresponding to regime boundaries.
Finally, it is worth mentioning that evidence such as observations of biphasic inhibitory post-synaptic potentials in FIGURE 4 | Characterization of the stationary dynamical regimes of the two-neuron system, for inward-rectifying inhibition with reversal potential −90 mV and A L = −4, in the space spanned by Γ C = (Γ 1 + Γ 2 )/2 and Γ D = (Γ 1 − Γ 2 )/2. The characteristics of the equilibria associated with the blue lenticular region and the green falciform regions are similar to those in the regions with the same color codes in Figure 3. Associated with the chevron-shaped red region are two stable equilibria in which one neuron is in one of two possible high states, and the second respectively in a high or low state. In the yellow region the system is monostable with both neurons in high states (in Figure 2F, the right-hand portion of the traces for largest Γ 2 are associated with this regime). Where the chevron and the lenticular region alone intersect, the network is tristable with one neuron residing in a high or one of two low states and the other neuron respectively in a low or one of two high states. Where the red, lenticular, and one falciform region intersect, the network has four stable states, and where they all intersect, six stable states. various CNS neurons (e.g., Alger and Nicoll, 1979;Davies et al., 1990) suggests that GABA-ergic pathways can activate both GABA A and GABA B receptors in target neurons. If this were to occur in the context of feedback inhibition in this model, it would lead to a mixture of inwardrectifying and non-rectifying inhibitory currents in the competitive neurons. Simulations show that such mixed inhibition results in dc behavior intermediate between the purely rectifying and non-rectifying cases, as might be expected. Interestingly, with a mixture of weakly-hyperpolarizing/ohmic and stronglyhyperpolarizing/inward-rectifying inhibition, a proportion of ohmic conductance as small as 25% of the total can eliminate the multistable solutions with multiple "high" states that are seen with inward-rectifying inhibition alone (e.g., the red regions in Figure 4).

DYNAMIC BEHAVIOR
Under non-stationary operating conditions, time-varying synaptic input signals play the role of driving or forcing functions in the dynamics of the network. Analysis of time-dependent behavior is complicated by the fact that not just neural state, but dynamical regimes themselves depend on synaptic inputs, and these regimes can change on a time scale comparable to changes in membrane potential. This can lead to complex time-domain behavior. Nevertheless, a number of conclusions can be reached by time-domain and ac simulations and analysis.
The kinetics of NMDARs associated with the input synapses, which do not play a direct role in the feedback loops, do affect network dynamics in the sense that they have a lowpass filtering effect on the input signals. Conventionally, NMDAR kinetics are regarded as exceptionally slow for an ionotropic receptor (Destexhe et al., 1994;Jahr, 1994;McBain and Mayer, 1994), although they may vary with subtype and with neuromodulatory state (with receptors incorporating the NR2A subunit, for example, having appreciably faster response times (Monyer et al., 1994;Flint et al., 1997;Cull-Candy et al., 2001) than other oligomeric combinations). The response of prototypical slow NMDARs to impulsive inputs (e.g., arrival of a presynaptic action potential or short burst thereof) has a time-to-peak of tens of milliseconds, and decay times on the order of 100 ms or more (Destexhe et al., 1994). Simulation of the ac response of a linearized NMDAR kinetic model (Destexhe et al., 1995) yields a lowpass −3 dB corner frequency at around 1 Hz. One result is that if the presynaptic signaling is spiking in nature, then relative variation in the state of a population of synaptic receptors due to each individual spike is insignificant when the firing frequency is more than a few tens of Herz. The implications of these slow kinetics for network function are that changes in network state may significantly lag the presynaptic input signals that drive them, and also that high-frequency changes in winning states would tend to be suppressed when the network is in a winner-takes-all regime.
An important question with respect to network dynamics is whether the fixed points identified by the stationary analysis are stable or unstable in nature. The answer is determined by the dynamical characteristics of the inhibitory feedback in conjunction with input state and other network parameters. The main focus of this study is on stable WTA behavior and the conditions that support or promote stable regimes in the network, and so Frontiers in Computational Neuroscience www.frontiersin.org February 2015 | Volume 9 | Article 12 | 5 loss of stability is considered mainly in this context. However, brief attention is given to operating characteristics in unstable regimes, particularly as pertains to WTA function, although a detailed treatment of unstable dynamics is beyond the scope of the paper.
To examine the issue of fixed-point stability of the feedback circuits, I include dynamic elements in the model network and apply a standard technique of linearization of the governing equations about the fixed points, as detailed in Methods. These elements include membrane capacitance in the competitive neuron, and selectively in the feedback path, a pole to model the effects of membrane conductance and capacitance in the inhibitory interneuron, and cascaded synaptic receptor models.
Consider first the simplest input scenario in which a single competitive neuron is activated. When membrane capacitance is added to the model neuron but the feedback remains instantaneous, a single pole residing in the left half of the complex plane is introduced in the dynamics and the fixed points of the circuit are stable. However, when independent rate equations are also present for one or more states in the feedback path, it is possible for the single-neuron system to lose stability on some subset of its fixed points. With this loss of stability, large-scale limit cycles-i.e., oscillations-occur. Unstable regimes coincide with small, or even negative, incremental membrane conductance in the competitive neuron-a condition that also prevails when a network with multiple active inputs enters a WTA mode of operation. (When the conductance is negative at a fixed point, instantaneous feedback stabilizes the network, but this effect is compromised with the introduction of sufficient feedback lag.) Conversely, when the system is not in a state that would support WTA behavior, it tends to be robustly stable. (This is the case for weakly hyperpolarizing ohmic inhibition, for which consistent fixed-point stability is seen for all feedback dynamics and input conditions considered in the study).
Destabilization of the single-neuron model system can be induced by the addition of even a single feedback pole, if the time constant and the loop gain are large enough. Such instability is enhanced by the introduction of multi-state dynamics to represent synaptic delays. As might be expected, when very slow (metabotropic) GABA B receptors are present, unstable behavior is more pronounced than with the much faster (ionotropic) GABA A receptors, in the sense that it occurs for smaller loop gains and a wider range of inputs. Oscillations in unstable regimes tend to be quite slow for GABA B inhibition (around 4 Hz for the models/parameters used herein), but considerably faster (20-30 Hz) for GABA A .
Analysis shows that when the single-neuron network enters an unstable regime, it does so via a Hopf bifurcation (Hopf, 1942;Andronov et al., 1966;Marsden and McCracken, 1976), i.e., the transition of a complex pole-pair between left-and right-halfplanes with variation of a parameter-in this case, the input strength. Furthermore, stability is typically lost and then regained via such bifurcations as the input increases in magnitude. This characteristic behavior is illustrated in Figure 5, in which stable and unstable regimes are illustrated for a particular network configuration (weakly hyperpolarizing, inward-rectifying inhibition) and choice of feedback dynamics. Qualitatively similar results are obtained for other network configurations that support WTA behavior. Figure 6A shows results of a time-domain simulation in which stability is lost and regained.
Analysis also shows that when stability is lost, it occurs via a supercritical Hopf bifurcation (in which a small limit cycle appears around the fixed point and increases in amplitude with increasing input). Stability is also regained via a supercritical Hopf bifurcation in the case of ohmic inhibition, but when the current is inward rectifying it occurs via a subcritical Hopf bifurcation (in which the fixed point becomes stable but a large-signal limit cycle remains outside a finite-sized basin of attraction about the fixed point). Thus, sustained oscillatory behavior can persist even after a fixed point has regained stability. The basin of attraction of the fixed point, however, grows rapidly with increasing input levels and eventually subsumes the entire state space.
What are the implications of this analysis for network operation when other inputs are activated? When an unstable regime is present in the single-neuron system, it defines a minimal range of input strengths over which the network will be unstable: even if the outputs of other competitive neurons were completely suppressed, sustained oscillations would take place if the largest input were in this range. However, simulations demonstrate that with multiple inputs, entrained, phase-locked oscillations can occur when input values are beyond this range, but within a regime that supports WTA behavior (e.g., within portions of the green region in Figure 5 that are not overlapped by the red). What would be damped oscillatory behavior in a single neuron can become unstable when additional inputs are activated, even though the fixed point solutions in this regime are themselves stable and WTA in character. With sufficient loop gain, this behavior can be seen in any of the model configurations supporting a WTA characteristic. It is associated with the summation of competitive outputs onto the inhibitory interneuron: the effect of phase-locked oscillations from the standpoint of an individual neuron is the same as an increase in its inhibitory loop gain, which is naturally destabilizing. This combination of large-signal limit cycles coexisting with fixed point stability is again characteristic of a system that undergoes a subcritical Hopf bifurcation.
The regime in which the network resides depends on the history of activation of the inputs: for example, entrained oscillations may grow when the onset of inputs to two or more competitive neurons is simultaneous, but if the same inputs were staggered in time, the initially-active neuron may settle to a degree that the entire network remains stable after all inputs have become active. Examples of such behavior are shown in the time domain in Figures 6B,C.
What are the implications of oscillatory behavior for the WTA function?-in particular, do peak membrane potentials in unstable regimes show any sort of WTA characteristic, when the corresponding fixed points for those inputs are WTA in nature? This question is germane due to the consistency of the present model with so-called PING models for gamma oscillations in vertebrate cortex (Whittington et al., 2000), along with the finding that such networks can implement an "E%-max winners take all" function (De Almeida et al., 2009)-which has led to the hypothesis that a multiple winners-take-all computation may in fact be a function of gamma oscillations (De Almeida et al., 2009). In the present Frontiers in Computational Neuroscience www.frontiersin.org February 2015 | Volume 9 | Article 12 | 6

FIGURE 5 | Graphs illustrating the fixed-point stability characteristics of one configuration of the WTA network when a single input is active.
Inhibitory synaptic conductance is inward-rectifying with reversal potential −70 mV. The resting membrane time constant τ R in the competitive neuron is set to 20 ms; membrane charging in the inhibitory interneuron is modeled with a pole with time constant 5 ms, and is cascaded with a four-state model for GABA B receptor dynamics (Destexhe et al., 1995) that introduces three additional poles in the feedback path. At left are depicted the real parts of the poles of the closed-loop system at its fixed points as functions of the input strength Γ , in this case for A L = −4. Each is coded by a different color. Where two curves overlie, the corresponding poles are complex conjugates. The arrows indicate the points at which a complex pole-pair transitions between left-and right-half-planes and bracket an input range over which the fixed points are unstable. In the middle part of this region, the poles are real and distinct, implying exponential divergence of trajectories away from the fixed point-but large-signal limit cycles (i.e., oscillations) are ultimately assumed by such trajectories. At right, the region over which the network is unstable in the space spanned by Γ and A L is shown in red (for purposes of comparison, the pink curve shows the boundary of this region when GABA A kinetics are substituted for GABA B ).
In the blue area the membrane conductance in the competitive neuron is negative for the input and inhibitory states at the fixed point, and in the green, WTA bistability would be possible if the input to a second neuron were active (and feedback instantaneous).
case, the answer is yes: simulations show that WTA behavior can be imputed to oscillatory regimes-but not invariably or strictly. At left in Figure 7 is shown an example of steady-state oscillations in a network of a particular configuration, in which the outputs of the neuron with the largest input exceeds the resting potential throughout its cycle, while the outputs of four other neurons with smaller inputs are uniformly suppressed below the resting potential. Conversely, Figure 7 also shows at right a result for another configuration in which the WTA characteristic is weakened by instability, in the sense that membrane potentials in non-winning neurons are not as strongly suppressed as they would be at the corresponding fixed points. However, although all of the neurons in this second case exceed the resting potential during their cycles, the dynamic gain during oscillation is still relatively high, since the smallest of the five inputs is only 20% smaller than the largest. Qualitatively, the behavior in these sorts of oscillatory regimes is seen to depend largely on the difference between the largest and second-largest inputs to the network. More generally, simulations show that unstable behavior can be relatively complex, for example showing evidence of long time-scale modes that mediate a switch between strict WTA behavior and lower-gain excursions as the system evolves in time for particular choices of configuration and input levels. Examination of such complexities are beyond the scope of this study. Interestingly, unstable behavior in this model may be suppressed by several mechanisms that might plausibly be found in biological neural circuits. One is saturation in the inhibitory interneuron, e.g., the limiting of its membrane potential by the finite reversal potential associated with its input synapses-a feature that is perhaps more realistic than the linear model used in the prior analyses. This sort of saturation was investigated in a limited number of simulations, and an example in which it suppresses entrained oscillations is depicted in Figure 6D. A second is a mixture of GABA A and GABA B -type inhibition in the competitive neurons, with the faster and more nearly ohmic GABA A channels exerting a stabilizing influence when inwardrectifying channels alone might result in unstable equilibria and limit cycle behavior. In qualitative terms, the fast activation of GABA A channels can allow the network dynamics to remain damped while slower changes in the states of the GABA B channels are ramping up or down-and even though the GABA A channels might not by themselves permit strong WTA behavior, their presence can allow the network to settle into stable WTA states. An example of such behavior is shown in Figure 6E.
Finally, simulations with more general time-varying inputs demonstrate that when the network is in a stable WTA mode (i.e., in the basin of attraction of a stable fixed point), it will subsequently remain stable for any limit-point transitions between winning neurons so long as the largest input remains outside the unconditionally unstable range. Figure 6 also serves to emphasize graphically the temporal sluggishness of the network response when both slow NMDAR kinetics and GABA B receptor dynamics are present in the model. In order for a network with these dynamics to act as a WTA in a biological context, the time scales on which the input signals operate (i.e., their persistence and the times between changes in their rank ordering) must be on the order of at least several 100 ms. This restriction would seem to preclude the network from functioning as a sequential WTA if its inputs were modulated at frequencies comparable, say, to that of the theta rhythm in mammals. Faster resolution of stable winning and losing states would naturally be the result, however, if faster NMDAR kinetics, and/or mixed GABA A and GABA B inhibition, were present. Membrane potential of the neuron in response to input steps with three different magnitudes of NMDAR conductance Γ is shown. Each step is applied at time t = 100 ms. For Γ = 4 the network is in a relatively low-gain stable state; for Γ = 8 it is in an unconditionally unstable state; and for Γ = 12 it is in a stable regime that would support WTA behavior.

Frontiers in Computational
(B,C) Responses when a first neuron N 1 receives a step input with Γ 1 = 12 at t = 100 ms, and a second neuron N 2 a step input whose amplitude is 90% that of the first (Γ 2 = 10.8). Membrane potential of N 1 is shown in black and that of N 2 in gray. In (B) the two steps are simultaneous, and the network breaks into entrained oscillations, even though it would be stable without the second input. In (C) the input to N 2 is delayed an additional 100 ms, and the network settles to a stable fixed point with a WTA characteristic. In (D) the network has the same configuration as in (A-C), except a saturating input-output characteristic has been attributed to the inhibitory interneuron N FB . The input steps are synchronous; but although there is a large transient depolarization in N 2 , the entrained oscillations that are present in (B) are squelched here. In (E) the inhibitory feedback is split between ohmic/weakly-hyperpolarizing and inward-rectifying/strongly-hyperpolarizing, with the former accounting for 25% of the total conductance. The inputs are stepped simultaneously to Γ 1 = 15, Γ 2 = 13.5 at t = 100 ms.

DISCUSSION
The results presented here demonstrate how a neuronal network with NMDA input synapses and lateral inhibition can act as a winner-takes-all circuit, and illuminate some of its significant operating characteristics. The possibility of WTA behavior is conditioned on the characteristics of the inhibitory feedback in the network-in particular, on the amount of loop gain and the electrical characteristics of synaptic ion channels associated with the inhibitory feedback. When WTA behavior is supported, it emerges as the magnitude of inputs to the network increase, with transitions between "winners" taking place with increasing gain, and then becoming abrupt as they are associated with limit points of multistable regimes (stable, that is, if feedback is instantaneous). The output of the winner-quantified herein as its level of membrane depolarization-carries analog information about the magnitude of its input when the network is in a WTA regime. As input strengths increase further, however, such regimes ultimately collapse, giving way again to finite-gain behavior, or to multistable regimes that do not have a WTA characteristic.
Thus, if the neuronal network is to operate as a WTA, the maximum NMDAR conductance must be bounded from above by some value that is determined by the strength and the characteristics of the inhibitory feedback. Biophysically, such a bound is a natural constraint that could be imposed by finite numbers of post-synaptic NMDA receptors available, and/or limits on presynaptic neurotransmitter release, at the inputs to the competitive neurons.
In a multistable regime in which transitions between winners are abrupt, if one neuron is established as a winner then the input to another neuron must exceed that of the first by some finite amount in order for such a transition to occur. This means that at any instant in time the "winning" neuron may not actually be the one with the largest input. Such hysteretic behavior, however, can serve the purpose of preventing transitions due to noise or to relatively small differences between the input levels of two contending neurons-in other words, it can prevent "jitter" between winners and could be interpreted as requiring some level of confidence in the result before a transition can occur. Analysis of time-domain behavior of the network with plausible dynamics in the inhibitory feedback loops show that it is possible for it to destabilize. With sufficient loop gain, as one or more inputs increase in magnitude, the network can enter a regime in which its fixed-point states (including some that would be multistable with instantaneous feedback) are unstable. When it receives sustained inputs in this regime, the network evolves into large-signal oscillations (entrained oscillations when there is more than one input). The peak membrane potentials during such oscillations may in some circumstances be regarded as WTA in character, and generally they reflect a high dynamic gain although the WTA function may be compromised relative to the corresponding fixed point solutions. As inputs increase further, the system exits this regime and enters a second in which there are stable fixed points with WTA characteristics-but these may coexist with large-signal limit cycles.
Destabilization of the network is most pronounced (i.e., occurs at lower loop gains and covers a larger range of input strengths) when the slow dynamics of (second-messenger) GABA B receptors are associated with the inhibitory synapses in the feedback loops. Very small or negative membrane conductance values in the competitive neurons, which are more prevalent when the inhibitory current is modeled as inward-rectifying, also promote instability. Thus, it might be expected that lateral inhibition implemented with GABA B synapses would compromise the ability of such a network to implement a stable WTA function. However, other biologically-plausible characteristics can mitigate against unstable behavior. When the simplifying assumption of linearity in the inhibitory interneuron is relaxed and saturating effects introduced, it can squelch large-signal oscillations and establish a stable WTA regime (although the network may pass through an unconditionally unstable regime at lower input values to reach it). In addition, the presence of GABA A as well as GABA B receptors in the inhibitory feedback pathway can have a significant stabilizing effect.
Although tangential to the main subject of the paper, it is of interest to note that when this model is in unstable regimes, it produces spontaneous oscillations whose frequencies are similar to the theta rhythm when GABA B kinetics are present in the feedback, and when GABA A kinetics are present, to the gamma rhythm (although somewhat slower for the particular parameter values used). It is well-known that the gamma as well as beta1 and beta2 bands are dependent on the time constants of GABA A receptor kinetics (Whittington et al., 2000). The results herein confirm that gamma-like oscillations can arise in a PING architecture with fixed or slowly-varying inputs, and in addition that behavior in these regimes is generally consistent with the hypothesis that a (single or multiple) winners-take-all function might be an inherent feature of gamma oscillations (De Almeida et al., 2009).
At this point, a word on the magnitudes of NMDAR conductance and the inhibitory loop gain constant in the model is in order. As defined in Methods, the loop gain parameter is simply the gain around any one of the feedback loops in the network when it is at rest; conceptually, this means that if a feedback loop in the quiescent network were broken at the level of a competitive neuron, then an increment of presynaptic depolarization at its output synapse would come back around the loop to effect a hyperpolarization |A L | times as large in the cell itself. The synaptic weightings implied by the single-digit values of A L considered in the study would therefore seem to be quite plausible. Analysis also suggests that the densities of active NMDARs required to achieve the reported network behaviors are biologically plausible. The numerical values assumed by the NMDAR conductance parameter may be misleading in this regard, because by convention it refers to the slope conductance at the channel reversal potential, where magnesium blockade is minimal. Where the blockade is more pronounced-as over most of the operating range of the network-the slope conductance of the NMDAR channels (relative to resting conductance) is typically much smaller than the Γ parameter itself.

Frontiers in Computational Neuroscience
www.frontiersin.org February 2015 | Volume 9 | Article 12 | 9 FIGURE 8 | Functions f N (solid black) and f I (solid gray) that respectively characterize the voltage-dependence of NMDAR ion channels and Kir channels (as defined by the second line in 5). The Kir channel reversal potential is set to −90 mV in this example. The dashed curve depicts the current-voltage relation of a membrane (with current normalized by a unit conductance to give a voltage) in which channels of both types are active, with the NMDAR conductance twice than of the Kir conductance, leading to electrical bistability. The middle zero of this function represents an unstable equilibrium and the other two, stable equilibria.
In summary, if such WTA microcircuits are in fact implemented in biological neuronal networks, they would endow nervous system with a building block of appreciable computation power (Maass, 2000), and could support functions ranging from the abstract (e.g., function approximation) to the psychophysicallyor ethologically-grounded (e.g., selective attention). But although the functional implications at the level of the primitive circuit are clear, any possible role in broader contexts (e.g., of the cortical circuits documented in the extensive literature on recurrent/reverberating networks) requires further consideration. (For example, observations of an "inverted U-shaped" dependence of persistent activity on external drive in assemblies of cortical cells (Williams and Goldman-Rakic, 1995;Brunel and Wang, 2001) could be related to the behavior of this circuit as it is driven up to, through, and then past a WTA regime by increasing input strength.) Finally, although the subject of learning (e.g., by longterm potentiation) is beyond the scope of this paper, it is clear that the model represents a biophysically-plausible component of a substrate for competitive learning. It would provide for Hebbian modification exclusive to a single winning neuron during any cycle of activation, and strengthening of excitatory synapses within the network would predispose neurons so affected to future wins.

MATHEMATICAL MODEL
In order to establish the fundamental behavior of the network, I consider a model which is intended to be as simple as practical while maintaining enough biological fidelity to characterize the possible range of behaviors of the neuronal networks it represents. The following specifications describe the framework for this model: • All neurons in the network shown in Figure 1 are electrically compact and representable as single electrical compartments; • Each neuron has a passive membrane conductance that can be approximated as ohmic, with an associated resting membrane potential (set to −60 mV in simulations), and these characteristics remain fixed in time; • The competitive neurons are biophysically homogeneous, i.e., they have the same resting membrane properties, and the synapses of each class are weighted identically from neuron to neuron; • Synapses onto the competitive neurons are represented as variable membrane conductances in the postsynaptic cell with values determined by presynaptic membrane potential. The nonlinear voltage dependence of NMDAR channel current is based on the model and parameters of Jahr and Stevens (1990). The reversal potential for the inhibitory synapses is assigned a value of either −70 mV or −90 mV in simulations, and inwardrectifying channels are modeled with the ad hoc current-voltage relationship used by Shoemaker (2011); • Neuronal activation involves a thresholding operation corresponding to half-wave rectification: output synapses are activated in proportion to depolarization but are inactive during hyperpolarization; • For the sake of simplicity, the inhibitory interneuron is approximated with a linear model unless otherwise indicated. It performs a simple summation of its inputs, and its behavior under non-stationary conditions is represented with first-order linear dynamics; • Postsynaptic conductance under non-stationary conditions is determined by convolution of the presynaptic activation with an impulse response function representing a linearized model of receptor kinetics (as given in Destexhe et al., 1995 for NMDARs, GABA B , and GABA A receptors).
This simple model neglects the generation of action potentials, and is thus applicable to cases in which the neurons either operate with graded potentials, or for spiking cells, in which outputs can be characterized by a spike rate that may in turn be related to mean membrane depolarization. In addition, although spiking input pathways could be represented with impulsive inputs in the time domain, continuous-time inputs are used throughout. This may be justified by the lowpass-filtering effects of slow NMDAR kinetics at the input synapses, as discussed in Section Dynamic Behavior. The model is not intended to represent networks in which spike timing is in some way critical, or in which individual spikes are otherwise significant computationally. In any of the competitive neurons N 1 , N 2 , N 3 , . . . , the membrane current balance equation may be written in the form where V m is membrane potential, and G N , G I , and G R are conductance parameters and the functions f N (V), f I (V; V r ), and f R (V; V r ) characterize the voltage dependence of channel currents for NMDARs, inhibitory receptors, and resting membrane conductance, respectively. The parameter V r in this notation indicates reversal potential; V rI and V rR are reversal potentials for inhibitory ion channels and resting conductance, respectively. The Frontiers in Computational Neuroscience www.frontiersin.org February 2015 | Volume 9 | Article 12 | 10 reversal potential for NMDAR channels is assumed to be 0V and is not included as a parameter. After the convention in (Shoemaker, 2011), the functions f N , f I , and f R are defined with dimensions of volts and are normalized to unity slope at their reversal potentials. G N is assumed to be governed by the activation of the NMDA synapse associated with the input to the neuron, and G I by the activation of the inhibitory feedback synapse. C m is membrane capacitance. Assuming that all competitive units in the network (i.e., neurons and synapses of each class that contact them) are biophysically identical, I append a numerical index to quantities that may vary between them, i.e., V m i and G Ni , i = 1... n, to indicate that they apply to the ith neuron N i among a set of n. In addition, it is convenient to eliminate a free parameter by multiplying (1) by R R = 1/G R , yielding the set of equations where Γ i = R R G n i and Γ I = R R G I are dimensionless synaptic conductances expressed relative to the resting conductance, and τ R = R R C m is the resting membrane time constant. The normalized dc transmembrane current I mi for the ith neuron is and collectively, the zeros of these currents define the fixed points of the system. The voltage dependence for NMDAR channel current is based on the model and parameters of Jahr and Stevens (1990), and the corresponding function f N may be written where b = 0.28 mM −1 · [Mg 2+ ] = 0.336 under an assumed extracellular magnesium concentration [Mg 2+ ] = 1.2 mM, and where k = 62V −1 . The voltage dependence of inhibitory synaptic ion channels takes one of the two forms where the form for Kir channels is the ad hoc function used in Shoemaker (2011), with d = 25 mV, e = 0.5, and c = −13.73 mV (which places the zero-crossing of the function at V rI ), and where tanh denotes hyperbolic tangent. The voltage-dependent functions for NDMAR and Kir channels are depicted in Figure 8. Under the assumption that the resting conductance is ohmic, The output of each competitive neuron is computed by applying a thresholding function h(V m ) ∼ = max[(V m − V rR ), 0] to the membrane potential to obtain the depolarization of the neuron above the resting potential. In practice, the function on the righthand side of this expression is splined with a quadratic function over the range V rR ± 1 mV to maintain continuity of both the thresholding function and its derivative. The output is identically zero for V m < V rR − 1 mV. I assume that the relationship between the NMDAR conductance Γ i and the input IN i to neuron N i is linear, and thus without loss of generality equate the two for purposes of stationary analysis. For dynamics, I assume that Γ i = h N * IN i , where * indicates temporal convolution, h N is a dimensionless unit impulse function based on the linearized 4-state NMDAR kinetic model of Destexhe et al. (1995), and IN i assumes the dimensions and scaling of the NMDAR conductance. The inhibitory conductance Γ I that is driven by the interneuron N FB takes on a common value for all competitive neurons, but in general is variable and dependent on the states of those units via a feedback equation. With the assumptions of electrical compactness and approximate linearity in N FB , under stationary conditions it may be written in the form where K is a constant of proportionality that includes the effects of the strength or weight of the input and output synapses of N FB . Now consider a single-neuron system with ohmic inhibitory conductance and with no input activation. Conceptually, break the feedback loop at the presynaptic side of the input to N FB . For an ideal thresholding function (i.e., ignoring the spline used in the numerical simulations), the right-handed derivative of Γ I with respect to the presynaptic membrane potential V m Pre is simply dΓ + I /dV m Pre = K for any V m Pre ≥ V rR . The derivative of the membrane potential of the competitive neuron with respect to Γ I is found to be dV m /dΓ I = (V rI − V m )/(Γ I + 1) by implicit differentiation, and evaluating this in the state that the closed-loop system assumes when it is not excited, i.e., V m = V rR and Γ I = 0, gives dV m /dΓ I = (V rI − V rR ). Thus, the loop gain under these conditions is I take the rightmost expression in (8) as the loop gain parameter: When the inhibitory channels are inward-rectifying, the quiescent loop gain is not precisely equal to this expression due to the nonlinearity of their channel current-voltage relation, but I maintain the same definition of A L so that it bears the same relation to the synaptic weights in the feedback paths. This allows straightforward comparison of the behavior of the two different types of inhibitory channels considered. With this definition, (7) becomes Frontiers in Computational Neuroscience www.frontiersin.org February 2015 | Volume 9 | Article 12 | 11 When considering the dynamic behavior introduced into the feedback loop by the kinetics of the inhibitory channels, the expression in (10) is replaced by where h I is a dimensionless unit impulse function based on one of the linearized four-state GABA kinetic models of (Destexhe et al., 1995). Finally, saturation of membrane potential in the inhibitory interneuron, when present, is modeled via the trans-synaptic response at its inputs. Identify (without loss of generality, since any scaling may be absorbed in the loop gain constant) the net input to N FB , n j=1 h(V m j − V rR ), with a current that is expressed across a unit membrane resistance to determine the depolarization ΔV mFB . When saturation is present, ΔV mFB is determined implicitly by the equation ΔV mFB = n j=1 h(V m j − V rR ) · (ΔV mFB + k sat · V rR )/(k sat · V rR ), which gives a hyperbolic input/output dependence as would be imposed by a limited reversal potential associated with the input synapses to N FB . The value k sat = 2 was used in simulations; with this value, full depolarization of a single competitive neuron (i.e., to a membrane potential of zero) drives the inhibitory interneuron halfway to saturation.
Stability analysis is carried out via linearized first-order state equations, written in terms of deviations of the states (indicated below by the prefix Δ) from their fixed-point values. The membrane potential V m of the single active competitive neuron is one state. The kinetics of the input synapses to the interneuron N FB are neglected, and the first state V FB in the feedback loop is assumed proportional to the membrane depolarization of N FB . (In practice, the scaling of the entire feedback loop is absorbed into this state.) The dynamics due to membrane resistance and capacitance contribute a single pole with time constant τ FB . When included, the dynamics of the output synapses, as represented by one or the other of the GABA receptor models of (Destexhe et al., 1995), are cascaded with this pole. These models include two "closed" states, C 1 and C 2 , where C 1 combines with a ligand L and may then transition toC 2 , from which an inactive state D and an open (conductive channel) state O are accessible. I associate C 1 + L with V FB , leaving three independent states that generate three additional poles. The state O is associated with the inhibitory feedback conductance Γ I in the target competitive neuron (the weight of the inhibitory synapse being, without loss of generality, absorbed into the loop gain parameter).
With these conventions, the linearized dynamic equations for the closed-loop system with feedback synapse dynamics can be written in vector-matrix form: where I m is the normalized dc transmembrane current in the competitive neuron as given in (3), k FB is a scaling constant yielding a unit step response from the receptor dynamics, and the partial derivatives are evaluated at the fixed point under consideration. Rate constants r 2 through r 6 are as defined by Destexhe et al. (1995) for the receptor type under consideration, with r Σ ≡ r 2 + r 3 + r 6 . The eigenvalues of the matrix on the right-hand side of (11) correspond to the poles of the system.

NUMERICAL ANALYSIS
Simulations were performed primarily with SPICE (Simulation Program with Integrated Circuit Emphasis), a tool that is intended for electronic circuit simulations and is optimized for solution of stiff nonlinear equations. Discrete elements are described by electrical constitutive relations involving their interconnection points or nodes, and governing equations for interconnected circuits are automatically generated by application of Kirchhoff 's current law (conservation of charge). Built-in linear elements or user-definable nonlinear voltage-controlled current sources were used for the implementation of nonlinear membrane conductances. SPICE can solve static equations to obtain dc/equilibrium states, and dc sweeps of inputs were performed to determine input/output relations (with sweeps in both directions used to identify limit point transitions). Membrane current-voltage relations could be tested by sweeping a virtual "voltage clamp" with the input and inhibitory conductances held at fixed-point values. SPICE also allows time-domain or transient simulations of dynamical systems, and was used to compute outputs in response to various input scenarios (such as the step functions in Section Dynamic Behavior). Analogical circuits were constructed to implement receptor dynamic models. The product T-SPICE (Tanner Research, Monrovia, CA) was used for the study. In addition, MATLAB (MathWorks, Natick, MA) was used to compute the poles of the system in the stability analysis. MATLAB and Excel (Microsoft, Redmond, WA) were used for supporting analysis and plot generation.