Kiril S Shterev
email: kshterev@imbm.bas.bg

Stefan K Stefanov
email: stefanov@imbm.bas.bg

A two-dimensional computational study of gas flow regimes past of square cylinder confined in a long microchannel

Keywords: transition, steady flow regime, unsteady flow regime, microflows, continuum model, SIMPLE-TS, molecular model, DSMC, SBT PACS: 47.15.Fe, 47.27.nd, 47.45.-n, 87.10.Rt 2000 MSC: 76F06, 76P05, 76M12, 82B80

In this paper, the transition between steady and unsteady regimes of the two-dimensional gas flow past a moving square cylinder confined between two parallel plates with equal temperatures was investigated for different subsonic and transonic speeds and blockage ratios B=3, 5, 10 and 20 (B is the ratio between microchannel height and square size). The influence of B on the drag coefficient of the square, pressure at stagnation point and transition between steady and unsteady regime were studied. A set of transition curves were obtained for different blockage ratios. The results showed a similar shape of the neutral curves, which shifts to lower Knudsen and Mach numbers when B decrease. A continuum approach based on the Navier-Stokes-Fourier equations was basically applied to all flow regimes. Navier-Stokes-Fourier equations was calculated numerically by using pressure based finite volume algorithm SIMPLE-TS. For the initial transition velocity slip regime, the results were compared to data obtained by using the particle Direct Simulation Monte Carlo (DSMC) method.

Introduction

Micromechanical devices are associated with rapidly emerging technologies, where new potential applications are continuously being developed, and fluid flow regime is important for their design [START_REF] Colin | Rarefaction and compressibility effects on steady and transient gas flows in microchannels[END_REF], [START_REF] Nabavi | Steady and unsteady flow analysis in microdiffusers and micropumps: a critical review[END_REF]. For example, a gaseous flow past a small particle confined in microchannel is a typical microfluidic configuration, which is considered here. The flow regime could different significantly under different conditions. Thus, in this paper a variety of established two-dimensional steady and unsteady regimes of flow past a cylinder with square cross-section confined in a microchannel is investigated as a function of gas rarefaction (Knudsen number), velocity of the moving obstacle (Mach number) and blockage ratio of the channel (ratio between microchannel height and square size). To simplify the problem, we consider the equivalent inverse problem, in which the obstacle is fixed in a local coordinate system and gas and channel walls are moving in the opposite direction with given velocity. Usually, a microchannel gas flow is modelled by using either continuum approach with slip boundary conditions or by using rarefied gas dynamics methods. The Knudsen number (Kn), defined as a ratio of the mean free path * ∞ to the macroscopic length scale of a physical system L * ,

Kn = * ∞ L * (1)
is a nondimensional parameter that determines the degree of gas rarefaction, and respectively, the degree of validity of the continuum model. For larger Knudsen numbers, where the non-equilibrium effects are significant, a molecular approach, based on kinetic theory models [START_REF] Cercignani | Rarefied Gas Dynamics. From Basic Concept to Actual Calculations[END_REF], [START_REF] Sone | Molecular Gas Dynamics: Theory, Techniques, and Applications[END_REF], or the particle DSMC (Direct Simulation Monte Carlo) method [START_REF] Bird | Molecular, Gas Dynamics and the Direct Simulation of Gas Flows[END_REF] has to be applied. In the present paper, we use both Navier-Stokes-Fourier continuum and molecular (DSMC) models to calculate the flow in the early transition regimes. The Navier-Stokes-Fourier continuum model [START_REF] Stefanov | Rayleigh-Bénard flow of a rarefied gas and its attractors. I. Convection regime[END_REF] with included statedepended transport coefficients determined by the first approximation of the Chapman-Enskog theory is used at low Knudsen numbers.

The peculiarity of fluid flows in and around micro-electro-mechanicalsystems (MEMS) could be related in most cases to low speed effects. The low Mach number typical for such flows makes pressure based methods very suitable to be used here. The pressure based algorithm SIMPLE-TS [START_REF] Shterev | Pressure based finite volume method for calculation of compressible viscous gas flows[END_REF] applied in the current study covers the entire slip regime.

The emergency of instability and prediction of a transition between different flow regimes is an important task in the study of gaseous flows around bodies with arbitrary shape. The laminar-turbulent properties of the flow field have a significant influence on skin friction and flow separation and, respectively on the lift and drag characteristics. Different transition prediction models have been developed to simplify both steady and unsteady flow regime predictions, see [START_REF] Patel | Turbulence models for near-wall and low Reynolds number flows-a review[END_REF], [START_REF] Johansen | Prediction of laminar/turbulent transition in airfoil flows[END_REF].

In particular, investigations of the final state of incompressible fluid past a square with non-slip boundary conditions are widely reported [START_REF] Sohankar | Numerical Simulation of Unsteady Flow Around a Square Two-Dimensional Cylinder[END_REF]. By contrast, there are a few investigations of the problem in slip and transition regimes. We refer to Cheng and Liao [START_REF] Cheng | DSMC Analysis of Rarefied Gas Flow Over a Rectangular Cylinder at All Knudsen Numbers[END_REF], who have analyzed a gaseous flow past a square in the transition regime ( moderate and high Knudsen numbers) for Mach number ranged from 0.85 to 8. Here we investigate the final state of compressible viscous flow with slip boundary conditions as a function of Knudsen and Mach numbers, and blockage ratio in the case of a two-dimensional flow past a square cylinder in a long microchannel. The investigated Knudsen numbers are in the range from 0.002 to 0.02. The free stream Mach numbers (M) are in the range from 0.1 to 0.6, while the local Mach number varies from 0.1 to ∼ 1.5. In order to be able to calculate all problems with a unified algorithm, the iterative pressure based finite volume method (FVM) with improved convergence with respect to time step (SIMPLE-TS) was used [START_REF] Shterev | Pressure based finite volume method for calculation of compressible viscous gas flows[END_REF]. A GPU (Graphical Processor Unit) version of algorithm SIMPLE-TS [START_REF] Shterev | GPU Calculations of Unsteady Viscous Compressible and Heat Conductive Gas Flow at Supersonic Speed[END_REF] with Initial Guess Extrapolation (IGE) [START_REF] Shterev | Iterative process acceleration of calculation of unsteady, viscous, compressible, and heat-conductive gas flows[END_REF] was used to speed-up the calculations.

A number of simulations by using the DSMC method has been performed to validate the SIMPLE-TS for moderate Knudsen numbers. The same DSMC algorithm has been applied previously in series of papers, [START_REF] Stefanov | Rayleigh-Bénard flow of a rarefied gas and its attractors. I. Convection regime[END_REF], [START_REF] Stefanov | Rayleigh-Bénard flow of a rarefied gas and its attractors. II. Chaotic and periodic convective regimes[END_REF], [START_REF] Stefanov | Rayleigh-Bénard flow of a rarefied gas and its attractors. III. Three-dimensional computer simulations[END_REF], for studying the flow instability of 2D and 3D Rayleigh-Bénard flow of a rarefied gas as a function of Knudsen and Froude numbers.

Continuum model equations

The time evolution of the continuum model is described on the basis of the Navier-Stokes-Fourier (NSF) equations [START_REF] Stefanov | Rayleigh-Bénard flow of a rarefied gas and its attractors. I. Convection regime[END_REF] for a compressible viscous gas with transport coefficients determined by the first approximation of the Chapman-Enskog theory for low Knudsen numbers.

The NSF equation system consists of the equations of continuity, momentums, energy and state for perfect gas. We introduce the nondimensional variables

u = u * V * th , v = v * V * th , p = p * p * ∞ , T = T * T * w , t = L * V * th , x = x * L * , y = y * L * , ( 2
)
where u is nondimensional horizontal velocity, v is nondimensional vertical velocity, p is nondimensional pressure, T is nondimensional temperature, t is nondimensional time, x and y are nondimensional coordinates of a Cartesian coordinate system, an asterisk denotes dimensional variables, infinity represent free stream conditions, L * is the dimensional reference length equal to the square size (L * = a * , see Figure 1), V * th = 2RT * w is the most probable molecular speed, R is the gas constant, T * w is dimensional temperature of channel and square walls, see Figure 1.

In our considerations, a hard-sphere gas description with viscosity coefficient µ and the heat conduction coefficient λ (first approximations are sufficient) reads:

µ = µ h √ T , µ h = (5/16)ρ * ∞ l * ∞ V * th √ π (3)
λ = λ h √ T , λ h = (15/32)c p ρ * ∞ l * ∞ V * th √ π (4)
where: ρ * ∞ is dimensional free stream density. Other used variables are: nondimensional density is ρ, Mach number of V * th is Ma = V * th / γRT * w = 6/5, the free stream Mach number is M = u * sq / γRT * w , dimensional square velocity is u * sq , the Prandtl number is P r = 2/3, C µ = (5/16) 2π/γ and γ = c p /c v = 5/3 (c p and c v are the heat capacities at constant pressure and constant volume respectively).

Then, the system of unsteady nondimensional continuum equations can be expressed in a conservative form as follows:

∂ρ ∂t + ∂(ρu) ∂x + ∂(ρv) ∂y = 0 (5) Du Dt = - 1 γρMa 2 ∂p ∂x + C µ 2Kn ρMa ∂ ∂x √ T ∂u ∂x + C µ Kn ρMa ∂ ∂y √ T ∂u ∂y + ∂v ∂x - 2 3
∂ ∂x The final state of the flow is determined for the case of a square cylinder with size a confined in a plane microchannel (height H ch ) as shown in Figure 1. The considered blockage ratios (B = H ch /a) are B = 3, 5, 10 and 20. The distance from square to the channel inlet flow is L a (inlet length).

√ T ∂u ∂x + ∂v ∂y (6)
The considered equation system represents two equivalent problem descriptions with respect to the coordinate system (depending on the observer viewpoint). In the first case, the Cartesian coordinate system (0xy) is fixed to the stationary channel. The square cylinder is moving with constant velocity -u sq along x-axes and the channel wall velocity is zero. The variables related to the coordinate system 0xy are denoted with overline. In the second case, a local Cartesian coordinate system (0xy) is fixed to the moving square cylinder. In the local coordinate system the velocity of the square cylinder is zero while the channel walls are moving with constant velocity u sq . In order to use a computational grid fixed to the moving square cylinder, this case is considered in the present paper. The relations between variables in both coordinate systems are given by the Galilean translation:

t = t, x = x -tu sq + constant, y = y, u = u -u sq , v = v, p = p, T = T (11)
In the local coordinate system, the equations ( 5) -( 9) and the boundary conditions at the inlet and outlet (see Figure 1) remain the same as given for the first coordinate system.

Velocity-slip and temperature-jump boundary conditions are imposed on the walls of the channel and the square. The transition between steady and unsteady final state occurs at Knudsen numbers lower than 0.02. Therefore, we impose first order velocity-slip and temperature-jump BCs as follows:

v s -v w =A σ Kn local ∂v ∂n s , T s -T w =ζ T Kn local ∂T ∂n s , (12)
where v s is velocity of the gas at the solid wall surface, v w is velocity of the wall, Kn local = Kn/ρ local is the local Knudsen number, ρ local is the local density, ∂v ∂n s is the derivative of velocity normal to the wall surface, T s is temperature of the gas at the wall surface, T w is temperature of the wall and ∂T ∂n s is the derivative of temperature normal to the wall surface. We have used the velocity-slip and temperature-jump coefficients A σ = 1.1466 and ζ T = 2.1904, calculated respectively in [START_REF] Albertoni | Numerical evaluation of the slip coefficient[END_REF], [START_REF] Bassanini | Comparison of kinetic theory analyses of linearized heat transfer between parallel plates[END_REF], from the BGK equation by using variational method (for details see [START_REF] Cercignani | The Boltzmann Equation and its Applications[END_REF]).

Computational considerations

In our considerations, we are interested mainly in the final flow state established after passing the transient period of flow evolution. Two final regimes that occur in the flow are distinguished and the transition between them is analyzed: Steady state regime: The flow regime after transient period of calculation is steady and macroscopic variable fields do not change in the course of time. Unsteady state regime: The flow regime after transient period of calculation is unsteady and macroscopic variable fields varies periodically in time (a non-periodical behavior is also possible but not observed in this investigation). A transition (neutral) curve delineated in the space of parameters M and Kn separates both steady and unsteady state regimes.

The final states of both regimes are investigated at Mach numbers from 0.1 to 0.6 for B = 3, 5, 10 and 20. Results from an earlier stage of this study for the case B = 3 at Mach numbers M = 0.1, 0.2, 0.3 and 0.4 were presented in [START_REF] Shterev | Determination of Zone of Flow Instability in a Gas Flow Past a Square Particle in a Narrow Microchannel[END_REF]. Here, a series of numerical calculations with different Knudsen number are performed for each considered Mach number. The simulation starts from a uniform state of the gas with temperature T w , pressure p ∞ and velocity equal to the channel walls. In SIMPLE-TS within each time step the iteration process is stopped, when the maximum residuals of the equations are less than 10 -9 . The final state is considered as unsteady when the latter criteria are not fulfilled in the course of time and non-symmetrical variations along y axis are developed in the macroscopic variable fields. The non-symmetrical oscillations in the investigated cases are time dependent and lead to corresponding variations of the drag coefficient. We consider that the unsteady flow reaches a final state when the drag coefficient exhibits a regular periodic time evolution behavior in any cross-section behind the square, see Figure 2. The drag coefficient C D is defined as:

C D = F * D 1 2 ρ * ∞ u * sq 2 A * , (13)
where F * D is dimensional drag force, u * sq is dimensional square velocity, and A * is dimensional cross section area of the square normal to the free stream flow. An important requirement for accurate calculations is to realize a flow regime when the influence of inflow and outflow boundaries becomes negligible. It is worth noting, that for the continuum approach the implementation of inflow and outflow boundary condition with no fixed pressure keeping mass conservation in computational domain is a task, yet not achieved completely [START_REF] Gadelhak | MEMS: Introduction and Fundamentals[END_REF]. Thus, in our considerations the pressure at the inflow boundary is fixed. The influence of this boundary condition can be neglected, when the disturbances arising around the moving square do not reach the channel inflow area, see Figure 3.

The system of equations ( 5) -( 9) is solved by using the algorithm SIMPLE-TS [START_REF] Shterev | Pressure based finite volume method for calculation of compressible viscous gas flows[END_REF]. Two variants of algorithm SIMPLE-TS were used. First one implies a Total Variation Diminishing (TVD) second order scheme, with SU-PERBEE limiter [START_REF] Roe | Some contributions to the modelling of discontinuous flows[END_REF] that approximates convective terms. The TVD scheme was applied as explicit scheme (Forward Euler) with regard to time. The second one, upwind first order scheme, approximates convective terms as an implicit scheme (Backward Euler) with regard to time. A second order central difference implicit scheme (Backward Euler) approximates the diffusion terms. Subsonic cases were calculated using SIMPLE-TS with second order TVD scheme for approximation of convective terms. A GPU (Graphical Processor Unit) version of algorithm SIMPLE-TS [START_REF] Shterev | GPU Calculations of Unsteady Viscous Compressible and Heat Conductive Gas Flow at Supersonic Speed[END_REF] with Initial Guess Extrapolation (IGE) [START_REF] Shterev | Iterative process acceleration of calculation of unsteady, viscous, compressible, and heat-conductive gas flows[END_REF] was used to speed-up the calculations.

The SIMPLE-TS algorithm uses a staggered-mesh system. Two meshes were used to ensure sufficient accuracy of results and reasonable computational time. For cases B=3 and Mach number M ≥ 0.4 and B=5 and M ≥ 0.55, spatial steps along x-axis ∆x min = 0.05 for x ∈ [-10; 20] (around the square) and ∆x max = 0.4 for x ≤ -30 and x ≥ 40 were used while step along y-axis is constant ∆y = ∆x min = 0.05. For smooth transition between ∆x min and ∆x max a linear polynomial relation was used . The uniform Cartesian mesh around the square, where gradients could be significant, is finer, while in zones with small gradients is coarser, see Figure 3. The staggered-mesh for all other considered cases was defined analogically with spatial steps: ∆x min = 0.1, ∆x max = 0.75 and ∆y = ∆x min = 0.1. The time step is sufficiently small to ensure convergence of numerical algorithm SIMPLE-TS and so that it is in a range from 0.01 to 0.04 for subsonic considered cases.

In the present paper, the particle DSMC method was used for validation of the continual approach in the early transition regime. The Stefanov's Simplified Bernoulli trials (SBT) (see for details [START_REF] Stefanov | On DSMC Calculations of Rarefied Gas Flows with Small Number of Particles in Cells[END_REF]) collision scheme was applied in the Transient Adaptive Sub-cells (TAS) [START_REF] Bird | Visual DSMC Program for Two-Dimensional and Axially Symmetric Flows, The DS2V Program Users Guide[END_REF] in order to improve the DSMC spatial accuracy. To be consistent with the continual approach, DSMC uses hard-sphere monoatomic gas collision model.

Results

We present mainly results for subsonic and transonic regimes because the calculations of a supersonic case required significant computational times because for transonic and supersonic regimes a flow transition occurs at very small Knudsen numbers between 0.001 and 0.0001. Additional complications raised for the continuum solver with second order TVD approximations of convective terms when non-physical oscillations appeared (especially at lower Knudsen numbers) in the numerical solution. This forced us to use first order upwind scheme requiring a very fine grid. As well-known the DSMC computational requirements also increase for lower Knudsen numbers. Nevertheless, the position of the shock wave in front of the square for B=3 turned out to be a very sensitive parameter requiring a detailed computational analysis. To illustrate the later case, mesh and particles number convergence of SIMPLE-TS and DSMC as shown in Figure 4 for Kn=0.02. It should be noted also that even the SIMPLE-TS was applied on fine grid the numerical solution was not sufficiently accurate in space to determine more sharply the position of the transition curve. To improve accuracy of the transition curve position, a new second order scheme and adaptive mesh was needed to determine place of the shock wave with reasonable computational resources. This circumstance explains why we restricted our detailed investigation up to M=0.6. The comparisons between DSMC and Navier-Stokes-Fourier continuum model [START_REF] Stefanov | Rayleigh-Bénard flow of a rarefied gas and its attractors. I. Convection regime[END_REF], [START_REF] Shterev | Pressure based finite volume method for calculation of compressible viscous gas flows[END_REF], [START_REF] Shterev | Strouhal number analysis for a Karman vortex gas flow past a square in a microchannel at low mach number[END_REF], show as expected a better agreement for smaller Knudsen numbers when the gas flow is close to continuum regime. Here we present a comparison between DSMC and SIMPLE-TS for B=3, M=0.6 and Kn=0.05. The channel length for this case is L ch = 300 while the inlet area length is L a = 200. The DSMC calculation with 7.5 million particles in computational domain, Figure 5 and Figure 6, shows an excellent agreement with the calculated by using SIMPLE-TS macroscopic fields of p, u, v, T and profiles of p and u. The blockage ratio B is bounded by two limits. One limit is B=1 that corresponds to movement of a piston in a channel. The fluid does not pass around the square and the length of disturbances in front of the square (correspondingly, drag coefficient) increases unlimited. The other case is when B tends to infinity. This situation corresponds to a free stream flow past a moving square cylinder in unconfined domain. For the unconfined flow we have used appropriate boundary conditions at the bottom and top boundaries of computational domain replacing the walls with free stream conditions as follows: p = 1, the other independent variables u, v and T were extrapolate linearly from the two closest non-ghost cells, or first-order extrapolation. The largest blockage ratio B=20 considered in this paper is sufficient to neglect the boundary effects and replace the wall boundaries with the described above free stream boundary conditions. Figure 7 shows a comparison of pressure profiles along x-axis between an unconfined flow past a square cylinder and the flow past a square confined in a microchannel for B=20. The length of disturbances in front of the square differs significantly. We consider that disturbance length in front of the square (L d ) is a distance from the front side of the square to the place where nondimensional pressure drop down to 1.001. length L d = 74 for the unconfined flow is significantly smaller then L d = 6138 for the flow past square in a microchannel. The difference between drag coefficients varies in wide range between 10 -20 percents for the considered range of Kn numbers, see Figure 8a. The results show that even for B=20 the flow "feels" the walls and the boundary conditions have a significant impact on the macroscopic quantities, especially, on the length of disturbed area in front of the square cylinder.

Figure 8 shows the drag coefficient as a function of Knudsen number for M = 0.2, 0.4 and 0.6 and B = 3, 5, 10 and 20. The Reynolds number is defined by the relation

Re = 16 5π γπ 2 M Kn ≈ γπ 2 M Kn , (14)
which is valid for compressible viscous model at small Knudsen numbers (see [START_REF] Stefanov | Comparison between Navier-Stokes and DSMC Calculations for Low Reynolds Number Slip Flow Past a Confined Microsphere[END_REF]). The local Mach number has been expressed using nondimensional variables as follows:

M local i,j = 2(u 2 i,j + v 2 i,j ) γT i,j
, for i = 1, 2...N x, and j = 1, 2...N y, The square size (equal to 1) compared to examined disturbance length (5000) is small and hardly visible. The pressure profile in front of the square is given from channel inlet to x = -0.5 while the pressure profile after the square is from x = 0.5 to the channel outlet.

where N x and N y are number of mesh nodes along x-axis and y-axis respectively.

The blockage ratio affects significantly flow properties such as drag coefficient, pressure at stagnation point of the square and the transition point between flow regimes. The dependence of the drag coefficient on the Knudsen number for different Mach numbers (in the same sub-figure) and different blockage ratios (in different sub-figures) is given in Figure 8. The range of values of each graphic is chosen large enough to capture the transition point between the flow regimes, where the drag coefficient curve undergoes a jump or a sudden change of its behavior when the flow transforms from one mode to another. According to eq. ( 14) the drag coefficient could be presented versus Reynolds number. However, as one can see Re, Mach and Kn are coupled and this does not allow to illustrate the effect of rarefaction on the transition point. In the vicinity of the transition point the flow behind the square in the steady state forms two symmetric vortices. With decrease of Kn at the same Mach the vortices become longer, respectively, more unstable and the flow looses the symmetry in this zone. The critical regime starts with a separation of the longer vortex followed by the separation of the second one. This mechanism is repeated periodically and the Karman vortex street is formed. In the zone, where the flow is stable the drag coefficient behavior is typical and it increases when the Knudsen number increases. On the other hand, the smaller blockage ration leads to a larger drag coefficient (C D ). With respect to the Knudsen number, the lower blockage ratio moves the critical transition point left, to lower Knudsen numbers, while the higher Mach number moves it right, to higher Kn values. Thus, the drag coefficient C D in the critical transition point for B=20 is around 1.65 to 1.72 for the considered Mach numbers. For B=10 the transition point moves to values of C D around 1.8-1.9, for B=5 to 2.44-2.55 and for B=3 the drag coefficient increase is up to 4.05-4.35 for considered Mach numbers. The drag coefficient C D at the transition point increases when Mach number increases. An exception case is C D for B=3 and M=0.6 because the gas flow is transonic (see Figure 10), while in all other cases considered in Figure 8 Lower blockage ratio decreases the cross section around the square, increases the velocity around the square and shift back the maximum velocity regions observed on both sides of the square, see Figure 10. The figures demonstrate the developed Karman vortex street past the square for this range of Knudsen numbers. The influence of the channel walls is also clearly seen when compared to the unconfined case (a). Figure 13 shows pressure variation at stagnation point for the closest to the transition point steady state regime. The pressure at stagnation point increases when B decreases as expected. The gas flow in a microchannel at low Mach number (M=0.1) leads to long propagation of small disturbances in front of the channel because the flow regime is close to incompressible limit. The precise calculations require extremely long microchannel and a lot of computational resources. A shorter channel did not contain the full length of disturbances but it determines the pressure at stagnation point within negligible differences compared to a much longer channel and one cannot see any difference in the position of the transition point between steady and unsteady state regime. On other hand, when the Mach number goes up the disturbance length decreases but the impact on the pressure at the stagnation point is bigger, see Figure 13. The inlet length cutoff for M=0.4 shows that a pressure difference at the stagnation point of 12% resulting in moving the trasnition point with respect to Knudsen number by about 20%. This short analysis demonstrates the importance of the pressure calculation at stagnation point within required accuracy. One can find a proportionality of the distances between the pressure curves and the distances between corresponding transition curves on Figure 14a.

The transition curve shape is similar for all considered B (Figure 14). When the blockage ratio is smaller the channel walls are closer to the square walls and stabilizes the gas flow (see Figure 10, Figure 11 and Figure 12). This effect shifts the corresopnding transition curve to lower Knudsen numbers. On the other hand, the smaller B leads to an increase of velocity near the square. This effect shifts the transition curve to lower Mach numbers.

Conclusions

We analyzed the influence of blockage ratio on the transition between steady and unsteady final state regimes of gas flow past a square cylinder in a microchannel. The results obtained by continuum model (SIMPLE-TS solver) were compared to the DSMC data for the case B=3, M=0.6 and Kn=0.05 showing an excellent agreement with the molecular model. The comparisons between unconfined and confined square in a microchannel illuminated the fact that even for large blockage ratio as B=20 the influence of channel walls on the flow around the square is significant. Presented are the transition curves for blockage ratios B=3, B=5, B=10, and B=20. The analysis showed that the transition curves for considered blockage ratios possess similar shape and the decrease of B leads to shifting of the transition curve to lower Knudsen and Mach numbers.

Figure 1 :

 1 Figure 1: Flow geometry for a square cylinder with size a confined in a channel with length L ch , height H ch and boundary conditions (BC) at inlet (BC in ) and outlet (BC out ).

Figure 2 :

 2 Figure 2: Time evolution of drag coefficient for B=3, M=0.6, Kn=0.008.

Figure 3 :

 3 Figure 3: Pressure along the channel central line (y = 0) for B=3, M=0.6 and Kn=0.008.

Figure 4 :

 4 Figure 4: Horizontal nondimensional velocity along the channel central line (y = 0), comparison between SIMPLE-TS (continuum approach) and DSMC (molecular approach) with Transient Adaptive Subcell techniques TAS for M=2.43, Kn=0.02.

Figure 5 :

 5 Figure 5: Comparison fields obtained by SIMPLE-TS and DSMC for B=3, M=0.6 and Kn=0.05 for: (a) p, (b) u, (c) v and (d) T .

Figure 6 :

 6 Figure 6: Comparison of profiles obtained by SIMPLE-TS (solid line) and DSMC (circles) for B=3, M=0.6 and Kn=0.05 along the channel central line (y = 0) for: (a) pressure and (b) horizontal component of velocity.

Figure 7 :

 7 Figure 7: Pressure along the channel central line (y = 0) for B=20, M=0.6 and Kn=0.02. Note: The square size (equal to 1) compared to examined disturbance length (5000) is small and hardly visible. The pressure profile in front of the square is given from channel inlet to x = -0.5 while the pressure profile after the square is from x = 0.5 to the channel outlet.

  the flow is subsonic. For B=3 and M=0.6 the drag coefficient C D have 3 parts: C A D , C B D and C C D . Part C A D corresponds to a steady state regime at larger Kn, where C D monotonically decreases when Kn decreases. For intermediate Kn, in the zone of the transition point from steady to unsteady-state regime the neutral curve changes its slope and starts increasing when Kn decreases (part C B D ). The drag in this part has small-amplitude but the energy conditions are changing sweepingly to jump on the third branch with C C D , see Figure 9a. Its slope changes once again for part C C D and C D decreases when Kn decreases in part C A D . C D oscillations in C B D increase significantly according to oscillations in part C B D for small change of Knudsen number, see Figure 9b.

Figure 8 :Figure 9 :

 89 Figure 8: Drag coefficient (C D ) as a function of Knudsen number (Kn) for the cases of M=0.2, 0.4 and 0.6 and unconfined flow and blockage ratio B=20 a), B=10 b), B=5 c), and B=3 d). Note: the drag coefficient curve is separated into two parts with respect to the transition point: the left part represents the drag coefficient in unsteady regimes while the right part -in steady-state regimes.

Figures 10 ,

 10 11 and 12 illustrate local Mach number, pressure and temperature flow fields for unconfined case (a) and different blockage ratios (b), (c), (d), and (e). The Knudsen varies in the range between 0.01 and 0.02.

  (a) Unconfined flow past a square and Kn=0.0190 (b) B=20 and Kn=0.01815 (c) B=10 and Kn=0.0145 (d) B=5 and Kn=0.0114 (e) B=3 and Kn=0.00885

Figure 10 :

 10 Figure 10: Unsteady state regime local Mach number fields and streamlines obtained closest to the transition curve for Mach number 0.6 and unconfined flow past a square and Kn=0.0190 a), B=20 and Kn=0.01815 b), B=10 and Kn=0.0145 c), B=5 and Kn=0.0114 d), and B=3 and Kn=0.00885 e).

  (a) Unconfined flow past a square and Kn=0.0190 (b) B=20 and Kn=0.01815 (c) B=10 and Kn=0.0145 (d) B=5 and Kn=0.0114 (e) B=3 and Kn=0.00885

Figure 11 :

 11 Figure 11: Unsteady state regime pressure field obtained closest to the transition curve for Mach number 0.6 and unconfined flow past a square and Kn=0.0190 a), B=20 and Kn=0.01815 b), B=10 and Kn=0.0145 c), B=5 and Kn=0.0114 d), and B=3 and Kn=0.00885 e).

  (a) Unconfined flow past a square and Kn=0.0190 (b) B=20 and Kn=0.01815 (c) B=10 and Kn=0.0145 (d) B=5 and Kn=0.0114 (e) B=3 and Kn=0.00885

Figure 12 :

 12 Figure 12: Unsteady state regime temperature field obtained closest to the transition curve for Mach number 0.6 and unconfined flow past a square and Kn=0.0190 a), B=20 and Kn=0.01815 b), B=10 and Kn=0.0145 c), B=5 and Kn=0.0114 d), and B=3 and Kn=0.00885 e).

Figure 13 :

 13 Figure 13: Pressure at stagnation point vs Mach number.

( a )Figure 14 :

 a14 Figure 14: The transition (neutral) curve as a function of Knudsen and Mach numbers.

Acknowledgements

We would like to acknowledge the financial support provided by the Bulgarian NSF under Grant DN-02/7-2016. This research has received funding from the European Union Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie MIGRATE grant agreement No. 643095.


Dominique G Béroule
email: dominique.beroule@limsi.fr

Pascale Gisquet-Verrier

Decision Making guided by Emotion A computational architecture

Keywords: retrieval cues, conditioning, neuromodulation, emotional valence, coincidence detection, Guided Propagation Network I

A computational architecture is presented, in which "swift and fuzzy" emotional channels guide a "slow and precise" decision-making channel. Reported neurobiological studies first provide hints on the representation of both emotional and cognitive dimensions across brain structures, mediated by the neuromodulation system. The related model is based on Guided Propagation Networks, the inner flows of which can be guided through modulation. A key-channel of this model grows from a few emotional cues, and is aimed at anticipating the consequences of on-going possible actions. Current experimental results of a computer simulation show the integrated contribution of several emotional influences, as well as issues of accidental all-out emotions.

INTRODUCTION

During its long-lasting evolution, the brain had to deal with a main physical constraint undergone by signal processing devices. Sensory systems are notably challenged by Uncertainty Principles in Fourier Analysis [START_REF] Cowling | Bandwidth versus time concentration: the Heisenberg-Pauli-Weyl inequality[END_REF], according to which: "the more a given analyzer locates a signal in the time domain, the less it can locate it in frequency/space and vice versa". In other words, setting up a sensory analyzer faces the following alternative: making it either "swift and somewhat ill-sighting" (like the peripheral vision), or "rather slow but sharp-eyed" (like the central vision). For surviving the primitive world predators, better being fast. This is why the ancient brain probably took up the first alternative, still operating today, on top of which in-depth but slower analyses have then developed as the environmental pressure decreased. The modern brain works with several parallel analyses of its environment, each tuned to a particular physical dimension, thus coping with the "uncertainty principle". The question of the interaction between these analyses then arises, as well as the possible generality of this perceptive architecture, e.g.: Could this sensory split in complementary channels extend to higher-order cognitive processes, as proposed in "Dualsystem" models of decision-making [START_REF] Franck | Multiple Systems in Decision Making: A Neurocomputational Perspective[END_REF]? This question will be addressed in the first section, with references to neurobiological data: evidence coming from the literature suggests that several brain regions, as well as neuromodulation systems, are busy selecting proper behaviors.

Beside the main stream of Artificial Neural Networks (ANN), Guided Propagation Networks (GPN) rely on inner flows of activity which are controlled by neuromodulation-like processes. This main feature guided our choice towards their use in the present study, given the recognized link between emotion and neuromodulation. Following a reminder of the GPNs main characteristics, a specific architecture is gradually introduced in the second section. A software implementation of this model has been run on the same training data with different types of internal modulation and input cues. Related results are reported in the last section.

II. CEREBRAL BASIS OF DECISION MAKING

A. Emotional marking of memory traces

For living organisms, dangerous situations must be identified as quickly as possible for an appropriate response to be elicited. In plants, a water deficit triggers chemical signals inducing the inhibition of photosynthesis, among other metabolic changes. For organisms owning brains, inner warning signals can even be delivered long before the occurrence of a threatening event. This happens thanks to the memory function, which keeps track of experiences leading to meaningful events. In order to anticipate dangerous situations, the uncertainty principle suggests that associated cues should better be memorized in their fuzzy rather than high-definition format. These so-called emotional cues seem more efficient, due to their swift recognition. Such a signal processing constraint is consistent with the fast warning circuit aroused by the amygdala, a tiny brain structure of the ancient limbic system [START_REF] Ledoux | The Emotional Brain, Fear, and the Amygdala[END_REF]. Fast decisions may involve emotions through physiological signals associated with their past outcomes, as proposed by the Emotional Markers hypothesis [START_REF] Damasio | Somatic Markers and the Guidance of Behavior: Theory and preliminary testing[END_REF]. The functional part of emotion would include the enhancement of the behaviors memorized as a rewarding experience. Conversely, actions having previously led to unpleasant events would likely be avoided. Summing up, the emotional system can be viewed as a background "swift channel" fed by fuzzy cues which have acquired an emotional salience, meant to influence the decisions made by slower and more precise channels.

Within a topological memory device aimed at decisionmaking, can we expect decisions to only rely on markers fixed at a particular acquisition time? In this case, among several patterns of activation then aroused, the negatively marked items would always be avoided in favor of those associated with a positive valence (emotional value). However, as the same cue may be followed by situations of opposite valences, decision-making should account for several on-going factors. Moreover, thanks to their ability to memorize events developing over large time-spans, human-beings can project in a distant future. If every negatively marked situation were systematically avoided, nobody would bear disagreements possibly leading to delayed benefits. Who would dare writing a scientific paper if not motivated by the perspective of either being nobelised one day, or more likely enjoying a conference backstage in Australia? The emotional value of the uncompleted paper cue may then take opposite valences, depending on the enhanced time scale, which may vary as a function of the writer's mood at the moment the draft enters her/his peripheral vision. Again, the valence of a memory trace cannot be accurately determined in advance, but has to be established at the time of memory retrieval, instead.

B. Emotional dimension across brain structures

Decision-making is a complex process supported by the prefrontal cortex in which different roles are played by subregions, namely -orbitofrontal (OFC), anterior cingulate (ACC) and dorsolateral prefrontal (DLPFC). The OFC would notably code for the reinforcing value of the diverse stimuli encountered by the organism, and be involved in guiding action [START_REF] Holland | Amygdala circuitry in attentional and representational processes[END_REF]. It would for instance be central for inhibiting responses that already led to undesirable consequences, and for enhancing an appropriate behavior despite the stronger activation provided by alternative ones [START_REF] Miller | An integrative theory of prefrontal cortex function[END_REF]. The guiding part played by this cortical structure, possibly in association with the ACC [START_REF] Kennerley | Double dissociation of value computations in orbitofrontal and anterior cingulate neurons[END_REF], would favor actions that revealed rewarding in the past experience.

From a more general point of view, decision-making involves different brain regions working together to evaluate environmental stimuli and transform this information into appropriate actions. Accordingly, this process requires skills in cognition, emotion, memory and motor control, all carried out by structures in interaction with reward circuits. For instance, it has been shown that the emotional reactivation provided by retrieval cues facilitates the retention performance, and hence decision-making. These cues strongly activate the amygdala, a brain structure coding for emotion, as well as the nucleus accumbens which is a part of the reward system [START_REF] Boujabit | Exposure to retrieval cues in rats induces changes in regional brain glucose metabolism: Role for the amygdala in the retrieval of memory?[END_REF]. Recent literature in rodents, monkeys and humans is consistent with the idea that a cortico-basal ganglia circuit constitutes the heart of the reward processing. This circuit works with the retrievaloriented DLPFC, and the emotion-oriented amygdala associated with regulating structures releasing neuromodulators such as dopamine (from the ventral tegmental area) and serotonin (from the dorsal raphe). Such cortico-striatothalamo-cortical loops have a dual organization that permits both parallel and integrative processing. More precisely, as depicted in Figure 1, fibers from different prefrontal areas converge within subregions of the striatum: through the organisation of the striato-nigrostriatal projections, the ventral striatum can influence the dorsal striatum. The non-reciprocal cortico-thalamic projection conveys information from rewardrelated regions, through cognitive and motor controls. The initial discovery that the ventral striatum funnels emotional information towards motor outcomes has been investigated in non-human primates for building the Ascending Spiral Model [START_REF] Haber | Striatonigrostriatal pathways in primates form an ascending spiral from the shell to the dorsolateral striatum[END_REF] [START_REF] Haber | The reward circuit: linking primate anatomy and human imaging[END_REF]. Such a view is close to the one supported by brain imaging studies in human, which postulate that cortical regions are functionally linked through a "cascade" of interactions [START_REF] Koechlin | The Architecture of Cognitive Control in the Human Prefrontal Cortex[END_REF]. According to this model, three main loops (limbic, associative and sensori-motor) form a spiral-like pattern between the cortex and the striatum. Each step would comprise two parallel projections from a given cortical region, respectively towards its dorsal and ventral striatum partner. While the dorsal striatum would send an inhibitory feedback response to its cortical afferent, the ventral striatum would "facilitate transfer of information to the next step in the spiral", by modulating the adjacent cortical region. The latter would follow the same track through the striatum to influence its neighbour cortical region, and so forth until the final motor cortex (Fig. 1).

Such a system allows the anticipation of both possible rewards and punishments to favour a proper behaviour. The related facilitation or inhibition can be obtained by different means. These effects may involve specialized cortical regions: within the OFC, rewards tend to activate the medial part, while punishments seem to activate more lateral regions, instead. The dopamine, which is largely involved in the reward system, may exert opposite effects through a possible connectivity pattern proposed in the Ascending Spiral Model (Fig. 1). Finally, dopamine release may also induce opposite effects depending on the type of binding receptors. Although the main function exerted by dopamine appears to be mediated by inhibitory D2 receptors, experiments also report the critical role of arousing D1 receptors in reward-related learning [START_REF] Beninger | Dopamine D1-like receptors and reward related incentive learning[END_REF].

C. Impairment of emotional circuits

Whereas most of our everyday choices do not cross our level of consciousness, patients with prefrontal damage lack this intuitive parallel evaluation of all possible options [START_REF] Bechara | Insensitivity to future consequences following damage to human prefrontal cortex[END_REF]. These patients are unable to perform optimal and rational decision-making. They miss the ability to anticipate the longterm consequences of their decisions, and preferably choose actions with the best short-term outcomes. Severe decisionmaking impairments are always associated with lesions of the prefrontal areas (OFC, ACC and DLPFC), an effect also observed in rat, probably due to a disruption of the basic modulating function [START_REF] Botreau | Effects of medial prefrontal cortex and dorsal striatum lesions on retrieval processes in rats[END_REF].

A unifying view of two pathologies linked with emotion has recently been proposed, according to which exposure to extreme reinforcers of behaviour such as a trauma (very negative) or a drug of abuse (very positive) may both induce hypersensitivity to retrieval cues. This change of sensitivity would be responsible for drug craving in addiction, intense reviviscences in Post Traumatic Stress Disorder (PTSD), as well as relapse in these two pathologies [START_REF] Gisquet-Verrier | Hypersensitivity to cue-elicited memory reactivation as a possible source for psychiatric pathologies such as relapse to drug addiction and post-traumatic stress disorder[END_REF].

III. NETWORK ACTIVITY GUIDED BY EMOTION

At the origin of ANNs, the Perceptron was already aimed at decision-making, namely for the categorization of visual patterns [START_REF] Rosenblatt | The perceptron: a probabilistic model for information storage and organization in the brain[END_REF]. As this research stream enlarged with other formalisms, a great demand arose for aiding professionals in their decision tasks. Satisfying this demand has first been addressed with computer tools close to the initial ANN architectures [START_REF] Hill | Artificial neural network models for forecasting and decision making[END_REF], then incorporating functional principles motivated by psychological data linked with research on emotion [START_REF] Levine | Neural network modeling of emotion[END_REF]. More recent developments have followed advances in the physiology of emotional circuits of the brain, including elements of neuromodulation [START_REF] Eckhoff | Optimality and robustness of a biophysical decision-making model under norepinephrine modulation[END_REF].

A. Sprouting of memory-pathways within modules

A Guided Propagation Network (GPN) is a modular architecture. Depending on its location within a matrix of channels and layers, each GPN module is focused on a particular representation of events (e.g.: typed words; Fig. 2). The same generic mechanisms are implemented in every module, all based on a main feature: an inner flow of activity. This spontaneous activity needs to be fed by a series of incoming stimuli for being able to cross the module and eventually reach its output where Feature-Detectors (FD) stand. Both rhythm and direction of the inner flow are given step by step while a given input pattern develops in time, leading to the full activation of the pattern-FD at the module output. The joint occurrence of the inner activity and input serial patterns implies that, at any moment, the front of the inner flow indicates the next stimuli to possibly occur, hence its predictive power.

In order to enable an unexpected stimulus and the inner flow to "hold future meetings", a new Elementary Processing Unit (EPU) can be brought into play in the course of processing (see Fig. 2). When a stimulus is found which did not participate in guiding the module inner flow, a EPU is taken from a pool, and assigned the role of coincidence detection between the stimulus and the current front. The memory location of a new EPU is thus given by the inner flow, provided the latter is focused enough. Little by little, chains of EPUs are created all along the network life by a non-supervised learning algorithm [START_REF] Béroule | The Never-ending Learning[END_REF], forming pathways comparable to the branches of a tree. FDs standing at the end of the branches may guide the activation of more central modules. In other words, GPN-tree seeds are all sowed at once from the end of a given branch towards the ramifications of another module below. This is a waste of energy, but there is no time-loss, hence the theoretical real-time processing of GPNs.

B. Guided propagation through modulation

In a former GPN-model of reading [START_REF] Béroule | Guided Propagation inside a Topographic Memory[END_REF], both precise perception of letters and scanner jumps were driven by the 'fuzzy and swift' perception of the text ahead of the letter currently fixed; the vocabulary subset sharing the same global shape could thus be enhanced. Despite its early response, this gross analysis was not to interfere with the on-going slow recognition process held in the neighboring channel. The solution taken on had been for a swift channel to decrease the decision thresholds of the EPUs participating in the representation of anticipated patterns. The upstream facilitation of EPU pathways has then been incorporated into the formalism, which incidentally opened the way to pattern production [START_REF] Toffano | A Functional Model of some Parkinson's Disease Symptoms using a Guided Propagation Network[END_REF].

Any given GPN module has therefore the potential for either: 1-acquiring, 2-identifying or: 3-generating time-space patterns in several application fields [START_REF] Béroule | An instance of coincidence detection architecture relying on temporal coding[END_REF].

A Modulation device is responsible for setting on-line the EPU parameters, which determine the module function among the aforementioned three ones. While a module is fed by a stimulus, its inner flow may thus be facilitated by upstream modulation signals (Fig. 3).

Figure 2. GPN module coding for words on the basis of one letter after the other (at the top). Two pathways (chains of EPUs: round cells, in dark) have already been built by the learning algorithm, when the following letters occurred for the first time, in succession:

-'s' at time t1, 'i' at t1+t, 'x' at t1+ 2t ; -'s' at t2, 'a' at t2+t, 'x' at t2+ 2t. Now, 'a' has just occurred after 's', but the 'd' Feature-Detector (FD, square cell) has just fired instead of the expected 'x'. A new pathway is thus being created (dotted items), ending with a new module output: the labeled 'sad' FD (at the bottom-right corner).

The inner flow originating from the top-left round cell can then be guided towards one of the current module FDs ('six', 'sax', 'sad'). The part of these word-FDs compares to the letter-FDs one, since the former may guide the inner flow of a more central module (not plotted here). The EPU standing at the "meeting spot" of the module incoming signals, is activated above its decision threshold n, expressed as a fraction of the EPU maximum input (where A is the maximum signal intensity, En the EPU Excitability and Rn the weight of the "Context input" fed by the inner flow). Above its threshold, the EPU output is proportional to its summed input, up to the A saturation value.

n n n E A R     ) 1 ( (1)
Although a module inner flow has a pace usually made by current incoming stimuli, it can also be induced to run ahead the present time. This is a precious asset in the framework of this study, for a module to be able to anticipate salient cues conveying an emotional valence. The inner flow is then reinforced by shifting the parameters introduced in formula (1) from the "restricted" to the "extended" mode (see Fig. 7), and this for every EPU of its become "far-sighted" module.

The logic behind the resulting multiple-channels architecture can be expressed in this way: the "aware" channel (Decision Maker, labeled III in Fig. 5) keeps pace with present stimuli, including proprioceptive ones from actions being performed (channel IV). In order to quickly select the future right action among possible ones, channels III and IV need to be guided by modulation signals (Fig. 4) issued from another channel (Cues Binder, II) in which the pathways of future events are activated in advance. The respective emotional values of these anticipated possible events need to be given by another channel (Emotion Detector, I).

C. Conditioning through coincidental module outputs

In a simulated "ascending spiral" architecture, channels I and II only respond to fuzzy cues, whereas channel III is activated by well-defined stimuli. Quite a few emotional cues initially feed the channel I. The recruitment of new emotional stimuli permits decision-making to be guided in a more subtle way. Through the sprouting of a new pathway in channel II, any neutral fuzzy stimulus can be paired with an emotional cue occurring shortly after it. In channel III, another pathway is simultaneously built, which receives the precise version of both stimuli (Fig. 5).

When the behavior begins with a neutral stimulus (e.g.: the bell in Pavlov's dog trial), and ends with an emotional cue (i.e.: the sight of meat, or Unconditioned Stimulus US), three parallel pathways are activated in their respective GPN channels, which results in cross-connections between them. This acquired connectivity pattern is aimed at conveying modulating signals from one module to its right-hand side neighbor, forming the ascending spiral. The channel I outflow selectively facilitates channel II pathways: each repetition of a sequence ending by an emotional cue increases the inner flow along a pathway of channel II, making it draw ahead the current input. In this way, a neutral stimulus becomes "conditioned" (CS), i.e.: able to fully activate a channel II pathway as if the US had actually occurred (e.g.: in the Pavlov trial, the dog salivates as soon as the bell rings). This schema applies to aversive conditioning as well, through which moving towards cues associated with stressful events is prevented. As far as the 'II=>III' cross-connection is concerned, the emotional value aroused by a channel II output induces the modulation of propagation within channel III (Fig. 5).

In order to implement 2 nd order conditioning (e.g.: a red light preceding the conditioned bell sound induces salivation),

Rn+1

Cha nnel III a stimulus which repeatedly triggers one of the US responses will join in the channel I input. This mechanism results in a "funnel effect", which tunes a subset of cues with the reward. Conversely, in aversive conditioning, "an umbrella is opened" under all the cues associated with a negative emotion, thus trying to avoid its revival. This capacity requires the parameters to be properly settled, so that a channel II pathway gets fully activated under the only influence of its initial stimulus.

In the GPN theory, the increment of the contextual weight (Rn parameter) along a pathway accounts for the integration of the previous stimuli, that is:

n n R R n     1 , 1 (2)
The output Cn of the n th EPU along a pathway (Fig. 6) is expressed above the EPU threshold as a linear function of the summed Stimuli Sn and the weighted output Cn-1 of the previous EPU.

) 1 )( 1 ( ) ( 1 1 2 1 1 1                 n n n n n n n n n n n n R R C R S R S R C R S C
In which:

0 1    n n S S and A C n  2 ) 1 )( 1 ( 1 1        n n n n n R R A R R C
For the internal flow to freely cross two EPUs without stimuli, the upper boundary of Excitability En+1 is such that the weighted output stays below the next EPU threshold (see (1)):

1 1 1 1 1 1 1 1 ) 1 ( ) 1 )( 1 ( ) 1 (                     n n n n n n n n n n n E A R R R A R R R E A R C R
Taking (2) into account, one gets: 2

1 1 1 1       n n n R R E (3)
For the inner flow to cross a single EPU, a similar rule can be obtained: 1

1 1 1 1       n n n R R E (4)
Returning to the assumption underlying a GPN model of Parkinson's disease symptoms [START_REF] Toffano | A Functional Model of some Parkinson's Disease Symptoms using a Guided Propagation Network[END_REF], the Rn parameter shown in Fig. 7 would compare to serotonin, whereas En would play a part akin to dopamine in the neuromodulation system.

D. Merged time-scales for Decision-making

Drugs of abuse are associated with immediate rewards, and cause serious trouble in the longer term. In that case, anticipating social and personal degradations may not be strong enough to counterbalance drug seeking. Figure 7 (adapted from [START_REF] Béroule | An instance of coincidence detection architecture relying on temporal coding[END_REF]). Couple of EPU control parameters. The Excitability En (vertical axis) and the inner flow input weight Rn (horizontal axis) determine the behavior of a given EPU. In the "restricted propagation" area (in white), a newborn EPU requires both its Stimulus and inner flow (Context) inputs to fire. The distribution of round dots follows the required arithmetic progression of Rn along a pathway of EPUs. In the (pebbly) "extended" propagation mode, En and Rn (square dots) allow older EPUs to fire thanks to the only inner flow. After having crossed the border between the "restricted" and the "extended" areas, several chained EPUs may get activated at once as the ( 4) and ( 3) conditions (see text) are satisfied (diamond shapes).

Conversely, an immediate effort may later on be rewarded. This dual time-scale influence can be modeled by taking advantage of the hierarchical nature of GPNs. Two layers with independent modulating circuits can be brought into play; one (central) layer integrating the stimuli of a more peripheral one. Thus, both modules of channel II may guide decision-making, sometimes with opposite modulation values (see Fig. 8).

IV. EXPERIMENTS

Experimental issues are investigated on the same instance of GPN, obtained through a preliminary training session. The network is first set according to the global architecture shown in both Fig. 8 and Fig. 9, where five modules initially own a single EPU at the source of their inner flow.

The network interface depends on patterns to be processed. Here, a compound behavior is represented by a sequence of short words, in which stimuli (consonants) alternate with actions (vowels). The following Letter-Detectors are involved: -A few unconditional stimuli or actions feed the channel I (i.e.: symbols 'd' of valence {-2}, 'r' {+2}, 'y' {+2}); -Twentyinitially neutralstimuli (consonants) in their fuzzy format (i.e.: lower-case) at the top of channel II; -Twenty stimuli in their high-definition format (i.e.: uppercase consonants) form the input of module IIIp; -Six action effectors (i.e.: upper-case vowels) are ready to deliver proprioceptive stimuli after having generated the action triggered by module IIIp.

A. Script of the experiments

An "elementary behavior" thus takes the shape of a short word. The software is reading a file, each line of which contains 2 or 3 such words (e.g.: 'BAD DOG GRUNT', 'BUZ ZARD'). While a file line is scanned, the concerned Letter-Detectors are activated one after the other at the network input.

Starting from nearly scratch, the software goes through the training of 100 combinations of 170 words which eventually brings 1200 EPUs and 4700 links in the current simulation. In the data, each first consonant is possibly followed by any of the vowels; in other words, after an initial cue, all available actions may be possible, among which the system will have to choose. This training period compares to the exploration of an environment comprising a few preset emotional cues; a series of actions in this environment may thus either be neutral or be concluded by a salient cue (i.e.: with valence {+/-2}).

When a new couple of words is read from the data file during training, memory pathways are sprouting in parallel inside every concerned module. Two elementary behaviors are memorized in IIIp, while IIIc grows a pathway combining them. The channel II only integrates a fuzzy version of the stimuli contained in these behaviors. For instance, the 'LOCK KEY' input results in the following pathways:

-in module Ip : one-EPU pathway labeled 'y' -in IIp : 'lck' and 'ky' / -in IIc : 'lckky' -in IIIp : 'LOCK' and 'KEY' / -in IIIc : 'LOCKKEY' An input ending with an emotional cue (i.e.: 'd', 'r', 'y') leads to new cross-modulation connections between modules outputs. As described in the previous section, these modulation links permit conditioning to operate in channel II with the help of channel I, as well decision-making in channel III, guided by channel II (Fig. 8). In this experiment, cues becoming "emotional" through 2 nd order conditioning are assigned half the valence of their conditioning stimulus, namely either {-1} or {+1}. At the end of training, depending on the number of repetitions required for a cue to be conditioned (parameter to be fixed beforehand), several new FDs have been created in the channel I input. The software is run on the training data (which takes 4 minutes on an Intel Pentium4 running at 1.5 GHz with 512 Mo of RAM). Given that new emotional cues appear in the course of training, some of them may have been created after the presentation of an input they could influence, hence a second presentation of the trial data.

During the following "test session", the system is given a cue (i.e.: a consonant), and is then free to decide what action to perform among the options previously met in the same context during training. The challenge for the IIIp inner flow is to be guided towards the most rewarding emotional target (Fig. 10).

It is worth noticing that less than 3% of the decision maker EPUs are firing simultaneously in the 'aware' channel III during this process, while around 6 times more EPUs are involved in the 'subconcious' channel II.

% B. An evaluation of each GPN module contribution

A channel II module can be withdrawn from the network either by switching off its inner flow generator, namely its root EPU, or by voiding its modulating output. Submitted to the same series of input stimuli, the output of the full architecture can thus be compared with the outcome of the system deprived of some modules. Among the alternative above, the weighting of modulation has been chosen for it offers the finest tuning of each module contribution to decision-making. Wp is the weight of modulation between IIp and IIIp; Wc holds the same part between IIc and IIIc. The results plotted in Fig. 11 are obtained with six different values of (Wp, Wc).

Tests are thus performed with these six types of modulation on 20 input cues (consonants). For each trial, this is the activity of the Effectors (channel IV) induced by a given initial cue which is considered significant. For example, 'B' may be followed by every vowel, including either 'A' or 'O' ('BAD DAY', 'BOY YARD'). The choice between the possible actions may be guided by both the peripheral modulation (IIp->IIIp) of the "word layer" and the central one (IIc->IIIc) of the "word-sequence layer". Given the opposite emotional valences of 'D'{-2} and 'Y'{+2}, the peripheral modulation favors the choice of 'O' because of the positive end of 'BOY' whereas 'A' is inhibited. On the contrary, the central modulation exerted via the last words cues ('DAY' and 'YARD') orientates decisionmaking towards the 'A' choice.

When both Wp and Wc are null, the consequent blocking of modulation drastically depresses the channel III. In absence of emotional information, the EPUs thresholds stay at their resting level, with only tiny differences reflecting the frequency information, namely the number of behaviors repetitions. This is why a general increase of the module EPUs Excitability can partly make up for the loss of action triggering, provided that one of the actions is more frequent than the others in the tested situation.

In current results, it appears that the anticipation of a close emotional cue has a modulating effect at least 3 times stronger than the prediction of a future one. This relative deficit of future cues can be attributed to the longer path taken by the central modulation (Fig. 10). This can also be interpreted in terms of "missing inhibition" from the deeper module IIIc towards module IIIp, given that in GPNs, only facilitating signals propagate upward within a given channel.

As far as channel I is concerned, switching it off after a training session does not impair the system global behavior, but prevents later emotional conditioning.

C. A simulation of hypersensitivity to retrieval cues

During the training session, the valence of the initial emotional cues is set to +/-2, while conditioned stimuli have gained a valence of +/-1. The system can then be submitted to stimuli of extreme valence, which implies the one-shot conditioning of the cues that precede them in the input. Furthermore, the modulation power of these all-out stimuli is so high that the pathways they touch are either fully inhibited or fully facilitated in comparison with conventional emotional modulations (Fig. 11).

Figure 11. Prototypical activity rate of GPN effectors as the result of emotional training. During training, every initial input stimulus was followed by any of the possible actions ('A', 'E', 'I', 'O', 'U'), followed itself by either immediate or future stimuli owning emotional valences. During the "decision making" test, the system had then to choose which action to perform in response to a given input stimulus. This figure shows the effectors activity rate induced by the same input stimulus ('C'), as a function of six different combinations of peripheral and central modulation weights (Wp, Wc) taking the following values along the front axis: 0:(0, 0); 1:(1, 0); 2:(1, 0.5); 3:(1, 1); 4:(0.5, 1); 5:(0, 1). As with any input stimulus of the corpus, no decision is made without modulation (null response of all the action effectors at the origin of the front horizontal axis). Here, the central modulation is not amplified by a factor of 3 (see IVB). Consequently, this is only when the peripheral influence is switched off (tag '5'), that the expectation of a future reward becomes significantly efficient (activation of the 'E' effector at 25% of the maximum intensity). The greater response (90%) is obtained in the current context by the 'I' action, associated in the training data with the highest immediate drug-like reward. When the data does not include this highly rewarding stimulus, the figure obtained is close to this one ('A': 43%; 'E': 25%), but free of the 'I' response (striped area). When no "drug" is experienced during the training session, an effector barely produces a strong response, but the model simulation remains selective, even without lateral inhibition between its action-effectors.

Let us assume that the 'X' character be chosen as a highly rewarding cue associated with drug consumption. After a single occurrence of the 'SIX' input, 'S' enters channel I with a positive valence. Then, any series of actions leading to 'S' (e.g.: I and A in 'BIAS') will be favored, as well as other stimuli leading to the highly facilitated 'SIX'. This way, 'S', as other possible hypersensitized cues, participates in recruiting several other stimuli, possibly in one-shot. This looks as if conditioning could in this case enlarge the "receptive field" of an all-out stimulus, making the system more likely to meet in its environment a cue driving it towards the abuse of drug.

Assuming limited modulation resources, the involved strongbut localfacilitation of some channel II pathways may affect the Modulation device global balance. This system impairment will be studied in connection with the brain neuromodulators imbalance caused by drug addiction [START_REF] Tassin | Uncoupling between noradrenergic and serotonergic neurons as a molecular basis of stable changes in behavior induced by repeated drugs of abuse[END_REF]. The treatment of compulsive drug-seeking and the avoidance behavior generated by a traumatism both evolve on a long-term scale, respectively leading to a possible relapse into addiction or revival of stressful situations. The assumption that similar mechanisms would cause the appearance of these different disorders will also be investigated with a new release of the computational model, in light of neurobiological data [START_REF] Gisquet-Verrier | Hypersensitivity to cue-elicited memory reactivation as a possible source for psychiatric pathologies such as relapse to drug addiction and post-traumatic stress disorder[END_REF].

V. CONCLUSION

Developed from a specific ANN, the model introduced here accounts for neurobiological evidence of the major part played by neuromodulation and emotional circuits in decision-making. The key component of the proposed architecture stands at the interface between the perception of emotional cues and modules oriented towards action. Referred to as "channel II", this component is able to chain fuzzy environmental cues, especially those which end by an emotional stimulus; it can then run ahead the present time and modulate the "decision maker" channel III so that its two internal flows are aspirated towards the most rewarding output, while triggering actions distributed along the preferred (path)way.

This architecture notably builds on the ascending spiral hypothesis linking the midbrain to the prefrontal cortex. A software implementation already confirmed the functionality of the involved mechanisms and allowed to start investigating their issues in the case of extreme behavioral reinforcers. As a representation of biological data, this model raises several expectations that could be considered later on, including: the hierarchical organization of emotional circuits, the greater (anticipatory) activation in subconscious emotional modules than in aware decision-making regions, the enlargement of emotional regions during a life-span, involving the creation of new cue-detectors, and the assumption that emotional consequences of a given situation be subconsciously replayed at every new instance of a similar situation. As a matter of fact, the reconstructive nature of memory is an important feature of retrieval processes, as often underlined since Bartlett [START_REF] Bartlett | Remembering: A Study in Experimental and Social Psychology[END_REF].

Over a life span, the memory landscape is gradually shaped in our mind, including easy routes that are frequently visited. Meanwhile, in the background, a swift system keeps vigil at specific cues of the environment, those which acquired salience in the past, because associated with an emotion. Whenever such retrieval cues are detected, the current landscape would be temporarily reshaped with a combination of hard-to-climb hills and promising valleys. If cues previously associated with an extreme emotion occur, either a high peak or a ravine adds a sharp item to the local relief: the memory of a traumatic episode grows an unreachable plateau, while a drug of abuse digs an inescapable deep hole. Decision-making is therefore restricted, together with the freedom to explore the full landscape. Fortunately, if the theory presented in this paper is right, these strong landscape modifications may not be permanent, which would allow their modulation to be locally controlled, and the causative conditioned cues to become extinct: an opportunity to be grasped in future work.

Figure 1 .

 1 Figure 1. Ascending spirals of connectivity of the striatum to midbrain (downwards) and prefrontal cortex (upwards) (Reproduced from [9]). Magnified oval region to the left shows hypothetical cortico-striatal connections. The upper projection activates directly a dopaminergic cell (a) in the Dorsal tier, resulting in inhibition. The lower projection terminates indirectly on a dopaminergic cell (b) of the Ventral tier via an interneuron, resulting in facilitation. (VTA: Ventral tegmental area, OMPFC: orbitomedial perfrontal cortex, DLPFC: dorsolateral prefrontal cortex).

Figure 3 .

 3 Figure3. Content of a GPN module. In the left image, a trained module contains 'inverted tree'-like pathways, fed by an inner flow. This awaiting stream is guided by 1/ series of stimuli from more peripheral modules and 2/ modulation from other modules (upward arrows). In the magnified window to the right, an Elementary Processing Unit (EPU) owns a threshold n for detecting the possible juxtaposition of its two inputs: stimulus and inner flow.

Figure 4 .

 4 Figure 4. Instance of EPU activity diagram. This EPU is primed by the inner flow at time tP. Its decision threshold (dotted line) then undergoes two modulation signals: an increase (inhibition) at time tI, a decrease (facilitation) at time tF, Thanks to this strong facilitation, the EPU treshold is more likely to be crossed at time tR, when a stimulus meets the inner flow.

Figure 5 .

 5 Figure 5. Representation of the potential interactions between GPN modules via a Modulation device (at the bottom). Both connectivity and color conventions partly come from Fig. 1. Downward straight arrows issued from the modules outputs stand for the activation of D and V tiers, the respective functions of which compare to the structures shown in Fig. 1: whereas D is in charge of inhibitory feedbacks (upward curved arrows), V sends modulating signals towards the next neighbor module to the right-hand side (upward leftto-right arrow). Modulation is either 'facilitating' (+) or 'inhibiting' (-). The modulating signals propagate upstream through channel III (striped areas), thus inducing one of the currently possible actions via the modulating device.

Figure 6 .

 6 Figure 6. Instance of a channel II GPN pathway, showing its associated symbols (Sn : Stimuli, Cn : Contextual input, Rn : contextual weight)

Figure 8 .Figure 9 .

 89 Figure 8. Full architecture of the model. The label of each module gives its location; e.g. IIIp corresponds to the 3 rd channel and peripheral layer, and IIc indicates the central module of channel II. The neuromodulation connectivity pattern plotted in Fig.5is reduced here to straight (blue) arrows at the modules bottom. The outflow of IIc is currently facilitating a pathway of IIIc (striped area), the EPU thresholds decrease of which is propagated upstream. As the IIIc inner flow follows this path, facilitation signals are sent towards a IIIp pathway (upward curved arrows). The latter is here inhibited at the same time by IIp, but not enough for avoiding the action selected by the facilitating flow issued from the central IIc.

Figure 10 .

 10 Figure 10. Activity histograms of pathways involved in guiding a decision through central modules (IIc and IIIc). An initial cue ('M') feeds channels II and III, channel II undergoes a swift cascade of activations across its two modules. Meanwhile, the EPUs linked with next possible actions are primed in IIIp. A facilitation wave (upward curved arrows) then selectively decreases channel III thresholds, making an effector ('I') fire in channel IV.

ACKNOWLEDGEMENT

The study reported in this paper has been made possible thanks to Vincent Bloch who 30 years ago got interested in the first release of GPNs, Jean-Marc Edeline who more recently initiated the connection between the authors, and Vincent Douchamps who put the Ascending Spiral Model to the front. This work has been supported by a "NeuroInformatique CNRS" grant.