Turing Computation with Recurrent Artificial Neural Networks

Abstract

We improve the results by Siegelmann & Sontag [1, 2] by providing a novel and parsimonious constructive mapping between Turing Machines and Recurrent Artificial Neural Networks, based on recent developments of Nonlinear Dynamical Automata. The architecture of the resulting R-ANNs is simple and elegant, stemming from its transparent relation with the underlying NDAs. These characteristics yield promise for developments in machine learning methods and symbolic computation with continuous time dynamical systems. A framework is provided to directly program the R-ANNs from Turing Machine descriptions, in absence of network training. At the same time, the network can potentially be trained to perform algorithmic tasks, with exciting possibilities in the integration of approaches akin to Google DeepMind's Neural Turing Machines.

Introduction

Methods

In this section we outline a mapping from Turing machines to R-ANNs. Our construction involves two stages. In the first stage a Generalized Shift [ 3 ] emulating a Turing Machine is built, and its dynamics encoded on the unit square via a procedure called Go¨delization, defining a piecewiseaffine linear map on the unit square, i.e. a NDA. In the second stage, the resulting NDA is mapped onto a first-order R-ANN. Next, the theoretical methods employed are discussed in detail. 2.1

Turing Machines

A Turing Machine [ 6 ] is a computing device endowed with a doubly-infinite one-dimensional tape (memory support with one symbol capacity at each memory location), a finite state controller and a read-write head that follows the instructions encoded by a transition function. At each step of the computation, given the current state and the current symbol read by the read-write head, the machine controller determines via the writing of a symbol on the current memory location, a shift of the read-write head to the memory location to the left (L) or to the right (R) of the current one, and the transition to a new state for the next computation step. At a computation step, the content of the tape together with the position of the read-write head and the current controller state define a machine configuration.

More formally, a Turing Machine is a 7-tuple MTM = (Q; N; T; q0; t; F; ), where Q is a finite set of control states, N is a finite set of tape symbols containing the blank symbol t, T N n ftg is the input alphabet, q0 is the starting state, F Q is a set of ‘halting’ states and is a partial transition function, determining the dynamics of the machine. In particular, is defined as follows: : Q

N ! Q

N fL; Rg: 2.2

Dotted sequences and Generalized Shifts

A Turing machine configuration can be described by a bi-infinite dotted sequence on some alphabet A; it can then be defined as:

s = : : : di 3 di 2 di 1 :di0 di1 di2 : : : ; where l = : : : di 3 di 2 describes the part of the tape on the left of the read-write head, r = di0 di1 di2 : : : describes the part on its right, q = di 1 describes the current state of the machine controller, and the dot denotes the current position of the read-write head, i.e. the symbol to its right. The central dot splits the tape into two one-sided infinite strings 0; , where 0 is the left part of the dotted sequence in reverse order. The first symbol in represents the current state of the Turing Machine, whereas the first symbol in represents the symbol currently under the controller’s head. The transition function can be straightforwardly extended to a function ^ operating on dotted sequences, so that ^ : AZ ! AZ.

A Generalized Shift acts on dotted sequences, and is defined as a pair MGS = (AZ; ), with AZ being the space of dotted sequences, : AZ ! AZ defined by where shifts the symbols to the left or to the right, or does not shift them at all, as determined by the function F (s). In addition, the Generalized Shift can operate a substitution, with G(s) being the function which substitutes a substring of length e in the Domain of Effect (DoE) of s with a new substring. Both the shift and the substitution are functions of the content of the Domain of Dependence (DoD), a substring of s of length `.

A Turing Machine can be emulated by a Generalized Shift with DoD = DoE = di 2 di 1 :di0 and the functions F; G appropriately chosen such that (s) = ^(s) for all s (see [ 7 ] for a detailed exposition). with (s) =

F (s)(s

G(s)) F : AZ ! Z G : AZ ! Ae (1) (2) (3) (4) (5) Go¨del codes (or Go¨delizations) [ 8 ] map strings to numbers and, in particular, allow the mapping of the space of one-sided infinite sequences to the real interval [0; 1]. Let AN be the space of one-sided infinite sequences over an alphabet A, s be an element of AN, rk the k-th symbol in s, : A ! N a one-to-one function associating each symbol in the alphabet A to a natural number, and g the number of symbols in A. Then a Go¨delization is a mapping from AN to [0; 1] R defined as: Conveniently, Go¨delization can be employed on a Turing machine configuration, represented as a dotted sequence : 2 AZ. The Go¨del encoding x and y of 0 and define a representation of s ( x( 0); y( )) known as symbol plane or symbologram representation, which is contained in the unit square [0; 1]2 R2. The choice of encoding x and y to use on the machine configurations is arbitrary. Therefore, to enable the construction of parsimonious Nonlinear Dynamical Automata our encoding will assume that always contains tape symbols only, and that the first symbol of 0 is always a state symbol, the rest being tape symbols only. Based on these assumptions, the particular encoding is defined as: (6) (7) 1 x( 0) = q(a1)nq 1 + X s(ak+1)ns knq 1;

k=1 1 y( ) = X s(bk)ns k;

k=1 with nq = jQj, i.e. the number of states in the Turing Machine, ns = jNj, i.e. the number of tape symbols in the Turing Machine, q and s enumerating Q and N respectively, and with ak and bk being the k-th symbol in 0 and respectively. 2.3.1

Encoded Generalized Shift and affine-linear transformations

The substitution and shift operated by a Generalized Shift on a dotted sequence s = : can be represented as an affine-linear transformation on ( x( 0); y( )), i.e. the symbologram representation of s. In particular, a substitution and shift on a dotted sequence can be broken down into substitutions and shifts on its one-sided components. In the following, we will show how substitutions and shifts on a one-sided infinite sequence can be represented as affine-linear transformations on its Go¨delization. These results will be useful in showing how the symbologram representation of a Generalized Shift leads to a piecewise affine-linear map on a rectangular partition of the unit square. Let s = d1d2d3 : : : be a one-side infinite sequence on some alphabet A. Substituting the n-th symbol in s with d^n yields s^ = d1 : : : dn 1d^ndn+1 : : :, so that (s) = (d1)g 1 + : : : (dn 1)g (n 1) + (dn)g n + (dn+1)gn+1 + : : : ; (s^) = (d1)g 1 + : : : (dn 1)g (n 1) + (d^n)g n + (dn+1)gn+1 + : : : ; = (s)

(dn)g n + (d^n)g n: As the previous example illustrates, Go¨delizing a sequence resulting from a symbol substitution is equivalent to applying an affine-linear transformation on the original Go¨delized sequence. In particular, the parameters of the affine-linear transformation only depend on the position and identities of the symbols involved in the substitution. Shifting s to the left by removing its first symbol or shifting it to the right by adding a new one yields respectively sl = d2d3d4 : : : and sr = b d1d2d3d4 : : :, where b is the newly added symbol. In this case

(sl) = (d2)g 1 + (d3)g 2 + (d4)g 3 + : : : and = g (s)

(d1); = g 1 (s) + (b)g 1: Again, the resulting Go¨delized shifted sequence can be obtained by applying an affine-linear transformation to the original Go¨delized sequence. 2.4

Nonlinear Dynamical Automata

A Nonlinear Dynamical Automaton (NDA) is a triple MNDA = (X; P; ), with P being a rectangular partition of the unit square, that is

P = fDi;j

Xj 1 i m; 1 j n; m; n 2 Ng; so that each cell Di;j is defined as the cartesian product Ii Jj , with Ii; Jj intervals for each bi-index (i; j), Di;j \ Dk;l = ; if (i; j) 6= (k; l), and Si;j Di;j = X. The couple (X; ) is a time-discrete dynamical system with phase space X = [0; 1]2 unit square) and with flow : X ! X, a piecewise affine-linear map such that Specifically, i;j takes the following form: [0; 1] being real jDi;j :=

R2 (i.e. the i;j . i;j (x) = ai;j

x ai;j y + i;j x 0 0 i;j y x y tThheeappipercoepwriisaeteabffirannec-lhin,eaanrdmthaups thaelsaoprperqoupirrieasteadsywniatcmhiicnsg, rauslea fu(nxc;tyio)n2ofJ1th;mecKurreJ1n;t nstKattoe.seTlhecatt is, (x; y) = i;j (x; y) () (x; y) = (i; j).

Each cell Di;j of the partition P of the unit square can be seen as comprising all the Go¨delized dotted sequences that contain the same symbols in the Domain of Dependence. That is, for a Generalized Shift simulating a Turing Machine, the first two symbols in 0 and the first symbol in . The unit square is thus partitioned in a number of I intervals equal to m = nqns, and one of J intervals equal to n = ns, with nq being the number of states in Q and ns the number of symbols in N, for a total of nqns2 cells. As each cell corresponds to a different Domain of Dependence of the underlying Generalized Shift in symbolic space, it is associated with a different affine-linear transformation representing the action of a substitution and shift in vector space. The transformation parameters (aix;j ; aiy;j ) and ( ix;j ; iy;j ) can be derived using the methods outlined in subsubsection 2.3.1. Thus, a Turing Machine can be represented as a Nonlinear Dynamical Automaton by means of its Go¨delized Generalized Shift representation. 3

NDAs to R-ANNs

The aim of the second stage of our methodology is to map the orbits of the NDA (i.e. i;j (x; y)) to orbits of the R-ANN, which we will denote by i;j (x; y).

Let ( ) denote the proposed map. Its role is to encode the affine-linear dynamics at each i;j branch in the architecture and weights of the network, and emulate the overall dynamics by suitably activating certain neural units within the R-ANN given the switching rule . Therefore, we generically define the proposed map as follows: (8) (9) = (I; A; ; ); (10) where I is the identity matrix mapping (identically) the initial conditions of the NDA to the RANN and A is the adjacency matrix specifying the network architecture and weights, which will be explained in subsequent sections. In addition, defines different neural dynamics for each type of the neural units, that is, = ( 1; 2; 3) corresponding to MCL, BSL and LTL, respectively (see below for the definitions of these acronyms). The details of the R-ANN architecture and its dynamics are subsequently discussed. 3.1

Network architecture and neural dynamics

The proposed map, , attempts to mirror the affine-linear dynamics (given by Equation 9) of an NDA on the partitioned unit square (see Equation 8) by endowing the R-ANN with a structure capturing the characteristic features of a piecewise-affine linear system, i.e. a state, a switching rule and a set of transformations. External Input

Machine Configuration

Layer

Branch Selection Layer

Linear Transformation

Layer (14) To achieve this, we propose a network architecture with three layers, namely a Machine Configuration Layer (MCL) encoding the state, a Branch Selection Layer (BSL) implementing the switching rule and a Linear Transformation Layer (LTL), as depicted in Figure 1.

The neural units within the various layers make use of either the Heaviside (H) or the Ramp (R) activation functions defined as follows:

H(x) = 0 if x < 0 1 if x 0 (11)

R(x) = 0 if x < 0 x if x 0 Since is a two-dimensional map, this suggests only two neural units (cx; cy) in the MCL layer encoding its state at every step. A set of BSL units functionally acts as a switching system that determines in which cell Di;j the current Turing machine configuration belongs to and then triggers the specific LTL unit emulating the application of an affine-linear transformation i;j on the current state of the system. The result of the transformation is then fed back to the MCL for the next iteration. On the symbolic level, one iteration of the emulated NDA corresponds to a tape and state update of the underlying Turing machine, which can be read out by decoding the activation of the MCL neurons. 3.1.1

Machine Configuration Layer

The role of the MCL is to store the current Go¨delized configuration of the simulated Turing Machine at each computation step, and to synaptically transmit it to the BSL and LTL layers. The layer comprises two neural units (cx and cy), as needed to store the Go¨delized dotted sequence representing a Turing Machine configuration (see Equation 7).

The R-ANNs is thus initialized by activating this layer, given the NDA initial conditions ( x( 0); x( )) which are identically transformed via I by the map ( ) as follows: (cx; cy) = ( x( 0); x( )) 1 = (I; ; ; )j( x( 0); x( )) (13) At each iteration, the units in this layer receive input from the LTL units, and are activated via the ramp activation function (Equation 12); in other words 1 (cx; cy) = (R(Pi tix); R(Pj tjy)). Finally, the MCL synaptically projects onto the BSL and LTL (refer to Figure 2 for details of the connectivity). 3.1.2

Branch Selection Layer

The BSL embodies the switching rule (x; y) and coordinates the dynamic switching between LTL units. In particular, if at the current step the MCL activation is (cx; cy) 2 Di;j = Ii Jj , with Ii = [ i; i+1) being the i-th interval on the x-axis and Jj = [ j ; j+1) being the j-th interval on the y-axis, the BSL units activate only the (tix;j ; tiy;j ) units in the LTL. In this way, only one couple of LTL units is active at each step. The switching rule is mapped by ( ) as follows: 2(x; y) = ( ; ; ; (x; y) = (i; j)): The BSL is composed of two groups of Heaviside (Equation 11) units, implementing respectively the x and the y component of the switching rule of the underlying piecewise affine-linear system, namely: i) the bx group receives input with weight 1 from the cx unit of the MCL layer, and comprises nqns units (i.e. bix; 1 i nqns); ii) the by group receives input with weight 1 from cy and comprises ns units (i.e. bjy; 1 i ns). The activation of the two groups of units is defined as: bix = H(cx bjy = H(cy i) j ) i = min(Ii); j = min(Ji): with with MCL BSL

Each bix and bjy BSL unit has an activation threshold, defined as the left boundary of the Ii and Jj intervals, respectively, and implemented as input from an always-active bias unit (with weight i for the bix unit and j for bjy). Therefore, an activation of (cx; cy) in the MCL corresponding to a point on the unit square belonging to cell Di;j , would trigger active all units bkx with k i. The same would occur for all neural units byk with k j.1 Each bix unit establishes synaptic excitatory connections (with weight h2 ) to all LTL units corresponding to cells Dk;i (i.e. (tkx;i; tyk;i)) and inhibitory connections (with weight 2h ) to all LTL units corresponding to cells Dk;i 1 (i.e. (tkx;i 1; tyk;i 1)), with k = 1; : : : ; ns; for a graphical representation see Figure 2. Similarly, each bjy unit establishes synaptic excitatory connections to all LTL units corresponding to cells Dj;k and inhibitory connections to all LTL units corresponding to cells Dj 1;k, with k = 1; : : : ; nqns. Together, the bix and bjy units completely counterbalance through their synaptic excitatory connections the natural inhibition (of bias h, which value and definition will be discussed in the following section) of the LTL units corresponding to cell Di;j (i.e. (tix;j ; tiy;j )). In other words each couple of LTL units (tix;i; tiy;j ) receives an input of Bxi + Byj , defined as follows: where the input sum

i Bx Byj = b = b x 2 i h + bix+1 j h + bjy+1 y 2 h 2 h 2 ; ; 8h Bxi + Byj = < h2 :0 if (cx; cy) 2 Di;j if cx 2 Ii; cy 62 Jj if (cx; cy) 62 Di;j or cx 62 Ii; cy 2 Jj only triggers the relevant LTL unit if it reaches the value h. That is, if (cx; cy) 2 Di;j then Bxi+Byj = h, and the pair (tix;i; tiy;j ) is selected by the BSL units. Otherwise (tix;i; tiy;j ) stays inactive as Bxi + Byj is either equal to h2 or 0, which is not enough to win the LTL pair natural inhibition. An example of this mechanism is shown in Figure 2 , where the LTL units in cell D1;2 are activated via mediation of bx = fb1x; b2x; b3xg and by = fb1y; b2yg. Here, both b3x and b2y are not excited since cx and cy, respectively, are not activated enough to drive them towards their threshold. However, b2x excites (with weights h2 ) the LTL units in cell D2;2 and D1;2 and inhibits (with weights 2h ) the LTL units in cell D2;1 and D1;1. Equally, b2y excites (with weights h2 ) the LTL units in cell D2;1, D2;2 and 1Note that the action of the BSL could be equivalently implemented by interval indicator functions represented as linear combinations of Heaviside functions. units excite cells fD2;1; D1;1g and fD1;1; D1;2; D1;3g, respectively, but these do not inhibit any cells (due to boundary conditions). 3.1.3

Linear Transformation Layer

The LTL layer can be functionally divided in sets of two units, where each couple applies two decoupled affine-linear transformations corresponding to one of the branches of the simulated NDA. On the symbolic level, this endows the LTL with the ability to generate an updated machine configuration from the previous one. In the LTL, a branch (i; j) of a NDA, i;j (x; y) = ( ix;j x + aix;j ; iy;j y + aiy;j ), is simulated by the LTL units (tix;j ; tiy;j ). Mathematically, this induces the following mapping: (tix;j ; tiy;j ) =

3i;j (x; y) = ( ; ; i;j (x; y); ): The affine-linear transformation is implemented synaptically, and it is only triggered when the BSL units provide enough excitation to enable (tix;j ; tiy;j ) to cross their threshold value and execute the operation. The read-out of this process corresponds to: A strong inhibition bias h (implemented as a synaptic projection from a bias unit) plays a key role in rendering the LTL units inactive in absence of sufficient excitation. The bias value is defined as follows ti;j x ti;j y h 2 = =

R( ix;j cx + ai;j

x R( iy;j cy + ai;j y h + Bxi + Byj ); h + Bxi + Byj ): max(aik;j + i;j;k ik;j ) with k = fx; yg: Hence, each of the BSL inputs Bxi and Byi contributes respectively to half of the necessary excitation ( h2 ) needed to counterbalance the LTL’s natural inhibition (refer to Equation 16 and Equation 17). The LTL units receive input from the two CSL units (cx; cy), with synaptic weights of ( ix;j ; iy;j ), and they are also endowed with an intrinsic constant LTL neural dynamics (aix;j ; aiy;j ). If the input from the BSL layer is enough for these neurons to cross the threshold mediated by the Ramp activation function, the desired affine-linear transformation is applied. The read-out is an updated encoded Turing machine configuration, which is then synaptically fed back to the CSL units (cx; cy), ready for the next iteration (or next Turing machine computation step on the symbolic level). 3.1.4

NDA-simulating first order R-ANN

The NDA simulation (and thus Turing machine simulation) by the R-ANN is achieved by a combination of synaptic and neural computation among the three neural types (MCL, BSL, and LTL) and with a total of nunits =

2 + ns + nsnq + 2ns2nq + |M{CzL} | {z BSL } | L{TzL }

1 b|ia{szu}nit neural units, where nq and ns are the number of states and the number of symbols in the Turing Machine to be simulated, respectively. These units are connected as specified by an adjacency matrix A of size nunits nunits, following the connectivity pattern described in Figure 1 and with synaptic weights as entries from the set

h h f0; 1; 2 ; 2 g [ faik;j

h j i = 1; : : : ; nqns ; j = 1; : : : ; ns ; k = x; yg; the second component being the set of biases.

An important modelling issue to consider is that of the halting conditions for the ANN, i.e. when to consider the computation completed. In the original formulation of the Generalized Shift, there is no explicit definition of halting condition. As our ANN model is based on this formulation, a deliberate choice has to be made in its implementation. Two choices seem to be the most reasonable. The first one involves the presence of an external controller halting the computation when some conditions are met, i.e. an homunculus [ 4 ]. The second one is the implementation of a fixed point condition, (18) (19) (20) (21) intrinsic to the dynamical system, representing a TM halting state as an Identity branch on the NDA. In this way a halting configuration will result in a fixed point on the NDA, and thus on the R-ANN. In other words, the network’s computation is considered completed if and only if 1(x0; y0) = (x0; y0): (22) In the present study we decided to use a fixed point halting condition, but the use of a homunculus would likely be more appropriate in other contexts such as interactive computation [ 9, 10, 11 ] or cognitive modelling, where different kinds of fixed points are required in order to describe sequential decision problems [ 12 ], such as linguistic garden paths [ 4, 10 ].

The implementation of the R-ANN defined like so simulates a NDA in real-time and, thus, it simulates a Turing Machine in real time. More formally, it can be shown that under the map ( ) the commutativity property = is satisfied, which extends the previously demonstrated commutativity property between Turing machines and NDAs [ 9, 13, 14 ]. 4

Discussion

In this study we described a novel approach to the mapping of Turing Machines to first-order RANNs. Interestingly, R-ANNs can be constructed to simulate any piecewise affine-linear system on a rectangular partition of the n-dimensional hypercube by extending the methods discussed The proposed mapping allows the construction, given any Turing Machine, of a R-ANN simulating it in real time. As an example of the parsimony we claim, a Universal Turing Machine can be simulated with a fraction of the units than previous approaches allowed for: the proposed mapping solution derives a R-ANN that can simulate Minsky’s 7-states 4-symbols UTM [ 15 ] in real-time with 259 units (as per Equation 21), approximately 1/3 of the 886 units needed in the solution proposed by Siegelmann and Sontag [ 1 ], and with a much simpler architecture.

In future work we plan to overcome some of the issues posed by the mapping and parts of its underlying theory, especially in relation to learning applications. Key issues to overcome are the missing end-to-end differentiability, and the need for a de-coupling of states and data in the encoding. A future development would see the integration of methods of data access and manipulation akin to that in Google DeepMind’s Neural Turing Machines [ 5 ]. A parallel direction of future work would see the mapping of Turing machines to continuous-time dynamical systems (an example with polynomial systems is provided in [ 16 ]). In particular, heteroclinic dynamics [ 12, 13, 17, 18 ] – with machine configurations seen as metastable states of a dynamical system – and slow-fast dynamics [ 19, 20 ] are promising new directions of research.

References

[1] H. T. Siegelmann and E. D. Sontag , “ On the computational power of neural nets,” Journal of computer and system sciences , vol. 50 , no. 1 , pp. 132 - 150 , 1995 .

[2] H. T. Siegelmann and E. D. Sontag , “ Turing computability with neural nets , ” Appl. Math. Lett , vol. 4 , no. 6 , pp. 77 - 80 , 1991 .

[3] C. Moore , “ Unpredictability and undecidability in dynamical systems,” Physical Review Letters , vol. 64 , no. 20 , p. 2354 , 1990 .

[4] P. beim Graben , B. Jurish , D. Saddy , and S. Frisch , “ Language processing by dynamical systems ,” International Journal of Bifurcation and Chaos , vol. 14 , no. 02 , pp. 599 - 621 , 2004 .

[5] A. Graves , G. Wayne, and I. Danihelka , “Neural turing machines, ” arXiv preprint arXiv:1410.5401 , 2014 .

[6] A. M. Turing , “ On computable numbers, with an application to the entscheidungsproblem , ” Proc. London Math. Soc , vol. 42 , 1937 .

[7] C. Moore , “ Generalized shifts: unpredictability and undecidability in dynamical systems ,” Nonlinearity, vol. 4 , no. 2 , p. 199 , 1991 .

[8] K. Go ¨del, “U¨ ber formal unentscheidbare sa¨tze der principia mathematica und verwandter systeme i,” Monatshefte fu¨r Mathematik und Physik , vol. 38 , pp. 173 - 198 , 1931 .

[9] P. beim Graben , “ Quantum Representation Theory for Nonlinear Dynamical Automata,” in Advances in Cognitive Neurodynamics ICCN 2007 , pp. 469 - 473 , Springer, 2008 .

[10] P. beim Graben , S. Gerth , and S. Vasishth , “ Towards dynamical system models of languagerelated brain potentials,” Cognitive neurodynamics , vol. 2 , no. 3 , pp. 229 - 255 , 2008 .

[11] P. Wegner , “ Interactive foundations of computing,” Theoretical Computer Science , vol. 192 , pp. 315 - 351 , 1998 .

[12] M. I. Rabinovich , R. Huerta , P. Varona , and V. S. Afraimovich , “ Transient cognitive dynamics, metastability, and decision making,” PLoS Computational Biology , vol. 4 , no. 5 , p. e1000072 , 2008 .

[13] P. beim Graben and R. Potthast , “ Inverse problems in dynamic cognitive modeling,” Chaos: An Interdisciplinary Journal of Nonlinear Science , vol. 19 , no. 1 , p. 015103 , 2009 .

[14] P. beim Graben and R. Potthast , “ Universal neural field computation ,” in Neural Fields , pp. 299 - 318 , Springer, 2014 .

[15] M. Minsky , “ Size and structure of universal turing machines using tag systems ,” in Recursive Function Theory: Proceedings, Symposium in Pure Mathematics , vol. 5 , pp. 229 - 238 , 1962 .

[16] D. S. Grac ¸a, M. L. Campagnolo , and J. Buescu , “ Computability with polynomial differential equations , ” Advances in Applied Mathematics , vol. 40 , no. 3 , pp. 330 - 349 , 2008 .

[17] I. Tsuda , “ Toward an interpretation of dynamic neural activity in terms of chaotic dynamical systems,” Behavioral and Brain Sciences , vol. 24 , pp. 793 - 810 , 2001 .

[18] M. Krupa , “Robust heteroclinic cycles, ” Journal of Nonlinear Science , vol. 7 , no. 2 , pp. 129 - 176 , 1997 .

[19] M. Desroches , M. Krupa , and S. Rodrigues , “ Inflection, canards and excitability threshold in neuronal models , ” Journal of mathematical biology , vol. 67 , no. 4 , pp. 989 - 1017 , 2013 .

[20] M. Desroches , A. Guillamon , R. Prohens , E. Ponce , S. Rodrigues , and A. E. Teruel , “ Canards, folded nodes and mixed-mode oscillations in piecewise-linear slow-fast systems,” SIAM Review, vol . in press, 2015 .