A dynamic neural field model of continuous input integration

Biological Cybernetics (2021) 115:451–471 https://doi.org/10.1007/s00422-021-00893-7 ORIGINAL ARTICLE A dynamic neural field model of continuous input integration Weronika Wojtak1,2 ·Stephen Coombes3·Daniele Avitabile4,5 ·Estela Bicho2·Wolfram Erlhagen1 Received: 1 April 2021 / Accepted: 1 August 2021 / Published online: 21 August 2021 © The Author(s), under exclusive licence to Springer-Verlag GmbH Germany, part of Springer Nature 2021 Abstract The ability of neural systems to turn transient inputs into persistent changes in activity is thought to be a fundamental requirement for higher cognitive functions. In continuous attractor networks frequently used to model working memory or decision making tasks, the persistent activity settles to a stable pattern with the stereotyped shape of a “bump” independent of integration time or input strength. Here, we investigate a new bump attractor model in which the bump width and amplitude not only reflect qualitative and quantitative characteristics of a preceding input but also the continuous integration of evidence over longer timescales. The model is formalized by two coupled dynamic field equations of Amari-type which combine recurrent interactions mediated by a Mexican-hat connectivity with local feedback mechanisms that balance excitation and inhibition. We analyze the existence, stability and bifurcation structure of single and multi-bump solutions and discuss the relevance of their input dependence to modeling cognitive functions. We then systematically compare the pattern formation process of the two-field model with the classical Amari model. The results reveal that the balanced local feedback mechanisms facilitate the encoding and maintenance of multi-item memories. The existence of stable subthreshold bumps suggests that different to the Amari model, the suppression effect of neighboring bumps in the range of lateral competition may not lead to a complete loss of information. Moreover, bumps with larger amplitude are less vulnerable to noise-induced drifts and distance-dependent interaction effects resulting in more faithful memory representations over time. Keywords dynamic neural field ·conservation law ·localized states ·stability ·input integration Communicated by Benjamin Lindner. BWeronika Wojtak w[email protected] Stephen Coombes [email protected] Daniele Avitabile d.a[email protected] Estela Bicho [email protected] Wolfram Erlhagen [email protected] 1Research Centre of Mathematics, University of Minho, Guimarães, Portugal 2Research Centre Algoritmi, University of Minho, Guimarães, Portugal 3Centre for Mathematical Medicine and Biology, School of Mathematical Sciences, University of Nottingham, Nottingham, UK 4Department of Mathematics, Vrije Universiteit, Amsterdam, The Netherlands 1 Introduction The ability of the brain to accumulate, combine and maintain information over behaviorally relevant time scales is fundamental for a wide variety of cognitive functions that bridge perception and action (Miller and Cohen 2001). This has been associated with persistent neural activity observed in many brain regions in working memory, decision making or planning tasks (Curtis and Lee 2010). The longstanding hypothesis that reverberations in recurrent neural networks can give rise to sustained activity throughout periods without external stimulation (Hebb 1949) has motivated an impressive amount of experimental and computational modeling work (for reviews see Wang 2003; Zylberberg and Strowbridge 2017). A particularly insightful theoretical framework has been the continuous bump attractor network (Amari 1977; Brody et al. 2003b; Compte et al. 2000) which has been widely used in the past to explain data from 5MathNeuro Team, Inria Sophia Antipolis Méditerranée Research Centre, Sophia Antipolis, France 123 452 Biological Cybernetics (2021) 115:451–471 behavioral experiments and to synthesize cognitive behavior in robotics systems (Erlhagen and Bicho 2006; Schöner 2016). It explains the emergence of stable representations of continuous-valued information, such as for example the direction of heading during navigation or the position of an object in space, by assuming a distance-dependent neuronal connectivity pattern in feature space. Typically, neurons tuned to similar values of a continuous variable excite each other, and those tuned to dissimilar values inhibit each other. When the neurons are ordered along a line by their selectivity, the self-stabilized activity pattern, which the intrinsic network dynamics develops in response to a transient input, is visualized as a spatially localized activity bump. Due to the assumed translation-invariant connection structure, the network can hold a continuous family of bumps, with each of the attractor states representing the memory of a specific input value. Continuous attractor networks provide some degree of robustness to perturbations away from the attractor manifold since perturbed network activity quickly decays back to the bump state. However, since bumps are neutrally stable, their position can be easily shifted along the manifold by weak external inputs (Amari 1977). This property can be exploited for instance to explain the capacity of a neural system to track the position of a moving object in real time (Wu et al. 2008), and for a robotics application (see also Bicho et al. 2000). Direct experimental evidence for a continuous bump attractor hypothesis comes from visuo-spatial working memory (WM) tasks in which a subject has to maintain the information about the spatial position of one or more objects over a short period of time. After the cessation of the tuned external input, the bump position is vulnerable to random activity fluctuations which may cause a drift of the bump over time (Camperi and Wang 1998; Kilpatrick and Ermentrout 2013). This model-derived predictive relationship between persistent neural population activity and the variability of spatial memory over time has been reported recently in a combined electrophysiological and behavioral task (Wimmer et al. 2014). Consistent with the neural dynamics in continuous attractor models of multi-item memory (Amari 1980; Ferreira et al. 2016; Krishnan et al. 2018), several behavioral experiments report distance-dependent interaction effects between memory traces which may be described as an attractive or repulsive bias of memorized nearby items (Almeida et al. 2015; Johnson et al. 2009, for a discussion of the neural underpinning of attraction versus repulsion see also Erlhagen et al. 1999). Despite their success in many applications, the explanatory power of classical continuous attractor models remains still limited since the memory bump has a stereotyped shape exclusively determined by the recurrent interactions within the network. A more sophisticated circuit-based model should regulate also the shape and rate of localized persistent activity depending on the qualitative and quantitative characteristics of the preceding inputs. It has been shown for instance that the level of persistent activity during the delay period of a spatial WM task correlates with stimulus contrast (Constantinidis et al. 2001) and can be modulated by additional spatially informative cues (Kuo et al. 2012; Wildegger et al. 2016). Importantly, the discharge rate of neural populations is predictive of the psychophysical task performance. The notion of a static bump as a neural correlate of WM is also inconsistent with recent evidence showing that persistent input-selective activity often changes systematically during the course of the delay period (e.g., ramping-type activity (Brody et al. 2003b)). One hypothesis is that the variability in activity is due to the continuous integration of other task relevant information such as for instance the time elapsed since stimulus presentation (Brody et al. 2003a). The capacity of neural circuits to rapidly change firing rates in response to transient inputs while maintaining firing at a constant level for extended periods of timeis fundamental for decision making processes (Curtis and Lee 2010; Sakai et al. 2006). This requires the accumulation and storage of evidence for competing options across time. Neural recordings from animals suggest that the level of activity of neural populations representing a certain option correlates with choice certainty (Kiani and Shadlen 2009). The observation that persistent activity related to the outcome of previous choices often continues during inter-trial intervals shows that the integration time can be quite long (Histed et al. 2009). Here, we analyze a new bump attractor model (Wojtak et al. 2016) based on the theoretical framework of dynamic neural fields that is able to hold in addition to a continuum of bump positions also a continuum of bump amplitudes. The construction of recurrent networks that can transform brief inputs into sustained firing rates at different levels has been proven to be a nontrivial problem (Koulakov et al. 2002; Sakai et al. 2006). Previous work was based on the assumption of networks with linear input-output functions or, when more realistic nonlinear neurons were used, heavily relied on numerical fine tuning of network parameters (Brody et al. 2003b; Seung 1998). A balanced relation of excitation and inhibition in the network must be kept so that the sustained activity at any time represents the time integral of past inputs. A recent dynamic field model addresses the existence of bumps with varying amplitudes in a lateral inhibition network with separate excitatory and inhibitory populations and a piecewise linear firing rate function (Carroll et al. 2014). To support a range of possible bump amplitudes, the recurrent excitation must be precisely chosen to be inversely proportional to the gain of the transfer function. To mitigate this fine tuning problem, additional mechanisms such as the use of bistable neurons with different activation thresholds have been proposed (Koulakov et al. 2002), for review of different computational approaches see (Sakai et al. 2006). Due to the bistability, the integrator network exhibits multi123 Biological Cybernetics (2021) 115:451–471 453 ple stable states even with imprecisely tuned feedback. The inherent trade-off is a loss of sensitivity to weaker input since a minimum input strength is necessary to transition between adjacent states. Our proposed model does not rely on a domain of linear input-output transfer and supports in principle a continuous integration of inputs of any strength and duration. It consists of two coupled dynamic fields of Amari type (Amari 1977). The field receiving external input has a classical Mexican-hat coupling with excitation dominating at short distances and surround inhibition. The second field integrates the activity from the first field with an inverted Mexican-hat connectivity and feeds its activity back locally. We first analyze the existence, stability and bifurcation structure of bumps with input-dependent shapes and discuss their relevance to modeling cognitive functions such as working memory and decision making. We then systematically compare the pattern formation process in the Amari and the two-field model in tasks with multiple localized inputs. Specifically, we are interested in evaluating how the local feedback mechanisms of the two-field model affect the interaction of nearby bumps and the robustness of memory maintenance in the presence of random activity fluctuations. 2 Definition of the two-field model We study a system of two coupled neural field equations defined on a one-dimensional finite domain Ω:[−L,L]: ∂u(x,t) ∂t=−u(x,t)+v(x,t)+I(x,t) +Ω w(|x−y|)f(u(y,t)−θ)dy, (1a) τv ∂v(x,t) ∂t=−v(x,t)+u(x,t) −Ω w(|x−y|)f(u(y,t)−θ)dy. (1b) Here the variables u(x,t)and v(x,t)represent the activity (e.g., membrane potential) of a neuron at the spatial location xat time tin two coupled populations, labeled uand v, respectively. The firing rate function f(x)is chosen as a Heaviside function with threshold θsuch that f(x)=1 for x≥0 and f(x)=0 otherwise. The weight function w(|x−y|)representing the synaptic connectivity between neurons at locations xand yhas a Mexican-hat shape given by the difference of two Gaussian functions: w(x)=Aexe−x2/2σ2 ex−Aine−x2/2σ2 in−winh,(2) with Aex >Ain >0 and σin >σ ex >0. The parameter winh >0 guarantees that spatially uniform inhibition exists at larger distances, |x|>xc. The time-dependent external input I(x,t)to the u-population is modeled as one or more Gaussians centered at positions xcj: I(x,t)=(Ht0(t)−Hte(t)) n  j=1 AIje−(x−xcj)2/2σ2 Ij,(3) where Ht0(t)represents the Heaviside step function with threshold t0≥0 controlling the start and the end of the input at times t0and te, respectively. In the following we use for the input duration the notation dI=te−t0. Input strength is controlled by AIj>0. The u-equation of system (1) without the term v(x,t)represents the field model originally introduced and analyzed by Amari (1977). He considers two separate inhibitory and excitatory subpopulations and derives the single population model with Mexican-hat connectivity by assuming that the wider ranging inhibitory feedback in the network is infinitely fast. This ensures that the spread of self-sustained excitation remains localized in a region activated by external input. Findings in neurophysiological and modeling studies of the cortical microcircuit of spatial WM stress the important role of local inhibitory feedback for generating persistent activity (Constantinidis and Wang 2004). A distinct subclass of interneurons has been described which integrates inputs from nearby excitatory cells and exhibits sustained responses with similar tuning properties during delay periods. As part of a feedback loop, these neurons might implement a tight balance of excitatory and inhibitory inputs to a target cell. Biologically plausible circuit models of WM have also addressed how cortical disinhibition of interneurons (e.g., by interneurons targeting interneurons) with different spatial tuning affect WM maintenance (Murray et al. 2014). The mechanistic two-field model implements a balanced local inhibitory feedback and a disinhibition of lateral inhibitory feedback by assuming that the v-population integrates the activity of the u-population with an inverted Mexican-hat connectivity. Note that a spatial integration with a Mexican-hat profile could be applied as well when −v(x,t) instead of +v(x,t)is used in the u-equation. To understand the role of the v-population in the pattern formation process, it is instructive to follow Amari’s strategy and also consider different time scales for the two populations. In the limiting case τv→0, the v-population is at quasi-equilibrium with v(x,t)=u(x,t)−Ω w(|x−y|)f(u(y,t)−θ)dy.(4) 123 454 Biological Cybernetics (2021) 115:451–471 Substituting this expression into the u-equation of (1) yields ∂u(x,t) ∂t=I(x,t). (5) The u-population thus implements a robust temporal integration of external inputs if fast inhibitory feedback loops counterbalance any input-induced changes in the recurrent interactions (Denève and Machens 2016; Lim and Goldman 2013). For the integration to be perfect, we have to assume a linear response function of the v-population, f(v) =v. The presence of nonlinearities in the neural dynamics such as saturation will make the integrator circuit insensitive to additional input when the saturation limit of the response function f(v) is reached (Lim and Goldman 2013). Notably, the integrator capacity does not depend on the precise choice of the connectivity function w(x)or the nonlinearity f(u). In the rest of the paper, we analyze and apply model (1) with identical timescales for the u-population and the vpopulation (τv=1). This allows us to directly compare the behavior of the two-field and the Amari models. 2.1 Existence and stability of one-bump solutions 2.1.1 Linear stability analysis à la Amari In the following we study the existence and stability of single bump steady state solutions of the model (1) in the case of a Heaviside firing rate function and initial condition u(x,0)+ v(x,0)=K, where K>0 is a constant. Let R[u(x,t)]be the region over which the field is excited R[u(x,t)>θ]=(x1(t), x2(t)) (6) and let the gradients of u(x,t)at the boundaries x1and x2of the excited region be c1=∂u(x1,t) ∂x,−c2=∂u(x2,t) ∂x.(7) After a short period of time dtthe excited region changes to R[u(x,t+dt)]=(x1(t+dt), x2(t+dt)). (8) The boundaries of the excited region satisfy at time t:u(xi,t)=θ, (9a) at time t+dt:u(xi+dxi,t+dt)=θ, (9b) where xi(t+dt)=xi+dxi,i=1,2. A Taylor expansion of the latter equation yields ∂u(xi,t) ∂xdxi+∂u(xi,t) ∂tdt=0,i=1,2,(10) Fig. 1 aThe plot of F(Δ) givenby(16)withK=1. For θ=0.2, there is one (stable) solution with Δ=4.14. For θ=0.9, there are two solutions: unstable solution with Δ=0.67 and stable solution with Δ=2.65. bThe plot of E(Δ) givenby(20). For θ=0.2, there is a local minimum with Δ=4.14. For θ=0.9, there is a local maximum with Δ=0.67 and a local minimum with Δ=2.65. The kernel w is given by (2)withAex =3, σex =1.4, Ainh =1.5, σinh =3and winh =0.2 for infinitesimally small dxiand dt. We can then rewrite (1)atx=xi(t)as ∂u(x,t) ∂t=−u(x,t)+v(x,t)+x2(t) x1(t) w(|x−y|)dy, (11a) ∂v(x,t) ∂t=−v(x,t)+u(x,t)−x2(t) x1(t) w(|x−y|)dy. (11b) Let W(x)=x 0 w(y)dy.(12) Since u(xi,t)=θ, we have from (11a) ∂u(xi,t) ∂t=−2θ+K+W(x2−x1). (13) 123 Biological Cybernetics (2021) 115:451–471 455 We have dx1 dt=−∂u/∂t ∂u/∂xx=x1 =−1 c1 (−2θ+K+W(x2−x1)), (14a) dx2 dt=∂u/∂t ∂u/∂xx=x2 =1 c2 (−2θ+K+W(x2−x1)). (14b) We can then describe the change of length of the excited region Δ(t)=x2(t)−x1(t) dΔ dt=1 c1 +1 c2(−2θ+K+W(Δ)).(15) The existence of solution of width Δ=x2−x1is determined by the roots of F(Δ) =−2θ+K+W(Δ) =0.(16) The stability condition is that a steady state of width Δis stable if dF(Δ) dΔ<0,(17) and unstable otherwise. We plot examples of function F(Δ) for different values of threshold θin Fig. 1a. 2.1.2 Lyapunov method The stability of bump solutions can be also analyzed by searching for local minima and local maxima of a Lyapunov functional associated with system (1) (French 2004; Kubota and Aihara 2005; Owen et al. 2007). For the two-field model with Heaviside firing rate function, the functional is given by E[u]=−1 2ΩΩ w(|x−y|)f(u(x,t)−θ) ×f(u(y,t)−θ)dxdy +(2θ−K)Ω f(u(x,t)−θ)dx. (18) For a bump of width Δ=(x2−x1)we have E(Δ) =−1 2x2 x1x2 x1 w(|x−y|)dxdy +(2θ−K)(x2−x1). (19) Using (12) we obtain E(Δ) =−Δ 0 W(x)dx+(2θ−K)Δ. (20) Figure 1b shows plots of E(Δ) for the example used in Fig. 1a. As indicated by the vertical lines, the local minima correspond to stable and the local maxima to unstable bumps. 3 Steady states of the model in the absence of input 3.1 Bumps with initial condition u(x,0)+v(x,0)=K We start the numerical investigation of one-bump solutions of the two-field model using the symmetric initial condition u(x,0)=Aue−x2/2σ2 u,v(x,0)=K−u(x,0), (21) where Au>0, σu>0 and K≥0. Stable bumps are found by numerically integrating (1) to a steady state. Following the numerical analysis of the Amari model (Amari 1980), we expect convergence if the width of the initial profile above threshold exceeds the width Δ1of the unstable pattern in Fig. 1. To find the unstable bump associated with a narrower initial profile, we apply numerical bifurcation techniques described in (Rankin et al. 2014). These techniques also allow us to track steady states of (1) when model parameters are varied. The bifurcation curve in Fig. 2a with the threshold θas continuation parameter shows that for K<2θbranches of stable (solid line) and unstable (dashed line) bumps coexist. For K≥2θ, the unique steady state is a stable bump. The panels on the right show particular solutions as indicated by the points P1,P2and P3on the bifurcation curve. Here and in all following figures, solid and dashed black lines represent u(x)and v(x), respectively. Like for the Amari model, the shape of the bump attractor does not depend on the shape of the initial profile (compare panels (d) and (e)). Note that for a subthreshold initial profile, Au<θ, the activity converges to the homogeneous solution (u,v)=(K/2,K/2). 3.2 Bumps with initial condition u(x,0)+v(x,0)=K(x) To work as a neural integrator, the bump shape of the twofield model (1) should depend on the initial conditions for the u-population and the v-population since they may be set by external input. Figure 3shows examples of bump solutions for the initial conditions u(x,0)=K(x), v(x,0)=0,K(x)=AKe−x2/2σ2 K, (22) which represent a homogeneous initial state for the vpopulation and a spatially structured state for the u-population. 123 456 Biological Cybernetics (2021) 115:451–471 Fig. 2 aBifurcation curve showing one-bump solutions of (1) with a Mexican-hat kernel (2) as the threshold parameter θis varied. Solid (dashed) lines represent stable (unstable) solutions. umax is the maximum of uover all x.(b-e) Examples of solutions at the points P1 (unstable) and P2(stable) for θ=0.8 (dotted line) and P3(stable) for θ=0.4 (dotted line) are shown. Solid and dashed black lines represent u(x)and v(x), respectively. For the solution at the point P3we show examples of different initial conditions, with solid (dashed) grey lines representing initial profiles u(x,0)(v(x,0)). They are given by (21) with Au=0.75, σu=1, K=1(d)andAu=1.5, σu=4, K=1 (e). For the remaining solutions, the initial condition is given by (21) with Au=1andσu=1, K=1. The kernel wis given by (2) with Aex =2, Ain =1, σex =1.25, σin =2.5andwinh =0.1 Fig. 3 One-bump solution of the two-field model (1) (black lines) for two different initial profiles (22) (grey lines) with AK=1, σK=1(a)and AK=2, σK=2(b). Parameters of the kernel as in Fig. 2. Threshold θ=0.5 (dotted line) As can be clearly seen when comparing the activation patterns in the two panels, the bump shape correlates perfectly with the shape of the initial activation profile of the upopulation. Figure 4shows results of numerical continuation of bump solutions of (1) with initial condition (22) and threshold θas a continuation parameter. The bifurcation curves for the smaller (a) and the larger (c) initial profiles are different but maintain the same qualitative features. Again we can see that in a certain range of threshold values stable and unstable bumps coexist. However, different to the case of the symmetric initial condition (21), there is an additional branch of stable subthreshold bumps given by u(x,t)=v(x,t)=K(x)/2. Panels (b) and (d) show examples of the three types of bump solutions at the points P1,P2and P3of the bifurcation curves. It can be seen that a larger initial profile leads to larger bumps and shifts their coexistence to higher threshold values. In the example, the width and the amplitude of the initial profile were manipulated affecting both the width and the amplitude of the evolving pattern. It is important to notice that stable bumps with a constant width and varying amplitude can be created for instance with a truncated Gaussian profile which limits the spatial range of the initial activation. We come back to this point when discussing the continuous integration of localized inputs of varying duration and strength in Sect. 4. 3.3 Multi-bump solutions 3.3.1 Numerical continuation of two-bump solutions Different from the Amari model (Laing and Troy 2003), the two-field model supports the existence of stable twobump solutions for which the inter-peak distance is within the competitive range of the Mexican-hat connectivity. We use numerical continuation to track two-bump solutions of (1)as the threshold parameter θvaries. Figure 5shows the resulting bifurcation curve with stable and unstable suprathreshold branches and a subthreshold branch. Different to the situation for single bumps, the subthreshold solution (P1) does not coexist with the suprathreshold two-bumps (P3 and P4). To evaluate the dependence of the inter-peak distance on model 123 Biological Cybernetics (2021) 115:451–471 457 Fig. 4 aand cBifurcation curves showing one-bump solutions of (1) with the initial condition (22) with AK=1, σK=1(a)andAK=2, σK=2(c) as the parameter θis varied. Solid (dashed) lines represent stable (unstable) solutions. Examples of solutions at the points P1,P3 (stable) and P2(unstable) for a narrower and a wider profile of K(x) are shown in panels (b)and(d), respectively. Threshold θ=0.7(b) and θ=1.25 (d) (dotted line). Parameters of the kernel as in Fig. 2 parameters, Fig. 6depicts the bifurcation curve when the lateral inhibition parameter σinh is varied. As illustrated by the two pairs of stable and unstable two-bumps, (P2, P3) and (P1, P4), respectively, the distance decreases with increasing σinh. The Amari model would predict for the same distances a pronounced repulsive interaction of the two localized activation patterns. 3.3.2 Existence and stability of N-bump solutions It is well known that classical continuous attractor networks have difficulties simultaneously stabilizing multiple bumps due to the constraints imposed by lateral inhibition and the interference between nearby patterns (Amari 1977; Wei et al. 2012). This limitation has been addressed in computational studies by introducing additional processing mechanism such as short-term synaptic facilitation (Rolls et al. 2013)orby using interaction kernels that are not of lateral-inhibition type (Ferreira et al. 2016; Laing and Troy 2003). However, already Amari (1980) pointed out that for the connectivity function (2) multiple stable bumps may exist at distances, |x|>xc, where the lateral inhibition is spatially uniform (for multiitem WM applications see e.g. Johnson et al. 2009; Wei et al. 2012). It is straightforward to derive a condition for the existence of an N-bump which is constrained by the summed inhibition that each bump receives from the N−1 neighboring patterns. The bump width depends on Nbut should be equal for each bump. For the numerical realization in a finite integration domain, one has to apply circular boundary conditions or open boundary conditions with a sufficiently wide lateral inhibition kernel. Without loss of generality, we assume that the bump is located in the interval [−Δ /2,Δ /2]and that the remaining N−1 pulses are located in the intervals [x1−Δ /2,x1+Δ /2],[x2−Δ /2,x2+Δ /2], …, xN−1−Δ /2,xN−1+Δ /2, with xi+1−xi>Δ+xs, i=1,2,...,N−1. That is, u(x)>θfor x∈[−Δ /2,Δ /2]∪ [x1−Δ /2,x1+Δ /2]∪...∪xN−1−Δ /2,xN−1+Δ /2, and u(x)≤θelsewhere. Equilibrium solutions of the two-population model are given by U(x)=V(x)+Δ /2 −Δ /2 w(y)dy+x1+Δ /2 x1−Δ /2 w(y)dy +...+xN−1+Δ /2 xN−1−Δ /2 w(y)dy, (23a) 123 458 Biological Cybernetics (2021) 115:451–471 Fig. 5 Bifurcation curve showing two-bump solutions of (1) with the initial condition (22)withAK=1.5, σK=1 as the parameter θis varied. Solid (dashed) lines represent stable (unstable) solutions. Examples of two-bump solutions at the points P1,P3,P4(stable) and P2(unstable) are shown. Threshold θ(dotted line): θ=1.9(P1), θ=0.82 (P2 and P3), θ=0.4(P4). Parameters of the kernel as in Fig. 2 Fig. 6 Bifurcation curve showing two-bump solutions of (1) with the initial condition (22)withAK=1.5, σK=1 as the parameter σinh is varied. Solid (dashed) lines represent stable (unstable) solutions. Examples of two-bump solutions at the points P1,P2(stable) and P3,P4 (unstable) are shown. σinh =3(P2and P3), σinh =4(P1and P4). Remaining parameters as in Fig. 2 V(x)=U(x)−Δ /2 −Δ /2 w(y)dy−x1+Δ /2 x1−Δ /2 w(y)dy −...−xN−1+Δ /2 xN−1−Δ /2 w(y)dy. (23b) Since at the boundaries of a bump U(−Δ /2)=U(Δ /2)=θ, U(x)+V(x)=Kand xi+Δ /2 xi−Δ /2 w(y)dy=−Δwinh,(24) we get the existence condition for N-bump solutions F(Δ) =−2θ+K+W(Δ) −(N−1)Δ winh =0.(25) The width of an individual bump of an N-bump solution can be thus found graphically by looking on the intersection of the plot of W(Δ) with the line 2θ−K+(N−1)Δwinh. Figure 7 shows an example for the case K=2θand N=2. The maximum number of stable bumps that the specific connection Fig. 7 Existence of None-bump solutions in the two-field model with the Mexican-hat kernel (2). Threshold θ=0.5andK=1. Straight lines show the condition (25)forN=2andN=6. Parameters of the kernel as in Fig. 2 function supports is N=6. A further increase of lateral inhibition would destabilizes the bumps since dW(Δ)/dΔ>0. The analogous condition for the existence of an N-bump solution in the Amari model is given by 123 Biological Cybernetics (2021) 115:451–471 459 Fig. 8 Solutions at time t=50 of the model (1) created with transient inputs I(x,t)given by (3) with variation of (b) input strength AI,(c) input width σI and (d) input duration dI. Parameters of the inputs: (a) AI1=1,σI1=1,dI=1(b) AI1=3,σI1=1,dI=1(c) AI1=1,σI1=6,dI=1(d) AI1=1,σI1=1,dI=3. The kernel wis given by (2)with Aex =2, Ain =1, σex =1.25, σin =2.5andwinh =0.1. Threshold θ=0.5 (dotted line), K=0 F(Δ) =−θ+W(Δ) −(N−1)Δ winh =0.(26) For the example of Fig. 7, the maximum number of bumps reduces in this case to N=3. In the next section we show that different to the Amari model with Mexican-hat connectivity, the two-field model supports multi-bump solutions created by inputs at any distance. 4 Input driven bump solutions In working memory applications, classical bump attractor models implement the encoding and maintenance of briefly presented visual input as an all-or-none process. The information is either stored in a self-stabilized bump with a shape defined by the recurrent interactions or the information is completely lost since neural activity decays back to resting state. However, a growing body of experimental evidence indicates that this binary classification may be insufficient as a description of WM storage, because the representational quality of WM items may differ depending on sensory attributes of the input or task demands. For instance, the salience of the visual input manipulated through varying stimulus strength (or contrast) is known to affect WM representations. Consistent with the notion of memory-associated persistent activity, a graded modulation of sustained population activity as a function of stimulus contrast has been described in a WM task (Constantinidis et al. 2001). Figure 8shows various bump solutions of the two-field integrator model in which input attributes are systematically varied. In panels (a) and (b), the attractor network converts a transient input into a persistent output which is proportional to the input amplitude. A comparison of panels (a) and (c) illustrates that the bump also represents faithfully input width. This property can be exploited for instance to encode uncertainty in the sense of an explicit probability code (Ma and Jazayeri 2014). Higher uncertainty is then represented by a wider activation pattern across the neural population tuned to a continuous input feature (e.g., representing ranges of movement direction (Erlhagen and Schöner 2002)). In panel (d), the localized input of panel (a) is applied but lasting three times as long. The amplitude of the evolving bump equals the bump amplitude in panel (b), as expected by perfect integration. The two-field model can thus be used to model the encoding of input duration (Wojtak et al. 2019)if one assumes that the input amplitude is first normalized by upstream sensory processing (Carandini and Heeger 2012). Figure 9shows the temporal evolution of a bump in two dynamic input scenarios which can be expected in any natural environment. In the first example (a), the strength of evidence continuously increases during the stimulation period. The bump amplitude reflects again the total external input applied to the population. The observed ramping-like activity is a hallmark of “drift-diffusion” models of decision making which assume the continuous accumulation of sensory and other evidence to a threshold (for review see Bogacz 2007). To investigate the dependence of the bump shape on input width, a truncated Gaussian input (I(x,t)=0,|x|>1.5σ) has been applied. As can be seen in panel (b), the bump width reaches its maximum well before the end of the stimulation period. For a read-out of the population activity by a downstream network this means that a gain in certainty about the most likely parameter value as expressed in the bump amplitude is not necessarily accompanied by more uncertainty 123 466 Biological Cybernetics (2021) 115:451–471 Fig. 17 Bump diffusion in the stochastic two-field model. (a)Bump solution at time t=60. The initial condition is a bump centered at x=0. (b) Space-time plot showing the bump drift. The magenta line indicates the bump center. (c) Variance of bump position computed for N=100 numerical simulations of the model with different input strengths AI. Remaining input parameters σI1=1, dI=1. The kernel wis given by (2) with Aex =2, Ain =1, σex =1.5, σin =2.5and winh =0.3. Threshold θ=0.25 (dotted line in panel (a)), K=0, =0.005, L=π. The spatial and temporal resolution is dx=0.01 and dt=0.01, respectively a small reduction of the input distance (2 dx, where dxis the numerical mesh size) causes attraction (Almeida et al. 2015). Figure 18 (a) shows overlaid activity profiles of an evolving two-bump solution taken at a fixed time in different model simulations. As can be clearly seen, the additive noise causes a switching between either attraction or repulsion in different trials. Panel (b) shows that the magnitude of the interference effect appears to be greatly reduced when stronger inputs are applied. Panel (c) shows this dependence on input strength/bump amplitude in a more quantitative manner by plotting the variance of the inter-peak distance as a function of time. To allow a direct comparison, Figure 19 depicts, for the same input distance, the activity pattern at a fixed time in two runs of the stochastic Amari model. It develops either a single bump (a) or two bumps with repulsion (b). The model simulations reveal that the single bump may be the result of a suppressive interaction or a fast merging of the input-induced activity patterns. 4.4 Gated integrator model The two-field model performs a continuous integration of input streams of any strength. This is a desirable feature for decision making since also weak evidence may bias decision processes involving multiple alternatives. One potential problem with the perfect integrator mechanism is that the pure accumulation of noise at field sites without specific input may in principle lead to the creation of bumps. Other neuro-inspired integrator models solve this problem by implementing a preset threshold mechanism that controls the minimum level of evidence entering the accumulation process (Cain et al. 2013; Koulakov et al. 2002). We propose a model extension inspired by experimental findings showing that neurons in local circuits may rapidly change their functional connectivity according to behavioral demands (Haider and McCormick 2009). In tasks requiring the accumulation of evidence, this modulation should affect the balanced excitation and inhibition in the local feedback loops responsible for the stabilization of population activity at any level. Only neurons that appear to be pre-activated by input above the noisy baseline should participate in these loops. The cellular and network mechanism by which functional connections are rapidly formed and broken are currently an active area of research (Haider and McCormick 2009). We apply the computationally simple idea of a state-dependent gating of the local coupling between the two fields by a threshold function g(u). A step function with threshold κ≤θis used for simplicity. The model reads ∂u(x,t) ∂t=−u(x,t)+v(x,t)g(u(x,t)−κ) +I(x,t)+1/2dW(x,t) +Ω w(|x−y|)f(u(y,t)−θ)dy, (32a) ∂v(x,t) ∂t=−v(x,t)+u(x,t)g(u(x,t)−κ) −Ω w(|x−y|)f(u(y,t)−θ)dy. (32b) The gating mechanism delays the perfect integration of evidence until a target location is selected by an initial input of strength AI>κ, rendering a bump creation by accumulated (weak) noise unlikely. Figure 20 compares the pattern formation in response to an input at x=0 in the stochastic model without (a) and with gating mechanism (c). Without gating, in addition to the input-driven bump at x=0 various noise-induced bumps at other locations develop. The impact on noise integration can be seen when comparing the sum of u(x)and v(x), representing the total of accumulated external input and noise (b). For the gating case, the sum represents the bell-shaped input slightly corrupted by noise and small random fluctuations at other field sites (d). It is important to notice that the gating mechanisms does not destroy the exis123 Biological Cybernetics (2021) 115:451–471 467 Fig. 18 Bump attraction/repulsion due to noise in the two-field model. Sample profiles of the u-population at time t=30 created with two inputs with AI1,2=1(a)andAI1,2=3(b). The percentage of “merged” solutions for inputs with amplitudes AI∈{1,2,3}in N=500 trials was respectively 60%, 98% and 100%. (c) Variance of the inter-peak distance as a function of time computed for N=500 numerical simulations of the model for different input strengths AI. The inputs with amplitudes AI1,2=1, AI1,2=2, AI1,2=3, are applied at positions xc1,2=±2.25. Remaining parameters of the inputs are σI1,2=1and dI=1. The kernel wis given by (2) with Aex =2, Ain =1, σex =1.5, σin =2.5andwinh =0.1. Threshold θ=0.25 (dotted line), K=0, =0.0025, L=3π,dx=0.01, dt=0.01 Fig. 19 Bump merging (a) and repelling (b) in the Amari model due to noise. The inputs with amplitudes AI1,2=1 are applied at positions xc1,2=±2.25. Remaining parameters of the inputs are σI1,2=1 and dI=1. Vertical lines indicate the sites xc1,2where the inputs were applied. The kernel wis given by (2)withAex =2, Ain =1, σex =1.5, σin =2.5andwinh =0.1. Threshold θ=0.25 (dotted line), =0.0025, L=3π,dx=0.01, dt=0.01 tence of subthreshold bumps since for sites xwith u(x)>κ a continuous integration of weak inputs is still supported. The gating mechanism is also important for resetting stimulusselective persistent activity to a homogeneous resting state. In WM applications of continuous attractor networks, it is typically assumed that a switch between a bump attractor and a stable resting state can be achieved by applying a transient inhibitory input to all neurons (Camperi and Wang 1998). Figure 21 illustrates this “forgetting” mechanism in a model simulation without noise. The neural integrator first develops a bump in response to a localized excitatory input which is followed by a strong decrease in activation when at time t=21 the homogeneous inhibition is applied to the u-population. Functionally, the two fields become decoupled since g(u)<κ,∀x, and the subthreshold activity in both layers converges to the homogeneous resting state, u=v=0, as predicted by the Amari model (Amari 1977). Without the gating mechanism, the field dynamics would stabilize a subthreshold bump. 5 Discussion Here we have presented a bump attractor model which is able to sustain localized activity patterns evoked by external inputs which differ in shape and amplitude. The model simulations demonstrate that the network dynamics not only supports the encoding of a single input value but also the storage of certain stimulus characteristics. This is in line with a growing body of experimental evidence showing a correlation between the level of persistent delay activity and stimulus attributes such as intensity/contrast (Constantinidis et al. 2001) or duration (Brody et al. 2003a). The possibility to stabilize localized activity patterns with a continuum of amplitudes can be in principle also used to model the experimentally observed dependence of persistent activity on task demands such as the amount of training (Qi et al. 2011; Histed et al. 2009) or varying choice certainty (Kiani and Shadlen 2009; Basso and Wurtz 1997). We have recently applied an extension of the two-field model to two spatial dimensions as part of a neuro-inspired robot control architecture for natural human-robot interactions. It implements the learning and 123 468 Biological Cybernetics (2021) 115:451–471 Fig. 20 Simulations of the two-field model (32)with different threshold values κand the Mexican-hat connectivity given by (2) with Aex =3, Ain =1.5, σex =1.4, σin =3 and winh =0.2, θ=0.75, K=1. The spatial fluctuations are given by a cosine-correlated noise with =0.005. The bump at x=0 is created with a localized input with AI1=2, σI1=1.5, dI=1. (a,b)no gating; (c,d)κ=0.15 Fig. 21 Space-time plot showing the forgetting mechanism. A bump is first triggered by a localized input applied at t=1. This is followed by a decay to resting state when at time t=21 a homogeneous inhibitory input is applied to the u-population of (32). The connectivity function isgivenby(2) with Aex =3, Ain =1.5, σex =1.4, σin =3and winh =0.2. The localized input is defined by AI1=1, σI1=1.5, dI=1 and the homogeneous inhibition by A=−1.5anddI=1. Thresholds κ=θ=0.5, K=1 adaptation of action planning and value-based decision making in a dynamic environment (Wojtak et al. 2021). The model works as a robust neural integrator since the bump amplitude faithfully reflects the accumulation of input signals over time. This contrasts with the behavior of most attractor networks applied in working memory and decision making tasks in which the stable state corresponds to a stereotyped pattern of activity exclusively determined by the recurrent interactions within the network. The input dependence is achieved by combining a lateral-inhibition type network connectivity with local inhibitory and excitatory feedback mechanisms that are able to stabilize input-induced changes of the population activity at any level. A biophysically mechanistic formulation of the processing principles based on two coupled field equations of Amari-type is amenable to mathematical analysis. Being able to rigorously analyze the existence and stability of bump solutions and their dependence on model parameters is a major advantage over more complex, biophysically realistic models when trying to explain the behavior of cognitive agents. Using an identical parameter set allowed us to directly compare the performance of the two-field model with the classical Amari model with Mexican-hat coupling. The network with the balanced local feedback exhibits several distinct features. First, two or more stable regions of excitation may exist at a distance where a gradient in the lateral inhibition profile predicts a merging or drifting of the bumps in the Amari model. Moreover, the lateral feedback excitation of the two-field model facilitates in general the formation of memory bumps. This is particularly evident in a sequential stimulus protocol in which lateral inhibition in the Amari model may prevent the encoding of an input in a stable activity pattern. This is true even for the case of a spatially uniform lateral inhibition for which a multi-bump solution exists. The two-field model thus predicts a higher working memory capacity without the need to refer to additional processing mechanism like topdown excitation (Edin et al. 2009) or a stimulus-selective synaptic facilitation (Mongillo et al. 2008). Second, the activation pattern induced by two nearby inputs do not merge completely into a single bump at an intermediate position but displays a wide plateau in its profile. This result is consistent with the prediction of biophysically realistic bump attractor networks(Weietal.2012) and has been described in electrophysiological stimulation studies (Vokoun et al. 2014). The higher uncertainty expressed by the wider activity distribution could be used by a downstream read-out system in cue integration and decision making tasks. Third, the two-field model supports the existence of stable subthreshold activity 123 Biological Cybernetics (2021) 115:451–471 469 patterns. In the context of the retrocueing paradigm, we have argued that such patterns might represent less-salient items that become suppressed below threshold by means of lateral inhibition from bumps representing cued items. The experimental observation that neural memory representations often degrade to a large extent during the delay period has been taken as evidence for an active-silent state of working memory encoded in synaptic changes (Mongillo et al. 2008). The model simulations however suggest that postulating such a latent memory state might not be necessary at least for the retrocueing task (Schneegans and Bays 2017). The information represented by weak sustained activity is not lost but can be restored with a spatially unspecific cue. Unlike in bistable attractor models of Amari type, the baseline activity of the neural integrator is not a stable state. A forgetting mechanism can thus not be implemented by simply applying a sufficiently strong inhibitory input that destabilizes existing bumps. For working memory applications, we have proposed a simple gating mechanism for the local feedback which restores the stable resting state of the Amari model without destroying the existence of subthreshold bump solutions. The possibility to represent stimulus attributes in the bump shape allows us to make predictions about the impact of for instance input strength (or contrast) on working memory performance. Simulations of the stochastic version of the two-field model reveal that the noise-induced drift of a single bump in the continuous attractor network decreases with increasing bump amplitude associated with stronger inputs. Larger bump amplitudes also greatly reduce the interaction effects of two input-induced bumps at a critical distance which can be described as repulsion or attraction when considering the relative peak positions. The neural integrator model shares the prediction that larger bumps provide more faithful input representations over time with the neural field model of lateral inhibition type introduced by Carroll et al. (2014). It consists of separate excitatory and inhibitory populations and implements a tight balance of excitation and inhibition within the network to support a continuum of possible bump amplitudes. Importantly, achieving this balance requires that the recurrent excitation is inversely proportional to the gain of a piecewise linear transfer function. Any change in the tuning of the form of the nonlinearity or the connectivity will destroy this balance. The two-field model is robust to such changes since the same connectivity and nonlinear transfer functions are used for the u-field and the v-field. An important consequence is that the assumption of a perfectly symmetric connectivity pattern in continuous bump attractor networks can be relaxed. Model simulations with a slightly asymmetric connectivity profile or a profile with random perturbations show that localized activity patterns with a biased or distorted shape are stabilized by the network dynamics. Mitigating the inherent structural instability of bump attractor models offers new perspective for applications of dynamic field theory that include for instance the learning of the connectivity pattern (Fotouhi et al. 2015). We will explore this line of research in future work. The two-field model shares the robustness to global changes in the neural gain and perturbations of the translationinvariant network structure with the firing rate model of working memory proposed by Lim and Goldman (2013). It consists of pools of interconnected excitatory and inhibitory neurons that are able to maintain graded levels of spatially tuned patterns of persistent activity without widespread lateral inhibition. It implements principles of negativederivative feedback control by assuming that positive feedback mediated by recurrent excitation and negative feedback mediated by recurrent inhibition have equal strength, but inhibition operates on a faster timescale. Perturbations that affect the tight excitation/inhibition balance disrupt the memory performance. Also the two-field model is sensitive to perturbations in the balance of recurrent inputs from neighboring neurons and the local inhibitory feedback that a neuron receives. For instance, applying a nonlinear response function with saturation for the v-population will cause a saturation limit of the input integration. However, it can be expected that lateral inhibition acts as a stabilizing factor if the imbalance is not too strong. Further work is needed to systematically investigate how changes applied to the gain, the amplitude or the spatial ranges of the recurrent interactions in just one of the fields affect the pattern formation process. Acknowledgements The work received financial support from FCT through the PhD fellowship PD/BD/128183/2016 the project “Neurofield” (PTDC/MAT-APL/31393/2017) and the research centre CMAT within the project UID/MAT/00013/2020. Availability of data and material Not applicable. Declarations Conflicts of interest The authors declare that they have no conflict of interest. Code availability Example codes implemented in MATLAB are available at https://github.com/w-wojtak/A-dynamic-neural-field-modelof-continuous-input-integration. A Numerical simulations Numerical simulations of the model were done in MATLAB using a forward Euler method with parameters L=30, N=12000, dx =2L/N=0.005, T=100, M=10000, dt =T/M=0.01, unless stated otherwise in specific examples. Numerical simulations of the stochastic model were done using the Euler-Maruyama method. To compute the 123 470 Biological Cybernetics (2021) 115:451–471 spatial convolution of wand fwe employ a fast Fourier transform (FFT), using MATLAB’s in-built functions fft and ifft to perform the Fourier transform and the inverse Fourier transform, respectively. Periodic boundary conditions are used. By choosing a sufficiently large domain size, we make sure that the localized patterns evolve sufficiently far from the boundaries. For performing numerical continuation, we use the method described in (Rankin et al. 2014) and adapt MATLAB code available in (Avitabile 2016). The main advantage of this method is that it can be applied directly to the full integral model. This is possible due to the usage of NewtonGMRES solvers combined with a fast Fourier transform (FFT) employed for computing the convolution term (Rankin et al. 2014). References Almeida R, Barbosa J, Compte A (2015) Neural circuit basis of visuospatial working memory precision: a computational and behavioral study. J Neurophysiol 114(3):1806–1818 Amari S (1977) Dynamics of pattern formation in lateral-inhibition type neural fields. Biol Cybern 27(2):77–87 Amari SI (1980) Topographic organization of nerve fields. Bull Math Biol 42(3):339–364 Avitabile D (2016) Numerical computation of coherent structures in spatially-extended systems. Second International Conference on Mathematical Neuroscience, Antibes Juan-les-Pins, 2016 Basso MA, Wurtz RH (1997) Modulation of neuronal activity by target uncertainty. Nat 389(6646):66–69 Bays PM, Taylor R (2018) A neural model of retrospective attention in visual working memory. Cogn Psychol 100:43–52 Bicho E, Mallet P, Schöner G (2000) Target representation on an autonomous vehicle with low-level sensors. Int J Robot Res 19(5):424–447 Bogacz R (2007) Optimal decision-making theories: linking neurobiology with behaviour. Trends Cogn Sci 11(3):118–125 Brody CD, Hernández A, Zainos A, Romo R (2003a) Timing and neural encoding of somatosensory parametric working memory in macaque prefrontal cortex. Cerebral Cortex 13(11):1196–1207 Brody CD, Romo R, Kepecs A (2003b) Basic mechanisms for graded persistent activity: discrete attractors, continuous attractors, and dynamic representations. Current Opin Neurobiol 13(2):204–211 Cain N, Barreiro AK, Shadlen M, Shea-Brown E (2013) Neural integrators for decision making: a favorable tradeoff between robustness and sensitivity. J Neurophysiol 109(10):2542–2559 Camperi M, Wang XJ (1998) A model of visuospatial working memory in prefrontal cortex: recurrent network and cellular bistability. J Comput Neurosci 5(4):383–405 Carandini M, Heeger DJ (2012) Normalization as a canonical neural computation. Nat Rev Neurosci 13(1):51–62 Carroll S, Josi´c K, Kilpatrick ZP (2014) Encoding certainty in bump attractors. J Comput Neurosci 37(1):29–48 Compte A, Brunel N, Goldman-Rakic PS, Wang XJ (2000) Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model. Cerebral Cortex 10(9):910– 923 Constantinidis C, Wang XJ (2004) A neural circuit basis for spatial working memory. Neurosci 10(6):553–565 Constantinidis C, Franowicz MN, Goldman-Rakic PS (2001) The sensory nature of mnemonic representation in the primate prefrontal cortex. Nat Neurosci 4(3):311–316 Curtis CE, Lee D (2010) Beyond working memory: the role of persistent activity in decision making. Trends Cogn Sci 14(5):216–222 Denève S, Machens CK (2016) Efficient codes and balanced networks. Nat Neurosci 19(3):375–382 Edin F, Klingberg T, Johansson P, McNab F, Tegnér J, Compte A (2009) Mechanism for top-down control of working memory capacity. Proc Nat Acad Sci 106(16):6802–6807 Erlhagen W, Bicho E (2006) The dynamic neural field approach to cognitive robotics. J Neural Eng 3:36–54 Erlhagen W, Schöner G (2002) Dynamic field theory of movement preparation. Psychol Rev 109(3):545 Erlhagen W, Bastian A, Jancke D, Riehle A, Schöner G (1999) The distribution of neuronal population activation (DPA) as a tool to study interaction and integration in cortical representations. J Neurosci Methods 94(1):53–66 Ferreira F, Erlhagen W, Bicho E (2016) Multi-bump solutions in a neural field model with external inputs. Physica D: Nonlinear Phenomena 326:32–51 Ferreira F, Wojtak W, Sousa E, Louro L, Bicho E, Erlhagen W (2020) Rapid learning of complex sequences with time constraints: A dynamic neural field model Fotouhi M, Heidari M, Sharifitabar M (2015) Continuous neural network with windowed Hebbian learning. Biol Cybern 109(3):321– 332 French DA (2004) Identification of a free energy functional in an integro-differential equation model for neuronal network activity. Appl Math Lett 17(9):1047–1051 Gazzaley A, Nobre AC (2012) Top-down modulation: bridging selective attention and working memory. Trends Cogn Sci 16(2):129–135 Haider B, McCormick DA (2009) Rapid neocortical dynamics: cellular and network mechanisms. Neuron 62(2):171–189 Hebb DO (1949) The organization of behavior: a neuropsychological theory. Chapman & Hall, J Wiley, New York Histed MH, Pasupathy A, Miller EK (2009) Learning substrates in the primate prefrontal cortex and striatum: sustained activity related to successful actions. Neuron 63(2):244–253 Jancke D, Erlhagen W, Dinse HR, Akhavan AC, Giese M, Steinhage A, Schöner G (1999) Parametric population representation of retinal location: Neuronal interaction dynamics in cat primary visual cortex. J Neurosci 19(20):9016–9028 Johnson JS, Spencer JP, Luck SJ, Schöner G (2009) A dynamic neural field model of visual working memory and change detection. Psychol Sci 20(5):568–577 Kiani R, Shadlen MN (2009) Representation of confidence associated with a decision by neurons in the parietal cortex. Sci 324(5928):759–764 Kilpatrick ZP, Ermentrout B (2013) Wandering bumps in stochastic neural fields. SIAM J Appl Dyn Syst 12(1):61–94 Koulakov AA, Raghavachari S, Kepecs A, Lisman JE (2002) Model for a robust neural integrator. Nat Neurosci 5(8):775–782 Krishnan N, Poll DB, Kilpatrick ZP (2018) Synaptic efficacy shapes resource limitations in working memory. J Comput Neurosci 44(3):273–295 Kubota S, Aihara K (2005) Analyzing global dynamics of a neural field model. Neural Process Lett 21(2):133–141 Kuo BC, Stokes MG, Nobre AC (2012) Attention modulates maintenance of representations in visual short-term memory. J Cogn Neurosci 24(1):51–60 Laing CR, Troy WC (2003) Two-bump solutions of Amari-type models of neuronal pattern formation. Physica D: Nonlinear Phenomena 178(3–4):190–218 123 Biological Cybernetics (2021) 115:451–471 471 Lim S, Goldman MS (2013) Balanced cortical microcircuitry for maintaining information in working memory. Nat Neurosci 16(9):1306– 1314 Lin PH, Luck SJ (2009) The influence of similarity on visual working memory representations. Vis Cogn 17(3):356–372 Ma WJ, Jazayeri M (2014) Neural coding of uncertainty and probability. Annual Rev Neurosci 37:205–220 Masse NY, Rosen MC, Freedman DJ (2020) Reevaluating the role of persistent neural activity in short-term memory. Trends Cogn Sci 24(3):242–258 Miller EK, Cohen JD (2001) An integrative theory of prefrontal cortex function. Annual Rev Neurosci 24(1):167–202 Mongillo G, Barak O, Tsodyks M (2008) Synaptic theory of working memory. Sci 319(5869):1543–1546 Murray JD, Anticevic A, Gancsos M, Ichinose M, Corlett PR, Krystal JH, Wang XJ (2014) Linking microcircuit dysfunction to cognitive impairment: effects of disinhibition associated with schizophrenia in a cortical working memory model. Cereb Cortex 24(4):859–872 Owen MR, Laing CR, Coombes S (2007) Bumps and rings in a twodimensional neural field: splitting and rotational instabilities. New J Phys 9(10):378 Qi XL, Meyer T, Stanford TR, Constantinidis C (2011) Changes in prefrontal neuronal activity after learning to perform a spatial working memory task. Cereb Cortex 21(12):2722–2732 Rankin J, Avitabile D, Baladron J, Faye G, Lloyd DJ (2014) Continuation of localized coherent structures in nonlocal neural field equations. SIAM J Sci Comput 36(1):B70–B93 Rhodes BJ, Bullock D, Verwey WB, Averbeck BB, Page MP (2004) Learning and production of movement sequences: Behavioral, neurophysiological, and modeling perspectives. Hum Mov Sci 23(5):699–746 Rolls ET, Dempere-Marco L, Deco G (2013) Holding multiple items in short term memory: a neural mechanism. PLOS ONE 8(4):e61078 Rose NS, LaRocque JJ, Riggall AC, Gosseries O, Starrett MJ, Meyering EE, Postle BR (2016) Reactivation of latent working memories with transcranial magnetic stimulation. Sci 354(6316):1136–1139 Sakai Y, Okamoto H, Fukai T (2006) Computational algorithms and neuronal network models underlying decision processes. Neural Netw 19(8):1091–1105 Schneegans S, Bays PM (2017) Restoration of fMRI decodability does not imply latent working memory states. J Cogn Neurosci 29(12):1977–1994 Schöner G (2016) Dynamic thinking: A primer on dynamic field theory. Oxford University Press, Oxford Seung HS (1998) Continuous attractors and oculomotor control. Neural Netw 11(7–8):1253–1258 Soto D, Mäntylä T, Silvanto J (2011) Working memory without consciousness. Curr Biol 21(22):R912–R913 Souza AS, Oberauer K (2016) In search of the focus of attention in working memory: 13 years of the retro-cue effect. Atten Percept Psychophys 78(7):1839–1860 Sprague TC, Ester EF, Serences JT (2016) Restoring latent visual working memory representations in human cortex. Neuron 91(3):694– 707 Sutterer DW, Foster JJ, Adam KC, Vogel EK, Awh E (2019) Itemspecific delay activity demonstrates concurrent storage of multiple active neural representations in working memory. PLOS Biol 17(4):e3000239 Trübutschek D, Marti S, Ojeda A, King JR, Mi Y, Tsodyks M, Dehaene S (2017) A theory of working memory without consciousness or sustained activity. Elife 6:e23871 Vokoun CR, Huang X, Jackson MB, Basso MA (2014) Response normalization in the superficial layers of the superior colliculus as a possible mechanism for saccadic averaging. J Neurosci 34(23):7976–7987 Wang XJ (2003) Persistent neural activity: experiments and theory. Cereb Cortex 13(11):1123–1123 Wang Z, Kruijne W, Theeuwes J (2012) Lateral interactions in the superior colliculus produce saccade deviation in a neural field model. Vis Res 62:66–74 Wei Z, Wang XJ, Wang DH (2012) From distributed resources to limited slots in multiple-item working memory: a spiking network model with normalization. J Neurosci 32(33):11228–11240 Wildegger T, Humphreys G, Nobre AC (2016) Retrospective attention interacts with stimulus strength to shape working memory performance. PLOS ONE 11(10):e0164174 Wilimzig C, Schneider S, Schöner G (2006) The time course of saccadic decision making: Dynamic field theory. Neural Netw 19(8):1059– 1074 Wimmer K, Nykamp DQ, Constantinidis C, Compte A (2014) Bump attractor dynamics in prefrontal cortex explains behavioral precision in spatial working memory. Nat Neurosci 17(3):431–439 Wojtak W, Coombes S, Bicho E, Erlhagen W (2016) Combining spatial and parametric working memory in a dynamic neural field model. In: Artificial Neural Networks and Machine Learning — ICANN 2016, Lecture Notes in Computer Science, Springer, vol 9886, pp 411–418 Wojtak W, Ferreira F, Bicho E, Erlhagen W (2019) Neural field model for measuring and reproducing time intervals. In: Tetko IV, K˚urková V, Karpov P, Theis F (eds) Artificial Neural Networks and Machine Learning - ICANN 2019: Theoretical Neural Computation. Springer International Publishing, Cham, pp 327–338 Wojtak W, Ferreira F, Vicente P, Louro L, Bicho E, Erlhagen W (2021) A neural integrator model for planning and value-based decision making of a robotics assistant. Neural Comput Appl 33:3737–3756 Wu S, Hamaguchi K, Amari Si (2008) Dynamics and computation of continuous attractors. Neural Comput 20(4):994–1025 Zylberberg J, Strowbridge BW (2017) Mechanisms of persistent activity in cortical circuits: possible neural substrates for working memory. Annual Rev Neurosci 40:603–627 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123