S Binczak

S Jacquir

J M Bilbault

V B Kazantsev

V I Nekorkin

Experimental study of electrical MFHN neurons

We present an electronical analogue circuit modelling a FitzHugh-Nagumo neuron with a modified excitability. To characterize this basic cell, the bifurcation curves between stability with excitation threshold, bistability and oscillations are investigated. An electrical circuit is then proposed to realize an unidirectional coupling between two cells, mimicking an inter-neuron synpatic coupling. In such a master-slave configuration, we show experimentally how the coupling strength control the dynamic of the slave neuron, leading to frequency locking, chaotic behavior and synchronization. These phenomena are then studied by phase map analysis. The architecture of a possible neural network is then described introducing different kind of coupling between neurons.

Introduction

As biological neurons have been well studied over the previous century, there is nowadays a strong interest in realizing electrical networks to reproduce the activity of biological neurons. Accurate models are usually used to describe biological neural network, such as Hodgkin and Huxley (HH) [1], and are subject to intensive studies [2,3] in engineering science. Never-theless, mathematicians often study the 2-component FitzHugh-Nagumo (FHN) model [START_REF] Fitzhugh | Impulse and physiological states in models of nerve membrane[END_REF][START_REF] Nagumo | An active impulse transmission line simulating nerve axon[END_REF], as a simplified version of the HH one. This set of equations leads indeed to analytical studies, enlightening the nonlinear dynamical mechanisms generating complex patterns in the neural information process [START_REF] Scott | Neuroscience: A mathematical primer[END_REF]. Although these differential equations are used to model the nerve membrane [START_REF] Scott | Neuroscience: A mathematical primer[END_REF][START_REF] Koch | Biophysics of computation: Information processing in single neurons[END_REF], we suggest to describe overall activity of the neuron by them. In this reductive approach, we propose an experimental electronic implementation of a neuron based on the FitzHugh-Nagumo (FHN) equation with modified excitability to conceive the electronic MFHN neuron [START_REF] Binczak | Experimental study of bifurcations in modified FitzHugh-Nagumo cell[END_REF]. In the first part, the electronic neuron is presented and characterized by its experimental bifurcation curves, showing that this MFHN circuit leads to complex dynamics of travelling waves [START_REF] Kazantsev | Selective communication and information processing by excitable systems[END_REF][START_REF] Dawson | Towards a global classification of excitable reaction-diffusion systems[END_REF] emerging from saddle homoclinic loop bifurcations. In the second part, we use the MFHN circuit as a basic cell to realize a master-slave configuration. Two cells are coupled in an unidirectional manner, which would correspond to two neurons coupled synaptically. After the presentation of the electronic circuit giving this coupling, we discuss the experimental conditions for which the master dynamics controls the excitability of the slave neuron leading to a shift of bifurcation curves, a variation of interspike frequencies or a phenomenon of intermittency route to chaos. Finally, we present the architecture of a possible experimental network including different electrical realizations of neuronal coupling.

2 Experimental description of one cell

Electrical circuit

The nonlinear circuit, as sketched in Fig. 1, can be described as follows: Part (A) is a parallel association of three different branches, two of them being resistive and commuted by silicium diodes (V d = 0.6 V ) while the third one is a negative resistor obtained with an operational amplifier. Due to diodes' com- muting behaviour, the resulting I-V characteristic is nonlinear and can be modelled by a cubic polynomial function for an appropriate set of parameters so that

I N L = f (U ) = 1 R 0 U - γ 2 U 3 3 , ( 1
)
where U and I N L are respectively the voltage and the corresponding current. The parameters R 0 and γ are obtained by a fitting approximation, e.g. by least mean square's method. As illustrated in Fig. 2, we obtain a good match between experimental results and equation ( 1) by setting R 0 = 1010 Ω and γ = 1.138 V -1 . This nonlinear resistor is in parallel

I N L = f (U ) from eq. (1) with R 0 = 1010 Ω and γ = 1.138 V -1 . Dashed line: Experimental I-V characteristic I 1 + I 2 = h(U ) with R 6 = 2021 Ω, R 7 = 690 Ω, L 1 = 10.2 mH, L 2 = 3.5 mH and E 1 = 0.4 V
with a capacitance and two branches in parallel including inductances, resistances and voltage sources, one of them being commuted by a silicium diode so that setting the conditions

R 6 L 1 = R 7 L 2 , E 2 = -V d,
and using a piecewise linear I-V description for diode D 7 , I 2 = 0 if U < 0. Therefore, using Kirchhoff's laws, the system of equations can be expressed in a normalized way by :

         dV dτ = V - V 3 3 -W dW dτ = ε g(V ) -W -η (2)
Where V = γU and W = γR 0 (I 1 + I 2 ) correspond, in biological terms, to the membrane voltage and the recovery variable;

τ = t R 0 C is a rescaled time, ε = R 0 R 6 C L 1 the recovery parameter and η = γ R 0 R 6 E 1 a bifurcation parameter. g(V ) is a piecewise linear function, g(V ) = αV if V ≤ 0 and g(V ) = βV if V > 0 where α = R 0 R 6 and β = L 1 + L 2 L 2 R 0 R 6 control
the shape and location of the W-nullcline [START_REF] Binczak | Experimental study of bifurcations in modified FitzHugh-Nagumo cell[END_REF].

Note also that the initial condition U ini can be loaded in the neuron via an analogue commutator controlled by a periodic signal V syn .

Experimental bifurcation curves of MFHN circuit

This section presents different dynamics of the MFHN neuron [START_REF] Binczak | Experimental study of bifurcations in modified FitzHugh-Nagumo cell[END_REF]. In the case, α = β, the system corresponds to the standard FitzHugh-Nagumo (FHN) equation where the recovery variable is driven by a linear function of the "membrane voltage". In the general case α = β, equation (2) exhibits much more complex dynamics. The phase portrait is very similar to the one occurring from the modified Morris-Lecar (ML) equations [START_REF] Rinzel | Analysis of neural excitability and oscillations[END_REF][START_REF] Koch | Biophysics of computation: Information processing in single neurons[END_REF] proposed to model barnacle muscle fibres and pyramidal cells. Considering the case α = 0.5 and β = 1.96, experimental bifurcation curves in the parameter plane (η,ε) is presented in Fig. 3. Several results are to be distinguished de- pending on the location of the nullcline points and on initial conditions. In domain (1), the systeme has one stable and unstable fixed point, while one point is a saddle. Then, if a perturbation of the rest state, O 1 , is large enough so that it lies between the points O 2 and O 3 , the system responds with an excitation pulse, as illustrated in Fig. 4. The upper inset shows an excitation pulse while the bottom one shows a perturbation reaching the resting state; the main figure gives the respective experimental phase portrait for the both cases. Domain (2) corresponds to the bistability case characterized by the existence of a stable fixed point and a stable limit cycle that has appeared from a big homoclinic loop bifurcation. Then, the model exhibits oscillations if the perturbation is large enough, as illustrated in Fig. 5. The inset shows a spiking train of pulses corresponding to the limit cycle in the experimental phase portrait. Otherwise, it decays to the rest state. In domain (3), the fixed point looses stability via a subcritical Andronov-Hopf bifurcation and only oscillations occur in the model, which are similar to the spiking train of pulses. Note that in the region (2), another saddle homoclinic loop bifurcation has taken place leading to a small unsta- ble limit cycle near the fixed point. Note also that, contrary to standard FHN, arbitrarily long interspike intervals can be found theoretically, as the two lower equilibrium points are merging. Finally, in domain (4), a single unstable fixed point exists leading to oscillations. In domains ( 2) and ( 3), the frequency of the oscillations is controlled by the parameters η and ε, leading to a wide range of values, as illustrated in Fig. 6.

In the following section, we study the dynamical behaviour of a network composed of two neurons connected unilaterally.

Unidirectional coupling of two cells

The neurons communicate mainly between them through specialized devices called synapses via chemical messages. The chemical synapse transmits the impulse unidirectionally. Therefore, it is interesting to conceive and to realize an electrical circuit in- ter (synaptic strength). Its circuit, as illustrated in the Fig. 8, includes an adder-inverter, an inverter and then a follower. U 1 (resp. U 2) is the voltage capacitor of the cell N 1 (resp. N 2). The value of the resistor R is fixed to 100 kΩ, which is large compared to the other components so that the current going though 2R is negligible, Ra = 10 kΩ, while Rc is a detuning parameter which allows to control the coupling parameter value. Note also that the initial condition can be loaded in the neuron via an analogue commutator controlled by V syn . Although it is usual to study a system with normalized variables, it

c - + - + R R R Ra Ra U 1 U 2 - + R Figure 8: Unidirectional coupling circuit.
is more convenient to describe electrical circuits with experimental variables. Therefore, we will keep both variables (as normalized V 1 and experimental U 1 ) in the following of this paper. Using Kirchoff's laws, the normalized equations corresponding to the coupling between two MFHN neurons in a master-slave configuration can be expressed by:

         dV i dτ = V i - V 3 i 3 -W i + DV 1 δ 2,i dW i dτ = ε i g(V i ) -W i -η i (3)
with i = {1, 2}, D = R 0 Rc , and where δ 2,i is a Kronecker symbol, so that δ 2,1 = 0 and δ 2,2 = 1. Therefore, the two neurons are coupled so that a part of current weighted by D via R, and generated by N 1 is included in N 2 . The two neurons are initially set to voltage U 1 ini and U 2 ini, due to the analogue commutators controlled by voltage V syn . When the initial conditions are loaded, these commutators are switched off while the two neurons are connected via a third commutator controlled by V syn . Note that the time delay between the two neurons has not been taken into account in this circuit, a master-slave configuration rendering it unnecessary.

The master in a resting state

When the voltage is so that V 1 is constant (the cell N 1 is in a resting state), the dynamics of the slave cell is given by eqs. (2) with effective bifurcation parameter

η2 = η 2 -DV * 1 , (4)
where V * 1 < 0 is the voltage at stable resting state of the master. Therefore, it implies a modification of the excitability of the cell N 2 corresponding to a shift in the (η, ε) plane illustrated on Fig. 9. The parameters are taken so that, when D = 0, the master neuron N 1 lies in domain (1), while the slave neuron N 2 is in domain (3) and generates a spiking train of pulses. When the unilateral coupling is increased and reaches a critical value, the neuron N 2 ceases to oscillate and stays in the resting state, meaning that the slave neuron has been moved from domain (3) to domain (1) of Fig. 3. This ability of neuron N 1 to inhibit neuron N 2 corresponds to the shift predicted by equation ( 4): As V * 1 < 0, increasing D implies to increase η2 and therefore the bifurcation curves of Fig. 3 are translated along abscissa, while the value of η 2 defined by the electrical parameters of neuron N 2 has not been changed. This result suggests that, for a defined activity of a slave neuron, the strength of this unilateral coupling should be above a critical value to give to the master neuron the control on the slave neuron. In Fig. 9, experimental values (D, η2 ) 3) of the neuron N 2 with ε 2 = 0.01 and when the master neuron N 1 lies in domain (1) in a resting state so that V * 1 = -1.05. Comparison shows a good match between experimental results (+) and equation (4) (continuous line), validating the unilateral coupling circuit.

The master in a spiking regime

In this section, we present some results when the master is in domain ( 2) or (3) and oscillates. As V 1 is varying in time, we cannot express a simple relationship between the parameters of neuron N 2 and V 1 , as in equation ( 4). Nevertheless, oscillations of neurons N 1 let V 1 be alternatively positive and negative, which implies that the bifurcation curves of neurons are translated along the abscissa in the plane (η, ) in a periodic manner (the position of saddle points of the cell N 2 is moved periodically). Thus, the slave neuron N 2 initially situated in the vicinity of a bifurcation curve may be able to cross sometimes this curve and develop a different dynamical behaviour.

Therefore, the master neuron can control the dy- namic of the slave one depending on either its intrinsic behaviour or on the coupling parameter D.

The influence of the frequency f m of the oscillating master neuron on the dynamical behaviour of the slave neuron is illustrated in Fig. 10, where the normalised eigenfrequencies f s /f m with f s , frequency of the oscillating slave neuron, are plotted in function of f m (given by η 1 as presented in Fig. 6) for two values of D.

Initially, the slave uncoupled neuron is so that it lies in domain ( 2) at the resting state, while the master neuron is oscillating.

In both cases, synchronization occurs for small values of f m , the locking range being wider for stronger D.

When f m is above a critical value depending on the coupling, chaotic behavior emerges, separating frequency locking plateaus. This result suggests that the excitability of the master neuron controls the dynamics of the slave neuron. It is interesting to note that chaotic response is observed in a wide range of frequencies, which in is accordance with numerous biological observations [START_REF] Korn | Is there chaos in the brain ? II. Experimental evidence and related models[END_REF]. It also suggests that the coupling parameter D changes the behaviour of the slave neuron. Note that D is the parameter controlling the inter- locking plateaus are separated by chaotic sequences, represented in dot lines in Fig. 11. An illustration of a chaotic signal is given in Fig. 13 for D = 0.02775. The corresponding probability of normalized interspike slave frequency f s/f m shows that in a chaotic regime, the interspike frequencies are widely distributed, indicating that the response spike timing can be very complex. In particular, the sequence has variable inter-spike periods which, however, should be associated with the initiating master spike sequence. Then, we define the relative spiking phase variable, φ n , characterizing relative time (phase) shift between a master spike and the following response spike [14]. The definition of φ n is illustrated in Fig. 14. We take it as

φ n = (t s n -t m n )/T, (5)
where t s n and t m n are the spike peak times of the slave and master unit, respectively, T is the master oscillation period and n is the number of successive slave spike in the response sequence. Figure 15 [START_REF] Korn | Is there chaos in the brain ? II. Experimental evidence and related models[END_REF]13) are plotted versus coupling coefficient D. It appears that for 1 : 1 frequency locking mode the phase tends to a single fixed point φ * . This value accounts the phase shift in the synchronized spiking sequence shown in Fig. 13 (f). Note that for more complicated frequency locking modes the phase shift is described either by stable fixed points with m < φ * < m + 1, where m is number of "integrated" master spikes to get a response, or by stable periodic orbits of a certain period. Note also that transition to chaos suddenly occurs with decreasing D. In this case the spiking phases "fill" some interval with arbitrary changing values from spike to spike. To characterize the phase sequences we calculate the phase map diagrams in the (φ n , φ n+1 )-plane (Fig. 16). For 1 : 1 locking mode the map has one stable fixed point φ * (for instance, φ * = 0.5 in Fig. 16 (c)). With crossing to the chaotic region the fixed point disappears and the phase starts to oscillate. Note the points are located near some curve in the (φ n , φ n+1 )-plane. Then the phase dynamics can be approximately described by a 1D phase map, φ n → φ n+1 . Note also that the map curve goes very close to the bisector line, φ n+1 = φ n . Then, the dynamics occurs as follows.

The phase shift between master and slave spikes is slowly increasing while evolving near the bisector line.

Then, the phase jumps out of the 1 : 1 response region. After one or a few oscillations it returns back to the φ n+1 ≈ φ n -region (Fig. 16 (b)). Such behavior indicates the intermittency chaotic dynamics of the phase variable. Note that the jumps corresponds to the appearance of "missed" spike in the response sequence relative to the master oscillations. More complex intermittency behavior with increasing number of jumps, indicating different missed spike numbers, are shown in Fig. 16 (a). Thus the response spike sequences obtained in the master-slave neuron model can be associated with spiking phase sequences which can be thought as a kind of neuronal "code". The "code" is tuned by the values of interaction coefficient, D. The relative spiking phase represents biologically relevant characteristics of real neuron dynamics. For example, oscillation phase processing plays a crucial role in object tracking by hippocampus neurons [13].

Architecture of the neural network

The previous study allowed us to validate the experimental setup where two neurons where connected unilaterally: The results have been indeed confirmed by numeral simulations of eq. ( 1) using a fourth order runge-kutta algorithm. Therefore, it is interesting to increase the size of the network to realize complex configuration mimicking small assemblies of neurons.

In this reductive point of view, the neurons can be connected via unilateral coupling, as presented previously or via bidirectional coupling, which would correspond to electrical synapses (gap junctions). The first mentioned coupling is improved to generate either inhibitory or excitatory connection (see Fig. 17) such as the current flowing from neuron n to neuron m can be expressed by I = D(V n -V c) where D can be either positive and negative and where V c is the voltage of the neuron n at its resting state. The bidirectional circuit is simply obtained by detuning resistors. Furthermore, an integration of such a circuit using inductor-like components could be realized, giving this opportunity of very large scale network.

Conclusion

We have presented analog electronic circuit implementing MFHN mathematical model. It reproduces the key features of neuron dynamics including excitability, oscillatory dynamics and bistability. We have experimentally investigated the main bifurcations of the model leading to oscillatory and excitable behavior. We have analysed a simple neuron network architecture with two cells in master-slave configuration, e. g. unidirectional coupling. We have shown that the intervals between successive spikes can be chaotic and depends on the intrinsic behaviour of the master neuron and on the coupling strength. The various dynamics of the slave neuron can be reached by detuning D, which could be useful in the point of view of neuronal plasticity.

Figure 1 :

 1 Figure 1: Diagram of the nonlinear circuit

3 Figure 2 :

 32 Figure 2: Experimental I-V characteristic (+ symbols) of part (A). Continuous line: I N L = f (U ) from eq. (1) with R 0 = 1010 Ω and γ = 1.138 V -1 . Dashed line: Experimental I-V characteristic I 1 + I 2 = h(U ) with R 6 = 2021 Ω, R 7 = 690 Ω, L 1 = 10.2 mH, L 2 = 3.5 mH and E 1 = 0.4 V

Figure 3 :

 3 Figure 3: Experimental bifurcation curves in the dimensionless plane (η, ε).

Figure 4 :

 4 Figure 4: Dimensionless phase plane (V,W) showing stability with excitation threshold corresponding to domain (1) of Fig.3 with η = 0.19 and ε = 0.2 (E 1 = 0.332 V and C = 1 nF ). Insets: Excitation pulse (upper inset); perturbation reaching resting state (lower inset). Parameters of insets: abscissa 20 µs/div; ordinate 225 mV /div.

Figure 5 :

 5 Figure 5: Dimensionless phase plane (V,W) showing bistability corresponding to domain (2) of Fig. 3 with η = 0.226 and ε = 0.2 (E 1 = 0.3981 V and C = 1 nF ). Inset: Spiking train of pulses corresponding to limit cycle in experimental phase portrait. Parameters of inset: abscissa 20 µs/div; ordinate 195 mV /div.

Figure 6 :Figure 7 :

 67 Figure 6: Frequency of an oscillating neuron in function of the parameters η and ε. The dot lines correspond to a jump from the resting state to the oscillating one.

2 DFigure 9 :

 29 Figure 9: Shifted bifurcation curve of the slave neuron N 2 between domains (1) and (3). Parameters: Master neuron N 1 : α 1 = 0.5, β 1 = 1.96, 1 = 0.01, η 1 = 0.108 leading to V 1 = -1.05 (i.e. U 1 = -921mV ). Slave neuron N 2 : α 2 = 0.5, β 2 = 1.96,

Figure 10 :

 10 Figure 10: Normalized eigen interspike frequency f s /f m in function of f m for two values of D. Parameters : ε 1 = 0.1934, ε 2 = 0.1938 and η 2 = 0.2163.

Figure 11 :Figure 12 :

 1112 Figure 11: Normalized eigen interspike slave frequency f s by interspike master frequency f m versus D with 1 = 0.1934, 2 = 0.1938, η 1 = 0.2527 and η 2 = 0.2262. action strength. Then changing D mimics variable "plasticity" of the synaptic interaction. We fix the

Figure 11 showsFigure 13 :

 1113 Figure11shows normalized spike frequency locking diagram. Spike frequency locking is defined as the ratio N : M where M is the number of master spikes in a certain time window and N is the number of successive response spikes of the slave neuron. Figure12 (a-f) illustrates the spike trains corresponding to plateaus (a-f) in Fig.11. With increasing value of D (from a to f) the number of response spikes in-

  illustrates the spiking phase bifurcation diagram obtained experimentally. Here the value of φ n , n = 1, 2, 3, . . .

Figure 14 :Figure 15 :

 1415 Figure 14: Qualitative view of superimposed master and slave time series shown by dashed and solid curves, respectively. The spiking phase, φ n , is defined as a time shift between master and slave oscillations peaks.

Figure 16 :

 16 Figure 16: Experimental (φ n , φ n+1 )-plane. The phase map are shown at the top of the figure. The corresponding temporal signals are given below. Zoom of each signal is presented at the bottom of the figure.

Figure 17 :

 17 Figure 17: Electrical excitatory or inhibitory coupling.

24 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 η ε (1) (2) (3) (4)

	0.15	0.16	0.17	0.18	0.19	0.2	0.21	0.22	0.23	0.

= 0.01

Acknowledgments

This research was supported in part by Russian Foundation for Basic Research (grants 03-02-17135), by grant of President of Russian Federation (MK 4586.2004.2). V.B.K. acknowledges Russian Science Support Foundation for financial support.


S Jacquir

S Binczak

P Athias

J M Bilbault
email: bilbault@u-bourgogne.fr

A theoretical approach of the propagation through geometrical constraints in cardiac tissue

Keywords: Travelling wavefront, Circular stationary wave, Blockade phenomenon, Cardiac tissue

The behaviour of impulse propagation in the presence of non-excitable scars and boundaries is a complex phenomenon and induces pathological consequences in cardiac tissue. In this article, a geometrical configuration is considered so that cardiac waves propagate through a thin strand, which is connected to a large mass of cells.

At this interface, waves can slow down or even be blocked depending on the width of the strand. We present an analytical approach leading to determine the blockade condition, by introducing planar travelling wavefront and circular stationary wave. Eventually, the influence of the tissue geometry is examined on the impulse propagation velocity.

Introduction

Reentry is the major mechanism of life-threatening ventricular arrhythmias associated with myocardial infarction scars [START_REF] Myerburg | Frequency of sudden cardiac death and profiles of risk[END_REF]. In abnormal conditions, electrical cardiac wavefront may be stopped by, for instance, electrically non-excitable scar or functional block, which may induce a reentry circuit in pathological endocardial tissue. Nevertheless, normal geometrical features may also imply variations of the conduction properties, as accessory atrioventricular pathways, consisting of narrow strands of myocytes coursing from atrium to ventricle [START_REF] Saffitz | The pathology of sudden cardiac death in patients with ischemic heart disease-arrhythmology for anatomic pathologists[END_REF]. Therefore, the architecture of the cellular network forming the myocardium [START_REF] Rohr | Form and function: Impulse propagation in designer cultures of cardiomyocytes[END_REF][START_REF] Wang | Action potential propagation in inhomogneous cardiac tissue: Safety factor consideration and ionic mechanism[END_REF] and the geometry of excitation wavefront [START_REF] Fast | Role of wavefront curvature in propagation of cardiac impulse[END_REF], 1995a,b] are important to characterize the impulse propagation. Mathematical models [START_REF] Kleber | Basic mechanisms of cardiac impulse propagation and associated arrhythmias[END_REF][START_REF] Noble | A modification of the Hodgkin-Huxley equation applicable to Purkinge fibre action and pacemaker potentials[END_REF] and nonlinear dynamics are widely used to study the impulse (action potential) propagation in cardiac cells [START_REF] Murray | Mathematical Biology[END_REF][START_REF] Keener | Mathematical Physiology[END_REF][START_REF] Scott | Nonlinear Science[END_REF]. Among them, one of the major model is described by the FitzHugh-Nagumo (FHN) equation [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF][START_REF] Nagumo | An active impulse transmission line simulating nerve axon[END_REF], which allows analytical approaches.

The goal of this paper is to study the propagation condition in a system composed of a thin strand of myocytes connected to a large mass of myocytes. This geometry is modeled using a modified version of the FHN 2D equation.

We show that, at the interface of the two parts of the system, the shape of the waves goes through a geometrical modification, which can lead to a decrease of the velocity up to a blockade phenomenon. Eventually, we determine the optimal width of the strand which minimizes the time delay at the interface.

Analytical determination of the propagation condition

We consider a thin sheet of myocardium such as a two-dimensional model is used to describe the impulse propagation. The medium is assumed to be isotropic and its schematic design is sketched in Fig. 1. In the bounded domain Ω, domain (1) is the corridor of width l 1 , domain (2) is the corridor of width l 2 with l 2 > l 1 and the corresponding equations modelling this system are:

               ∂V ∂t = D∆V -f (V ) -W ∂W ∂t = ε(V -γW ) , ( 1
)
where D is the diffusion parameter, t is the time, ∆ is the continuous Laplacian operator, V is the transmembrane voltage. W is the recovery variable, which indicates the capacity of the medium to revert to its resting state after the propagation of impulses due to potassium current and the nonlinear function f (V ) represents the behaviour of the sodium current (it also includes a leak potassium current [START_REF] Scott | Neuroscience[END_REF]). Usually, f (V ) is a cubic polynomial function, but a common simplification is to approximate this function by a piecewise linear expression [START_REF] Fath | Propagation failure of traveling waves in a discrete bistable medium[END_REF]. Furthermore, characterizing the velocity of the action potential can be performed by determining the velocity of the leading edge of this wave, which is given by the condition W (t) = 0 ∀t. The system

(1) becomes

∂V ∂t = D∆V -[V -H(V -α)], ( 2
)
where α is a threshold between the passive and the active role of the sodium conductance (0 < α < 1/2) and H is the Heaviside step function. This system is completed by the Neumann boundary conditions, so that

∂V ∂n δΩ = 0 , (3)
where ∂ ∂n δΩ denotes the outer normal derivative at the boundary δΩ of the bounded domain Ω. To investigate a propagation at the interface, we initiate a wavefront in domain (1). Due to the width of the strand, we assume this wave to be planar and propagating towards domain (2). The symmetry induced by this choice of configuration leads to reduce the system in domain

(1) to a continuous one-dimensional medium so that it can be described by the following bistable equation,

∂V ∂t = D ∂ 2 V ∂x 2 -[V -H(V -α)] . (4)
Introducing the travelling frame coordinate ξ = x -ct, where c is the front velocity and x is the spatial parameter representing the longitudinal axis of the strand, Eq. ( 4) can be written so that

D ∂ 2 V ∂ξ 2 + c ∂V ∂ξ -[V -H(V -α)] = 0 , ( 5
)
which yields, due to C 0 continuity, the propagating solution in domain (1):

               V (ξ) = αe λ 2,1 ξ , if V < α V (ξ) = (α -1)e λ 1,2 ξ + 1, if V ≥ α , (6)
where

λ 1,2 = - c 2D ± 1 2 ( c D ) 2 + 4 D .
Note that in [START_REF] Zemskov | Stability of travelling fronts in a piecewise-linear reaction-diffusion system[END_REF], the stability of this kind of solution in an infinite medium has been analyzed. Note also that this solution has been observed and confirmed by numerical simulations.

When this wave reaches the interface, a symmetry breaking occurs, depending on the width l 1 (resp. l 2 ) of the domain (1) (resp. domain ( 2)) leading to different scenarios:

• if l 1 ≤ l c l 2
, the incoming wave is pinned at the interface (see Fig ( 8)),

where l c is a critical width to be determined,

• if l 2 > l 1 > l c
, the wave is transformed in either circular (l 1 small) or elliptic one (l 1 large) (see Fig ( 9)),

• if l 1 = l 2 , we observe a planar propagation in domain ( 2), according to Eq. ( 6).

In the blockade situation, numerical simulations show that the shape of the stationary wave is circular. Therefore, it corresponds to the existence of a circular stationary solution in domain (2), precluding the possibility of a propagating front.

In order to determine the critical value l c , we need to express the stationary solution in the case of a circular wave. Because of this circular symmetry, we look for a time independent solution V (r) with r radial coordinate, leading to express Eq. ( 2), so that

D( ∂ 2 V ∂r 2 + 1 r ∂V ∂r ) = V -H(V -α). ( 7
)
Introducing ρ = r 1 D , Eq. ( 7) becomes

( ∂ 2 V ∂ρ 2 + 1 ρ ∂V ∂ρ ) = V -H(V -α) , (8)
leading to two versions of the modified Bessel equation:

• if V < α, Eq. (8) becomes ∂ 2 V ∂ρ 2 + 1 ρ ∂V ∂ρ -V = 0 and the stationary solution is V (ρ) = A 1 I 0 (ρ) + B 1 K 0 (ρ) , ( 9
)
where I 0 and K 0 represent the modified Bessel functions of the first and second kind respectively.

• if V ≥ α, Eq. (8) becomes ∂ 2 V ∂ρ 2 + 1 ρ ∂V ∂ρ -V + 1 = 0 and the stationary solution is V (ρ) = 1 + A 2 I 0 (ρ) + B 2 K 0 (ρ). ( 10
)
We impose the following conditions leading to find the constants A 1 , B 1 , A 2 and B 2 :

• V (r → ∞) = 0 leading to A 1 = 0. • V (r = 0) = V 0 > α leading to B 2 = 0.
Eventually, due to C 0 continuity at r = r a where V (r a ) = α, we get for the stationary wave solution

                 V (r) = αK 0 ( r √ D ) K 0 ( ra √ D ) , if r ≥ r a V (r) = 1 + (α -1) I 0 ( r √ D ) I 0 ( r a √ D ) , if r < r a . (11)
Moreover, r a is determined by:

r a = √ DI -1 0 ( α -1 V 0 -1 ). (12)
In Eq. 11, V 0 is the maximal amplitude of an initial condition corresponding to r = 0 and is expressed by

V 0 = 1 + α -1 I 0 ( 1 D r a ) (13)
The value of V 0 only depends on the nonlinear threshold α and it is not influenced by the diffusion properties of the tissue (see Fig 2).

An illustration of a circular stationary wave is presented in Fig. 3. The stationary solution given by Eq. ( 11) allows the determination of the critical width l c . Indeed, the blockade situation corresponds to the case where an incoming planar propagating wave in domain (1) gives birth to the circular stationary wave in domain (2). In this configuration, and due to the boundary conditions, only the half plane solution of Eq. ( 11) is taken into account.

The numerical simulations (see Figs ( 8) and ( 9)) show that at first approximation, it is necessary that at least a half-disc of radius r c is under excited conditions, i.e. with V > α. Thus, for calculating a minimal area generating a propagation, r a is assumed to be equal to r c .

In order to determine r c , we use an analogy of Eq. ( 7) with the Poisson equation used in the field of electrostatic phenomena. The theorem of the divergence is given by

surf --→ gradV . -→ dS = vol ∆V.dτ ( 14
)
where surf is the closed surface directed towards outside and vol is the interior volume.

In our case, the study is restricted to a thin sheet of cells, i.e a bidimensional medium. Therefore, the theorem ( 14) becomes

lig --→ gradV . --→ dM = surf ∆V.dS, ( 15
)
where lig is the closed line directed towards outside and surf is the interior surface.

We calculate the two terms of Eq. (??). J 1 is the integral corresponding to the closed line and J 2 is the integral of the interior surface.

Referring to Fig. ( 4), J 1 = φ 1 + φ 2 where φ 1 and φ 2 are flows.

Thus, the flow coming from the domain ( 1) is given by

φ 1 = -2r c (V i -1) √ D , ( 16
)
with V i corresponds to the value of the wave to the interface of the two domains given by Eq (13). In addition, the outgoing flow of the half-disc, calculated from Eq. ( 11), is given by

φ 2 = -πα r c √ D K 1 ( r a √ D ) K 0 ( ra √ D ) . ( 17
) Then J 2 = surf ∆V. dS = f (V ) D dS = r c r=0 V -1 D πrdr, as when V > α, f (V ) = V -1, with V = 1 + (α -1) I 0 ( r √ D ) I 0 ( ra √ D ) .
The integral J 2 is thus expressed by

J 2 = r c 0 (α -1) D I 0 ( r √ D ) I 0 ( ra √ D ) π r dr. ( 18
)
Setting z = r √ D and using the fact that I 0 (Z) is very appreciably constant and equal to 1 as long as |z| 1, J 2 can be evaluated and expressed by:

J 2 = πr 2 c 2D (α -1) I 0 ( r c √ D ) . ( 19
)
Thus the critical width l c of the domain (1), corresponding to the blockade of the wave, is obtained implicitly from r c , knowing that l c = 2r c . The value of r c is given by J 1 -J 2 = 0:

( 4 π + r c √ D ) (1 -α) I 0 ( r c √ D ) -2α K 1 ( rc √ D ) K 0 ( r c √ D ) = 0. (20)
A comparison between this theoretical prediction and numerical simulations shows a good match between them as illustrated in Fig. 5 (based on a finite difference scheme using a 4 th order Runge-Kutta on a 100 * 100 adaptative mesh). This good correlation is in fact a justification a posteriori of the assumption that a half disc of ray r c = l c 2 must have its cells in excited state (i.e. V > α) so that the propagation initiated in the corridor can continue in the wide second medium. Our study gives in fact a little severe condition on l c

as one can see it on the figure (5). It would be interesting of to continue this study by considering a propagation going from a broad domain to a narrow domain (for example in one kind of bottleneck).

In addition to the blockade phenomenon confirmed numerically when l 1 ≤ l c , the numerical simulations indicate that the wave propagation in the domain

(2) induced a delay because of the transformation of the planar wave in a curved wavefront. In fact, this delay is necessary so that the cells reach their excitation threshold, such as is observed in experiments [START_REF] Saffitz | The pathology of sudden cardiac death in patients with ischemic heart disease-arrhythmology for anatomic pathologists[END_REF]. In the following part, numerical simulations showing the relationship between the velocity of the wave, the intrinsic parameters of the model and the geometrical characteristics of the cardiac tissue are presented.

Wavefront velocities in the medium

Besides the numerical confirmation of the blockade phenomenon when l 1 ≤ l c , numerical simulations indicate that the emergence of propagating wave in domain (2) induces a time delay due to the transformation of the planar travelling wave to a non-zero curvature wavefront. Actually, this delay is required so that cells reach the threshold excitation, as observed experimentally [START_REF] Saffitz | The pathology of sudden cardiac death in patients with ischemic heart disease-arrhythmology for anatomic pathologists[END_REF].

In this section, we present a numerical study showing the relationship between velocity, intrinsic parameters of the model and geometrical features. The basic morphological features of the transition region between a domain (1) and a domain (2) are depicted on Fig. 6, where l is the width of the domain (1), v c is the wavefront velocity between two positions along the longitudinal axis in the domain (1), v o is the wavefront velocity between two positions along the longitudinal axis in the domain (2) and v t is the wavefront velocity between one position in the domain (1) and a second position along the longitudinal axis in the domain (2). As previously described, we initiate a planar wavefront in the corridor travelling towards the large area. Figure 7 shows the different velocities of the wavefront when α = 0.4 and D = 0.5, in function of l (expressed in cell number). The results indicate that the velocity of the circular wave in domain (2) depends on the width l of the strand, which changes the curvature of the transformed wave, as suggested in [START_REF] Kogan | Excitation wave propagation within narrow pathways: Geometric configurations facilitating unidirectional block and reentry[END_REF]. The numerical results illustrate the fact that there exists a critical value l c of l below which the wavefront fails to propagate through the subdomains interface, as shown in Fig. 8. When l is large enough, velocities are equal and the wavefront is no more influenced by the interface. Between these two cases, the symmetry breaking implies a time delay corresponding to the geometrical transformation of the wave from a planar to a non planar wavefront (see Fig. 9 for instance, where a circular wave appears just after the transition region).

Let t a be the necessary normalized time needed to the wavefront to cross the interface so, that

t a = d t v t - d o v o - d c v c d o v o + d c v c ( 21
)
with d i , the distance between the different parts of the medium for which the velocity v i (i={t,o,c}) is measured. Figure 10 shows that this overall delay can be either positive or negative depending on the value of l. Note that delays may imply reentrant arrhythmias due to a lack of synchronicity between wavefront coming from different pathways. Our results indicate that, besides the case of large corridors, there exists an optimum width l opt around which overall delays are negligible (see Fig. 10). Therefore, there is a possible relationship between length and width of corridors preventing dramatic changes of cardiac activity.

Discussion and conclusions

In this study, an analytical solution of travelling front in the case of the planar propagation has been determined. The criterion for the initiation of the cir-cular waves has been investigated, based on the existence of stationary waves.

We showed that the existence of the stationary wave prevents the propagation.

Again, this breakdown depends not only on the nonlinearity threshold and on the diffusion coefficient but also on the geometrical morphology of the cardiac syncytium and on the form of the excitation wavefront. However, our results concerning the condition of propagation are valid in the case of a very thin strand of cells (domain 1). To generalize and to affine our approach, it would be necessary to use the spheroidal functions using the elliptic coordinates.

There exists specific geometric configuration of two-dimensional narrow pathways, which allow unidirectional conduction blocks and reentry even when membrane properties and intercollated resistance are uniform [START_REF] Kogan | Excitation wave propagation within narrow pathways: Geometric configurations facilitating unidirectional block and reentry[END_REF]. We suggest that this work could be applied to avoid the reentry phenomena. In other words, our results could modify the strategies used during arrhythmia ablation procedures. Indeed, precise radiofrequency applications inside scar tissues could improve conduction velocity in other areas and prevent reentries [START_REF] Chillou | Isthmus characteristics of reentrant ventricular tachycardia after myocardial infarction[END_REF]. It would correspond to fix the optimized width so that the delay becomes negligible, which would assures synchronicity between waves coming from different pathways.

Captions

Fig1: Schematic diagram of the geometry of the tissue. l 1 (resp. l 2 ) is the width of the domain (1) (resp. domain ( 2)). Arrow indicates direction of propagation.

Fig2: Initial condition V 0 versus α for all diffusion parameter D.

Fig3: Cross-section of the standing wave with D = 1, α = 0.2, r a = 0.8186 leading to V 0 = 0.319.

Fig4:

  Fig4: Diagram illustrating flows entering and outgoing of the half disc of ray r c . Fig5: l c versus α (D = 1). Theoretical results (cf Eq.(??)) in continuous line (crosses), numerical results in dashed line (stars) with error bars ±0.05. Fig6: Schematic diagram of the geometry of the tissue using during simulation.l is the width of the thin strand and d i is the distance between the different parts of the medium for which the velocity v i (i={t,o,c}) measured.

Fig7:

  Fig7: Wavefront normalized velocity versus the width of the strand expressed in cell numbers with parameters D = 0.5 and α = 0.4. Fig8: Simulation of the blockade phenomenon with parameters l < l c , D = 1 and α = 0.2. Fig9: Simulation of the propagation phenomenon with parameters l > l c , D = 1 and α = 0.2. Fig10: Time t a versus the width of the strand expressed in cell numbers with parameters D = 0.5 and α = 0.4.

Fig

  Fig. 4.


Sabir Jacquir
email: sjacquir@u-bourgogne.fr

Stephane Binczak

Jean-Marie Bilbault

NDES2005 Submission Style ANALYTICAL DETERMINATION OF INITIAL CONDITIONS LEADING TO FIRING IN NERVE FIBERS

Keywords: Stationarity, reaction-diffusion, nerve fibers

An analytical solution characterizing initial conditions leading to action potential firing in smooth nerve fibers is determined, using the bistable equation. In the first place, we present a non-trivial stationary solution wave, then, using the perturbative method, we analyze the stability of this stationary wave. We show that it corresponds to a frontier between the initiation of the travelling waves and a decay to the resting state. Eventually, this analytical approach is extended to FitzHugh-Nagumo model.

Introduction

In 1952, Hodgkin and Huxley (HH) [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]] described the behaviour of electrical impulses in the axon of the squid by a set of partial differential equations. Their model produces a realistic description of the action potential propagation along a smooth axon. There are other models which arise from the simplification of the HH model, such as the FitzHugh-Nagumo (FHN) model [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF]; [START_REF] Nagumo | An active impulse transmission line simulating nerve axon[END_REF]]. This one produces the qualitative dynamics of the nerve membrane and it is largely used to investigate the mechanism of the propagation or its failure in a one-dimensional excitable medium [START_REF] Neu | Initiation of propagation in a one-dimensional excitable medium[END_REF]; [START_REF] Lewis | Wave-block in excitable media due to regions of depressed excitability[END_REF]], in inhomogeneous media [START_REF] Prat | Stability of front solutions in inhomogeneous media[END_REF]], and in discrete reaction-diffusion systems [START_REF] Binczak | Experimental propagation failure in a nonlinear elctrical lattice[END_REF]; [START_REF] Keener | The effects of discrete gap junctions coupling on propagation in myocardium[END_REF]; [START_REF] Bell | Excitability behavior of myelinated axon models, in reaction-diffusion equations[END_REF]], for instance. In a piecewise form, it has also been used to study threshold phenomena for the reaction-diffusion systems [START_REF] Terman | Threshold Phenomena for a Reaction-Diffusion System[END_REF]]. When focussing on the leading edge of impulses, this FHN model reduces to the bistable equation. The aim of this paper is to determine the conditions allowing the excitation of nerve fibers (initiation of the travelling pulse) or its setting in latency. The paper is organized as follows: In section 2, the analytical shape of the stationary wave which corresponds to initial conditions, is calculated and its stability is studied for the bistable equation. The generalization of the stationary solution to FHN model is presented in section 3. Finally, the section 4 is devoted to concluding remarks.

Stationarity in the Bistable Model

In a first place, let us consider the bistable equation in an infinite continuous onedimensional medium such as:

∂V ∂t = D ∂ 2 V ∂x 2 -f (V ) . ( 1
)
In this model of the smooth nerve fiber, V denotes the voltage accross the membrane, D represents the diffusion coefficient, t is the time and x is the longitudinal axis of the fiber. The nonlinear function f (V ) = V (V -α)(V -1) represents the behaviour of the sodium current with α acting as a threshold between the passive and the active role of the sodium conductance (0 < α < 1), so that the steady state V = 0 is the resting state and the steady state V = 1 is the excited one. The system is completed by Neumann conditions, so that lim x→±∞ dV dx = 0 . In this context, we look for non uniform stationary waves so that

lim x→±∞ V (x) = V ∞ and V (x = 0) = V ∞
, where V ∞ is the steady state 0 or 1. Due to symmetry, we can reduce this problem to a semi-infinite medium by setting dV dx (x = 0) = 0 . Depending on the choice of V ∞ , several results are to be distinguished to find a stationary wave.

Case where

V ∞ = 0
We look for a time-independent solution V = V s (x) of eq. (1):

D d 2 V dx 2 = V (V -α)(V -1) . ( 2
)
Multiplying both sides by dV dx and taking into account the Neumann conditions requirements, eq. ( 2) yields, after integration, to

dV dx = 2 D ( V 4 4 -(α + 1) V 3 3 + α V 2 2 ) = 2 D F(V ) . ( 3
)
Where F(V ) is the antiderivative function of f (V ). In the case where α ≤ 1/2, F(V ) has three real roots and can be expressed by

F(V ) = 1 4 V 2 (V -V 1 )(V -V 2 ) with 0 < V 1 < V 2 , V 1,2 = 2 3 (α + 1) ∓ 1 3 √ 4α 2 -10α + 4, ( 4
) so that V 1 V 2 = 2α.
As illustrated in Fig. 1, V = 0 and V = 1 are the minima of F(V ) corresponding to stable states and V = α is a local maximum corresponding to an unstable state. To integrate eq. ( 3), F(V ) must be positive, therefore excursion of V must lie in the range [0,V 1 ]. Thus, eq. ( 3) leads to

V s Z V 1 dV V (V -V 1 )(V -V 2 ) = ± 1 2D x Z 0 dx . ( 5
)
which yields to Eq. ( 6) corresponds to the stationary wave whose width depends on the diffusion parameter D and on the threshold parameter α. The amplitude of V s is maximum when x = 0 so that V s (0) = V 1 and is only a function of α. Because of symmetry, this result can be extended to an infinite medium, as illustrated in Fig. 2. When α 1, eq. ( 6) becomes

V s (x) = V 2 1 + ( V 2 V 1 -1) cosh 2 ( x 2 α D ) (6)
V s (x) = 3 2 α sech 2 ( x 2 α D ) , (7)
which corresponds, when D = 1, to the results obtained in [START_REF] Neu | Initiation of propagation in a one-dimensional excitable medium[END_REF]]. The existence of the stationary wave given by Eq. ( 6) (in accordance with the results obtained in [START_REF] Bountis | Stationary Pulses and Waves Front Formation in an Excitable Medium[END_REF]]) is assured only when α ∈ [0, 1/2]. As we shall see, the condition α ∈ [1/2, 1] implies another solution.

Case where V ∞ = 1

In this case, eq. ( 2) yields to

dV dx = 2 D F(V ) . ( 8
)
The antiderivative function F(V ) has three real zeros only if α ≥ 1/2. In this case, it can be expressed by

F(V ) = 1 4 (V -1) 2 (V -V 1 )(V -V 2 ) where V 1 < V 2 < 1, V 1,2 = 2α 3 -1 3 ∓ 1 3 √ 4α 2 + 2α -2 (9)
As previously, to integrate eq. ( 8), F(V ) must be positive, therefore excursion of V must lie in the range [V 2 , 1], as illustrated in Fig. 3. Thus, eq. ( 8) leads to

V s Z V 2 dV (V -1) (V -V 1 )(V -V 2 ) = ± 1 2D x Z 0 dx , ( 10
)
which yields to with Γ = D(1 -α). Here again, the width of this stationary wave depends on D and α, while the minimum value of the amplitude V min = V 2 , located at x = 0 is only function of α. An illustration of such a wave is given in Fig. 4, where the medium is extended from a semi-infinite one to an infinite one. Note that the existence of this stationary wave is assured when α ∈ [1/2, 1], complementing the first result.

V s (x) = V 1 sinh 2 ( Γx 2 ) -V 2 cosh 2 ( Γx 2 ) + (1 -2α)/3 V 1 cosh 2 ( Γx 2 ) -V 2 sinh 2 ( Γx 2 ) -1 , ( 11)

Stability of the stationary solution

To study the stability of the stationary wave, we examine the evolution of an infinitesimal disturbance v(x,t) superimposed on the stationary wave V s (x). To determine its evolution equation, we introduce a scaling parameter Ψ 1 so that V (x,t) can be expressed by

V (x,t) = V s (x) + Ψv(x,t) . ( 12
)
Equation ( 12) is inserted in eq. ( 1) leading to express the evolution of V s (x) and v(x,t) separately. Using this perturbative method to the power of Ψ inferior to 2, eq. ( 1) is developped into two equations depending on the power of Ψ:

     Ψ 0 : D d 2 V s dx 2 = V s(V s -α)(V s -1) Ψ 1 : ∂v ∂t = D ∂ 2 v ∂x 2 -h(x)v , ( 13
) with h(x) = 3V s (x) 2 -2(α + 1)V s (x) + α.
The first equation (Ψ 0 ) corresponds to the stationary case, while the second one (Ψ 1 ) gives the evolution of the perturbation v(x,t). Let us focus on a dirac perturbation located at x = 0, as v(x,t = 0) = v 0 δ(x) . In this case, one can approximate h(x) so that h(x) = h(0) = h 0 becomes uniform in the region of x = 0. Considering this approximation and using the Fourier transform and the Inverse Fourier transform, the solution of eq. ( 13) is given by:

v(x,t) = v 0 e -h 0 t √ 4πDt e -x 2 4Dt . ( 14
)
If h 0 > 0 then v(x,t) → 0 therefore V (x,t) = V s (x,t), which indicates that the stationary wave is stable. If h 0 < 0 then v(x,t) diverges, which is the case around x = 0 ∀α. This result indicates that this stationary wave is unstable and acts as a frontier between two different evolutions: From numerical simulations and according to [START_REF] Neu | Initiation of propagation in a one-dimensional excitable medium[END_REF]], if the initial signal is beneath V s (x), it will converge to the resting state V (x,t) = 0, otherwise it generates a travelling wave, as illustrated in Fig 5, where, in the top (resp. bottom) part, a negative (resp. positive) perturbation is superimposed to the stationary solution, leading to a decay to zero (resp. a travelling wave). In the following section, the precedent results are generalized to FitzHugh-Nagumo model.

Generalization to the FHN Model

The FHN model includes a recovery variable variable W indicating the capacity of the medium to revert to its resting state [START_REF] Scott | Nonlinear Science[END_REF]] , and is expressed by

∂V ∂t = D ∂ 2 V ∂x 2 -f (V ) -W ∂W ∂t = ε(V -γW ) . ( 15
)
The stationary solution of this equation must satisfy

D ∂ 2 V ∂x 2 -f (V ) - V γ = 0 and W = V /γ . ( 16
)
Applying the previous methodology, the antiderivative function eq. ( 16) leads to

dV dx = 2 D F(V ) , ( 17
)
where

F(V ) = 1 4 V 2 (V -V 1 )(V -V 2 ) , V 1,2 = 1 3γ (2αγ + 2γ ∓ 4α 2 γ 2 -10αγ 2 + 4γ 2 -18γ)
A nontrivial stationary wave only exists if V 1 and V 2 exist and are real, i.e. if α < 1 2 and γ > γ c > 0 with Under these conditions, integrating eq. ( 14) yields to express the stationary wave, so that

γ c = 9 (1 -2α)(2 -α) > 0 . ( 18
     V s = V 2 1 + ( V 2 V 1 -1) cosh 2 ( x 2 V 1 V 2 2D ) W s = V s γ -1 (19)
An illustration of such a wave is presented in Fig. 6, while Fig. 7 shows the relationship between α and γ leading to the existence of stationary waves, according to eq. ( 18). Numerical simulations have shown that a positive (resp. negative) perturbation superimposed to V s in the region of x = 0 leads to the emergence of a travelling impulse (resp. a decay to zero), while a perturbation superimposed to W s implies an opposite evolution, which indicates that this stationary wave corresponds to minimum spatial initial conditions from which a travelling wave can emerge.

Conclusion

In this study, the criterion for the initiation of the travelling waves has been investigated, based on the existence of stationary waves. We suggest that these results could be applied to determine the minimum punctual electrode stimulation needed to generate an action potential. Referring to [START_REF] Neu | Initiation of propagation in a one-dimensional excitable medium[END_REF]], these stationary waves can be seen as some critical nucleus, which can be useful to determine the relationship between amplitude and width of conditions initiating travelling waves and to optimize the size of electrodes or stimulation areas.

February

  Fig. 1. F(V) in the case where α = 0.3.

Fig. 4 .

 4 Fig. 3. F(V) in the case where α = 0.7.

Fig. 5 .

 5 Fig. 5. Evolution of the perturbed stationary wave in the case where α = 0.2 and D = 1.

Fig. 6 .Fig

 6 Fig. 6. Stationary wave of the FHN model (α = 0.2 and γ = 10).