Benoît Fabrèges

A PARADIGM FOR WELL-BALANCED SCHEMES FOR TRAVELING WAVES EMERGING IN PARABOLIC BIOLOGICAL MODELS

We propose a methodology for designing well-balanced numerical schemes to investigate traveling waves in parabolic models from mathematical biology. We combine well-balanced techniques for parabolic models known in the literature with the so-called LeVeque-Yee formula as a dynamic estimate for the spreading speed. This latter formula is used to consider the evolution problem in a moving frame at each time step, where the equations admit stationary solutions, for which well-balanced techniques are suitable. Then, the solution is shifted back to the stationary frame in a well-balanced manner. We illustrate this methodology on parabolic reaction-diffusion equations, such as the Fisher/Kolmogorov-Petrovsky-Piskunov Equation, and a class of equations with a cubic reaction term that exhibit a transition from pulled to pushed waves. We show that the numerical schemes capture in a consistent way simultaneously the wave speed and, to an extent, the so-called Bramson delay.

Introduction

Many biological systems may give rise to spatial spreading phenomena, such as for instance the invasion of rodents [START_REF] Skellam | Random dispersal in theoretical populations[END_REF] or the collective displacement of a population of Escherichia coli cells in a micro-channel [START_REF] Adler | Chemotaxis in Bacteria[END_REF]. Mathematical modelling of this phenomena has contributed to a better understanding of these spreading phenomena and gives also rise to a quantitative framework, which may for instance account for the spreading speed. On the mathematical side, we may cite as models, that have attracted much attention, reaction-diffusion equations (e.g. [START_REF] Fisher | The wave of advance of advantageous genes[END_REF][START_REF] Kolmogorov | Etude de l équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique[END_REF][START_REF] Aronson | Nonlinear diffusion in population genetics, combustion, and nerve pulse propagation[END_REF][START_REF] Hadeler | Travelling fronts in nonlinear diffusion equations[END_REF][START_REF] Fife | The approach of solutions of nonlinear diffusion equations to travelling front solutions[END_REF]) and parabolic models of chemotaxis (e.g. [START_REF] Keller | Traveling bands of chemotactic bacteria: A theoretical analysis[END_REF][START_REF] Saragosti | Directional persistence of chemotactic bacteria in a traveling concentration wave[END_REF]). The biological question of the spreading of a population translates then into the mathematical question of the existence of a traveling wave, i.e. solutions that are stationary in a moving frame of reference.

The mathematical investigation of traveling waves for biological models naturally raises the question of the numerical investigation of these traveling wave solutions. Numerically, this question is difficult, as it requires an accurate scheme over large domains in order to capture precisely the behavior of the solution over large time. Moreover, the numerical scheme should meet the following specifications: (i) find the accurate spreading speed; (ii) find the accurate traveling wave profile. In fact, these two questions are related, since mathematically the wave speed is intimately connected to the wave profile: wave profiles tend to decrease exponentially at their leading edge and their exponential decay parameter can be tied to the wave speed through a dispersion relation. Hence, we wish to find numerical schemes that capture at the same time accurately a wave speed, as well as an exponential decay, which requires great precision even at very low orders of magnitude, because of the nature of exponential decays.

In this article, we will use so-called well-balanced (WB) schemes in order to investigate these questions. The notion of WB schemes has been introduced by [START_REF] Greenberg | A Well-Balanced Scheme for the Numerical Processing of Source Terms in Hyperbolic Equations[END_REF] and consists roughly speaking in numerical schemes that preserve the steady state of the continuous model, i.e. the steady state of the numerical scheme matches the discretization of the steady state of the continuous model. Because of this property, these schemes compute faithfully the stationary profiles. Many different approaches to such schemes exist and we refer the reader for instance to [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF] for a general discussion on these techniques. However in their essence, WB schemes accurately balance the flux terms and the source terms. WB schemes have been investigated for parabolic equations [START_REF]Viscous Equations Treated with $\mathcalL$-Splines and Steklov-Poincaré Operator in Two Dimensions[END_REF][START_REF] Gosse | NUMERICAL HIGH-FIELD LIMITS IN TWO-STREAM KINETIC MODELS AND 1D AGGREGATION EQUATIONS[END_REF]. In these works, the strategy consists in interpolating the function in each cell with L -splines (see [START_REF]L-Splines and Viscosity Limits for Well-Balanced Schemes Acting on Linear Parabolic Equations[END_REF] for an overview on the concept of L -splines), which correspond to the stationary solutions of the problem. The interpolation is then continous inside each cell, but at the nodes of the grid it may have C 1 -discontinuities. The time integration then is computed by using the C 1 defect at each node. Of note, most referenced models are one-dimensional models, but in recent years the leap to higher dimensions has attracted interest [START_REF] Després | The structure of well-balanced schemes for Friedrichs systems with linear relaxation[END_REF][START_REF]A two-dimensional version of the Godunov scheme for scalar balance laws[END_REF][START_REF] Bretti | Diffusive limit of a two-dimensional well-balanced approximation to a kinetic model of chemotaxis[END_REF].

Yet, in the aformentioned studies, most often the schemes were designed to describe a stationary state, not a traveling wave. Whilst in some cases the propagation speed can be computed explicitly and thus be used in the implementation of the scheme, in most cases the propagation speed is not known a priori. This adds to the difficulty of designing a WB scheme for traveling waves. For instance in the work [START_REF] Calvez | Travelling Chemotactic Aggregates at Mesoscopic Scale and BiStability[END_REF], the authors have proposed a numerical scheme, which is well-balanced for a wave with speed σ = 0, and have used it in order to capture traveling waves with nonzero speed. Hence, the proposed scheme is stricto sensu not WB for the observed traveling wave.

In the present article, in order to overcome this difficulty, we implement a strategy, which measures dynamically the spreading speed. We use an estimate of this spreading speed, which has been proposed in [START_REF] Leveque | A study of numerical methods for hyperbolic conservation laws with stiff source terms[END_REF] and has recently been refered to as LeVeque-Yee formula [START_REF] Mascia | Phase transitions of biological phenotypes by means of a prototypical PDE model[END_REF]. This formula is based on the following observation: suppose that u(t, x) = u(xσt) is a traveling wave profile, which admits at ±∞ fixed limit values u(±∞), then u satisfies:

∂ t u + σ∂ x u = 0.
Integrating over R and rearranging the terms, we obtain the following identity:

(1)

σ = - R ∂ t udx u(+∞) -u(∞)
The LeVeque-Yee formula is then merely a discretized version of Identity (1):

(2)

σn LY = ∆x t n -t n-1 i u n-1 i -u n i u n I -u n 0 ,
where u is computed on a Cartesian grid (t n ) n∈N × (i∆x) i=0,...,I . For solutions that are not traveling waves, the right hand-side of Identity (1) remains well defined and can be interpreted as an estimate of their spreading speed.

The contribution of this article consists in combining the LeVeque-Yee formula as a dynamic speed estimate with the methodology of WB schemes for parabolic equations. The schemes presented in this article share the following common structure:

(1) Estimate at each time step the spreading speed σn via the LeVeque-Yee Formula (2).

(2) Consider over the time step the problem in the frame of reference moving at speed σn and use a WB scheme to integrate numerically the solution over the time step. (3) Given the solution integrated in time in the moving frame, shift back the solution to the stationary frame. To do so in a WB manner, interpolate in each cell the stationary solution, with the boundary conditions prescribed by the values obtained at the preceding step. Finally, in order to obtain the values on the grid of the stationary frame, use the value of the obtained stationary solution at the corresponding position. In order to illustrate this methodology, we apply it to reaction-diffusion equations of the type:

(3)

∂ t u -∂ xx u = f (u).
Many theoretical results are known about these types of equations, which will be used as a base of comparison for our numerical results. First, we study the case f (u) = u(1u) of the Fisher/Kolmogorov-Petrovsky-Piskunov (F/KPP) equation [START_REF] Fisher | The wave of advance of advantageous genes[END_REF][START_REF] Kolmogorov | Etude de l équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique[END_REF], which has been a prototype of traveling wave phenomena in parabolic equations. Then, we consider a cubic reaction function of the form f (u) = u(1u)(1 + au) for a parameter a > 0 introduced in [START_REF] Hadeler | Travelling fronts in nonlinear diffusion equations[END_REF]. This class is relevant, because it exhibits a transition from pulled to pushed waves (see Subsection 1.2), which illustrates the behaviors of our numerical scheme on these regimes. Finally, we stress that this methodology is not restricted to reaction-diffusion equations and can be applied to parabolic equations with an advection term, such as for instance models of bacterial chemotaxis (e.g. [START_REF] Keller | Traveling bands of chemotactic bacteria: A theoretical analysis[END_REF][START_REF] Saragosti | Directional persistence of chemotactic bacteria in a traveling concentration wave[END_REF]).

1.1. Fisher/Kolmogorov-Petrovsky-Piskunov Equation. The Fisher/Kolmogorov-Petrovsky-Piskunov (F/KPP) equation ( 4)

∂ t u -∂ xx u = u(1 -u),
introduced independently in [START_REF] Fisher | The wave of advance of advantageous genes[END_REF] and [START_REF] Kolmogorov | Etude de l équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique[END_REF], is a prototype of reaction-diffusion equations. It describes a population u(t, x) of individuals undergoing unbiased motion, modeled through a diffusion operator, as well as competition among the population (e.g. for ressources) through a logistic growth term u(1u): the higher the density u becomes, the lower the growth rate 1u becomes, with a saturation when the population reaches the maximal density u = 1. It admits two steady states u = 0 and u = 1: the first one is unstable, whilst the second is stable. Interestingly, for velocity σ Equation (4) admits nonnegative traveling wave solutions of the form:

(5) u σ (-∞) = 1, u σ (+∞) = 0 and u σ (t, x) = u σ (x -σt)
if and only if σ ≥ σ * = σ F/KP P := 2 (see [START_REF] Fisher | The wave of advance of advantageous genes[END_REF][START_REF] Kolmogorov | Etude de l équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique[END_REF][START_REF] Aronson | Nonlinear diffusion in population genetics, combustion, and nerve pulse propagation[END_REF]). These solutions are in fact invariant by translation and we will fix a specific solution by the convention that u σ (0) = 1 2 . Hence the equation describes a linear spatial invasion of the state u = 0 by the state u = 1. Furthermore for a wide variety of initial data u 0 , one can show that in a certain sense, which we will specify, u the solution of Equation ( 4) converges to the traveling wave with minimal velocity σ * = σ F/KP P . Take for example an initial datum u 0 , which satisfies:

(6) lim inf x→-∞ u 0 (x) > 0 and lim sup x→+∞ ∂ x u 0 (x) u 0 (x) ≤ -1.
Then, define the position of the level set x c (t) for c ∈ (0, 1):

(7) x c (t) = sup{x|u(t, x) = c}
Through a probabilistic interpretation of Equation ( 4) (see [START_REF] Hamel | A short proof of the logarithmic Bramson correction in Fisher-KPP equations[END_REF][START_REF] Nolen | Convergence to a single wave in the Fisher-KPP equation[END_REF] for a proof based on PDE argument), first Bramson [START_REF] Bramson | Maximal displacement of branching brownian motion[END_REF], and later Uchiyama [START_REF] Uchiyama | Brownian first exit from and sojourn over one sided moving boundary and application[END_REF] for inital conditions satisyfing [START_REF] Calvez | Travelling Chemotactic Aggregates at Mesoscopic Scale and BiStability[END_REF], proved that there exists a constant x ∞ depending on the initial datum and c, such that:

(8) x c (t) = 2t - 3 2 ln(t) + x ∞ + o(1) and lim t→+∞ u(t, x + x c (t)) = u σ * (x).
The logarithmic delay in the asymptotic expansion of x c (t) has ever since been referred to as the Bramson delay.

The result is in a sense very surprising, because it shows convergence to the traveling wave with speed σ F/KP P but in the frame shifted by a logarithmic correction term. As it has been observed in [START_REF] Hamel | A short proof of the logarithmic Bramson correction in Fisher-KPP equations[END_REF], this logarithmic correction term is due to subtle interactions between the diffusion and the reaction term, when u ∼ 0. Hence, from a numerical point of view, capturing the Bramson delay is likely a very delicate task, since it requires an accurate computation of the profile of u at very low orders of magnitude. As far as we know, we are not aware of any work on a numerical scheme for Equation ( 4), which addresses the question of the Bramson delay.

1.2. Reaction-Diffusion Equations with a Transition from Pulled to Pushed Waves. In [START_REF] Hadeler | Travelling fronts in nonlinear diffusion equations[END_REF], the authors have introduced the reaction-diffusion equations with a cubic term:

(9) ∂ t u -∂ xx u = u(1 -u)(1 + au).
For this model, the minimal wave speed solutions can be computed explicitly and it can be shown the following formula for the wave speed:

σ * = 2 if a ≤ 2 a 2 + 2 a if a > 2 .
The dichotomy given by this formula extends into a dichotomy on the qualitative behavior of the waves, which is referred to as pulled or pushed regime [START_REF] Stokes | On two types of moving front in quasilinear diffusion[END_REF]. For a ≤ 2, the wave is in the pulled regime: qualitatively the spreading is driven by the dynamics at the very leading edge of the front, where u ∼ 0. For a > 2, the wave is in the pushed regime: qualitatively, the spreading is subject to a contribution from the whole population. These regimes are also reflected in the asymptotic expansion of the level sets x c (t). For a > 2, as a consequence of the exponential convergence to the traveling wave proven in [START_REF] Rothe | Convergence to pushed fronts[END_REF], we have: [START_REF] Fife | The approach of solutions of nonlinear diffusion equations to travelling front solutions[END_REF] x

c (t) = σ * t + x ∞ + o e -ωt ,
for a constant ω > 0. For a < 2, the level sets are again corrected by the logarithmic Bramson delay [START_REF] Hamel | A short proof of the logarithmic Bramson correction in Fisher-KPP equations[END_REF]:

(11) x c (t) = 2t - 3 2 ln(t) + x ∞ + o(1).
The critical case a = 2 has recently been investigated in [START_REF] Giletti | Monostable pulled fronts and logarithmic drifts[END_REF][START_REF] Derrida | Cross-overs of bramson's shift at the transition between pulled and pushed fronts[END_REF] and following the terminology introduced in [START_REF] An | Pushed, pulled and pushmi-pullyu fronts of the Burgers-FKPP equation[END_REF] is referred to as the pushmi-pullyu regime. The level sets are corrected by a logarighmic delay with a constant -1 2 :

(12)

x c (t) = 2t - 1 2 ln(t) + x ∞ + o(1).
1.3. Outline of the paper. In this article, we combine the LeVeque-Yee formula with the methodology of WB schemes for parabolic equations in order to design numerical schemes that are WB for traveling waves.

In Section 2, we present our methodology for WB numerical schemes applied to reaction-diffusion equations.

In Section 3, we numerically assess the WB schemes by applying them to the F/KPP Equation ( 4) and then to the cubic case [START_REF] Fabrèges | WB Scheme for F/KPP Equation[END_REF]. The scheme is compared to an operator splitting approach. We show that our scheme captures accurately the propagation speed and we compare the asymptotic spreading with the logarithmic Bramson correction [START_REF] Després | The structure of well-balanced schemes for Friedrichs systems with linear relaxation[END_REF][START_REF] Fisher | The wave of advance of advantageous genes[END_REF][START_REF] Giletti | Monostable pulled fronts and logarithmic drifts[END_REF] to the asymptotic expansion in the pulled case. As a result we will show that the presented methodology works particularily well in the pulled case, but does not outperform standard schemes in the pushed case.

A WB Scheme for Reaction-Diffusion Equations

In order to solve numerically reaction-diffusion equations (3) in a WB manner, we propose a numerical scheme based on [START_REF]Viscous Equations Treated with $\mathcalL$-Splines and Steklov-Poincaré Operator in Two Dimensions[END_REF] with the difference that we also investigate an implicit Euler method for time integration instead of the standard explicit Euler method. In fact, as it has been pointed out in [START_REF]Viscous Equations Treated with $\mathcalL$-Splines and Steklov-Poincaré Operator in Two Dimensions[END_REF], for the explicit Euler method, one requires essentially a parabolic CFL condition in order to guarantee stability. For the implicit Euler method, through numerical investigation, we observe that a hyperbolic CFL condition is sufficient. Moreover, resolving the time integration implicitly comes computationally at not too high of a cost, because of the special tridiagonal structure of the problem.

We consider a Cartesian grid (t n ) × (x i ), with x i = i∆x, for i ∈ Z for a given ∆x. The time points (t n ) are set dynamically because of the lack of an a priori bound to satisfy a CFL condition. We set ∆t n = t n+1t n . We introduce the parallel grid with points z i = i∆z, ∆z = ∆x, Ci+ 1 2 = (z i , z i+1 ) and consider u n i ≈ u(t n , x i ). The procedure of the numerical scheme can be summarized as follows:

(1) We use LeVeque-Yee Formula (2) on the profile u in order to obtain an estimate σn of the propagation speed.

(2) Given σn , on the time interval (t n , t n+1 ) we consider u in the moving frame

(t, z) = (t, x -σn (t -t n ))
and denote it by ū, which leads to Equation:

(13)

∂ t ū -σn ∂ z ū -∂ zz ū = f (ū).
In order to deal with the nonlinear term in the right hand-side of Equation (3), we freeze the non-linear contribution to the term on each mesh (t n , t n+1 ) × Ci+ 1 2 by considering that f (ū) ≈ ū F n i+ 1 2 , where ūn

i+ 1 2 = u n i +u n i+1 2 , F n i+ 1 2 = F (ū n i+ 1
2

) and F (u) = f (u) u . As noted in the Introduction, the dynamics of the F/KPP Equation (4) or Equation ( 9) with a ≤ 2 is predominantly set by the dynamics at the leading edge of the profile, where u ∼ 0. There, the linearization of Equation (3) approximates extremely well the nonlinear Equation. Hence it is reasonnable to expect that this way of handling numerically the nonlinearity will not influence much the dynamics of the propagation. Thanks to the freezing of the nonlinearity, we obtain a linear problem for which we can apply the methodology of L -spline interpolation, as proposed in [START_REF]Viscous Equations Treated with $\mathcalL$-Splines and Steklov-Poincaré Operator in Two Dimensions[END_REF], leading to the values ūn+1 i , which approximate ū(t n+1 , z i ).

(3) It remains to shift back to the stationary frame (t, x). The values ūn+1 i are used to extrapolate a value for u(t n+1 , x i ), i.e. ū(t n+1 , z i + σn ∆t n ). To do so we consider the solution û of the following stationary problem in each cell Ci+ 1 2 :

(14)    -σn ∂ z û -∂ zz û = û F n i+ 1 2 , for z ∈ Ci+ 1 2 û(z i ) = ūn+1 i and û(z i+1 ) = ūn+1 i+1 , with ûn+1
i+ 1 2 = ūn+1 i +ū n+1 i+1 2 and F n+1 i+ 1 2 = F ûn+1 i+ 1 2
. Finally, we set:

u n+1 i = û(z i + σn ∆t n ).
In order to avoid taking values from cells which would not be direct neighbors, we enforce the following CFL condition:

(15) σn ∆t n ≤ ∆x.

2.1. LeVeque-Yee Formula for F/KPP Equation. We apply the LeVeque-Yee formula to the profile u in order to obtain an estimate for the propagation speed:

(16) σn = ∆x ∆t n-1 i u n-1 i -u n i u n I -u n 0 .
The leftmost (resp. rightmost) value u n 0 (resp. u n I ) will in practice (possibly after a short transitory phase) be very close to the steady state 1 (resp. the steady state 0). Moreover in order to satisfy CFL condition (15), we set dynamically the time points of the grid, in order to avoid imposing a restriction on the speed estimate.

2.2.

A WB Scheme in the Moving Frame. As explained above, we consider the evolution of f in the moving frame (t, z) = (t, xσn ∆t n ), i.e. ū(t, z) = u(t, z + σn (tt n )), where ū solves Equation [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF]. In order to approximate numerically this solution, we use the method of L -spline interpolation proposed in [START_REF]Viscous Equations Treated with $\mathcalL$-Splines and Steklov-Poincaré Operator in Two Dimensions[END_REF]. The scheme is based either on the explicit integration in time rule: [START_REF] Gosse | NUMERICAL HIGH-FIELD LIMITS IN TWO-STREAM KINETIC MODELS AND 1D AGGREGATION EQUATIONS[END_REF]

ūn+1 i = ūn - ∆t n ∆x L n i+ 1 2 -R n i-1 2 ,
or on the implicit integration in time rule: [START_REF] Graham | Precise asymptotics for fisher-kpp fronts[END_REF]

ūn+1 i = ūn - ∆t n ∆x L n+1 i+ 1 2 -R n+1 i-1 2 , where R n+1 i-1 2 , L n+1 i+ 1 2 , R n i-1 2 , L n
-σn ∂ z w -∂ zz w = w F n i+ 1 2 , for z ∈ Ci+ 1 2 w(z i ) = ūn+1 i and w(z i+1 ) = ūn+1 i+1 . Given the solution w, we then set R n+1 i+ 1 2 := w (z - i+1 ) and L n+1 i+ 1 2 := w (z + i ). One can view L n+1 i+ 1 2 -R n+1 i-1 2 (19) (20)
as a measure of the defect of C 1 -smoothness at z i . The time integration (18) will therefore be stationary if ū solves Equation ( 19) on the whole space, which in particular implies that it is C 1 as a solution to a second-order elliptic equation. Thus this procedure is well-balanced in the sense that the numerical scheme admits as stationary state the discretization of the continuous stationary state. For the explicit case, one merely changes the boundary conditions in [START_REF] Greenberg | A Well-Balanced Scheme for the Numerical Processing of Source Terms in Hyperbolic Equations[END_REF] to w(z i ) = ūn i and w(z i+1 ) = ūn i+1 and follows the same procedure. Numerically, we have observed that the explicit time integration is stable under the parabolic CFL condition that:

∆t n ≤ 1 2 ∆x 2
The coefficient 1 2 can be slightly higher, but for 0.6, we have observed instability. This empirical CFL condition is in the same spirit as an explicit CFL condition in [START_REF]Viscous Equations Treated with $\mathcalL$-Splines and Steklov-Poincaré Operator in Two Dimensions[END_REF]. For the implicit time integration, we have observed that under the hyperbolic CFL condition [START_REF]Viscous Equations Treated with $\mathcalL$-Splines and Steklov-Poincaré Operator in Two Dimensions[END_REF], the scheme remains stable.

Let us now describe the implicit (linear) relation between

(L n+1 i+ 1 2 , R n+1 i+ 1 2
) and (ū n+1 i , ūn+1 i+1 ). The same approach applies mutatis mutandis for the explicit case. We solve Problem [START_REF] Greenberg | A Well-Balanced Scheme for the Numerical Processing of Source Terms in Hyperbolic Equations[END_REF] with respect to (ū n+1 i , ūn+1 i+1 ). On the cell Ci+ 1 2 , it is a simple second-order differential equation with constant coefficients. We compute the roots of its characteristic polynomial:

(21) µ 2 + σn µ + F n i+ 1 2 = 0 The discriminant is ∆ n i+ 1 2 = (σ n ) 2 -4 F n i+ 1 2
and three cases exist:

(1) Suppose

∆ n i+ 1 2 > 0, then set µ n i+ 1 2 ,± = -σ n ± ∆ n i+ 1 2

2

. Then for z ∈ Ci+ 1 2 , we have:

w(z) = a -e µ n i+ 1 2 ,- z + a + e µ n i+ 1 2 ,+ z
Therefore by setting the following matrices:

(22) M = e µ n i+ 1 2 ,- zi e µ n i+ 1 2 ,+ zi e µ n i+ 1 2 ,- zi+1 e µ n i+ 1 2 ,+ zi+1 and M =    µ n i+ 1 2 ,-e µ n- i+ 1 2 zi µ n i+ 1 2 ,+ e µ n+ i+ 1 2 zi µ n i+ 1 2 ,-e µ n- i+ 1 2 zi+1 µ n i+ 1 2 ,+ e µ n+ i+ 1 2 zi+1    ,
we obtain the following two systems:

ūn+1 i ūn+1 i+1 = M a - a + and L n+1 i+ 1 2 R n+1 i+ 1 2 = M a - a + .
Hence, by setting S i+ 1 2 = M M -1 , we obtain the following identity:

(23)

L n+1 i+ 1 2 R n+1 i+ 1 2 = S i+ 1 2 ūn+1 i ūn+1 i+1 (2) Suppose ∆ n i+ 1 2 = 0, then we set µ n i+ 1 2 := - σn 2
and the solution is of the shape:

w(z) = (az + b)e µ n i+ 1 2 z .
The definition of the matrix S i+ 1 2 follows mutatis mutandis as precedingly.

(3) In the case ∆ n i+ 1 2 < 0, the solution w is of the shape:

w(z) = a cos -∆ n i+ 1 2 z + b sin -∆ n i+ 1 2 z e -σn z 2 ,
and again the definition of the matrix S i+ 1 2 then follows mutatis mutandis as precedingly. Relation ( 23) can be expressed explicitly, but we will not do so for the sake of concision. Then, combining Relation [START_REF] Kolmogorov | Etude de l équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique[END_REF] with the time integration (18) leads to a tridiagonal system, which can be inverted via Thomas' algorithm.

2.3.

A WB Shift to the Stationary Frame. In the next step, we need to compute the values u n+1 i (in the stationary frame) from the values ūn+1 i (in the moving frame (t, z) = (t, xσn (tt n ))). To do so, we solve Equation ( 14), which we recall:

   -σn ∂ z û -∂ zz û = û F n+1 i+ 1 2 , for z ∈ Ci+ 1 2 û(z i ) = ūn+1 i and û(z i+1 ) = ūn+1 i+1 ,
where we recall that F n+1

i+ 1 2 = F ûn+1 i+ 1 2
is an updated version of the nonlinearity with ûn+1

i+ 1 2 = ūn+1 i +ū n+1 i+1 2
. Then, we set:

(24) u n+1 i = û(z i + σn ∆t n ).
Enforcing the CFL condition (15) avoids taking values from cells which would not be direct neighbors. Equation ( 14) can be solved in the same manner as previously. One takes N i+ 1 2 = M i+ 1 2 and for (φ 1 , φ 2 ) the fundamental system of solutions, that we have used above:

(25) Ñi+ 1 2 = φ 1 (z i + σn ∆t n ) φ 2 (z i + σn ∆t n )
This then leads to the following 1 × 2 matrix:

T i+ 1 2 := Ñi+ 1 2 N -1 i+ 1 2 , and: (26)
u n+1 i = T i+ 1 2 ūn+1 i ūn+1 i+1 .

Numerical Assessements

We start by comparing the preceding WB scheme for F/KPP Equation (4) with two other numerical schemes, which are described just below. A similar comparison is carried out for the Equation [START_REF] Fabrèges | WB Scheme for F/KPP Equation[END_REF]. The code for the numerical schemes in this Section can be found in [START_REF] Fabrèges | WB Scheme for F/KPP Equation[END_REF].

3.1. Alternative Schemes.

2 0 2 1 2 2 2 3 2 4 2 5 2 6 1 ∆ x -4 -3 -2 -1 0 log 10 (|σ -2|)
WB scheme Implicit WB scheme "0-Wave" WB scheme "0-Wave" Implicit WB scheme OS Scheme (Strang splitting)

u n+1 = R ∆t 2 e -∆t∂xx R ∆t 2 u n
For the reaction term, in general we may use an implicit Euler integration, but in the special case for the F/KPP equation, we can get an exact inegration as follows:

R ∆t 2 v = e ∆t 2 e ∆t 2 -1 + 1 v .
The heat operator is integrated via the Crank-Nicolson method with v n+1 = e -∆t∂xx v n satisfying for all i:

v n+1 i -v n i ∆t = v n+1 i-1 -2v n+1 i + v n+1 i+1 2∆x 2 + v n i-1 -2v n i + v n i+1 2∆x 2 .
3.1.2. WB Scheme for the "0-Wave". The second scheme is based on the same method as the WB scheme mentioned before, but we impose that σn = 0. This choice corresponds to the choice that has been investigated in [START_REF]Viscous Equations Treated with $\mathcalL$-Splines and Steklov-Poincaré Operator in Two Dimensions[END_REF]. As mentioned in the Introduction, this approach is only WB for a stationary state, or equivalently a traveling wave with speed 0. For the implementation of this scheme, one does not need to shift from the moving to the stationary frame, as both coincide. Or, with the notations from above, we have:

T i+ 1 2 = 1 0 .

3.2.

The Asymptotic Propagation Speed. The F/KPP Equation ( 4) has an asymptotic progation speed σ of 2. We compare the ability of different schemes to capture this asymptotic velocity: the WB scheme of Section 2 in its explicit and implicit version, the OS scheme of Section 3.1.1 and the WB scheme for the "0-Wave" of Section 3.1.2 also in its explicit and implicit version.

In order to capture the asymptotic propagation speed, the F/KPP equation is solved over the time interval [0, 1500]. Because the expected spreading speed is σ = 2, the domain considered should be at least twice bigger: here, we take the domain [0, 3080]. This large domain avoids roughly speaking boundary effects. The initial condition is the sigmoid function u 0 (x) = 1 -1 1+e -3(x-40) . Figure 1 represents the velocity σn obtained with the LeVeque-Yee formula at the final time t = 1500 for the schemes considered here. With a coarse mesh, the asymptotic propagation speed obtained is far from 2 for some schemes and the spatial domain had to be extended to [0, 6080] in order to avoid boundary effects for these schemes. The same time step is chosen for all the considered schemes at each iteration. It is the one satisfying

2 1 2 3 2 5 2 7 2 9 1 ∆ x -1 0 1 2 log 10 (|Shift -3/2|)
Implicit WB scheme (Hyper.)

Implicit WB scheme (Para.)

OS Scheme (Strang splitting) all the CFL conditions which is the parabolic one, ∆ n t ≤ 1 2 ∆ 2 x , needed for the explicit WB schemes. As the meshes become finer, all the numerical schemes seem to converge to the same result, which is extremely close to the asymptotic propagation speed of 2. Yet, one can observe that the speed does not quite converge to 2. Heuristically, this discrepancy may be a consequence of the Bramson delay. Indeed, differentiating formally the asymptotic expansion (8) leads ẋc (t) = 2 -3 2t + o 1 t . Then plugging in the value t = 1500 would explain the observed discrepancy of the order 10 -3 . Moreover, it can be noted that the WB scheme, both in its explicit and implicit version, captures extremely well even for very coarse meshes the asymptotic behavior that is captured by the other schemes for much finer meshes.

3.3.

Capturing the Bramson Delay. We recall that the Bramson delay is the logarithmic correction that appears in the expression of the position of a level set x c (t) for c ∈ (0, 1) given by Equation [START_REF] Després | The structure of well-balanced schemes for Friedrichs systems with linear relaxation[END_REF]. In the simulations, we track the position of the level set x 0.5 (t) by linearly interpolating the solution in the mesh cell where the level set is. The coefficient 3 2 of the Bramson delay is approximated by fitting the model m α,β,γ (t) = α log(t) + β + γ √ t to x 0.5 (t) -2t over the last half of the simulation. As the next term in the asymptotic expansion (8) is -3 √ π √ t [START_REF]Refined long-time asymptotics for fisher-kpp fronts[END_REF], the γ √ t is added in order to improve the fitting of the model. We try to capture this Bramson delay by refining the discretization of the spatial domain. Using schemes that depend on a parabolic CFL condition results in a huge number of iterations in time which becomes rapidly intractable in a reasonable amount of time. Moreover, the "0-Wave" implicit WB scheme does not capture the asymptotic propagation speed as well as the implicit WB scheme. Thus, we only fit the Bramson delay for the implicit WB scheme and the OS scheme (see Figure 2). The WB scheme and the OS scheme are of order 1 and 2 in time respectively. In order to compare these schemes, we also plot the fitting of the WB scheme using small enough time step to mimic a second order in time scheme. In this case, the time step is ∆ n t = min(0.05, 1 2 ∆ 2 x ). Yet, this parabolic CFL condition again leads to computational costs that are too high, when ∆ x is too small. Thus we have computed the implicit WB scheme with parabolic CFL condition only for ∆ x ≥ 2 -6 . For all computed schemes, the shift obtained seems to converge to a value slightly different from 3/2. One possible explanation for this discrepancy might be higher order terms in the expansion (8) (see for instance [START_REF] Graham | Precise asymptotics for fisher-kpp fronts[END_REF]). Here again, the implicit WB scheme captures the Bramson delay even for coarse meshes. Moreover, for the implicit WB scheme, there is no significant difference between a time step ∆ n t = O(∆ x ) and a time step ∆ n t = O(∆ 2 x ).

3.4. Asymptotic Propagation Speed and Bramson Delay for the Cubic Case. For Equation ( 9), we compare the solutions obtained with the WB and the OS schemes in these three different regimes • pulled regime with a = 1,

• pushmi-pullyu transition regime with a = 2,

2 1 2 3 2 5 2 7 2 9 1 ∆ x -3.0 -2.5 -2.0 -1.5 log 10 (|σ -2|)
Implicit WB scheme (Hyper.)

Implicit WB scheme (Para.)

OS Scheme (Strang splitting)

Figure 3. Asymptotic propagation speed for the cubic case in a pulled regime (a = 1). For the OS scheme, we observe at ∆x = 2 -3 a change in sign of σ -2. 4, and 5, illustrate the convergence of the numerical asymptotic speed in the three regimes. In all cases, the velocities are well captured by the schemes, but the implicit WB scheme with a parabolic CFL performs especially better than the other ones. However, as already mentioned, the computational cost of the parabolic CFL is too high to use it with fine meshes. In the pushed regime (a = 3), it can be noted that for all schemes as the meshes get finer, the speed gets closer and closer to the predicted speed (Figure 5). This is likely due to the exponential convergence in time given by asymptotic expansion [START_REF] Fife | The approach of solutions of nonlinear diffusion equations to travelling front solutions[END_REF], which states that the expansion speed at t = 1500 of the continuous solution should be extremely close to the predicted speed, contrary to what happens in the pulled case due to the Bramson delay. Moreover, both schemes under hyperbolic CFL condition perform roughly equally well in the pushed case. Concerning the Bramson delay, the results for the pulled regime (see Figure 6) are very similar to the ones for the F/KPP equation. In particular, the WB scheme approximates well the shift even for coarse meshes. The Bramson delay of 1/2 in the pushmi-pullyu transition regime (see Figure 7) is captured by all schemes. Yet again, the WB scheme recover the 1/2 shift even for coarse meshes, whilst the OS scheme catches up only for finer meshes.

i+ 1 2 1 2, L n+1 i+ 1 2

 111 correspond to numerical fluxes at the point z = z i of two interpolated functions in Ci-1 2 and Ci+ 1 2 . With the implicit time integration, the values R n+1 i-are obtained by solving the linear stationary problem:

Figure 1 . 3 . 1 . 1 .

 1311 Figure 1. Evolution of LeVeque-Yee velocity σ for the equation at time t = 1500 for different schemes with respect to the space discretization size ∆ x .

Figure 2 .

 2 Figure 2. Fitting of the coefficient 3/2 of the Bramson delay with the implicit WB scheme and the OS scheme for the F/KPP Equation.

Figure 4 .

 4 Figure 4. Asymptotic propagation speed for the cubic case in the pushmi-pullyu transition regime (a = 2). For the OS scheme, we observe at ∆x = 2 -4 a change in sign of σ -2.

Figures 3 ,

 3 Figures 3,[START_REF] Bramson | Maximal displacement of branching brownian motion[END_REF], and 5, illustrate the convergence of the numerical asymptotic speed in the three regimes. In all cases, the velocities are well captured by the schemes, but the implicit WB scheme with a parabolic CFL performs especially better than the other ones. However, as already mentioned, the computational cost of the parabolic CFL is too high to use it with fine meshes. In the pushed regime (a = 3), it can be noted that for all schemes as the meshes get finer, the speed gets closer and closer to the predicted speed (Figure5). This is likely due to the exponential convergence in time given by asymptotic expansion[START_REF] Fife | The approach of solutions of nonlinear diffusion equations to travelling front solutions[END_REF], which states that the expansion speed at t = 1500 of the continuous solution should be extremely close to the predicted speed, contrary to what happens in the pulled case due to the Bramson delay. Moreover, both schemes under hyperbolic CFL condition perform roughly equally well in the pushed case.

Figure 5 .

 5 Figure 5. Asymptotic propagation speed for the cubic case in a pushed regime (a = 3).

Figure 6 .

 6 Figure 6. Bramson delay for the cubic case in a pulled regime (a = 1).

Figure 7 .

 7 Figure 7. Bramson delay for the cubic case in the pushmi-pullyu transition regime (a = 2).

Acknowledgements

The authors are extremely grateful to Vincent Calvez and Thierry Dumont, who have been essential to the present work. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 865711).

Faculty of Science, Koç University, Istanbul, Turkey Email address: mete.demircigil@math.cnrs.fr Institut Camille Jordan (ICJ), Université Claude Bernard Lyon 1, Lyon, France Email address: fabreges@math.univ-lyon1.fr


Guido Schillaci

Alejandra Ciria

Egidio Falotico

Bruno Lara

Cecilia Laschi

Modulating Multi-Modal Integration in a Robot Forward Model for Sensory Enhancement and Self-Perception

Keywords: Internal models, multi-modal integration, perceptual optimisation, sensory enhancement, self-perception

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Prediction has been considered as having a fundamental role in human cognition and in perception, at least since Helmholzt proposal of unconscious inferences [START_REF] Helmholtz | Concerning the perceptions in general[END_REF]. In cognitive robotics, prediction and prediction error are seen traditionally as the key tool for learning and model updating. However, the effects of prediction error on perception has not been widely modeled [START_REF] Lara | Embodied cognitive robotics and the learning of sensorimotor schemes[END_REF]. Few examples link prediction error directly with perception [START_REF] Bechtle | On the sense of agency and of object permanence in robots[END_REF][START_REF] Lang | A deep convolutional neural network model for sense of agency and object permanence in robots[END_REF][START_REF] Schillaci | Body representations for robot ego-noise modelling and prediction. towards the development of a sense of agency in artificial agents[END_REF][START_REF] Cauli | Correcting for changes: Expected perception-based control for reaching a moving target[END_REF]. GS, EF and CL are affiliated with the BioRobotics Institute and the Dept. of Excellence in Robotics & AI, Scuola Superiore Sant'Anna, Pisa, Italy.

AC is affiliated with Facultad de Psicologia, Universidad Nacional Autónoma de México, México City, Mexico.

BL is affiliated with Centro de Investigación en Ciencias, Universidad Autónoma del Estado de Morelos, Cuernavaca, Mexico CL is affiliated with Dept. of Mechanical Engineering, National University of Singapore, Singapore This work presents a multi-modal forward model for a simulated iCub humanoid robot.This model allows predicting and attenuating visually perceived movements, either generated by the robot or produced by other objects in the environment. Moreover, sensory attenuation processes are implemented to enable the robot to maintain the sight of task-relevant objects that could otherwise be occluded by other moving entities. In fact, predictive processes can enable a low-cost enhancement of the perceptual capabilities of the agent. Previous studies [START_REF] Bechtle | On the sense of agency and of object permanence in robots[END_REF][START_REF] Lang | A deep convolutional neural network model for sense of agency and object permanence in robots[END_REF] are extended by including external moving objects during training and testing, and by analysing the impact of these events on sensory attenuation. In addition, the proposed model can modulate the integration of different input modalities (i.e., visual, proprioceptive, and motor command). Inspired by the work of Shim et al. [START_REF] Shim | Robust deep multi-modal sensor fusion using fusion weight regularization and target learning[END_REF] and Patel et al. [START_REF] Patel | Sensor modality fusion with cnns for ugv autonomous driving in indoor environments[END_REF], the proposed model modulates sensory integration by a set of values that dynamically weight the importance of the different modalities for their contribution to visual prediction. The impact of different modulations in the predictive process during an object detection task is analysed using a mechanism for extraction and manual modulation of multi-modal integration weights.

The paper is organized as follows. Sections I-A and I-B frame out this work within the research on sense of agency in human and robots, and on perceptual enhancement in robots. Section II introduces the proposed computational model and the experiments. Results are discussed in section III. Finally, section IV draws the conclusions and outlines venues for further research.

A. Sense of agency and multi-modal integration

During their ontogenetic development, starting in the uterus, foetuses acquire their own body schema, fundamental for their efficient action planning and interaction with the environment [START_REF] Zoia | Evidence of early development of action planning in the human foetus: a kinematic study[END_REF][START_REF] Reissland | The development of anticipation in the fetus: A longitudinal account of human fetal mouth movements in reaction to and anticipation of touch[END_REF]. These schemas have been proposed to be composed of associations of different sensory and motor regions in the brain. In the scientific literature of the last decades, inverse and forward models have emerged as the standard computational candidate to encode them. In the literature on sense of agency, the forward model is also known as the comparator model [START_REF] Haggard | Sense of agency in the human brain[END_REF].Theories claim that, when a motor command is executed, an efference copy is generated to predict the consequences of such action on what is sensed by the organism. The reafference feedback allows the comparison between the predicted sensory consequences and the actual sensory consequences of an action. This comparison provides the capability to distinguish sensory changes consequence of voluntary actions from those generated by external causes. The comparator model theory claims that this comparison is the necessary and sufficient condition for the sense of agency. Under this theory, the resulting magnitude of the prediction error is directly linked to the sense of agency. Prediction error is also thought to be behind sensory attenuation and enhancement depending on the context and the task at hand [START_REF] Juravle | Contextdependent changes in tactile perception during movement execution[END_REF][START_REF] Cardoso-Leite | A new look at sensory attenuation: Actioneffect anticipation affects sensitivity, not response bias[END_REF][START_REF] De Lange | How do expectations shape perception[END_REF]. Furthermore, sensory attenuation has been argued by scholars to be one of the markers of pre-reflective subjective experience [START_REF] Bays | Attenuation of self-generated tactile sensations is predictive, not postdictive[END_REF][START_REF] Schillaci | Exploration be-haviors, body representations, and simulation processes for the development of cognition in artificial agents[END_REF][START_REF] Schillaci | Body representations for robot ego-noise modelling and prediction. towards the development of a sense of agency in artificial agents[END_REF][START_REF] Pyasik | Shared neurocognitive mechanisms of attenuating self-touch and illusory self-touch[END_REF][START_REF] Georgie | An interdisciplinary overview of developmental indices and behavioral measures of the minimal self[END_REF][START_REF] Hafner | Prerequisites for an artificial self[END_REF] and as one of the prerequisites for a sense of object permanence [START_REF] Bremner | Perception of object persistence: The origins of object permanence in infancy[END_REF][START_REF] Bechtle | On the sense of agency and of object permanence in robots[END_REF][START_REF] Lang | A deep convolutional neural network model for sense of agency and object permanence in robots[END_REF].

It has been suggested that reafference can help to understand the origin of a 'self'. Reafference is a concept introduced by von Holst and Mittelstaedt [START_REF] Holst | Das reafferenzprinzip[END_REF] to make a distinction between the sensed consequences of actions and externally caused sensations, or exafference. Reafference is a characteristic of living organisms, offering a mechanism to predict sensory changes consequence of their actions, such as optical flow during motion. Proprioception is a basic form of reafference, which, together with efferent motor signals and visual input, contributes to the visual self-recognition and, thus, to the sense of agency [START_REF] Balslev | Proprioception contributes to the sense of agency during visual observation of hand movements: Evidence from temporal judgments of action[END_REF]. Daprati et al. [START_REF] Daprati | Remembering actions without proprioception[END_REF] suggest that the contribution of proprioceptive feedback for the sense of agency is timeconstrained: while movement is ongoing, a constant monitoring of movement, as well as an actualization of the sense of agency, are needed. It has been proposed that the central nervous system adjusts the weights given to the reafference inflow (e.g. visual and proprioceptive information) and efferent outflow (i.e. motor commands) to estimate the position and velocity of the hand [START_REF] Wolpert | An internal model for sensorimotor integration[END_REF]. Thus, optimal sensorimotor integration for state estimation is based on previous experiences and on the noise or reliability of the signals. This is in line with the idea of a higher weighting of efferent information to rapidly differentiate ourselves from objects and others when visual information is not reliable [START_REF] Salomon | self pop-out": agency enhances selfrecognition in visual search[END_REF]. Here, it is suggested that the sense of agency can be partially explained by the capability to accurately predict self-produced sensory changes. A further clue to the sense of agency can be given by the capability to monitor prediction errors to dynamically assign weights to motor, proprioceptive, and visual signals for optimal state estimation during movement in a task within a context.

In artificial agents, the sense of agency has been modeled mostly using the magnitude of the prediction error as an indicator. In Schillaci et al. [START_REF] Schillaci | Is that me? sensorimotor learning and self-other distinction in robotics[END_REF], an artificial agent observes trajectories performed by itself or by another agent. Given that the agent is better at predicting the consequences of its own movements, the magnitude of the prediction error on the velocity profile indicates the authorship of the trajectories. In line with this idea, in Schillaci et al. [START_REF] Schillaci | Body representations for robot ego-noise modelling and prediction. towards the development of a sense of agency in artificial agents[END_REF], Bechtle et al. [START_REF] Bechtle | On the sense of agency and of object permanence in robots[END_REF], Lang et al. [START_REF] Lang | A deep convolutional neural network model for sense of agency and object permanence in robots[END_REF], prediction error has been used as the main marker for distinguishing between sensory data generated by a humanoid robot and by another entity. A similar paradigm is adopted here, although tested in a more complex experimental setting, where multiple moving objects are present in the scene. In addition, we analyse the dynamic modulation of multi-modal integration weights, and propose this to be also an essential prerequisite for the development of subjective experience in artificial system.

B. Enhancing perceptual capabilities in artificial systems

This work focuses on visual perception in robotics and on the particular issue of visual occlusions. Despite decades of development, occlusion is still a challenge, especially for accurate motion estimation [START_REF] Liu | Selflow: Selfsupervised learning of optical flow[END_REF]. In computer vision, motion is estimated through optical flow calculations. Optical flow can be defined as the flow of visual changes resulting from the relative motion between an observer and the scene, or objects in the scene. Detecting and estimating movements is of crucial importance for any living agent, as well as for robots. Movements can be produced by different sources, including the agent itself, and anticipating them is fundamental for interaction.

Machine learning and computer vision literature on optical flow estimation is substantial. Robotics studies can be found addressing different applications, including navigation [START_REF] Zhang | Flowfusion: Dynamic dense rgb-d slam based on optical flow[END_REF] and object manipulation [START_REF] Arriola-Rios | Modeling of deformable objects for robotic manipulation: A tutorial and review[END_REF]. Visual occlusions have been the subject of interest also for non-engineering sciences. Several related studies can be found in the developmental psychology literature, in particular on how infants deal with visual occlusions and how object permanence is perceived. In fact, the sense of object permanence is considered to be a milestone in infant development [START_REF] Bremner | Perception of object persistence: The origins of object permanence in infancy[END_REF]. Likely occurring already during the firsts months of life, it is through the acquisition of the knowledge that objects persist in the environment, even if they are occluded, that infants start to form internal representations of objects and to handle them. Very early stages of development are already characterised by predictive capabilities, which infants typically exhibit through anticipatory gaze behaviours [START_REF] Applin | Six-month-old infants predict agents' goal-directed actions on occluded objects[END_REF] Studies show that infants' knowledge of the permanence or non-permanence of objects is embodied in their predictive tracking [START_REF] Bertenthal | Phenomenal permanence and the development of predictive tracking in infancy[END_REF].

Few robotics studies have linked the investigation on movement estimations with the cognitive development of a sense of object permanence. Here and in previous work [START_REF] Lang | A deep convolutional neural network model for sense of agency and object permanence in robots[END_REF][START_REF] Bechtle | On the sense of agency and of object permanence in robots[END_REF], it is argued that such a capability may be in part the result of predictive and sensory attenuation processes. Attenuating expected information -such as visually detected self-generated movements -could help maintaining the focus of attention on an object of interest. In this work, a new end-to-end multisensory model for enhancing robot perpeption and that deals with self-generated body movements occlusions or body self-occlusions, as well as with movements generated by other objects in the environment, is presented.

II. METHODOLOGY

We propose a multi-modal forward model for predicting changes in the visual field. The model is learnt with sensory data coming from a humanoid robot and is used for sensory enhancement and self-perception. The proposed framework has two distinctive features: first, a mechanism for attenuating predicted visual changes that are either generated by the robot itself or by other objects in the scene; second, the capability to learn the influence that each sensory and motor input has in the prediction and how this is modulated to perform selective sampling of sensory information.

A similar mechanism -i.e., precision weighting -has been proposed in the predictive processing literature (see [START_REF] Ciria | Predictive processing in cognitive robotics: a review[END_REF] for a review). In precision weighting, information incoming from different modalities is weighted according to the expected confidence given a certain task in a specific context [START_REF] Parr | Working memory, attention, and salience in active inference[END_REF]. The architecture proposed here is capable of learning such weights and of modulating them. This is performed through a mechanism inspired on features extraction processes and typical of autoencoder neural networks.

The model architecture is inspired by the works of Shim et al. [START_REF] Shim | Robust deep multi-modal sensor fusion using fusion weight regularization and target learning[END_REF] and Patel et al. [START_REF] Patel | Sensor modality fusion with cnns for ugv autonomous driving in indoor environments[END_REF], which combine convolutional networks with sensor fusion techniques and fusion weights regularization through auxiliary uni-modal branches for classification purposes. In their NetGated architecture [START_REF] Patel | Sensor modality fusion with cnns for ugv autonomous driving in indoor environments[END_REF], authors employ a gating and sensor fusion approach to integrate different modalities.The authors aim at preventing sensor failures during a classification task, as interference of a noisy channel can be gated off when the corresponding fusion weight is suppressed. A late fusion approach is adopted, where low-dimensional features extracted through uni-modal convolutional branches are integrated. Fusing the different modalities consists of multiplying the uni-modal features with the corresponding fusion weights. Such fusion weights are learned through a further extraction process, which is applied on each uni-modal low-dimensional feature.

However, as noted by Shim et al. [START_REF] Shim | Robust deep multi-modal sensor fusion using fusion weight regularization and target learning[END_REF], similar gating architectures can lead to inconsistency in the sensor fusion weights. NetGated may develop unstable fusion weights, which do not well reflect the relevance of the corresponding modalities in the specific task. To overcome this issue, Shim et al. [START_REF] Shim | Robust deep multi-modal sensor fusion using fusion weight regularization and target learning[END_REF] propose a series of variations to the original NetGated architecture.Among the proposed extensions, there is the introduction of auxiliary uni-sensory processing paths during training, and of fusion weight regularisation using the losses generated by such auxiliary branches.

We adopt a similar approach and extend it. First, we develop a forward model which integrates different sensory inputs and motor commands using a gating system as proposed by [START_REF] Shim | Robust deep multi-modal sensor fusion using fusion weight regularization and target learning[END_REF], and that produces predictions about changes in the visual input. Second, we introduce the possibility to modulate the fusion weights, in fact regulating the integration of the different inputs during a prediction step. This is achieved by overriding the values of the multi-modal integration weights resulting from the model training with those of additional inputs passed to the model, when needed. Moreover, model predictions are combined with upcoming visual inputs, and used in sensory attenuation processes. The forward model is characterised by the following inputs and outputs.

Inputs:

• visual input V recorded at time (t), consisting of a 32×32 pixels grayscale image captured from the left camera of a simulated iCub robot; • proprioceptive input P recorded at time (t), consisting of the absolute positions of the 16 joints of the left arm and hand (including fingers) of the iCub robot; • motor command M applied at time (t), consisting of a 16-D vector representing the target positions for the same left arm and hand joints. Outputs:

• the magnitude of the optical flow (OF) computed between the visual inputs at time (t) and at time (t+1). The resulting OF image consists of a 32 × 32 pixels matrix. In the model structure, depicted in Figure 1, each input signal flows through a sequence of layers in the respective uni-modal branch. The three uni-modal branches differ from the type of processed input. The visual branch consists of a sequence of convolutional, dropout (rate: 40%) and maxpooling layers terminating on a feature layer of 256 neurons. The proprioceptive and motor branches consist of sequences of dense and dropout (rate: 40%) layers, terminating in two feature layers of the same dimension as the visual one.

Input signals flow through different paths in the feature layers. The path extracting the multi-modal integration weights concatenates, first, the three uni-modal 256-D feature layers. The concatenated tensor is reduced through a series of dense layers to a 3-dimensional tensor, which is then processed by a SoftMax layer. This ensures that the output activation of the three neurons sum up to one. The resulting layer is split into three single neurons, each representing the weight of the corresponding modality for the multi-modal integration. Finally, regularisation of the weight activity is performed. Multi-modal integration weights are multiplied to the unimodal 256-D features extracted earlier, modulating, in fact, the output of each uni-modal branch. At this point, weighted multi-modal integration is performed by summing up the incoming signals. The result of the multi-modal integration is thus passed through a sequence of dense, reshape and upsampling layers, until the desired shape [START_REF] Holst | Das reafferenzprinzip[END_REF][START_REF] Holst | Das reafferenzprinzip[END_REF] of the model output is achieved.

Importantly, the three uni-modal branches (from the input layers to the 256-D feature layers) follow also a parallel, auxiliary path. In particular, each uni-modal input is mapped to an exact copy of the model output layer (layers and connections between them are shared with the main model). This allows to calculate auxiliary uni-modal losses and thus to estimate how efficient is every single modality in generating the desired output. Hence, the full architecture is characterised by four outputs: the main model output and three uni-modal branch outputs. During training, the same optical flow information is passed to the four outputs of the model.

Model optimisation is performed through an Adam optimiser [START_REF] Kingma | A method for stochastic optimization[END_REF]. A loss function that combines the main model loss with a weighted average of the auxiliary losses is adopted. In particular, the following loss function has been implemented:

Loss = Loss main + α • N n=1 w n • Loss n aux (1)
where Loss main is the loss of the main model, α is a constant (set to 0.5 in the experiments presented here), N is the number of modalities (i.e., visual, proprioception and motor), w is the weight used for multi-modal integration and Loss n aux is the auxiliary loss of the modality n. Each loss is computed as the mean squared error between the targets and the predicted outputs. In addition, activities of the three multi-modal integration weights are regularised as follows:

w n = w n -β • (w n -e -(Loss n aux ) 2 ) ( 2
)
where β is a constant set to 0.05 in the experiments presented here, and the exponential is passed through a sigmoid function. These equations are adapted from [START_REF] Shim | Robust deep multi-modal sensor fusion using fusion weight regularization and target learning[END_REF].

A. Sensory attenuation

We propose a mechanism that combines forward model predictions and upcoming visual inputs with the aim of attenuating visually perceived movements.Given the visual and proprioceptive inputs recorded at time (t) and a motor command, the forward model generates an optical flow (OF) prediction, i.e., an prediction of the movements to be perceived in the upcoming visual input. The predicted and the observed magnitudes of the OF images are then binarised applying a threshold to these values. The binarised OF prediction is used as a mask for the upcoming visual input, V (t+1). In particular, pixels in V (t + 1) whose coordinates correspond to those OF pixels that have maximum magnitude values are attenuated.

Here, attenuation is performed simply by substituting such pixels with those of a background image that has been stored at the beginning of the experiment.

In the experiments presented below, we measure sensory attenuation in an object detection task. In the experimental setup, a set of nine objects are placed in front of the robot. These objects are tagged with fiducial markers to ease their detection, as illustrated in Figure 2. We measure the quality of sensory attenuation by comparing the number of objects detected in the observed image -where occlusions due to the movements of the robot or of other objects may occur -with those detected in the image where predicted changes in the image are attenuated. Performances under different training and test configurations are examined, as well as under different modulations of the multi-modal integration weights.

B. Design of the experiments

A series of experiments has been carried out in which we analyse different aspects of the proposed model. First, we trained and tested the model on different datasets, analysing its learning performance under the different configurations. Second, we tested the capability of the proposed system to enhance robot perception in an object detection task (as described in Section I-B), as well as its capability to characterise self-generated movements by means of prediction errors. Finally, we measured sensory attenuation and object detection performance under different modulations of the multi-modal integration weights.

Three datasets have been collected for training and testing the model. Each sample of the datasets is characterised by visual, proprioceptive and motor information grabbed at a time (t), and the optical flow estimated between the visual inputs at time (t) and (t + 1). Datasets have been collected under three different experimental conditions simulated in the iCub humanoid robot simulator:

• DS1: the iCub robot, alone in the scene, randomly explores the movements of its left arm, hand and fingers. • DS2: the robot is moving as in DS1, while coloured balls fall from the sky from random positions. A new ball is created every 20 seconds. • DS3: same conditions as in DS2, but with more simulated objects (a new ball is created every second). A set of nine static objects was also present in the scene, as depicted in Figure 2, during all the data collection processes. Objects were tagged with fiducial markers to ease their detection. Each dataset consisted of 25000 samples, gathered with a sampling rate of 5 samples per second.

We designed a set of experiments combining the use of the described sample collections either as training or test datasets, as depicted in Table I Experiments 1, 2 and 3 have been carried out to test whether prediction errors can be used as a means to detect changes in sensory information that are not produced by the robot itself. In particular, three models are trained with the same dataset (DS1: iCub alone in the scene), however each one of them is tested with a different test dataset (DS1 in experiment 1; DS2 in experiment 2; DS3 in experiment 3). We expect prediction errors to be, in average, higher in exp. 2 and 3where moving objects are present in the scene -than in exp. 1. This study extends the experiments presented by Bechtle et al. [START_REF] Bechtle | On the sense of agency and of object permanence in robots[END_REF] and Lang et al. [START_REF] Lang | A deep convolutional neural network model for sense of agency and object permanence in robots[END_REF] by adopting a more complex computational model and experimental setup. In addition, a sensory attenuation and object detection experiment has been carried out for each of the models trained as in Table I. We test the contribution of sensory attenuation processes in enhancing robot perception: attenuating self-generated movements, or those generated by other entities, allows detecting objects even if they are occluded in the raw visual input. In experiments 4 and 5 (see Table I), the proposed forward model has been exposed, sequentially, to sensorimotor data belonging to different training datasets. Each of these experiments consisted of six training phases. The training dataset was changed for each phase. The aim of these experiments is to analyse the re-modulation of the multi-modal integration weights of the model during the learning process. We expect the model to adapt to such dynamic circumstances, and to re-modulate its MIW when visual information becomes too unreliable to be used for predicting upcoming sensory input, as occurring in DS3. Finally, experiments 6 and 7 have been carried out using both DS1 and DS3 as training datasets.

Each experiment has been run 10 times. A run consisted of a 10-epochs offline training session, where batches of 32 samples were processed in each training iteration. Moreover, in each experiment we quantify the impact of multi-modal integration weights modulation onto the attenuation process. In particular, we analysed the object detection performance of each model on images of the corresponding test dataset, under the following configurations: (i) no sensory attenuation is carried out; (ii) sensory attenuation is performed using the model's prediction -i.e, using the multi-modal integration weights resulting from the model training (MIW from now); (iii) sensory attenuation is performed using the model's prediction resulting from manually modulating the multi-modal integration weights. In this last configuration, six different sets of weights have been tested, as described in Table II. Object detection performance has been measured as the average number of objects detected in the given image dataset. The results of the experiments are described in the following section. Datasets, scripts for setting up the experimental environment and source code are openly available 2

III. RESULTS

A. Experiments 1, 2 and 3

Experiments 1, 2 and 3 have been carried out to test whether sensory attenuation is better when objects are occluded by selfgenerated changes in the visual input than when they are due to other external moving objects (falling balls). Three models are trained with the same dataset (DS1: iCub alone in the scene), and tested with DS1, DS2 and DS3, respectively. The mean validation loss for each of the 10 epochs of the 10 runs has been compared over the experiments (see Fig. 3).

We used a one-way Welch's Heteroscedastic F Test with Trimmed Means and Winsorized Variances for independent samples to test for equal means (Welch's Test in what follows). The Welch's Test showed a significant difference in the validation loss across the experiments 1, 2 and 3, F (2, 145.438) = 3846.481, p < 0.0001. The post-hoc paired comparisons using a Holm correction showed a significant difference in validation loss between exp. 1, M = 4.69e -05 , SD = 1.15e -05 , and exp. 2 (DS2, M = 8.31e -05 , SD = 1.47e -05 , p < 0.0001), between exp. 1 (DS1) and exp. 3 (DS3, M = 0.0003769, SD = 2.87e -05 , p < 0.0001), and between exp. 2 (DS2) and exp. 3 (DS3, p < 0.0001). In line with the theories presented in section I and the results presented by Bechtle et al. [START_REF] Bechtle | On the sense of agency and of object permanence in robots[END_REF] and Lang et al. [START_REF] Lang | A deep convolutional neural network model for sense of agency and object permanence in robots[END_REF], prediction errors were significantly lower in exp. 1 where the iCub was alone in the scene than exp. 2, where there were also other moving objects, and even to a greater extent when the rate of the presence of other moving objects increased, as occurred in exp. 3. In other words, average prediction error was significantly lower when the iCub was observing only its own movements.This supports the idea that prediction errors of a multi-modal forward model can be used as a cue for detecting activities performed by other individuals and, ultimately, for distinguishing self-generated sensory consequences from those generated by others. Figure 4, shows the results of the object detection tests for exp. 1, 2 and 3. Each sub-figure shows a series of plots illustrating the mean number of objects detected in the original images and in those processed after sensory attenuation for each WS. In general, as can be seen in Figure 4, a perceptual enhancement of the objects occurred using the MIW resulting from the model training, as well as with the six custom multimodal integration weights sets (WS0-5), in comparison to the number of objects detected where no sensory attenuation occurred (i.e, baseline condition). In addition, an object detection statistical analysis has been carried out, where optical flow predictions were generated using the MIW, as well as using the WS0-5 (Table I). We tested the contribution of the sensory attenuation processes in enhancing the robot object perception. Specifically, we tested if attenuating visually perceived movements, either selfgenerated or generated by other entities, allows detecting objects even if they are occluded in the raw visual input. The quality of sensory attenuation and the object detection performance under different configurations and multi-modal integration weight modulations was measured and analyzed. The mean of the total 'object detected' without any sensory attenuation was used as a baseline, and the differences on the 'objects detected' was calculated for each experiment. Thus, the data used for the statistical analysis was the result of subtracting the observed mean of each of the ten epochs of the ten total runs minus the corresponding baseline mean of each experiment. A total of 100 mean differences were calculated for each experimental condition, the MIW, and each one of the six WS0-5 (Table II), in each one of the three experiments. II). The blue line over training epochs, represents the baseline of objects detected, which is the mean of objects detected in the test image dataset for each experiment. The baseline is constant, as no sensory attenuation occurred and occlusions, due to robot movements or objects motion, may have occurred.

In exp. 1, the Welch's test showed a significant difference in the objects detected between the conditions F(6, 245.343) = 9.652, p < 0.0001. The post-hoc paired comparisons using a Holm correction showed a significant difference in the objects detected between the MIW (M = 0.434, SD = 0.222) and the WS5 weights set (M = 0.276, SD = 0.169). This suggests that the dynamic MIW outperformed the WS5 weights set (vision: 0.1, proprioception: 0.3, motor: 0.6). There were no significant differences in the objects detected between the MIW and the other remaining five different custom multi-modal integration weights sets. 3 Contrary to exp. 1, the Welch's test in exp. 2 (F (6, 245.350) = 1.524, p = 0.17), as well as in exp. 3 (F (6, 245.677) = 0.114, p = 0.99), showed no significant differences in the objects detected between the conditions.

From the trends depicted in Figure 4, it seems that the more the visual modality becomes noisy the more the model tends to rely to a greater extent on the motor modality. This trend could explain why, in exp. 1, a statistically significant difference was observed in the comparison between the object detection performance of the MIW and that of the model using WS5. This suggests that, when the visual modality is reliable, up-weighting the motor modality to predict the optical flow produces a worse sensory attenuation -and thus worse object detection -performance. Although the rest of the comparisons were not statistically significant, promising trends can be observed in Figure 4, where in exp. 2 and 3, the performances of WS5 becomes closer to those of the main model, suggesting that the MIW increased to some extent the weights assigned to the motor modality.

Finally, two qualitative analyses on how the system modulates, during testing, the dynamic multi-modal integration weights were performed. The first analysis used a sample sequence of the iCub moving its arm occluding the objects (Figure 5), and the second, using a sequence of the iCub hand leaving the field of view of the robot camera (Figure 6). The aim of these qualitative analyses was to elucidate how the visual, proprioceptive, and motor weights were dynamically adjusted as a result of an incremental occlusion of the objects caused by a self-generated movement, as well as the opposite situation. First, as the robot arm incrementally occluded the objects throughout the movement, the proprioceptive and motor modalities weights were dynamically adjusted to be more influential for predicting the optical flow (see third row of Figure 5). It is plausible to think that, given the ample movement of the robot arm, proprioceptive and motor information needed to be highly weighted to correctly predict the OF product of self-generated movements. In the second test, the iCub hand is leaving the field of view of the camera. At the beginning of the sample sequence (see Figure 6), the visual weights increased because the hand is almost out of the field of view. Then, the hand appeared again in the field of view and as a consequence the multi-modal weights were almost equally weighted. However, while the hand disappeared from the visual field throughout the sample sequence, the visual weights remained high.

B. Experiments 4 and 5

The aim of experiments 4 and 5 was to analyse the potential re-modulation of the MIW during the learning process, changing training dataset over six phases during the experiment. DS3, andDS3). Additionally, in Figure 7 the validation loss is shown for the dataset tests DS1 and DS3 used in exp. 4 and in exp. 5, respectively. Interestingly, even when both experiments were trained using the same sequences of datasets within phases, in exp. 4 the validation loss outstandingly outperforms the validation loss obtained in experiment 5. These findings can be interpreted as the model being able to learn and dissociate the visual consequences of self-generated movements and thus accurately predict the OF of these changes in all the test phases. In exp. 5, the validation loss curve shows that the model was not as good in predicting the optical flow caused by other sources than the robot. Although in phases 1, 4, and 5, the validation loss is lower than the loss during training, the prediction error is considerably higher in these phases than in exp. 4. Moreover, in exp. 5, the validation loss curve bumps within phases and only showed a tendency to decrease in the first two phases. This further supports the idea that prediction errors of a multi-modal forward model can be used as a cue for distinguishing self-generated sensory consequences from those generated by others, also in dynamic contexts. As expected, the model was able to adapt to dynamic circumstances, for instance, to re-modulate its weights when visual information became too unreliable to be used for predicting upcoming sensory input, as when being trained with DS3.

Figure 8 shows the average object detection rates in exp. 4 and 5. In comparison with the baseline, an enhancement of the robot visual perceptual capabilities when applying predictive and sensory attenuation processes can be seen in all the experimental conditions. Additionally, a tendency for switching the MIW can be seen depending on the dataset used. In both tests using the DS1 (exp. 4) and DS3 (exp. 5), the WS with high vision weights, such as WS0 and WS1, tended to decrease their performance in phases 1, 4, and 5, where DS3 (whose samples have noisy visual inputs) was used during training. On the contrary, WS4 and WS5, i.e. the sets with high proprioceptive and motor weights, tended to perform better in those phases. As in previous experiments, a statistical analysis on the mean object detection has been carried out. For each phase, a total of 100 mean differences were calculated for each experimental condition, the MIW and each WS. In exp. 4, in its phase 0, the Welch's Test showed a significant difference in the objects detected between experimental conditions F(6, 244.968) = 3.039, p = 0.009. Although there were significant differences between some of the post-hoc paired comparisons, none of them were between the dynamic MIW and the six WS shown in Table III. In phase 1, there were significant differences F(6, 242.956) = 32.390, p < 0.0001. The post-hoc paired comparisons showed that the dynamic MIW significantly outperformed all the six WS. In phase 2, there were also significant differences F(6, 243.565) = 31.337, p < 0.0001. However, contrary to phase 0 and 1, in the post-hoc paired comparisons, the WS2 weights set (M = 0.574, SD = 0.073) performed equally well as the dynamic MIW (M = 0.571, SD = 0.072; p = 1.0). The same pattern of results were observed in phase 3 (F(6, 243.508) = 26.241, p < 0.0001), the WS2 weights set (M = 0.584, SD = 0.090) performed equally well as the dynamic MIW (M = 0.579, SD = 0.090; p = 1.0), but here, also the WS4 weights set (M = 0.547, SD = 0.109; p = 0.14). Finally, in phase 4 and 5 significantly differences were found F(6, 239.687) = 40.078, and F(6, 238.417) = 36.647, respectively. In both phases, the post-hoc paired comparisons showed that the WS5 weights set (phase 4: M = 0.606, SD = 0.086; phase 5: M = 0.606, SD = 0.092) performed equally well as the dynamic MIW (phase 4: M = 0.571, SD = 0.072; phase 5: M = 0.594, SD = 0.071; p = 1.0).

Table III relates to exp. 4 and shows that the model readapted its MIW when switching training dataset. In phase 2 and 3, the model adjusted its MIW to perform as in those weights sets where proprioceptive and motor modalities have more influence than the visual one. Yet, the visual modality still has an influence in predicting on DS1. When the visual modality becomes too noisy during training (DS3), as in phases 4 and 5, the model seems to rely more on the motor modality. This trend can explain why, in phase 4 and 5, the WS5 weights set performed equally well as the dynamic MIW to predict on DS1. In exp. 5, the Welch's test showed that the difference in the objects detected between experimental conditions in phase 0 was not statistically significant (F(6, 245.702) = 0.946, p = 0.46). Therefore, post-hoc paired comparisons were not performed, nor reported in Table IV. In phase 1, there were statistically significant differences (F(6, 238.699) = 20.872, p < 0.0001). The post-hoc paired comparisons showed that the dynamic MIW (M = 1.399, SD = 0.473) significantly outperformed all the WS, except for the WS5 weights set (M = 1.43, SD = 0.482; p = 0.17). In phase 2, there were also statistically significant differences (F(6, 245.274) = 2.804, p < 0.0001). Although there were significant differences between some of the post-hoc paired comparisons, none of them were between the dynamic MIW and the custom WS as shown in Table IV. In phase 3, the difference was not statistically significant (F(6, 245.422) = 0.343, p = 0.91), as in phase 0, hence post-hoc paired comparisons were not performed. Finally, in phase 4 and 5, statistically significantly differences where found F(6, 237.872) = 18.343, p < 0.0001, and F(6, 235.715) = 25.749, p < 0.0001, respectively. As in phase 1, in phase 4, the dynamic MIW (M = 1.417, SD = 0.477) significantly outperformed all the WS, except for WS5 weights set (M = 1.445, SD = 0.488; p = 1.0). However, in phase 5, although the dynamic MIW (M = 1.418, SD = 0.478) significantly outperformed all the WS, the WS5 weights set was outperformed the dynamic MIW (M = 1.453, SD = 0.489; p = 0.01).

Exp 5

Ph The aim of experiments 6 and 7 was to test if the dissociation between visual changes produced by self-generated movements and changes that are externally produced in changing contexts can occur when different datasets are combined in the same training phase. This experimental manipulation can be understood as being exposed to a noisy context, where abrupt visual external changes can suddenly occur, and the occurrence of these changes is not predictable. For the simulation of this type of context, DS1 and DS3 were combined and randomly presented during training. DS1 was adopted as test dataset for exp. 6 and DS3 as test dataset for exp. 7. Figure 9 shows the trends of the loss and the validation loss of the model under the different experimental conditions, as depicted in Table I. The Welch's test showed that the difference in the validation loss between exp. 6 and 7 was statistically significant (F(1, 79) = 209.645, p < 0.0001). As expected, the prediction error was significantly lower when the iCub was observing its own movements (DS1) than when there were other objects moving in the environment (DS4), despite during training the model was exposed to both datasets equally likely. This finding can be interpreted as the model being able to dissociate the selfgenerated visual changes from the externally visual changes and to accurately predict the resulting optical flow of its movements. Figure 10 shows the results of the mean object detection tests for exp. 6 and 7, depicted in Table I. As in the previous experiments, an enhancement of the robot visual perceptual capabilities for object detection can be seen in all the experimental conditions of both experiments in comparison with the baseline. Results of the statistical analyses of the comparisons of object detection under the different configurations are reported in Table V In exp. 6, the Welch's test showed a statistical significant difference in the objects detected between experimental conditions (F(6, 241.557) = 4.372, p = 0.0003). The posthoc paired comparisons indicate that the WS2 weights set (M = 0.495, SD = 0.115; p = 0.51), the WS4 weights set (M = 0.511, SD = 0.097; p = 0.91), and the WS5 weights set (M = 0.511, SD = 0.101; p = 0.31), performed equally well as the dynamic MIW (M = 0.534, SD = 0.071). In exp. 7, a similar pattern of results was observed 4 .

Therefore, in exp. 6, the dynamic MIW tended to rely mainly on proprioceptive and motor modalities, but also recruited the visual modality, as for the WS2 weights set. Interestingly, in exp. 7, where DS3 was used as test dataset, the model tended to rely on proprioceptive and motor modalities to a greater extent, likely due to the fact that the visual modality was not as reliable as in exp. 6. Therefore, when the model was trained under a noisy context, where abrupt visual external changes occurred suddenly, and the occurrence of these changes were not predictable, the multi-modal trained weights tended to be higher on the proprioceptive and motor modalities. However, it is important to highlight that during training there were also data examples that were not noisy and only self-generated visual changes occurred, and even in 4 A statistical significant difference occurred among experimental conditions (F(6, 243.472) = 19.166, p < 0.0001), but here, only the WS4 weights set (M = 1.388, SD = 0.285; p = 0.095), and the WS5 weights set (M = 1.436, SD = 0.292; p = 1.0), performed equally well as the dynamic MIW (M = 1.441, SD = 0.282). some examples non visual changes occurred (DS1). In this examples, the model should up-weight the visual modality as can be seen in Figure 6, which can explain why the model, when being tested with DS1 performed equally well as the WS2 weights set.

IV. CONCLUSIONS AND FUTURE WORKS

This work has presented a computational model for multimodal integration in a simulated humanoid robot. The model allows anticipating and attenuating movements from the visual input that are either self-generated by the robot or generated by other objects. Thanks to these processes, the visual perceptual capabilities of the robot can be enhanced. Here, sensory enhancement has been tested on a experiment on object detection and object permanence: through anticipation and attenuation processes, objects that are originally occluded by self-generated movements or those generated by other entities can be maintained in the field of view of the robot. The results presented in this study are in line with the idea that a sense of object permanence, an important capability acquired by infants during early developmental stages, may rely in part on predictive processes. Moreover, the prediction errors generated by the computational model under study can be used to characterise parts of the visual inputs as produced by the robot itself or by other objects as demonstrated by exp. 1, 2, and 3. This is a fundamental capability for self-perception and for self-other distinction.

An important feature of the computational model proposed here is the capability to extract multi-modal integration weights. This is the capability to learn how influential each input modality is in the process of predicting an upcoming visual input. We extended a model proposed by Shim et al. [START_REF] Shim | Robust deep multi-modal sensor fusion using fusion weight regularization and target learning[END_REF] and Patel et al. [START_REF] Patel | Sensor modality fusion with cnns for ugv autonomous driving in indoor environments[END_REF], providing the possibility to manually modulate the multi-modal integration weights, and analysing the model's capability to learn these features and to adapt them to dynamic environmental circumstances.

The observed trends in the modulation of multi-modal integration suggest different interesting venues for further research. The current process for attenuating expected movements from the visual input is based on the binarisation of the magnitude of the optical flow. This introduces an additional error that may bias the attenuation and object detection performances -which depend on the chosen threshold. Further work should test a more robust and comprehensive algorithm for combining predicted and observed optical flows and for allowing also partial attenuation. The direction of the movements could also be taken into account in the predictive model. Moreover, a better mechanism for encoding the background images should be considered. Convolutional networks could be used to encode the current background and context into a short-term and working memory representation. This would allow also extending the algorithm to support head movements, which are likely to produce dense optical flow information over the whole visual input. Moreover, measures to allow an incremental learning of the model and to prevent catastrophic forgetting issues in the dynamic context tests should be adopted. More complex object detection scenarios could be also developed. Nonetheless, multi-modal integration should be also validated on an experiment that is agnostic to the object detection task. For instance, the comparisons between the multi-modal integration using weights resulting from the model training and the custom multi-modal integration weights could be also tested on the raw prediction error of the model. Moreover, custom multi-modal integration weights, not fixed for the entire duration of the test sequence, should be tested.

Fig. 1 .

 1 Fig. 1. An illustration of the implemented forward model and a screenshot of the humanoid robot iCub.

Fig. 2 .

 2 Fig. 2. The set of objects positioned in front of the robot.

  DS3 DS3 TABLE I: Design of experiments. Legend: DS1: iCub alone in the scene; DS2: iCub and balls falling in the scene; DS3: iCub and many balls falling in the scene.

Fig. 3 .

 3 Fig. 3. Loss (blue plots) and validation loss (orange plots) of exp. 1, 2 and 3.

Fig. 4 .

 4 Fig. 4. Mean and standard deviation of objects detected in experiments 1, 2 and 3 over training epochs and applying sensory attenuation using the MIW and W0-5 (TableII). The blue line over training epochs, represents the baseline of objects detected, which is the mean of objects detected in the test image dataset for each experiment. The baseline is constant, as no sensory attenuation occurred and occlusions, due to robot movements or objects motion, may have occurred.

Fig. 5 .

 5 Fig. 5. A sample sequence showing the dynamics of the MIW while the robot hand is entering the field of view of the robot camera. The first row shows the observed images v(t + 1); the second row shows the observed OF; the third row shows the MIW; the fourth row shows the OF predicted by the model; finally, the fifth row shows the image resulting from the sensory attenuation.

Fig. 6 .

 6 Fig. 6. A sample sequence showing the dynamics of the MIW while the robot hand is leaving the field of view of the robot camera.

Fig. 7 .

 7 Fig. 7. Loss and validation loss of experiments 4 and 5.

Figure 7

 7 Figure 7 shows the trends of the loss of the model using MIW in exp. 4 and 5, in each one of the six training phases with their respective training dataset (DS1, DS3, DS1, DS1, DS3, and DS3). Additionally, in Figure7the validation loss is shown for the dataset tests DS1 and DS3 used in exp. 4 and in exp. 5, respectively. Interestingly, even when both experiments were trained using the same sequences of datasets within phases, in exp. 4 the validation loss outstandingly outperforms the validation loss obtained in experiment 5. These findings can be interpreted as the model being able to learn and dissociate the visual consequences of self-generated movements and thus accurately predict the OF of these changes in all the test phases. In exp. 5, the validation loss curve shows that the model was not as good in predicting the optical flow caused by other sources than the robot. Although in phases 1, 4, and 5, the validation loss is lower than the loss during training, the prediction error is considerably higher in these phases than in exp. 4. Moreover, in exp. 5, the validation loss curve bumps within phases and only showed a tendency to decrease in the

Fig. 8 .

 8 Fig. 8. Means of object detection in experiment 4 and 5 using DS1 and DS3 respectively.

Fig. 9 .

 9 Fig. 9. Loss and validation loss of experiments 6 and 7.

Fig. 10 .

 10 Fig. 10. Comparison of the object detection tests between experiments 6 and 7.

  .

		visual	proprio	motor
	WS0	0.6	0.3	0.1
	WS1	0.5	0.3	0.2
	WS2	0.3	0.4	0.3
	WS3	0.45	0.1	0.45
	WS4	0.2	0.3	0.5
	WS5	0.1	0.3	0.6

TABLE II :

 II Custom multi-modal integration weights sets (WS0-6 from now) tested in the object detection experiments.

TABLE III :

 III Comparisons of means of object detection with DS1 as test set in experiment 4. (Ph.: phase)

TABLE IV :

 IV Comparisons of means of object detection with DS4 as test set in experiment 5. (Ph.: phase)

	C. Experiments 6 and 7

  .

		Exp. 6	Exp. 7
	Model vs WS0	*	***
	Model vs WS1	*	***
	Model vs WS2	n.s	***
	Model vs WS3	*	***
	Model vs WS4	n.s	n.s
	Model vs WS5	n.s	n.s

TABLE V :

 V Comparisons exp. 6 and 7

When adopted as training dataset, 95% of the total

samples in the corresponding collection were used, whereas only 5% of the total 25000 samples was used if adopted as testing dataset. Datasets were shuffled

Source code is available here: https://github.com/guidoschillaci/icub perceptual optimisation.git. Tensorflow 2 has been used. Datasets are available in the Zenodo platform[START_REF] Schillaci | Dataset from simulated iCub humanoid robot for sensory enhancement and perceptual optimisation[END_REF].Experiments have been tested on different platforms, including an NVidia Jetson Nano and cloud computing resources -gently offered by the OCRE (EU-H2020 Open Clouds for Research Environment) project. Scripts for setting up Docker environments for the different platforms are also available at the provided links.

The post-hoc paired comparisons of the MIW performance on the mean 'objects detected' between the six different custom multi-modal integration weights sets were of special interest for the present work. Therefore, only the post-hoc paired comparisons between the main model and the six different custom weights sets are reported.

ACKNOWLEDGEMENTS

GS has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 838861 (Predictive Robots). Predictive Robots is an associated project of the Priority Programme "The Active Self" funded by the DFG (German Research Foundation).