T Raharijaona
email: thibaut.raharijaona@univ-amu.fr

J Serres

E Vanhoutte

F Ruffier

Toward an insect-inspired event-based autopilot combining both visual and control events

This paper presents the autopilot and the behavior of a "simulated bee" traveling along two different tunnels using both visual and control events. The computational gain of an event-based PID controller compared to its time-based version is usually discussed because the event detector is computationally expensive. By combining visual and control events, the newly suggested event-based autopilot requires very low computational resources. In particular, the event detector which computes the control error and tests its magnitude is activated only when a new contrast is detected by the optic motion detectors that assesses the optic flow, i.e. only when the magnitude of the optic flow error could have changed. This new event-based control strategy used faithfully the visual information already available in the optic flow sensor to reduce even further the computational cost. The "simulated bee" was equipped with: (i) a minimalistic compound eye comprising 10 or 8 local motion sensors (depending of the tunnel configuration) measuring the optic flow magnitude, (ii) two optic flow regulators updating the control signals whenever specific optic flow criteria changed and (iii) three event-based controllers taking into account both error signals and visual events, each one in charge of its own translational dynamics. The "simulated bee" managed to travel safely along the tunnels without requiring any speed or distance measurements, using very low computational resources, by (i) concomitantly adjusting the side thrust, vertical lift and forward thrust only when both a visual contrast and a change of optic flow control error were detected, and (ii) avoiding collisions with the surface of the tunnels and decreasing or increasing its speed, depending on the clutter rate perceived by motion sensors.

I. INTRODUCTION

The performance of an event-based controller in terms of computational load is usually controversial since the eventdetector criterion is periodically calculated. In a classical event-based PID controller, the event detector requires many computations such as to calculate error and to compare its magnitude at each time step. The idea is to take benefit of event-based nature of the bio-inspired optic flow sensor to perform the autopilot computation only when it is needed. Insect-size Micro Aerial Vehicles (MAV) are increasingly becoming a reality ( [START_REF] Kushleyev | Towards a swarm of agile micro quadrotors[END_REF], [START_REF] Ma | Controlled flight of a biologically inspired, insect-scale robot[END_REF], [START_REF] Dunkley | Visual-inertial navigation for a camera-equipped 25g nano-quadrotor[END_REF], [START_REF] Floreano | Science, technology and the future of small autonomous drones[END_REF]) and have to be fitted with sensors and flight control systems enabling them to perform any kind of aerial maneuvers extremely close to obstacles both indoor and outdoor. Consequently, moving jointly very close to obstacles and at a high flight speed represents a big challenge for MAV applications. For the purpose of MAV applications with scarce resources, computationally-frugal methods would improve the biomimicry of the optic flowbased autopilot. Winged insects do not require any sensory system such as a GPS or an Inertial Measurement Units (IMU), nor any emissive proximity sensors such as ultrasonic or laser range finders. On the other hand, they rely heavily on Optic Flow (OF), which is the angular velocity of contrasted points, edges, or surfaces caused by the relative motion between the observer and the surrounding objects ( [START_REF] Gibson | The perception of the visual world[END_REF], [START_REF] Koenderink | Facts on optic flow[END_REF]). They make use of the OF to jointly avoid walls and adjust their speed in tunnel-like environments with either variable section ( [START_REF] Portelli | Honeybees' speed depends on dorsal as well as lateral, ventral and frontal optic flows[END_REF], [START_REF] Srinivasan | Visual control of navigation in insects and its relevance for robotics[END_REF]), an abrupt change in tunnel width ( [START_REF] Baird | Minimum viewing angle for visually guided ground speed control in bumblebees[END_REF]), a sudden change in unilateral optic flow by manipulating the wall texture from a checkerboard to axial stripes [START_REF] Linander | Bumblebees measure optic flow for position and speed control flexibly within the frontal visual field[END_REF]. We can therefore take inspiration from their OF-based strategies to develop smart MAVs autopilots ( [START_REF] Beyeler | Vision-based control of near-obstacle flight[END_REF], [START_REF] Shoemaker | Optic flow estimation on trajectories generated by bio-inspired closed-loop flight[END_REF], [START_REF] Srinivasan | Visual control of navigation in insects and its relevance for robotics[END_REF], [START_REF] Roubieu | A biomimetic vision-based hovercraft accounts for bees complex behaviour in various corridors[END_REF], [START_REF] Expert | Flying over uneven moving terrain based on optic-flow cues without any need for reference frames or accelerometers[END_REF]). Electrophysiological studies on flies motion sensitive neurons have shown that their OF sensitivity is texture-and contrastindependent as well as it can be explained by a "time-oftravel" scheme ( [START_REF] Franceschini | Directionally selective motion detection by insect neurons[END_REF]).

A B

Surge axis

Sway axis

Heave axis The output of "time-of-travel" scheme can be only updated when a change in the optical contrast with time is detected ( [START_REF] Franceschini | Directionally selective motion detection by insect neurons[END_REF], [START_REF] Vanhoutte | Time-of-travel methods for measuring optical flow on board a micro flying robot[END_REF]). An asynchronous signal is therefore generated in the output of a "time-of-travel" algorithm, which can be conveniently used in event-based control systems ( [START_REF] Durand | Further results on event-based PID controller[END_REF]). Because of the slight delay in the camera read-out, event-based collision avoidance algorithms have been only developed so far under open loop conditions in slightly cluttered corridors, but they were expected to be ready for use in the near future in closed loop on-board robotic platforms ( [START_REF] Clady | Asynchronous visual event-based time-to-contact[END_REF], [START_REF] Milde | Bioinspired event-driven collision avoidance algorithm based on optic flow[END_REF]). Event-based control systems have also been used for visual odometric purposes in the field of robotics ( [START_REF] Censi | Low-latency event-based visual odometry[END_REF]). Visual event: Contrast Detection Forward OF set-point Fig. 2. A visual event which consists of a contrast detection activates the event detector of each visual feedback loop. A visual feedback loop owns its OF set-point: a speed control loop (in blue) and a positioning control loop (in red). The event-based surge controller adjusts the pitch angle θ pitch as presented in Fig. 1A (that determines Vx via the bees' surge dynamics) on the basis of whichever sum of the two coplanar (horizontal or vertical) OFs measured is the larger. This value is compared with the forward OF set-point ω setF wd . The event-based surge controller adjusts the forward speed to minimize the error F wd . The event-based positioning controller adjusts the roll angle θ roll (or the wing stroke amplitude ∆Φ), which determines the distances to the walls, the ground, or the roof, depending on the sway (or heave) dynamics on the basis of whichever of the measured OFs is the largest. The latter value is compared with the positioning OF set-point ω setP os . The direction of avoidance is given by a Control direction Selector that multiplies the control signal by a direction factor depending on the direction of the maximum OF signal (see in Fig. 3).

T LMS 1 (45°) = 2.82°S L LMS 4 (-45°) = 2.82°L MS 6 (-135°) = 2.82°L MS 3 (135°) = 2.82°L MS 2 (90°) = 4°L MS 5 (-90°) = 4°L MS 8 (-90°) = 4°L MS 7 (-45°) = 2.82°L MS 10 (90°) = 4°L MS 9 (135°) = 2.82°T
In the present paper, we describe in detail an eventbased version of the ALIS (ALIS stands for "AutopiLot using an Insect-based vision System") which was originally designed with fixed-step digital controllers to perform tunnelfollowing tasks ( [START_REF] Portelli | Modelling honeybee visual guidance in a 3-D environment[END_REF]). Preliminary results were presented in a work-in-progress paper ( [START_REF] Serres | Event-based visual guidance inspired by honeybees in a 3d tapered tunnel[END_REF]), however, in this new version, the event-based ALIS autopilot is described in detail including new simulation results showing its ability to fit with bees' behavior and bees' compound eye by combining both visual and control events.

In section 2, the simulation set-up and the simulated agent mimicking a honeybee fitted with a minimalistic visual system including only 10 local OF measurements are described. In section 3, the ALIS autopilot and the event-based discrete time controllers are described in detail. In section 4, we describe the trajectories of the simulated agent fitted with the event-based ALIS autopilot in a doubly tapering tunnel (with a tapering angle of 7 • in both planes) or in straight tunnel with an abrupt change in width. The results of the simulations performed show that the event-based ALIS autopilot generates a small number of events making the "simulated bee" travel along tunnels endowed with crosssection changes, thus greatly decreasing the number of control signal computations required while mimicking bees' behavior during tunnel-following task.

II. SIMULATION SET-UP

A. Software set-up

In-silico experiments were carried out on two standard computers running: (i) a first one equipped with the MATLAB/Simulink c software program at a sampling rate of 1 kHz for OF sensing purposes and feeding the eventbased control systems, and (ii) a second one equipped with the MORSE software program1 (MORSE stands for "Modular OpenRobots Simulation Engine") to simulate the visual interactions with the virtual 3D environment. The two computers were linked together by a TCP/IP protocol in which the maximum delay was limited to 3 ms. The virtual 3D environment was built with the Blender software program 2 . The four surfaces of the tunnel were lined with high resolution photographs of natural panoramic scenes ( [START_REF] Brinkworth | Biomimetic motion detection[END_REF]). The 20 scalar photoreceptor signals were transmitted via TCP/IP to a second computer processing the 10 local OF computations (Fig. 1B-C), the 3 event-based controllers and the 3 translational dynamics of the "simulated bee" using the Matlab/Simulink software program (Fig. 1A). The Matlab/Simulink program then sent the 3D position (x, y, z) of the "simulated bee" back to the MORSE simulator.

B. Optic flow

The OF vector field perceived by a simulated agent depends in particular on the structure of the 3D environment ( [START_REF] Gibson | The perception of the visual world[END_REF], [START_REF] Koenderink | Facts on optic flow[END_REF]). The OF can be defined by a vector field of the apparent motion of objects, surfaces, and edges in a visual scene generated by the relative motion between an agent and the scene. The OF field ω (1) is a combination between two components: a translational ω T and a rotational OF ω R ( [START_REF] Koenderink | Facts on optic flow[END_REF]).

ω = ω T + ω R (1)
It has been established that hymenopterans stabilize their gaze by compensating for any body rotations ( [START_REF] Viollet | Feed-forward and visual feedback control of head roll orientation in wasps (polistes humilis, vespidae, hymenoptera)[END_REF]). It is assumed that the head of the artificial flying insect is perfectly stabilized canceling all rotation due to body pitch and roll with respect the inertial frame and hence stabilized along the tunnel axis. Therefore, the head of the artificial flying insect experienced only translational OF and each OF sensor will receive a purely translational OF ( ω R = 0). The translational OF (expressed in rad/s) can be defined as follows:

ω T = - V -( V • d) • d D ( 2
)
where d is a unit vector describing the viewing direction, V is the translational velocity vector, and D is the distance from the object seen by an OF sensor.

C. Simulated flying agent

The behaviour of the "simulated bee" introduced into the in-silico experiments was based on findings which may explain how a flying honeybee controls its speed and avoids walls in tunnel-like environments ([9], [START_REF] Portelli | Honeybees' speed depends on dorsal as well as lateral, ventral and frontal optic flows[END_REF], [START_REF] Srinivasan | Visual control of navigation in insects and its relevance for robotics[END_REF]). Each translational axis of the bee's dynamic model was uncoupled, as occurs in a quadrotor (Fig. 1A). Each of the bee's translational dynamics can be defined by a first order transfer function with a time constant of 0.22 s (see [START_REF] Portelli | Modelling honeybee visual guidance in a 3-D environment[END_REF] for details).

Pitch angle, roll angle, and wing stroke amplitude were bounded on the basis of data previously published on bees ( [START_REF] Portelli | Modelling honeybee visual guidance in a 3-D environment[END_REF]). In this indoor study, the simulated agent was not subjected to any wind and could only move in translation.

The present agent was equipped with a set of 10 OF sensors for the doubly tapered tunnel (LMS 1 to LMS 10 in Fig. 1) and with a set of 8 OF sensors (LMS 1 to LMS 8 in Fig. 1) for the tunnel with an abrupt change in width without ceiling: each of these sensors consisted of just two photoreceptors driving a Local Motion Sensor (LMS) based on a "time-of-travel" scheme ( [START_REF] Franceschini | Directionally selective motion detection by insect neurons[END_REF], [START_REF] Vanhoutte | Time-of-travel methods for measuring optical flow on board a micro flying robot[END_REF]). The visual axes of the two adjacent photoreceptors were assumed to be separated by the angle ∆ϕ, and each photoreceptor's angular sensitivity was assumed to be a Gaussoid function with an angular width at half height of ∆ρ = ∆ϕ. The inter-receptor angles ∆ϕ were therefore smaller in the frontal or rear parts of the visual system (∆ϕ 45 • = ∆ϕ • × sin(45 • ) = 2.82 • ), than in the lateral, ventral, or dorsal parts (∆ϕ 90 • = 4 • ) as observed in the bee's compound eye ( [START_REF] Seidl | Die Sehfelder und Ommatidien-Divergenzwinkel von Arbeiterin, Königin und Drohn der Honigbiene (Apis mellifica)[END_REF]). A sine-law gradient can be used either in the horizontal or vertical planes to compute the corresponding inter-receptor angle ∆ϕ. 6 OF sensors were located in the horizontal plane oriented at azimuthal angles ϕ of ±45 • , ±90 • , and ±135 • (Fig. 1A), and the other 4 were oriented in the vertical plane at elevation angles θ of -45 • , ±90 • and +135 • (Fig. 1B).

III. VISUAL EVENT-BASED DISCRETE TIME

CONTROLLERS

In this section, we describe how we adapted the eventbased control approach introduced in [START_REF] Durand | Further results on event-based PID controller[END_REF] for visual eventbased OF regulation purposes. As presented in the schematic of the event-based autopilot in Fig. 2, three control input signals were computed when a visual event occurs, corresponding to the three uncoupled translational degrees of freedom on the surge (x), sway (y) and heave axes (z). We consider in this work (i) as a "visual event", an event-related 3. The control direction selector automatically selects the tunnel surface to be followed (wall, ground or roof) by multiplying the control signal (the output from the Positioning controller) by a direction factor u(3) that depends on the direction of the largest OF signal. Note that the sway and heave dynamics can be driven alternately, depending on which OF (side or vertical) is the largest at any given time. Adapted from [START_REF] Portelli | Modelling honeybee visual guidance in a 3-D environment[END_REF]. condition based on the visual contrast detection inside the bio-inspired optic flow sensors and (ii) as a "control event", an event-related condition based on optic-flow error signals.

Lead controllers are introduced into the heave and sway control systems to improve their stability. The natural integrator that relates the simulated agent's position to its speed makes both the heave and sway control systems reach a zero steady state error. However, a PI controller is introduced into the surge control system because OF is proportional to speed, but inversely proportional to distance (see [START_REF] Ma | Controlled flight of a biologically inspired, insect-scale robot[END_REF]). Consequently, we have to add an integral action to cancel the steady state error.

A. Visual event-based PI controller on the surge axis

A discrete time PI controller was designed. The modelling of the proportional part was quite straightforward, and the backward difference approximation method was used to model the integral part. The resulting code is: The event-based controller is tuned with the values K p = 20.10 -3 , K i = 6.10 -2 and T s = 1ms.

if newVisualEvent ==

B. Visual event-based PD controllers on the sway and the heave axes

Two discrete time PD controllers were also designed to deal with the sway (y) and heave axes (z). Modelling the proportional part was quite straightforward, and the backward difference approximation method was used to model the derivative part. The following code was thus obtained for each controller: The event-based controllers on the sway axis were tuned with the values K p = 20.10 -4 , T d = 0.1 s, and on the heave axis, they were tuned with the values K p = 20.10 -3 , T d = 0.1 s. The input sign was on the basis of Fig. 3 adapted from [START_REF] Portelli | Modelling honeybee visual guidance in a 3-D environment[END_REF]. The parameter elim, was chosen in order to significantly reduce the control updates, while making it possible to perform reference tracking of the OF set-points.

IV. RESULTS

Frontal and fronto-lateral spatial configurations of the optic flow sensors (±30 • and ±45 • , respectively) were tested in the following simulated experiments.

A. Doubly tapered tunnel

In this simulation, as presented in Fig. 2, the following parameters are chosen such as:

k 1 ω 1 = max ω 45 • L sin 2 (45 • ) , ω 90 • L sin 2 (90 • ) , ω 135 • L sin 2 (135 • ) (3)
k 2 ω 2 = max ω -45 • R sin 2 (45 • ) , ω -90 • R sin 2 (90 • ) , ω -135 • R sin 2 (135 • ) (4)
k 3 ω 3 = max ω -45 • V sin 2 (45 • ) , ω -90 • V sin 2 (90 • ) (5)
k 4 ω 4 = max ω 90 • D sin 2 (90 • ) , ω 135 • D sin 2 (135 • ) (6)
Each OF output signal ω i (i = 1, 2, 3, 4) is weighted as explained in [START_REF] Roubieu | A biomimetic vision-based hovercraft accounts for bees complex behaviour in various corridors[END_REF]. The subscripts in ω L , ω R , ω V and ω D stand for the OF output signal in the left, right, ventral and dorsal directions. In figure 4, the simulated environment was a doubly tapering tunnel 6m long, 1m wide, 1m high with a slope angle of 7 • . Figure 4A shows a perspective view of the tunnel. The simulated bee entered the tunnel at zero speed, with the initial coordinates x 0 = 0.1m and various couples/pairs of y 0 and z 0 (Fig. 4B). Within the first 800 ms, the simulated bee was controlled in the open loop mode to enable all the sensors measure the OFs. From the time 800 ms onwards, the OF was regulated in the closed loop mode via the event-based controllers. Figure 4C shows three trajectories in the vertical plane (x, z) and Fig. 4D shows these trajectories in the horizontal plane (x, y), plotted every 300ms. The simulated bee can be seen to have gradually increased both its ground clearance (Fig. 4C) and its right and left clearance (Fig. 4D), while the forward speed (Fig. 4E) increased automatically up to 1.5 m/s. Figure 5 shows the control input signals with a starting point x 0 = 0.1m, y 0 = 0.5m and z 0 = 0.3m. The control input signals composed of the pitch, roll and wing stroke angles are also plotted in Fig. 5A, 5B and 5C, respectively. During the horizon of simulation of six seconds at the sampling time T s = 1ms, as shown in Fig. 5D, the computation of the condition if (abs(e-e old)>elim) is activated only 799 times by the visual events. Therefore, it can be seen that the number of control updates is greatly reduced to 154, 83 and 10, for the pitch, the roll and the wing stroke control signals, respectively. This be compared with fixed-step PID controllers computed at 1kHz. Figure 6A plots the OF with respect to the set-point of 540 • /s driving the surge dynamics. In figure 6B, the OF is plotted with respect the set-point of 315 • /s driving the sway and heave dynamics. Both systems of OF regulation gave satisfactory results.

A B

Fig. 6. Simulated systems of OF regulation in the case of a 3D trajectory starting at the initial coordinates x 0 = 0.1m; y 0 = 0.5m; z 0 = 0.3m. (A) Forward feedback signal determined by taking the largest sum of two diametrically opposed OF sensors (horizontal OF sensors, vertical OF sensors). (B) Positioning feedback signal determined by taking the largest output from the OF sensors (right OF sensors; left OF sensors, ventral OF sensors, dorsal OF sensors).

B. Tunnel with an abrupt change in width

In figure 7A, the simulated environment was a tunnel with an abrupt change in width whose dimensions are 12m long, 1m wide and 1m high. The abrupt change is located at midway. Within the first 1.5s, the simulated bee was controlled in the open loop mode to enable all the sensors Fig. 7.

Comparison between two spatial configurations of optic flow sensors: frontal configuration (±30 • , in red) and fronto-lateral configuration (±45 • , in black) (A) Perspective view of the tunnel with an abrupt change in width. (B) Two simulated bees 3D trajectories starting at the same initial coordinates x 0 = 0.1m; y 0 = 0.7m; z 0 = 0.2m plotted every 300 ms. In black, 6 OF sensors were located in the horizontal plane oriented at azimuthal angles of ±45 • , ±90 • , and ±135 • (Fig. 1B), and the other 4 were oriented in the vertical plane at elevation angles of -45 • , ±90 • and +135 • (Fig. 1C). In red, the sensors are located in the horizontal plane oriented at azimuthal angles of ±30 measure the OF. From the time 1.5s onwards, the OF was regulated in the closed loop mode combining both visual and control events. Figure 7B shows two simulated bees 3D trajectories starting at the same initial coordinates x 0 = 0.1m; y 0 = 0.7m; z 0 = 0.2m plotted every 300 ms for two different sets of parameters, in black using (3)-( 5), in red using ( 7)- [START_REF] Baird | Minimum viewing angle for visually guided ground speed control in bumblebees[END_REF].

k 1 ω 1 = max ω 30 • L sin 2 (30 • ) , ω 90 • L sin 2 (90 • ) , ω 150 • L sin 2 (150 • ) (7)
k 2 ω 2 = max ω -30 • R sin 2 (30 • ) , ω -90 • R sin 2 (90 • ) , ω -150 • R sin 2 (150 • ) (8)
k 3 ω 3 = max ω -30 • V sin 2 (30 • ) , ω -90 • V sin 2 (90 • ) (9)
k 4 ω 4 = 0 since the dorsal OF is equal to zero (the tunnel has no ceiling). While figures 7C and 7D exhibit similar results in the vertical and horizontal planes, Fig. 7E shows a decrease in the forward speed V x before the abrupt change with azimuthal angles equal to ±30 • . The constriction point is anticipated sooner than for azimuthal angles equal to ±45 • . Moreover one can note that the forward control loop is much more faster than the position control loop as depicted in Figs. 7D and7E. Consequently, in the future simulations, the parameters given by ( 7)-( 9) and k 4 ω 4 = 0 are chosen. Figure 8C shows 3 trajectories in the vertical plane (x, z) and Fig. 8D shows these trajectories in the horizontal plane (x, y), plotted every 300 ms. The simulated bee can be seen to have gradually increased the forward speed (Fig. 8E) up to 2m/s and decreased the forward speed as soon as the abrupt change was detected thanks to the OF sensors oriented at the azimuthal angles ϕ of ±30 • . Figure 9 shows the control input signals, the control updates and the visual events with x 0 = 0.1m, y 0 = 0.7m and z 0 = 0.2m. One can note that the computation of the condition on the error leading the eventbased controllers is activated 1030 times by the visual events (see Fig. 9D). It results low numbers of control updates, 236 for the pitch, 130 for the roll and 25 for the wing stroke angles. Figure 10A plots the OF with respect to the setpoint of 540 • /s driving the surge dynamics. In figure 10B, the OF is plotted with respect the set-point of 315 • /s driving the sway and heave dynamics. The event-based controllers cancel the disturbance introduced by the abrupt change in tunnel width. Using OF measurement coming from 30 • to 150 • to control the speed of the "simulated bee" allows it to reduce its speed before the constriction point (blue shaded zone in Fig. 8E), as bees did in similar tunnel configuration ( [START_REF] Baird | Minimum viewing angle for visually guided ground speed control in bumblebees[END_REF], [START_REF] Linander | Bumblebees measure optic flow for position and speed control flexibly within the frontal visual field[END_REF]).

V. CONCLUSIONS

In this paper, we present an innovative visual event-based autopilot mimicking the direct Optic Flow (OF) feedback A B Fig. 10. Simulated systems of OF regulation in the case of a 3D trajectory starting at the initial coordinates x 0 = 0.1m; y 0 = 0.7m; z 0 = 0.2m control loops which may occur in honeybees. Likewise, honeybees tiny brains contain smart OF sensors, called Velocity-Tuned motion-sensitive descending neurons (or VT neurons) which are spiking neurons sensitive to the OF ( [START_REF] Ibbotson | Evidence for velocity-tuned motion-sensitive descending neurons in the honeybee[END_REF]): each spike can be taken to constitute a visual event occurring in the insects neuronal system. Our proposed visual event-based control system combines both error signals and visual events. In a near future, it will be likely to become a major tool which can be used to understand more clearly how insects process the event-based visual information they receive in order to react appropriately: this information ranges contrasts and local motion to wing movements generating insects 3D movements in any tunnel configuration. This is the first time we tended to explain using a simulation that a decrease in the forward speed before the abrupt change in tunnel width can be observed only by considering OF measurements coming from frontal azimuthal angles ±30 • instead of fronto-lateral azimuthal angles ±45 • (Fig. 7E). These results are particularly consistent with bumblebees behavior in similar situations ( [START_REF] Baird | Minimum viewing angle for visually guided ground speed control in bumblebees[END_REF], [START_REF] Linander | Bumblebees measure optic flow for position and speed control flexibly within the frontal visual field[END_REF]). However, no difference in the simulation is observed in the lateral displacement with both spatial configurations (Fig. 7D) before reaching the constriction point, which is also similarly observed in bumblebees ( [START_REF] Linander | Bumblebees measure optic flow for position and speed control flexibly within the frontal visual field[END_REF]). The bumblebees behavior across an abrupt change in tunnel width may be therefore explained jointly by their dynamic capabilities and their viewing angles. The results of the present in-silico experiments show that the visual event-based control system based on OF consisting in updating the control signals solely when visual contrasts are detected provides a suitable basis for traveling along unknown tunnels endowed with cross-section changes. It is now proposed to perform further studies involving a more sophisticated OF-based autopilot traveling along more complex tunnels, rooms, or mazes, as well as to implement the visual event-based ALIS autopilot in the near future onboard a 25-gram nano-quadrotor ( [START_REF] Dunkley | Visual-inertial navigation for a camera-equipped 25g nano-quadrotor[END_REF]).

Fig. 1 .

 1 Fig. 1. (A) Top view describing the orientation of 6 Local Motion Sensors (LMSs) in the horizontal plane corresponding to 12 simulated photoreceptors. (B) Side view showing the orientation of 4 LMSs in the vertical plane corresponding to 8 simulated photoreceptors.

  Fig.3. The control direction selector automatically selects the tunnel surface to be followed (wall, ground or roof) by multiplying the control signal (the output from the Positioning controller) by a direction factor u(3) that depends on the direction of the largest OF signal. Note that the sway and heave dynamics can be driven alternately, depending on which OF (side or vertical) is the largest at any given time. Adapted from[START_REF] Portelli | Modelling honeybee visual guidance in a 3-D environment[END_REF].

Fig. 4 .Fig. 5 .

 45 Fig. 4. (A)Perspective view of the doubly tapering tunnel. (B) Simulated bees 3D trajectories starting at the initial coordinates x 0 = 0.1m; y 0 = 0.5m; z 0 = 0.5m (in black), x 0 = 0.1m; y 0 = 0.4m; z 0 = 0.5m (in blue) and x 0 = 0.1m; y 0 = 0.5m; z 0 = 0.3m (in red) plotted every 300ms, (C) Trajectories in the vertical plane (x, z), and in (D) in the horizontal plane (x, y), plotted every 300 ms. (E) Forward speed Vx profiles.

  Fig. 7. Comparison between two spatial configurations of optic flow sensors: frontal configuration (±30 • , in red) and fronto-lateral configuration (±45 • , in black) (A) Perspective view of the tunnel with an abrupt change in width. (B) Two simulated bees 3D trajectories starting at the same initial coordinates x 0 = 0.1m; y 0 = 0.7m; z 0 = 0.2m plotted every 300 ms. In black, 6 OF sensors were located in the horizontal plane oriented at azimuthal angles of ±45 • , ±90 • , and ±135 • (Fig. 1B), and the other 4 were oriented in the vertical plane at elevation angles of -45 • , ±90 • and +135 • (Fig. 1C). In red, the sensors are located in the horizontal plane oriented at azimuthal angles of ±30 • , ±90 • , and ±150 • , and in the vertical plane at elevation angles of -30 • , ±90 • and +150 • . (C) Trajectories in the vertical plane (x, z), and in (D) in the horizontal plane (x, y). (E) Forward speed Vx profiles.

Fig. 8 .

 8 Fig. 8. Comparison between three sets of initial conditions six sensors located in the plane oriented at azimuthal angles of ±30 • , ±90 • , and ±150 • , and the other 4 were oriented in the vertical plane at elevation angles of -30 • , ±90 • and +150 • . (A) Perspective view of the tunnel with an abrupt change in width. (B) Simulated bees 3D trajectories starting at the initial coordinates x 0 = 0.1m; y 0 = 0.7m; z 0 = 0.5m (in black), x 0 = 0.1m; y 0 = 0.4m; z 0 = 0.5m (in blue) and x 0 = 0.1m; y 0 = 0.7m; z 0 = 0.2m (in red) plotted every 300 ms. (C) Trajectories in the vertical plane (x, z), and in (D) in the horizontal plane (x, y), plotted every 300 ms. (E) Forward speed Vx profiles.

Fig. 9 .

 9 Fig. 9. Simulated control signals, control updates and visual events during a 3D trajectory starting at the initial coordinates x 0 = 0.1m; y 0 = 0.7m; z 0 = 0.2m. (A) The pitch angle is updated 236 times, (B) The roll angle is updated 130 times, (C) The wing stroke angle is updated 25 times, (D) 1030 constrast detections stand for the visual events over the simulation in the tunnel with an abrupt change in width.

19 May

2015: MORSE release 1.3 under GNU GPLv2, https://www.openrobots.org/wiki/morse/ 2 15 Nov. 2015: Blender release 2.76 under GNU GPLv2, https://www.blender.org/

ACKNOWLEDGMENT

We thank R. Brinkworth and D. O'Carroll (Adelaide Uni.) for kindly making their panoramic images available to us.

This work was supported partly by CNRS (Life Science;


Pramod S Kaushik

Maxime Carrere

Frédéric Alexandre

Surampudi Bapi Raju

A biologically inspired neuronal model of reward prediction error computation

The neurocomputational model described here proposes that two dimensions involved in computation of reward prediction errors i.e magnitude and time could be computed separately and later combined unlike traditional reinforcement learning models. The model is built on biological evidences and is able to reproduce various aspects of classical conditioning, namely, the progressive cancellation of the predicted reward, the predictive firing from conditioned stimuli, and delineation of early rewards by showing firing for sooner early rewards and not for early rewards that occur with a longer latency in accordance with biological data.

I. INTRODUCTION

One of the main motivations to design neural networks is to get adaptive functions mimicking natural learning. In addition, inspiration from biology can be deep and exploit the paradigm of neural computation to propose also a model of the underlying brain circuitry. One of the earliest attempts to understand how animals learn involved pairing an unconditioned stimulus (US) with a cue or conditioned stimulus (CS) and observing that animals start responding to the CS after some point in time [START_REF] Pavlov | Conditioned Reflexes (V.Anrep, trans.)[END_REF]. This is the basis of pavlovian learning, a fundamental learning mechanism in animals, which has been addressed by several models of neural networks [START_REF] Schmajuk | Stimulus configuration, classical conditioning and the hippocampus[END_REF], [START_REF] Le Pelley | The role of associative history in models of associative learning: a selective review and a hybrid model[END_REF]. The model described here focuses on the mechanism of reward prediction error within pavlovian learning and does not deal with other conditioning phenomena.

Concerning the link to the brain architecture, when neurons in a cerebral structure called VTA were recorded on a similar procedure in primates [START_REF] Schultz | A neural substrate of prediction and reward[END_REF], dopaminergic neurons in VTA were found to shift their firing from the US to the CS allowing a link to the Temporal Difference (TD) learning algorithm put forward by Sutton and Barto as a possible mechanism that animals might use to learn [START_REF] Sutton | Reinforcement Learning: An introduction[END_REF]. In short, the algorithm says that the dopaminergic firing that the US causes is an error signal that the brain uses for learning. It is proposed to correspond to the Reward Prediction Error (RPE), comparing predicted and actual rewards. This error signal flows back in a recursive manner from the US to the CS, canceling the peak of dopamine at the time of the US and creating one at the time of the CS.

The TD framework has a high explanatory value while remaining simple. But, despite its usefulness, it remains a rather high level model that does not precisely account for the knowledge accumulated on the brain mechanisms associated to pavlovian learning.

The RPE computation involves a high number of neural structures which makes explaining it in a biologically plausible manner difficult. Still, considerable progress has been achieved [START_REF] Brown | How the basal ganglia use parallel excitatory and inhibitory learning pathways to selectively respond to unexpected rewarding cues[END_REF], [START_REF] Suri | Temporal difference model reproduces anticipatory neural activity[END_REF], [START_REF] O'reilly | Pvlv: the primary value and learned value pavlovian learning algorithm[END_REF], [START_REF] Tan | A dopamine-acetylcholine cascade: simulating learned and lesion-induced behavior of striatal cholinergic interneurons[END_REF], [START_REF] Hazy | Neural mechanisms of acquired phasic dopamine responses in learning[END_REF], [START_REF] Vitay | Timing and expectation of reward: a neurocomputational model of the afferents to the ventral tegmental area[END_REF].

The available models possess certain mechanisms in common. One of them is the dual pathway mechanism proposing that distinct circuits cancel the peak of dopamine at the US and create another one at the CS. The models usually differ in the structures implicated and the origin of the timing signals. Some vary by the signal that cuts off the arrival of the reward: the O'Reilly PVLV model [START_REF] O'reilly | Pvlv: the primary value and learned value pavlovian learning algorithm[END_REF] uses a ramping expectation similar to the one in the TD algorithm while the Vitay model [START_REF] Vitay | Timing and expectation of reward: a neurocomputational model of the afferents to the ventral tegmental area[END_REF] proposes an oscillatory mechanism which peaks at the expected time through the ventral striatum (VS), inhibiting the reward signal. Analyzing the performances as well as the properties of these models can be interesting to decipher cerebral mechanisms but also to develop more powerful algorithms in Machine Learning.

To acquire a deeper understanding of the roles of the cerebral structures involved in the RPE computation, we need a neurocomputational model implicating them more faithfully in Pavlovian learning, The approach described here is a model that proposes a circuit involving a dissociation between magnitude expectation and timing expectation and explaining more precisely how the computation of the reward prediction error happens inside the VTA.

Our model dissociates the processing of reward magnitude and reward timing and delegates it to two structures of the Medial Temporal Lobe, the Basolateral Amygdala (BLA) and the Ventral Striatum (VS) respectively. Unlike other models, our expectation signal is based on VTA GABA neurons that ramps at the expected time to inhibit the US dopamine reward signal. These VTA GABA neurons receive their input from one of the sub-populations in the Peduncolopontine Nucleus (PPN) called PPN FT(Fixation Target), which exhibits persistent activity between the CS and US [START_REF] Okada | Different pedunculopontine tegmental neurons signal predicted and actual task rewards[END_REF].

II. MODEL OVERVIEW

The model attempts to explain how the dopamine reward prediction error is computed in appetitive conditioning in the VTA. The model is shown in Figure 1. The functioning of the model can be explained in four phases of functioning, described in the following paragraphs, together with references to the main biological evidences supporting them:

A. US Firing and Learning

When reward is delivered, it is reported to fire the Lateral Hypothalamus (LH) and activate the LH → PPN RD (Reward delivery) → VTA Dopamine pathway resulting in US dopamine firing prior to any sort of learning [START_REF] Semba | Afferent connections of the laterodorsal and the pedunculopontine tegmental nuclei in the rat: A retro-and antero-grade transport and immunohistochemical study[END_REF] [START_REF] Lokwan | Stimulation of the pedunculopontine tegmental nucleus in the rat produces burst firing in a9 dopaminergic neurons[END_REF]. A pathway from LH to BLA ensures that BLA firing for CS has the exact amplitude as the US firing [START_REF] Sah | The amygdaloid complex: anatomy and physiology[END_REF]. The VTA US dopamine firing alerts the BLA to recognize there is a reward and it progressively learns to associate the CS (reaching BLA through the inferotemporal cortex, IT) to the US gated by the US VTA dopamine [START_REF] Bissière | Dopamine gates ltp induction in lateral amygdala by suppressing feedforward inhibition[END_REF]. There is concurrent learning in the VS for the time duration of the cue and reward delivery [START_REF] Takahashi | Temporal specificity of reward prediction errors signaled by putative dopamine neurons in rat vta depends on ventral striatum[END_REF].

B. CS Firing

Another nucleus of the Amygdala, the Central Nucleus (CE) appears to get activated during this learning by the BLA [START_REF] Ledoux | Emotion circuits in the brain[END_REF], enabling the VTA dopamine to undergo phasic bursts of the same amplitude at the presentation of the CS through the IT → BLA → CE → PPN RD → VTA Dopamine pathway [START_REF] Cheng | Organization of corticostriatal and corticoamygdalar projections arising from the anterior inferotemporal area te of the macaque monkey: a phaseolus vulgaris leucoagglutinin study[END_REF] at the end of conditioning. During conditioning, part of the reward magnitude is learnt resulting in partial conditioning causing partial cancellation of the US and also showing a partial firing at the arrival of the CS like other dual pathway models.

C. Expectation

The expectation signal is what ultimately cancels the predicted reward enabling the dopamine to fire its background rate at the time of reward delivery. The magnitude and timing of the reward are handled by BLA and VS respectively. The presentation of the CS, which reportedly fires the IT and thereby the Orbitofrontal Cortex (OFC) [START_REF] Carmichael | Sensory and premotor connections of the orbital and medial prefrontal cortex of macaque monkeys[END_REF], activates the VS and its neurons learn the interval timing and it acts similar to a negative integrator and progressively lowers the inhibition that VS exerts on PPN FT, as reward delivery is approached. Thus, it conveys the precise time where the PPN FT can increase the inhibition through VTA GABA and cancel the dopamine.

The magnitude of expectation originates from the CS firing in the Central Amygdala (CE) and maintained in the PPN FT through a self sustaining mechanism [START_REF] Okada | Reward prediction-related increases and decreases in tonic neuronal activity of the pedunculopontine tegmental nucleus[END_REF] [START_REF] Kobayashi | Reward prediction error computation in the pedunculopontine tegmental nucleus neurons[END_REF]. The GABA firing in the VTA is reflective of this [START_REF] Yau | Pontomesencephalic tegmental afferents to vta non-dopamine neurons are necessary for appetitive pavlovian learning[END_REF] and the PPN FT integrates the magnitude from the Central Amygdala (CE) and timing information from the VS to achieve the ramping signal that encodes both time and magnitude of the reward delivery.

D. Early reward

It has been reported in the literature [START_REF] Schultz | A neural substrate of prediction and reward[END_REF] as a hallmark of reinforcement learning in VTA that an omitted reward causes a dip of dopamine at the time of the expected reward, which can be easily explained by the dual pathway mechanism, dissociating the creation of a dopaminergic peak at the time of the CS and the cancellation of the peak at the time of the US. More difficult to explain is the fact that an early reward does not cause dips at the time of the expected US. Here it is explained using a different mechanism due to the sustained nature of the expectation signal. It posits that there is an inhibition from PPN RD to PPN FT and the early reward that flows through PPN RD resets the expectation of PPN FT.

III. MODEL DESCRIPTION A. Computational principles

The proposed model is composed of computational units where each unit represents a population and computes the mean activity of the population. A time-dependent firing rate describes the dynamics across time for each population. V (t) represents the membrane potential of the unit and the firing rate is a positive scalar of V (t) given by U (t). Each unit is represented by the following equations:

τ. dV (t) dt = (-V (t) + g exc (t) -g inh (t) + B + η(t)) (1)
U (t) = (V (t)) + ( 2
)
where τ is the time constant of the cell, B is the baseline firing rate and η(t) is the additive noise term chosen randomly at each time step from an uniform distribution between -0.01 and 0.01. The incoming afferent currents g exc and g inh represent the weighted sum of excitatory and inhibitory firing rates respectively, the weight representing the synaptic weights between the populations. Some of the populations require an incoming tonic component converted to a short phasic transformation. This is done by the following equations

τ. dx(t) dt = (-x(t) + x(t)) (3)
φ τ,k (x(t)) = (x(t) -k.x(t)) + (4)
where x(t) integrates the incoming input x(t) with a time constant τ , while φ τ,k (x(t)) represents the positive part of the difference between x(t) and x(t). The constant k is a constant that controls how much of the tonic component is kept, a k value of 0 indicates the entire tonic component to be preserved and a k value of 1 outputs the entire phasic component from the tonic input.

A Bound function is used when the firing of a population is described with an upper and a lower limit in certain populations

ψ(x) =    0 if x < 0 x if 0 < x < 1 1 if x > 1 (5)
There is also a threshold function used in some populations and it outputs 1 when the input exceeds a threshold Γ, 0 otherwise:

∆ Γ (x) = 0 if x < Γ 1 otherwise (6)
The learning rule defined in the model is based on the Hebbian learning rule. The evolution over time of the weight w(t) of a synapse between the neuronal population pre (presynaptic neuron) and the neuronal population post (postsynaptic neuron) is governed by:

dw(t) dt = (α.U pre (t).U post (t)) ( 7
)
where w is the weight term, α the learning rate and U (t) is indicating the firing rate of the presynaptic and postsynaptic neuronal populations.

B. Population definitions

1) Representations of inputs: IT and LH are the populations used for inputs for the CS and the US respectively and they are represented simplistically by a square wave signal as given below:

U (t) = I(t) + (8)
where I(t) is an external input resulting either from a stimulus or from a reward.

2) Basolateral Amygdala: The BLA learns to associate the CS with the US and provides the magnitude expectation that eventually cancels the US dopamine. The BLA receives inputs from the IT.

τ. dV (t) dt = (-V (t) + φ τ exc,k (g exc (t)) + η(t)) (9)
U (t) = (V (t)) + (2)
with τ = 10ms, τ exc = 10ms, k= 1.

The CS is learnt by updating the synaptic weights between IT and BLA and the learning rule is given by:

dw(t) dt = D.α.U pre (t).(U mag -U post (t)) + (10)
where D indicates the presence of the US corresponding to the dopaminergic neuronal modulation from the VTA, α is the learning rate equal to 0.003, U mag is the magnitude of LH firing, U pre and U post are the firing rates of presynaptic and postsynpatic neurons respectively.

3) Central Amygdala: The CE is the output nuclei of the amygdala in this model and it projects to both the PPN nuclei, relaying information from the BLA. The CE projects to the PPN RD neurons that convey US and CS firing to the VTA dopamine neurons and PPN FT neurons that convey expectation.

The equations for the membrane potential and the firing rate are the same as Equation 9and Equation 2 respectively, with τ = 20ms, τ exc = 5ms, k = 1.

4) Peduncolopontine nucleus:

The PPN has two distinct populations in this model for reward and expectation. a) PPN RD: The PPN Reward Delivery neurons signal occurrence of the CS and the US from the CE and the LH respectively. It also has a sub-population of inhibitory neurons that inhibit the PPN FT expectation neurons.

The equations for the membrane potential and the firing rate are the same as Equation 9and Equation 2 respectively, with τ = 5ms, τ exc = 5ms, k = 1.

b) PPN FT:

The PPN FT neurons encode the expectation and are subdivided into two populations, one holding the magnitude and the other delivering the expectation to the VTA GABA neurons.

PPN FT Magnitude: The PPN FT Magnitude neurons receive information from the CE and are inhibited by the PPN RD neurons. They serve to maintain a constant magnitude that is conveyed to the other population of PPN FT neurons (PPN FT Relay).

τ. dV (t) dt = (-V (t) + (g exc (t)) -g inh (t) + η(t)) (11)
U (t) = (V (t)) + (2)
with τ = 5ms.

PPN FT Relay: The PPN FT Relay population receives information from the PPN FT Magnitude population and is inhibited by the VS that conveys the timing signal and the output of these neurons is passed to the VTA GABA neurons enabling final cancellation.

The equations for the membrane potential and the firing rate are the same as Equation 11and Equation 2 respectively, with τ = 5ms.

5) Ventral Striatum and OFC:

The Ventral Striatum handles the timing by reducing its inhibition at the required moment of reward delivery thereby conveying the precise moment to cancel the predicted reward. The OFC in this case, indicates the presence of the stimulus and has an excitatory effect on the Ventral Striatum. The timing model of VS described here is a simplified timing model comprising a negative integrator similar to the timing algorithm in [START_REF] Rivest | Adaptive drift-diffusion process to learn time intervals[END_REF]. The integrator here has an amplitude of 1 at the beginning of the trial and after weight updating, decreases its firing to 0 at the precise time of reward delivery adjusting its slope.

Mechanism of timing: The timing mechanism described here is an abstract method describing the time for a fixed interval with the weights encoding the duration of the interval.

τ. dV (t) dt = (g exc (t)) -V.∆ Γ (φ τ mod,k (g mod (t)) -B) + η(t)) (12)
U (t) = (g exc (t) -∆ Γ (φ τ mod,k (g mod (t)) -B) -ψ(V (t)) + (13)
with τ = 1ms, τ mod = 5ms, k = 1, Γ = 6 and B is the baseline firing rate from VTA dopamine to VS. Γ ensures a minimum threshold to be achieved for the VTA dopamine phasic firing to enable modulation. ψ() is a bound function. where α is the learning rate equal to 0.4 The first term decreases the weights based on α and the weights keep decreasing until the bound is reached when ∆ Γ (U (t)) becomes greater than 0 at the time of the reward. The correcting update is the second term of the weight updating and the slope is increased with a weight increase encoding the duration of the interval.

It should be noted that the model postulates the learning of time to be happening before the learning of value of the stimulus i.e. its magnitude.

6) VTA:

The VTA in this model is divided into two populations based on the type of neurons as found in [START_REF] Cohen | Neurontype-specific signals for reward and punishment in the ventral tegmental area[END_REF] and the neurobiological assumptions are derived from [START_REF] Eshel | Arithmetic and local circuitry underlying dopamine prediction errors[END_REF] where a ramping VTA GABA signal did not have a significant influence on tonic dopamine firing and only affected phasic dopamine. Refer [START_REF] Eshel | Arithmetic and local circuitry underlying dopamine prediction errors[END_REF] for details.

VTA Dopamine: The VTA dopaminergic neurons convey the final reward prediction error of the system. The VTA Dopamine neurons initially fire for the US reward which progressively gets canceled and at the same time predict the US through the phasic firing it undergoes upon arrival of a CS.

τ. dV (t) dt = (-V (t) + φ τ exc,k (g exc (t)) + η(t)) (9)
U (t) = (V (t) + B) + (15)
with τ = 5ms, τ exc = 5ms, k = 1 and B is the baseline firing rate of the VTA Dopamine equal to 0.2 VTA GABA: The VTA GABA neurons encode expectation and receive their inputs from the PPN FT neurons.

τ.

dV (t) dt = (-V (t) + (g exc (t)) + η(t)) ( 16)

U (t) = (V (t)) + (2)
with τ = 20ms

This model is implemented in Python, and is using the DANA library for neuronal computation [START_REF] Rougier | DANA: Distributed (asynchronous) Numerical and Adaptive modelling framework[END_REF].

IV. EVALUATION OF THE MODEL

The paradigm used to evaluate the model is a simple CS-US associative learning task and considers also how the expectation cancels out the dopamine peak at the time of the reward. The trial duration is 500 time steps with each time step corresponding to 1ms. The stimulus is presented at the 10th time step and is kept switched on till the arrival of the reward at the 400th time step (400ms). The reward and the stimulus have by default a magnitude of 1. The number of trials for the entire conditioning to happen was 9 trials.

A. Initial Trial

During the first trial of the conditioning (Figure 4), the BLA hasn't yet learnt to associate the CS with the US. The BLA recognizes the reward signal through LH and the VTA. The CS firing gradually builds up and encodes the final magnitude of the US at the end of conditioning. The VTA fires on the delivery of the reward as shown in Figure 4. The Ventral Striatum is yet to learn the timing of the interval duration and the VTA dopamine US firing enables the Ventral Striatum to learn the duration of the interval in subsequent trials and this learning of the time happens before the learning of the magnitude. Since the CS has not been recognized as rewarding, there is no expectation at the arrival of the CS.

B. Partial Conditioning

After a few trials (four in our simulations), the magnitude of the reward is partially encoded in the BLA and the BLA fires upon the arrival of the CS. The synaptic weights between IT and BLA are updated after each rewarding trial. The learning of time of the US happens before the learning of magnitude and at this stage, the interval time has been completely learnt. The Ventral Striatum has no inhibition at the end of the interval time and allows the entire expectation to inhibit the VTA dopamine neurons. This results in partial CS firing and partial cancellation of dopamine. This partial cancellation is achieved through partial expectation developed as a result of CS firing through PPN FT neurons activating the VTA GABA. Partial conditioning results in both the CS and the US showing some firing (Figure 5).

C. Complete Reward cancellation

At the end of conditioning, the reward has been fully learnt and the VTA GABA neurons ramps to its maximum, canceling the entire US signal coming from the LH. The BLA neurons that drive firing for the CS has now encoded the magnitude of the US and drives the expectation that finally maintains the background firing of the dopamine at the point of reward delivery as shown in Figure 6.

D. Early Reward

The model is consistent with physiological data that does not treat all early rewards as the same. The expectation is substantial even at the half way point and dopamine firing is not observed for those "early" rewards that have a longer latency (that come after the half way mark). The earlier rewards fire more than the ones with longer latency and the earlier rewards that fire do not fire as much as the "unpredicted" reward. [START_REF] Fiorillo | The temporal precision of reward prediction in dopamine neurons[END_REF] The Figure 7, shows a reward delivered before the half way point (at the 100th time step) invokes a dopamine firing but less than the firing showed for an "unpredicted" reward while Figure 8, shows an early reward delivered after the half way mark (at the 300th time step), does not invoke any firing.

V. DISCUSSION

The model adds to the literature of computing the reward prediction error and is different from the other dual pathway models owing to its dissociation of magnitude and timing signals [START_REF] Brown | How the basal ganglia use parallel excitatory and inhibitory learning pathways to selectively respond to unexpected rewarding cues[END_REF], [START_REF] O'reilly | Pvlv: the primary value and learned value pavlovian learning algorithm[END_REF]. This dissociation results in a distributed manner of processing and would enable the system to be more robust and maintain the information of one dimension even if the other is changed, say even if the previous time of the interval duration is changed for the same stimulus, the value of the stimulus need not be relearned again enabling faster transitions compared to the Temporal Difference (TD) algorithm. The simulated model has a continuous representation of time and evaluates how the inhibitory signal is modulated throughout the duration of CS and US differing it from other models such as PVLV [START_REF] O'reilly | Pvlv: the primary value and learned value pavlovian learning algorithm[END_REF] which has a single point of inhibition at the expected time of the US. The model also attempts to explain the early reward scenario in a biologically plausible manner by showing that not all early rewards are the same and the ones with longer latency (after the halfway mark of the interval) don't cause any firing at all.

Since the model is also an attempt to understand the circuits behind reward processing in the brain, it makes the following predictions:

1) CE and PPN FT encode magnitude of expectation This hypothesis does not have a separate CS-US processing and suggests the canceling of the US reward has its origins in the CS firing, one of the predictions is that inhibiting or partial lesions in the Central Amygdala (CE) in a conditioning task at the time of the CS would decrease the PPN FT expectation and possibly result in a positive prediction error in VTA dopamine instead of complete cancellation.

2) PPN through VTA GABA cancels dopamine

The expectation that is encoded in PPN FT firing through VTA GABA provides the final cancellation signal required for maintaining the baseline firing of dopamine. There is a specific projection from PPN to VTA GABA that should implement this.

3) VS encodes timing

The model hypothesizes VS to learn the timing and progressively lower the inhibition that it exerts on PPN FT causing the VTA GABA to ramp and cancel the final reward.

4) Early Reward

Upon receiving early reward through the LH → PPN RD pathway, these PPN RD neurons should inhibit and reset the PPN FT neurons canceling the expectation. PPN FT neurons wouldn't show any activity after the arrival The model postulates that interval timing is learnt before magnitude of the stimulus.

The effects of the VS lesion experiments described in [START_REF] Takahashi | Temporal specificity of reward prediction errors signaled by putative dopamine neurons in rat vta depends on ventral striatum[END_REF] could be reproduced by this model. A VS lesion in this model would make the system lose its ability to track time but magnitude expectation will be preserved, resulting in the VTA GABA neurons undergoing elevated firing throughout the duration of the CS and US. Thus, the VTA dopamine neurons would only fire for those rewards that have magnitude greater than the expected reward and not for early rewards as described in the study.

The model also contends that the computation of reward processing is highly complex in the brain and could involve synchronous circuits that excite VTA dopamine neurons and simultaneously inhibit the VTA GABA neurons to cause phasic firing in VTA dopamine neurons [START_REF] Aggarwal | Neural control of dopamine neurotransmission: implications for reinforcement learning[END_REF]. Such a synchrony is not in this model's current scope. Recent studies also show LH has direct afferents to VTA GABA and reward firing could be caused by the GABAergic neurons in the LH inhibiting the VTA GABA neurons [START_REF] Nieh | Inhibitory input from the lateral hypothalamus to the ventral tegmental area disinhibits dopamine neurons and promotes behavioral activation[END_REF] and thus disinhibiting VTA dopamine neurons. The model speculates the same mecha-nisms of this model could be used in such a case where the integration of time and magnitude could converge on VTA GABA instead of PPN FT which would then only have the magnitude component and the VS still providing the timing information through direct afferents to VTA GABA. Early rewards with longer latencies would still be unable to cause VTA dopamine firing since the LH GABA might not be able to fully suppress the VTA GABA to cause a disinhibition in VTA dopamine.

VI. IMPLICATIONS AND FUTURE WORK

This model is a biologically inspired way of looking at reward prediction errors. It tries to approximate behavior in a computational manner but is also an attempt to explain the circuit and the computations involved in reward prediction error processing. It aims to describe the precise nature of how the circuits of the brain solve the phenomenon of classical conditioning. It could play the "critic" in the actorcritic reinforcement learning paradigm and could also be used to extend to those cases where the reward prediction error computation does not happen properly (for e.g., in addiction or anhedonia) [START_REF] Redish | Addiction as a computational process gone awry[END_REF]. The model currently does not account for aversive conditioning and hopes to include it in the future.

VII. CONCLUSION

The neurocomputational model described here represents a model-free reinforcement learning system and learns the CS-US association in classical conditioning. It achieves this by proposing a dissociation between magnitude and timing expectation separating it from traditional reinforcement learning approaches. The model posits the brain could be solving the dimensions involved in classical conditioning separately in such a distributed manner. Interestingly, in a modular view, such a system, having the components required to process a given natural phenomena broken down into its elemental dimensions, could enable the same dimensions to be combined with new elemental dimensions to process other natural phenomena.

Whereas the model remains simple in its nature and description and is based on homogenous units, it can reproduce a variety of aspects in classical conditioning in accordance with physiological data, namely, predicting the US by CS firing, canceling the expected US by a ramping expectation, the twin peaks of the CS firing and US firing during partial conditioning, sooner early rewards causing firing but not long latency early rewards and sooner early rewards not firing as much as "unpredicted" reward.

Fig. 1 .

 1 Fig. 1. Model diagram illustrating the neuronal structures and their connections involved in RPE computation. Pointed arrows represent excitatory connections, while rounded arrows represent inhibitory projections. Dashed lines represent learnable connections, while solid represent fixed connections in the model.

Fig. 2 .

 2 Fig. 2. The slope is decreased at every iteration until it exceeds the duration (the red line) enabling exact correction of the weight encoding the duration to be found (the black line). The colors indicate the progressive iterations As described in figure 2, weight is updated after each iteration according to the following rule: dw(t) dt = (-α.w + ∆ Γ (U (t)).w.(U (t)/(1 -U (t))) (14)

Fig. 4 .

 4 Fig. 4. Initial Trial shows arrival of the reward in LH (US) and the VTA dopamine neurons firing as a result. There isn't any firing in the PPN FT or VTA GABA since the CS is yet to be recognized as rewarding and no expectation has developed as a result. The VTA US firing also shows in the VS which enables the VS to finally learn the interval duration. All the firing rates of the populations are scaled down to fall between 0 and 1.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Partial Conditioning shows the magnitude partially learnt but the timing fully learnt. The partial magnitude learning is reflected in an expectation that shows in the firing of the CS and in partial cancellation at the US resulting in the twin peaks as observed in physiological data

Fig. 7 .

 7 Fig. 7. Sooner Early Reward invokes a firing but less than the unpredicted reward in the VTA dopamine neurons

VIII. ACKNOWLEDGEMENTS

The authors would like to acknowledge the following grants which have been a major support to this research.

• Indo-French CEFIPRA Grant for the project Basal Ganglia at Large (No. DST-INRIA 2013-02/Basal Ganglia dated 13-09-2014) • Internships programme at INRIA, 6 month Internship with Team Mnemosyne at INRIA Bordeaux -Sud-Ouest


Bogdan Aman
email: bogdan.aman@gmail.com

Gabriel Ciobanu

Maciej Koutny
email: maciej.koutny@newcastle.ac.uk

Behavioural Equivalences over Migrating Processes with Timers

Keywords: mobility, timer, process algebra, bisimulation, behaviour, equivalence

The temporal evolution of mobile processes is governed by independently operating local clocks and their migration timeouts. We define a formalism modelling such distributed systems allowing (maximal) parallel execution at each location. Taking into account explicit timing constraints based on migration and interprocess communication, we introduce and study a number of timed behavioural equivalences, aiming to provide theoretical underpinnings of verification methods. We provide several relationships between such behavioural equivalences.

Introduction

Process calculi are a family of formalisms used to model distributed systems. They provide algebraic laws allowing a high-level description and analysis of concurrent processes, behavioural equivalences (e.g., bisimulations) between processes, and automated tools for the verification of interaction, communication, and synchronization between processes. During the past couple of decades, a number of calculi supporting process mobility were defined and studied; in particular, π-calculus [START_REF] Milner | Communicating and mobile systems -the π-calculus[END_REF], fusion calculus [START_REF] Parrow | The fusion calculus: Expressiveness and symmetry in mobile processes[END_REF], join calculus [START_REF] Fournet | The join calculus: A language for distributed mobile programming[END_REF], and mobile ambients [START_REF] Cardelli | Mobile ambients[END_REF]. Various specific devices were introduced to obtain such formalisms, including explicit locations in distributed π-calculus [START_REF] Hennessy | A distributed π-calculus[END_REF], explicit migration and timers in timed distributed π-calculus [START_REF] Ciobanu | Timers for distributed systems[END_REF], and timed mobile ambients [START_REF] Aman | Timed mobile ambients for network protocols[END_REF]. Time is an important aspect of distributed computing systems, and can play a key role in their formal description. Since time is a complex subject, its introduction to the domain of process calculi has received a lot of attention in, e.g., [START_REF] Baeten | Discrete time process algebra: Absolute time, relative time and parametric time[END_REF][START_REF] Berger | Towards Abstractions for Distributed Systems[END_REF][START_REF] Chen | Timed Processes: Models, Axioms and Decidability[END_REF][START_REF] Hennessy | A process algebra for timed systems[END_REF][START_REF] Nicollin | The algebra of timed processes, atp: Theory and application[END_REF][START_REF] Yi | A Calculus of Real-Time Systems[END_REF]. Papers like these assume the existence of a global clock which is usually required in the description of complex systems. However, there are several applications and systems for which considering a global clock would be inappropriate. This paper presents such an approach, in which local clocks operate independently, and so processes at different locations evolve asynchronously. On the other hand, processes operating at the same location evolve synchronously. Overall, temporal evolution of mobile processes is governed by independent local clocks and their migration timeouts.

The ever increasing complexity of mobile processes calls for the development of effective techniques and tools for the automated analysis and verification of their properties including, in particular, behavioural equivalences between systems. Bisimulation is the most common mathematical concept used to capture behavioural equivalence between processes. The corresponding equivalence relation, called bisimilarity, is used to abstract from certain details of the systems, and is widely accepted as a standard behavioural equivalence over processes. Several kinds of bisimulations had been defined (e.g, strong or weak bisimulation for π-calculus [START_REF] Milner | Communicating and mobile systems -the π-calculus[END_REF]).

In the paper [START_REF] Ciobanu | Modelling and verification of timed interaction and migration[END_REF], we defined TiMo, a basic language for mobile systems in which it is possible to add timers to control process mobility and interaction. After that, in [START_REF] Ciobanu | Timed mobility in process algebra and petri nets[END_REF], a local clock was assigned to each location of a system modelled in TiMo. Each such clock determines the timing of actions executed at the corresponding location. Then, starting from TiMo, we created a flexible software platform supporting the specification of agents and their physical distribution, allowing also a timed migration in a distributed environment [START_REF] Ciobanu | Flexible software architecture and language for mobile agents[END_REF]. We obtained this implementation by using an advanced software technology, creating a platform for mobile agents with time constraints. In this paper, we consider the same specification language for mobile agents with timeouts, i.e., TiMo, and define various behavioural equivalences taking into account both timers and locations (for migration). We study both the strong bounded timed bisimulation and strong open bounded timed bisimulation. We then relax the bisimulations being considered, by applying the so-called "up-to" technique. In [START_REF] Ciobanu | Timed migration and interaction with access permissions[END_REF] we discussed the TravelShop running example, in which clients buy tickets to predefined destinations from travel agents. One can use the bisimilarities defined in this paper to differentiate, for example, between two travel agents using the same databases for prices, but having different delays for providing the answer (in an urgent situation, the faster one would be preferred). On the other hand, it is possible to define equivalence classes of agents offering similar services with respect to the waiting time (possibly up to an acceptable time difference). As TravelShop is a non-trivial example, in this paper we use a simpler one in order to illustrate how the newly defined bisimulations can be applied.

The paper is organised as follows. We start in Section 2 with a brief presentation of TiMo, including its syntax and operational semantics, and introduce an example to illustrate the basic features of TiMo. In Section 3, we formally define in the framework of TiMo a number timed bisimulations, and discuss some of their properties. Conclusion and references end the paper. Due to space limitations, proofs of formal results are omitted.

Syntax and Semantics of TiMo

Timing constraints for migration allow one to specify what is the time instant for a process to move from one location to another. A timer denoted by ∆3 associated to a migration action go ∆3 work indicates that the process moves to location work after 3 time units. It is also possible to constrain the waiting for a communication on a channel; if a communication action does not happen before a deadline, the waiting process gives up and switches its operation to an alternative. A timer ∆5 associated to an output action a ∆5 ! 10 makes the channel available for communication only for the period of 5 time units. We assume suitable data sets including a set Loc of locations and a set Chan of communication channels. We use a set Id of process identifiers, and each id ∈ Id has the arity m id .

The syntax of TiMo, an extension of TiMo from [START_REF] Ciobanu | Flexible software architecture and language for mobile agents[END_REF], is given in Table 1, where P are processes and N networks. Moreover, for each id ∈ Id there is a unique process definition of the form: id(u 1 , . . . , u m id ) df = P id , where the u i 's are distinct variables playing the role of parameters. The process 0 is called primitive. In Table 1, it is assumed that:

• a ∈ Chan is a channel; • lt ∈ N ∪ {∞}
is a deadline, where lt stands for local time;

• v is a tuple of expressions built from values, variables and allowed operations; • u is a tuple of variables; • l is a location or a location variable; • P is a purely technical device which is used in the subsequent formalization of the structural operational semantics of TiMo; intuitively, P specifies a process P which is temporarily (i.e., until a clock tick) stalled and so cannot execute any action.

Processes P ::= a ∆lt ! v then P else P (output) a ∆lt ?(u) then P else P (input)

go ∆lt l then P else P (move)

P | P (parallel) 0 (termination) id(v) (recursion)
P (stalling) Located processes L ::= l[[P ]] Networks N ::= L L | N Table 1. TiMo Syntax
The only variable binding constructor is "a ∆lt ?(u) then P else P " which binds the variable u within P (but not within P ). We use f v(P ) to denote the free variables of a process P (and similarly for networks). For a process definition, we assume that f v(P id ) ⊆ {u 1 , . . . , u m id } and so the free variables of P id are parameter bound. Processes are defined up to alpha-conversion, and {v/u, . . .}P is obtained from P by replacing all free occurrences of a variable u by v, possible after alpha-converting P in order to avoid clashes.

Intuitively, a process "a ∆lt ! v then P else P " attempts to send a tuple of values v over channel a for lt time units. If successful, it continues as process P ; otherwise it continues as process P . Similarly, "a ∆lt ?(u) then P else P " is a process that attempts for lt time units to input a tuple of values and substitute them for the variables u. Mobility is implemented by a process "go ∆lt l then P else P " which moves from the current location to the location given by l after lt time units. If the movement does not take place in lt time units, because location l does not exist (l ∈ Loc), the process continues at the current location with P . Note that since l can be a variable, and so its value is assigned dynamically through the communication with other processes, migration actions support a flexible scheme for the movement of processes from one location to another. Processes are further constructed from the (terminated) process 0 and parallel composition P |P . A located process l[[P ]] specifies a process P running at location l, and a network is composed out of its components N | N . A network N is well-formed if the following hold: there are no free variables in N , and there are no occurrences of the special symbol in N . The set of processes is denoted by P, the set of located processes by L, and the set of networks by N .

By delaying the migration to another location, we can model in a simple way the movement time of processes along the network which is, in general, outside the control of a system designer.

The first component of the operational semantics of TiMo is the structural equivalence ≡ on networks; it is the smallest congruence such that the first four equalities in Table 2 hold. Its role is to rearrange a network in order to apply the action rules which are also given in Table 2. Using the first four equalities in Table 2, one can always transform a given network N into a finite parallel composition of located processes of the form l

1 [[P 1 ]] | . . . | l n [[P n ]
] such that no process P i has the parallel composition operator at its topmost level. Each located process

l i [[P i ]] is called a component of N , and the parallel composition is called a component decomposition of the network N . Table 2 introduces two kinds of rules, N λ -→ N and N √ l --→ N .
The former is an execution of an action λ, and the latter a time step at location l. In the rule (Move), if location l ∈ Loc, the process from location l migrates to l . In the rule (Time), N -→ l means that no action at location l (i.e, no action λ of the form l l or a v @l or id@l) can be applied to N . It can be noticed that in rule (Time) we use negative premises, i.e., an activity is performed in the absence of other actions. This is due to the fact that sequencing the evolution in time units can only be defined using negative premises, as done for sequencing

(NComm) N | N ≡ N | N (NAssoc) (N | N ) | N ≡ N | (N | N ) (NSplit) l[[P | P ]] ≡ l[[P ]] | l[[P ]] (LComm) l[[P | Q]] ≡ l[[Q | P ]] (Move) l ∈ Loc l[[go ∆0 l then P else Q]] l l --→ l [[ P ]] (Com) l[[a ∆lt ! v then P else Q | a ∆lt ?(u) then P else Q ]] a v @l ----→ l[[ P | {v/u}P ]] (Call) l[[id(v)]] id@l ---→ l[[ {v/u}P id ]] (Par) N λ -→ N N | N λ -→ N | N (Equiv) N ≡ N N λ -→ N N ≡ N N λ -→ N (Time) N -→ l N √ l --→ φ l (N ) Table 2. TiMo Operational Semantics
processes in [START_REF] Bloom | Bisimulation can't be traced[END_REF][START_REF] Groote | Transition system specifications with negative premises[END_REF][START_REF] Moller | Axioms for Concurrency[END_REF]. Moreover, φ l (N ) is obtained by taking the component decomposition of N and simultaneously replacing all components:

l[[a ∆lt ω then P else Q]] by l[[Q]] if lt = 0; l[[a ∆lt-1 ω then P else Q]] otherwise. l[[go ∆lt l then P else Q]] by l[[Q]] if lt = 0; l[[go ∆lt-1 l then P else Q]] otherwise.
where ω stands for ! v or ?(u). After that, all the occurrences of the symbol in N are erased (this is done since processes that were unfolded or interacted with other processes or migrated need to be activated).

The rules of Table 2 express executions of individual actions. A complete computational step in distributed systems is captured by a derivation of the form N Λ@l = = =⇒ N , where Λ = {λ 1 , . . . , λ m } (m ≥ 0) is a finite multiset of actions for some location l (i.e., actions λ i of the form l l or a v @l or id@l) such that

N λ1 -→ N 1 . . . N m-1 λm --→ N m √ l --→ N .
That means that a derivation is a condensed representation of a sequence of individual actions followed by a clock tick, all happening at the same location. Intuitively, we capture the cumulative effect of the concurrent execution of the multiset of actions Λ at location l, and say that N is directly reachable from N . Note that if there is only a time progression at a location l, we write N ∅@l ==⇒ N . We now give an example that shows the modelling features of our formalism. The example will be used in the next section to illustrate various timed bisimulations defined in this paper.

Example 1. To exemplify how TiMo works, we use an UrgentTravel running example in which a client process attempt to initiate a travel process as soon as possible it receives an offer. The scenario involves three locations with their local clocks and three processes. The role of each location is as follows: (i) office is a location where the client process starts and ends its journey, (ii) agency i , for 1 ≤ i ≤ 2 are the travel agencies from which the client can found out about the tickets prices. The role of each of the processes is:

• client is a process that initially resides the office location, and is determined to pay for a flight as soon as it receive one offer from the two existing agencies. After sending an email to each agency, it awaits for the quickest response to initiate the travel;

• agent i (1 ≤ i ≤ 2
) is a process that resides inside the agency i location, answers to the emails received from clients and sends them the agency offers.

In this formalism we use timers in order to impose deadlines on the execution of communications and migrations. Each location has its local clock which determines the timing of actions executed at that location. The timeout of a migration action indicates the network time limit for that action . We assume local clock for the three locations such that:

• office and agency 1 have the same time unit;

• agency 2 has a double time unit (i.e., equals with 2 time units of the local clock of the locations office and agency 1 ).

The specifications that capture the essential features of the scenario are: agent i = a ∆5 ! offer i then agent i else agent i client = go ∆2 agency 1 then a ∆1 ?(x) then go ∆2 office then d ∆1 ! x else 0 else 0 else 0 | go ∆2 agency 2 then a ∆1 ?(x) then go ∆2 office then d ∆1 ! x else 0 else 0 else 0 | d ∆6 ?(y) then travel else client In order to simplify the writing we have not added any alternatives for the communication action d ∆t ! x . Next we illustrate the execution of U rgentT ravel.

The expired timers cause an immediate migration or the selection of an offer received in the communication action.

The initial network is

U rgentT ravel = office[[client]] | agency 1 [[agent 1 ]] | agency 2 [[agent 2 ]] ∅@office =====⇒ ∅@office =====⇒ ∅@agency 1 = ===== =⇒ ∅@agency 1 = ===== =⇒ ∅@agency 2 = ===== =⇒ 5× (Time) office[[go ∆0 agency 1 then a ∆1 ?(x) then go ∆2 office then d ∆1 ! x else 0 else 0 else 0 | go ∆0 agency 2 then a ∆1 ?(x) then go ∆2 office then d ∆1 ! x else 0 else 0 else 0 | d ∆4 ?(y) then travel else client]] | agency 1 [[a ∆3 ! offer 1 then agent 1 else agent 1 ]] | agency 2 [[a ∆4 ! offer 2 then agent 2 else agent 2 ]] {office agency 1 , office agency 2 }@office = ======================== =⇒ 2× (Move) {a offer 1 @agency 1 }@agency 1 = ================== =⇒ {a of f er2 @agency 2 }@agency 2 = =================== =⇒ 2× (Com) office [[d ∆4 ?(y) then travel else client]] | agency 1 [[agent 1 | go ∆2 office then d ∆1 ! offer 1 else 0]] | agency 2 [[agent 2 | go ∆2 office then d ∆1 ! offer 2 else 0]] ∅@office =====⇒ ∅@office =====⇒ ∅@agency 1 = ===== =⇒ ∅@agency 1 = ===== =⇒ ∅@agency 2 = ===== =⇒ 5× (Time) office [[d ∆1 ?(y) then travel else client]] | agency 1 [[a ∆2 ! offer 1 then agent 1 else agent 1 | go ∆0 office then d ∆1 ! offer 1 else 0]] | agency 2 [[a ∆4 ! offer 2 then agent 2 else agent 2 | go ∆1 office then d ∆1 ! offer 2 else 0]] {agency 1 office}@agency 1 =================⇒ 1× (Move) {d offer 1 @office}@office ================⇒ 1× (Com) office [[travel{offer 1 /y}]] | agency 1 [[a ∆2 ! offer 1 then agent 1 else agent 1 ]] | agency 2 [[a ∆4 ! offer 2 then agent 2 else agent 2 | go ∆1 office then d ∆1 ! offer 2 else 0]]

Timed Bisimulation in TiMo

In what follows, we define various equivalences for processes and networks by considering their temporal behaviour. In a similar way as in timed distributed π-calculus of [START_REF] Ciobanu | Behaviour equivalences in timed distributed π-calculus[END_REF], we start by extending the standard notion of strong bisimilarity to take into account timed transitions.

Definition 1 (strong timed simulation and bisimulation).

Let R ⊆ N × N be a binary relation on networks of processes.

1. R is a strong timed simulation (ST simulation) if

(N 1 , N 2 ) ∈ R ∧ N 1 λ -→ N 1 =⇒ ∃N 2 ∈ N : N 2 λ -→ N 2 ∧ (N 1 , N 2 ) ∈ R .
2. R is a strong timed bisimulation (ST bisimulation) if both R and R -1 are strong timed simulations. 3. The strong timed bisimilarity is the union ∼ of all ST bisimulations. Essentially, the above definition treats timed transitions just as any other transitions, and therefore coincides with the original notion of bisimilarity for labelled transition systems. It is easy to check that ∼ is an equivalence relation, and also the largest strong timed bisimulation. From the point of view of the evolutionary behaviour of TiMo networks, a crucial result is that strong timed bisimulation can be used to compare the complete computational steps of two systems.

Theorem 1. Let N 1 , N 2 be two networks such that N 1 ∼ N 2 . If N 1 Λ@l = = =⇒ N 1 then there exists N 2 ∈ N such that N 2 Λ@l = = =⇒ N 2 and N 1 ∼ N 2 .
Example 2. Inspired by the UrgentTravel of Example 1, consider the following two located processes:

L 1 = agency 1 [[a ∆1 ?(x) then go ∆4 office then d ∆1 ! x else 0 else 0 | a ∆3 ! offer 1 then agent 1 else agent 1 ]] L 2 = agency 1 [[a ∆2 ?(x) then go ∆3 office then d ∆1 ! x else 0 else 0 | a ∆4 ! offer 1 then agent 1 else agent 1 ]]
This means that all the processes having the same timer for go are equivalent, even the timers of the communication channels differ. However, if the client reaches office after different number of time units, the two located processes are not strong timed bisimilar, i.e., (L 1 ∼ L 2 ), since they have different evolutions in time after 3 units of time.

L 1 {a offer 1 @agency 1 }@agency 1 = ================== =⇒ (Com) agency 1 [[go ∆4 office then d ∆1 ! offer 1 else 0 | agent 1 ]] ∅@agency 1 = ===== =⇒ ∅@agency 1 = ===== =⇒ ∅@agency 1 = ===== =⇒ 3× (Time) agency 1 [[go ∆1 office then d ∆1 ! offer 1 else 0 | a ∆2 ! offer 1 then agent 1 else agent 1 ]] ∅@agency 1 = ===== =⇒ (Time) agency 1 [[go ∆0 office then d ∆1 ! offer 1 else 0 | a ∆1 ! offer 1 then agent 1 else agent 1 ]] . . . L 2 {a offer 1 @agency 1 }@agency1 = ================== =⇒ (Com) agency 1 [[go ∆3 office then d ∆1 ! offer 1 else 0 | agent 1 ]] ∅@agency 1 = ===== =⇒ ∅@agency 1 = ===== =⇒ ∅@agency 1 = ===== =⇒ 3× (Time) agency 1 [[go ∆0 office then d ∆1 ! offer 1 else 0 (HPar) | a ∆2 ! offer 1 then agent 1 else agent 1 ]] {agency 1 office}@agency1 =================⇒ (Move) office[[d ∆1 ! offer 1 ]] . . .
The above definition of equivalence compares the evolution of whole networks, but does not provide means for reasoning about equivalence of compositionally defined networks. To be able to reason in this way, one may augment (only for the purpose of dealing with equivalences) the operational semantics of processes with two additional rules relating to communication, which intuitively represent individual evolutions of interacting processes:

(Snd) l[[a ∆lt ! v then P else Q]] a! v @l -----→ l[[ P ]] (Rcv) l[[a ∆lt ?(u) then P else Q]] a? v @l -----→ l[[ {v/u}P ]]
Even when these two new rules are added, rule (Com) remains unchanged. However, in this paper, we concentrate on the network-level evolutions, leaving the detailed treatment of compositionally defined network to a future work.

Strong Bounded Timed Equivalences

We first introduce some additional terminology:

• If |Loc| = n, then by t = (t 1 @l 1 , . . . , t n @l n ) ∈ (N@Loc) n we denote a tuple in which each location l k ∈ Loc has attached a natural number t k . • for t, u ∈ (N@Loc) n , we say that t ≤ u if ∀l ∈ Loc, t l @l and u l @l we have that t l ≤ u l . • σ is an function used to substitute free names, with other names.

• A timed relation over N is any relation R ⊆ N × (N@Loc) n × N . • The identity timed relation is ident t df = {(N, t, N ) | N ∈ N , t ∈ (N@Loc) n }. • The inverse of a timed relation R is R -1 df = {(N 2 , t, N 1 ) | (N 1 , t, N 2 ) ∈ R} .
• The composition of timed relations R 1 and R 2 is

R 1 R 2 df = {(N, t, N ) | ∃N ∈ N : (N, t, N ) ∈ R 1 ∧(N , t, N ) ∈ R 2 } . • If R is a timed relation and t ∈ N then R t df = {(N 1 , N 2 ) | (N 1 , t, N 2 ) ∈ R} is R's t-projection. We also denote R ∞ df = t∈(N@Loc) n R t . • R is a timed equivalence if R ∞ is an equivalence relation, and an equivalence up-to time u ∈ N if 0≤t<u R t is an equivalence relation.
The first notion of equivalence takes into account the timed behaviour requiring an exact match of transitions of two networks, for their entire evolution. Sometimes these requirements are too strong. According to [START_REF] Posse | Theory and implementation of a real-time extension to the pi-calculus[END_REF] where a similar approach is presented, real-time distributed systems usually require a certain behaviour within a given amount t of time units. That is why we restrict our equivalences up-to a time t, namely to a local time t l for each location l ∈ Loc, defining bounded timed equivalences.

Definition 2 (strong bounded timed simulation and bisimulation).

Let R ⊆ N × (N@Loc) n × N be a timed relation over N .

1. R is a strong timed simulation (SBT simulation) if (N 1 , (. . . , t l @l, . . .), N 2 ) ∈ R N 1 √ l --→ N 1 =⇒ ∃N 2 ∈ N : N 2 √ l --→ N 2 (N 1 , (. . . , (t l -1)@l, . . .), N 2 ) ∈ R and (N 1 , t, N 2 ) ∈ R N 1 λ -→ N 1 =⇒ ∃N 2 ∈ N : N 2 λ -→ N 2 (N 1 , t, N 2 ) ∈ R ,
where λ is not of the form √ l . 2. R is a strong timed bisimulation (SBT bisimulation) if both R and R -1 are strong bounded timed simulations. 3. The strong bounded timed bisimilarity is the union of all SBT bisimulations.

Theorem 2. Let N 1 , N 2 be two networks s.t. N 1 t N 2 , with t = (. . . , t l @l, . . .).

If N 1 Λ@l = = =⇒ N 1 then there exists N 2 ∈ N such that N 2 Λ@l = = =⇒ N 2 and N 1 t N 2 , with t = (. . . , (t l -1)@l, . . .). Corollary 1. Let N 1 , N 2 be two networks such that N 1 N 2 . If N 1 Λ@l = = =⇒ N 1 then there exists N 2 ∈ N such that N 2 Λ@l = = =⇒ N 2 and N 1 N 2 .
Example 3. Consider the two networks of Example 2. Even if those networks have different definitions, they are strong bounded timed bisimilar before the 3rd time unit at location agency 1 (N 1 (3@agency1,0@office) N 2 ) since they have the same evolutions, at location agency 1 , during this period. Hence N 1 and N 2 cannot be identified by strong timed bisimulation, but this is possible by using strong bounded timed bisimulation up-to time 3 at location agency 1 . However if t > (3@agency 1 , 0@office) we have that N 1 t N 2 .

Proposition 1. 1. Identity, inverse, composition and union of SBT bisimulations are SBT bisimulations. 2.

is a timed equivalence. 3.

is an SBT bisimulation. 4.

is the largest SBT bisimulation.

Strong bounded timed bisimulation satisfies the following properties showing that the equivalence up-to a certain deadline includes the equivalence up-to any bound below that deadline.

Proposition 2. If N t N then N u N , for every u ≤ t.

Strong Open Bounded Timed Equivalences

In this section we are looking to define a timed equivalence which is closed under arbitrary substitutions. According to the following example, the SBT equivalence is not closed under substitution.

Example 4. Consider the client from Example 1 decides to query someone from the office to which agency to go. We have two situations: the client decides to use the receives information (process P 1 ) or not (process P 2 ). In order to simplify the notations we use some generic processes P and Q to illustrate the evolution of the client after moving to a agency 1 or agency 2 .

N 1 = l[[P 1 ]] | agency 1 [[R 1 ]] | agency 2 [[R 2 ]] N 2 = l[[P 2 ]] | agency 1 [[R 1 ]] | agency 2 [[R 2 ]
] P 1 = query ∆1 ! agency 1 then P else Q | query ∆5 ?(u) then (go ∆3 u then P else Q ) else 0, P 2 = query ∆1 ! agency 1 then P else Q | query ∆5 ?(u) then (go ∆3 agency 1 then P else Q ) else 0. We have N 1 t N 2 , for any t = (t 1 @l, t 2 @agency 1 , t 3 @agency 2 ), since applying the rules for communication on channel query: {query agency 1 @l}@l = ============= =⇒ l[[go ∆3 agency 1 then P else Q ]] (Com), we obtain the same network N that has the same evolution (N t N ). However, if for these two processes P 1 and P 2 we consider the substitution σ = {agency 2 /agency 1 } (the client receives the answer agency 2 ), after applying (N 1 , (. . . , t l @l, . . .), N 2 ) ∈ R The intuition is that two located processes are equivalent whenever all possible instantiations (substitutions of their free names) have matching transitions.

l[[P 1 σ]]
N 1 σ √ lσ --→ N 1 =⇒ ∃N 2 ∈ N : N 2 σ √ lσ --→ N 2 (N 1 , (. . . , (t l -1)@l, . . .)σ, N 2 ) ∈ R and (N 1 , t , N 2 ) ∈ R N 1 σ λ -→ N 1 =⇒ ∃N 2 ∈ N : N 2 σ λ -→ N 2 (N 1 , t σ, N 2 ) ∈ R ,
Theorem 3. Let N 1 , N 2 be two networks s.t. N 1 o t N 2 , with t = (. . . , t l @l, . . .). If N 1 Λ@l = = =⇒ N 1 then there exists N 2 ∈ N such that N 2 Λ@l = = =⇒ N 2 and N 1 o t N 2 , , with t = (. . . , (t l -1)@l, . . .)σ. Corollary 2. Let N 1 , N 2 be two networks such that N 1 o N 2 . If N 1 Λ@l = = =⇒ N 1 then there exists N 2 ∈ N such that N 2 Λ@l = = =⇒ N 2 and N 1 o N 2 .

Relaxing Timed Bisimulation by Using "Up-To" Technique

In what follows we use the "up-to" technique presented in [START_REF] Sangiorgi | A theory of bisimulation for the π-calculus[END_REF][START_REF] Sangiorgi | The π-Calculus -a theory of mobile processes[END_REF] in the context of bounded timed bisimulations. The standard proof technique to establish that N 1 and N 2 are bisimilar is to find a bisimulation R s.t. (N 1 , N 2 ) ∈ R and R is closed under transitions, namely the derivatives (N 1 , N 2 ) of (N 1 , N 2 ) are also in R. Sometimes it is difficult to find directly such a relation R. Instead, there is an useful alternative technique, the so-called bisimulation "up-to" some relation R : for a relation R, which is not a bisimulation, if (N 1 , N 2 ) ∈ R, then the derivatives (N 1 , N 2 ) are in R . Under certain conditions we can establish that N 1 and N 2 are bisimilar. For this technique, in [START_REF] Sangiorgi | On the proof method for bisimulation (extended abstract)[END_REF] is presented a general framework that works on untimed labelled transition systems. We cannot make direct use of that framework, but we can extend it in a straightforward manner to timed labelled transition systems. We begin by introducing a notion of "progressing" a timed relation towards another timed relation:

Definition 5. Let R, R be any timed relations. We say that R strongly progresses to R , written R R , if

(N 1 , (. . . , t l @l, . . .), N 2 ) ∈ R N 1 √ l --→ N 1 =⇒ ∃N 2 ∈ N : N 2 √ l --→ N 2 (N 1 , (. . . , (t l -1)@l, . . .), N 2 ) ∈ R and (N 1 , t, N 2 ) ∈ R N 1 λ -→ N 1 =⇒ ∃N 2 ∈ N : N 2 λ -→ N 2 (N 1 , t, N 2 ) ∈ R .
where λ is not of the form √ l . This definition is similar to that of SBT bisimulation, except that the derivatives (N 1 , t, N 2 ) must be in R rather than in R.

Remark 1. R is a SBT bisimulation if and only if R R.

Proposition 6. If R R and R is a SBT bisimulation, then R is also a SBT bisimulation.

Therefore, to establish that N 1 t N 2 it is enough to find a relation R with (N 1 , t, N 2 ) ∈ R which strongly progresses to a known SBT bisimulation R . The choice of R depends on the particular equivalence we are trying to establish. One of the most common cases is when R = . However, in general we may not have a relation R known to be a bisimulation. Nevertheless we may find that R progresses to a relation R = F(R) for some function F over relations. The idea is that if R progresses to F(R) and F satisfies certain conditions, then R is included in . Thus, to establish N 1 t N 2 we need to find such an F whenever R contains (N 1 , t, N 2 ).

These functions F over relations are characterized in [START_REF] Sangiorgi | The π-Calculus -a theory of mobile processes[END_REF] as follows: Definition 6. A function F is strongly safe if for any timed relations R, R , if R ⊆ R and R R , then F(R) ⊆ F(R ) and F(R) F(R ).

More details about the "up-to" techniques and safe functions can be found in [START_REF] Sangiorgi | The π-Calculus -a theory of mobile processes[END_REF].

Conclusion

This paper presents an approach in which local clocks operate independently, and so processes at different locations evolve asynchronously; on the other hand, processes operating at the same location evolve synchronously, and temporal evolution of mobile processes is governed by independent local clocks and their migration timeouts. A computational step captures the cumulative effect of the concurrent execution of a multiset of actions executed at a certain location.

In process calculi like distributed π-calculus, timed distributed π-calculus and other formalisms with explicit migration operators, bisimulations are used to compare behaviours of mobile processes evolving in distributed systems with explicit locations. Bisimulations are behavioural equivalences used to study the properties of a concurrent system by verifying its bisimilarity with a system known to enjoy those properties. Moreover, given the model of a system, bisimulations can be used to consider equivalent simplified models.

In this paper we define three behavioural equivalences between migrating process in distributed systems in terms of local time and locations. We prove that strong timed bisimulation can be used to compare the complete computational steps of two networks. We also prove that two networks that are strong bounded timed bisimilar up to a certain deadline t remain equivalent at any time u before t. Finally, we define a strong open bounded timed equivalence that is closed under arbitrary substitutions; to establish that two networks are strong bounded timed equivalent, it is enough to find a relation that strongly progresses to a known SBT relation.

Among the next steps to follow, we intend to work with these behavioural equivalences over mobile agents evolving in distributed systems by using a software platform presented extensively in [START_REF] Ciobanu | Flexible software architecture and language for mobile agents[END_REF].

  l[[P 1 ]] {query agency 1 @l}@l = ============= =⇒ l[[go ∆3 agency 1 then P else Q ]] (Com), and l[[P 2 ]]

{query agency 2

 2 @l}@l = ============= =⇒ l[[go ∆t agency 2 then P else Q ]] (Com), andl[[P 2 σ]] {query agency 1 @l}@l = ============= =⇒ l[[go ∆t agency 1 then P else Q ]] (Com).we obtain that the networks evolve in a different way, and thusN 1 t N 2 .Thus we define the following equivalence. Definition 3 (strong open bounded timed simulation and bisimulation). Let R ⊆ N × (N@Loc) n × N be a timed relation over N . 1. R is a strong open bounded timed simulation (SOBT simulation) if, for every substitution σ,

  where λ is not of the form √ l . 2. R is a strong open bounded timed bisimulation (SBT bisimulation) if both Rt and R -1 are strong open bounded timed simulations. 3. The strong open bounded timed bisimilarity is the union o of all SOBT bisimulations.

Acknowledgement

The work of Bogdan Aman and Gabriel Ciobanu was supported by a grant of the Romanian National Authority for Scientific Research, CNCS-UEFISCDI, project number PN-II-ID-PCE-2011-3-0919, while the work of Bogdan Aman was also supported by POSDRU/89/1.5/S/49944.