Bahram Houchmandzadeh

Marcel Vallade

Exact results for fixation probability of bithermal evolutionary graphs

Keywords: Evolutionary graphs, fixation probability, fitness, probability generating functions

One of the most fundamental concepts of evolutionary dynamics is the "fixation" probability, i.e. the probability that a mutant spreads through the whole population. Most natural communities are geographically structured into habitats exchanging individuals among each other and can be modeled by an evolutionary graph (EG), where directed links weight the probability for the offspring of one individual to replace another individual in the community. Very few exact analytical results are known for EGs. We show here how by using the techniques of the fixed point of Probability Generating Function, we can uncover a large class of of graphs, which we term bithermal, for which the exact fixation probability can be simply computed.

Introduction.

Evolutionary dynamics is a stochastic process due to competition between deterministic selection pressure and stochastic events due to random sampling from one generation to the other. One of the most fundamental concepts of evolutionary dynamics is the fixation probability, i.e. the probability that a mutant spreads and takes over the whole community(Patwa and Wahl [START_REF] Patwa | The fixation probability of beneficial mutations[END_REF]). In the framework of the Moran model (Moran [START_REF] Moran | The Statistical processes of of evolutionary theory[END_REF]) for a well mixed population, where an individual's offspring can replace any other one in the community, the fixation probability is

π f = 1 -r -m0 1 -r -M (1)
where M is the size of the community, m 0 the original number of mutants and r the relative fitness of the mutants. A similar result was reached by Kimura (Kimura [START_REF] Kimura | On the probability of fixation of mutant genes in a population[END_REF]) for the Fisher-Wright model under the diffusion approximation. The idea of well mixed population is however far from realistic. Natural communities are geographically extended and subdivided into patches that exchange individuals (figure 1a). Maruyama (Maruyama [9,[START_REF] Maruyama | A simple proof that certain quantities are independent of the geographical structure of population[END_REF]) was the first to cast the The star configuration, which can be proved to differ from a non-structured population for its fixation probability.

problem of evolutionary dynamics into a Moran process on graphs (or islands, in his terminology) and under harsh approximations, concluded that the fixation probability does not depend on the population structure. The first formal proof that Maruyama's results are not always correct was given by Lieberman et al (Lieberman et al. [8]) who showed that for a Moran process on a star graph (Figure 1b), the effective fitness of the mutant, in the large population limit, can be enhanced to r 2 due to topological effects. In order to do so, Lieberman et al. considered evolutionary graphs (EG) where nodes contain exactly one individual, whether wild type or mutant, connected by directed links representing the geographical (or social) connectivity. Lieberman et al. also extended their results to the funnel topology with K layers where the effective fitness, in the limit of large population, can be amplified to r K , but provided only a sketch of the proof. Beyond the special cases considered by Lieberman et al, very few exact analytical results are known. A review of the present state of known results is given by Shakarian et al (Shakarian et al. [13]).

Most of the results of the EG are obtained through Monte Carlo numerical simulations. These simulations however scale as 2 M where M is the number of nodes. A new numerical scheme has been proposed (Barbosa et al. [START_REF] Barbosa | Early appraisal of the fixation probability in directed networks[END_REF]) to accelerate the speed of these simulations, but the computation of the fixation probability of large graphs is still very time consuming.

It would therefore be important to know the exact fixation probability of a large class of graphs that can be used as an approximation of closely related graphs or as a benchmark for assessing the progress in numerical simulation schemes. This is the aim of the present article.

We recently proposed a new method (Houchmandzadeh and Vallade [START_REF] Houchmandzadeh | The fixation probability of a beneficial mutation in a geographically structured population[END_REF]), based on the fixed points of the time dependent Probability Generating Function (fp-PGF) which can efficiently approximate the fixation probability of large, arbitrary graphs by solving only a system of M second order algebraic equations. We show in the present article that the fp-PGF method can also be used to derive exact results for a large class of graphs that we call bithermal, which are an extension of the isothermal graphs considered by Lieberman et al.

The temperature of a node is related to the imbalance between the sum of the weights of the incoming and outgoing links. In isothermal graphs, all nodes are balanced and have the same temperature T = 1. Bithermal graphs are bipartite, with two kinds of nodes at either temperature T A or T B . The star topology (figure 1b) is one particular example of such graphs, some other particular examples are shown in Figure 2 . We show here that the fixation probability of these graphs is a simple rational function of the fitness r and of the number of nodes M A and M B in each class

π f = f (M A , M B , r)
The exact expression of this function is given by equation ( 21) and its plot by figure 4a. When there is the same number of nodes in each class (M A = M B ), bithermal graphs become isothermal and the function f above is equal to the Moran expression [START_REF] Antal | Evolutionary dynamics on degreeheterogeneous graphs[END_REF]. When the imbalance between the number of nodes in each class is large (M A ≫ M B or M A ≪ M B ), the fixation probability tends toward 0 or 1-1/r 2 , depending on the nature of the Moran process (birth-death or death-birth).

This article is organized as follow. In the next section, we recall the continuous time stochastic process of Moran on graph and its associated Master equation ; the third section is devoted to the Probability Generating Function method ; In the fourth section, we apply these results to the bithermal graphs and give their exact fixation probability. The last section is devoted to some generalizations and conclusions.

Continuous time Moran process on graph.

Consider a community of M individuals, which can either be wild type (WT) with fitness 1 or mutant with relative fitness r = 1+s. The individuals are spread spatially and the progeny of an individual i can replace individual j according to a connectivity map. The connectivity map can be envisioned as a graph, where each node i contains exactly one individual, either mutant or wild type ; the weight of a link m ij specifies the probability for the progeny of an individual at node i to replace an individual at node j (Figure 1a). The coefficients m ij are collected into a connectivity matrix m. As the number of individual is fixed, it is sufficient to specify the number of mutants (0 or 1) on each node n = (n 1 , n 2 , ...n M ) at a given time to have complete information about the system at this time. We consider here a continuous time model where birth (or death) events occur randomly with rate µ (Houchmandzadeh and Vallade [START_REF] Houchmandzadeh | Alternative to the diffusion equation in population genetics[END_REF]). The probability density for a node i to decrease or increase its number of mutants by one unit during a time interval dt is (Houchmandzadeh and Vallade [6])

W - i (n i ) = (µ/M )n i M k=1 m ki (1 -n k ) (2)
W + i (n i ) = (µ/M )r (1 -n i ) M k=1 m ki n k (3)
It is crucial at this step to distinguish between two kinds of Moran processes (Antal et al. [START_REF] Antal | Evolutionary dynamics on degreeheterogeneous graphs[END_REF]). In the first case which we call D-B (Death first and then replacement, also called Voter Model), a death occurs first on a node, then immediately one connected node duplicates and send its progeny to this node. Equation ( 2) is therefore the probability density that a mutant dies at node i during dt, and is replaced by the progeny of a WT on a connected node k. Equation ( 3) is the probability density that a WT dies at node i and is replaced by a mutant on a connected node. Without loss of generality, the mutant's advantage r is included in this line, either as a decreased mortality or a better replacement success once a death event has occurred.

In the case of B-D processes (also called Invasion Process), a birth occurs first on a node k and the progeny is then sent to a connected node i to replace the local resident. The transition probabilities are still expressed by the same equations (2,3), but the quantity µ now denotes the birth rate.

Although the rate equations (2,3) are similar for these two processes, the normalization conditions of m ki coefficients are different :

M i=1 m ik = 1 (D-B); M i=1 m ki = 1 (B-D) ∀k (4)
The temperature of a node is defined for these processes as

M i=1 m ki = T k (D-B) ; M i=1 m ik = T k (B-D) ∀k (5)
Because of the normalization constraint (4), we must have M j=1 T j = M for both processes. This means that if some nodes are cold (T < 1), others must be hot (T > 1).

The Moran process is a one-step stochastic one, where during an infinitesimal interval dt, only one birth or death event can occur. Equations [START_REF] Barbosa | Early appraisal of the fixation probability in directed networks[END_REF][START_REF] Ewens | Mathematical Population Genetics[END_REF] are transition probabilities between states n on the one hand and states a i n = (n 1 , ...n i -1, ...n M ) and a † i n = (n 1 , ...n i +1, ...n M ) on the other. The probability P (n, t) of observing state n at time t obeys the Master equation

∂P (n, t) ∂t = {m} W (m → n)P (m, t) -W (n → m)P (n, t) (6)
The EG process we are considering has two absorbing states 1 = (1, 1, ...1) and 0 = (0, 0, ...0). Once a mutant has invaded all the nodes or has been eliminated from all of them, it is fixed or lost and there is no further evolution (until a new mutant appears by random mutation) : W (1 → m) = 0 and W (0 → m) = 0 as can be deduced from equations [START_REF] Barbosa | Early appraisal of the fixation probability in directed networks[END_REF][START_REF] Ewens | Mathematical Population Genetics[END_REF]. Note that the probability of reaching state 1 from an initial state n, i.e. the fixation probability π(n) can in principle be found using the Kolmogorov's backward equation (Ewens [3]) :

m (π(n) -π(m)) W (n → m) = 0 (7) π(0) = 0 ; π(1) = 1 (8)
which is a set of linear equations in the unknowns π(n). The direct resolution of the above set of equation however can be attempted only in special cases. The best example of a direct solution is the unstructured population, where the graph is fully connected and m ij = 1/M ; the EG dynamics can then be mapped into a biased one-dimensional Brownian motion and solved by standard techniques (Ewens [START_REF] Ewens | Mathematical Population Genetics[END_REF]), which yields the well known result [START_REF] Antal | Evolutionary dynamics on degreeheterogeneous graphs[END_REF] . Other cases, such as the star topology in the limit of large population considered by Lieberman et al (Lieberman et al. [START_REF] Lieberman | Evolutionary dynamics on graphs[END_REF]), use such careful mapping. The mapping method however is hard to generalize.

The fp-PGF method.

Computation of the fixation probability can be simplified if instead of the linear system [START_REF] Kimura | On the probability of fixation of mutant genes in a population[END_REF], we use the dynamics of the Probability Generating Function φ(z, t) = {n} P (n, t)z n where the variable z = (z 1 , z 2 , ...z M ) is conjugate to n = (n 1 , n 2 , ..., n M ) and z n = z n1 1 z n2 2 ...z nM M . Time is measured in generation time units M/µ. Note that φ(0, t) = 0 and φ(1, t) = 1. From the Master Equation ( 6), we can derive the dynamics of the PGF (Houchmandzadeh and Vallade [START_REF] Houchmandzadeh | The fixation probability of a beneficial mutation in a geographically structured population[END_REF]) which reads:

∂φ ∂t = M k=1 f k (z) ∂φ ∂z k - M i,k=1 g i,k (z) ∂ 2 φ ∂z i ∂z k (9)
where for D-B processes,

f k (z) = z k M i=1 m ki (z i -1) -r -1 (z k -1) (10)
g i,k (z) = m ki (z i -1)(z i -r -1 )z k (11)
For the B-D process, the first order term is slightly different and reads

f k (z) = z k M i=1 m ki (z i -1) -(T k /r)(z k -1)
Solving the PGF equation ( 9) would seem at least as formidable as solving directly the Master Equation ( 6). However, for the computation of the fixation probabilities, we are only interested in the large time limit t → ∞. At large time, the mutant is either fixed or lost, therefore the stationary solution φ s (z) to which the PGF converges is simply

φ s (z) = π 0 + π f M i=1 z i (12)
where π 0 and π f are the loss and fixation probabilities and implicit functions of the initial conditions. It can also be checked manually that ( 12) is indeed a solution of [START_REF] Maruyama | A markov process of gene frequency change in a geographically structured population[END_REF]. The problem of finding π 0 and π f becomes trivial if the PGF possesses a fixed point

ζ = (ζ 1 , ...ζ M ) such that ∂φ ∂t z=ζ = 0
In this case, we have

φ(ζ, 0) = φ(ζ, ∞) = π 0 + π f M i=1 ζ i
and, as π f + π 0 = 1,

π f = 1 -φ(ζ, 0) 1 - M i=1 ζ i
As the quantity φ(ζ, 0) is known from the initial conditions, finding a fixed point of the PGF determines entirely the fixation probability. Note that once a fixed point has been found, the fixation probability for any initial condition can be trivially computed.

For the initial condition of the mutant appearing at random with probability 1/M on one node, φ(z 1 , ...z M ; t = 0) = (1/M ) i z i and therefore the fixation probability is given by

π f = 1 -(1/M ) i ζ i 1 -i ζ i (13)
The condition for a point ζ to be a fixed point is

f k (ζ) = 0 ∀k (14) g i,k (ζ) + g k,i (ζ) = 0 ∀i, k (15)
Whether such a fixed point exists or not depends on the connectivity matrix m ik . For an isothermal graph where all the nodes k have temperature T k = 1, it is easy to check that ζ = r -1 1 is a fixed point : and the fixation probability (13) of the isothermal graph is equal to the result (1) for unstructured populations. We see here how easily this theorem can be obtained from the fixed point of the PGF.

f k (r -1 1) = r -1 (r -1 -1) -r -1 (r -1 -1) = 0 g i,k (r -1 1) = m ik (r -1 -1)(r -1 -r -1 )r -1 = 0

bithermal graphs.

We now consider a subset of bipartite graphs that we call bithermal [START_REF] Barbosa | Early appraisal of the fixation probability in directed networks[END_REF] for which the fixation probability can be determined in algebraic closed form. In these graphs, there are M A nodes of type A at temperature T A and M B nodes of type B at temperature T B . More over, we require that two nodes can be connected only if they are at different temperatures, and if a node i is connected to a node j, then j is also connected to i. Of course, we suppose that the graph is connected, i.e. there is always a path from any node i to any node j. With appropriate numbering of the nodes, the connectivity of such graph is a block matrix of the form m = 0 α β 0 (figure 3). The star topology (figure 1b) is such a graph with M A = 1 central and M B peripheral nodes. For a D-B process, all the elements of the star's β block are equal to 1/M B and all the elements of α block are equal to 1. Therefore

T A = M B and T B = 1/M B .
The star graph can be generalized to the case where the number of central nodes M A > 1 (figure 2a). For a general bithermal graph the number of B nodes to which an A node connects and the weights of these links can be arbitrary, as long as the constraints on temperatures are respected. Note that T A and T B are not independent. By summing over all the elements m ki of the β or α block of the connectivity matrix m, we have for a D-B process

M A = M B T B ; M B = M A T A
for a B-D process, the role of M A and M B are exchanged, but in both cases, we have

T A T B = 1 (16)
We now search for the bithermal graphs which have an exact fixed point. Following the example of the isothermal graphs, we look for a fixed point

ζ = (ζ 1 , ...ζ M ) where ζ i = ζ A if i ∈ A and ζ i = ζ B if i ∈ B .
Let us first consider the case of D-B processes. Using condition (14) for k ∈ A or B we obtain a set of two algebraic equations :

f A (ζ) = ζ A (ζ B -1)T A -r -1 (ζ A -1) = 0 f B (ζ) = ζ B (ζ A -1)T B -r -1 (ζ B -1) = 0
the solution of which is given by

ζ A = u(r, M A /M B ) (17) ζ B = u(r, M B /M A ) (18) u(r, T ) = 1/r + T r + T (19)
This solution also satisfies condition (15) if for k ∈ A and i ∈ B, we have the following relation between the coefficient of the connectivity matrix m :

m ki m ik = - (ζ B -1)(ζ B -r -1 )ζ A (ζ A -1)(ζ A -r -1 )ζ B = 1 T A = M A M B
which we can express as a relation between the two blocks β and α of the connectivity matrix m:

α = T A β ⊺ (20)
We observe here that this is the sufficient condition to form exactly solvable bithermal graphs : form an M B × M A matrix β where the sum of elements in each column is 1 and the sum of elements in each row is 1/T A = M A /M B ; form the block matrix m from β and α = T A β ⊺ . The β block contains M A M B coefficients subject to M A + M B summation rules, so the number of bithermal graphs with exact solutions is indeed very large when M A or M B are large. An important subset of bithermal graphs that always has an exact fixed point is a set we call symmetric bithermal graphs. For these graphs, all the existing links from a node A to a node B (respect. B to A ) have the same weight m AB (resp. m BA ). The generalized star graph (fig 2a) belongs to this subset. For these sets, the weight of a link is determined only by the number of connected nodes, and each A (resp. B) node has always the same number of neighbors. The symmetric subset automatically satisfies all the constraints and always has exact fixed points. individual. This restriction can be relaxed and we can let each node contain a number of individuals N ≥ 1 (figure 5) . This is equivalent to the original island model of Maruyama (Maruyama [9]) or can alternatively be envisioned as a bi-level graph as defined by Shakarian et al (Shakarian et al. [START_REF] Shakarian | A review of evolutionary graph theory with applications to game theory[END_REF]), where each node is mapped into N sub-nodes conserving the original topology. As we have shown earlier (Houchmandzadeh and Vallade [6]), this parameter does not alter the expression of the PGF and the fixed points are still computed by the same equations. The fixation probability is slightly modified and reads

π f = 1 -(M A ζ A + M B ζ B )/(M A + M B ) 1 -ζ N MA A ζ N MB B
The approach can be extended even further and allows for different numbers N k of individuals for each node k (Houchmandzadeh and Vallade [START_REF] Houchmandzadeh | The fixation probability of a beneficial mutation in a geographically structured population[END_REF]).

Another direction toward which this work can be extended is the n-thermal graphs. Here we have considered only bithermal graphs composed of two types of nodes which we can represent by an A -B topology. In principle, we can generalize the method to graphs containing P types of nodes O 1 , ..., O P : The M I nodes belonging to type O I have the temperature T I and a hot class can only be connected to cold classes and vice versa. We could in principle form a polymeric topology such as O 1 -O 2 ... -O P , branched systems, closed rings and so on. The exploration of these topologies implies an analytic study of the roots of algebraic equations of order 2P and is beyond the scope of the present article.

To summarize, we have obtained exact analytical results for a wide class of graphs that we have called bithermal in the field of Evolutionary Graph Theory. EGT is a cornerstone for our understanding of evolution, because natural population are always geographically extended and cannot be a priori approximated as "well mixed". Exact results in EGT however have been hard to obtain because in each case, a mapping into a one-dimensional Brownian motion has to be constructed ; whether such a mapping exists or not for a particular problem is not a trivial problem. The method we develop is radically different : by using the continuous time version of the Moran model and the dynamics of the PGF, we reduce the problem of finding an exactly solvable model into finding the roots of algebraic equations. We have illustrated the power of this dynamical method through our study of bithermal graphs. We believe that the method we have presented can be a powerful tool to get exact results for the fixation probability of more complex evolutionary graphs.

Figure 1 :

 1 Figure 1: Evolutionary Graphs. (a) individuals are spread in space. A node i can send its offspring only to connected nodes j with a probability m ij . A non-structure population can be considered as a fully connected graph with uniform migration probability m ij = 1/N . (b)The star configuration, which can be proved to differ from a non-structured population for its fixation probability.

Figure 2 :

 2 Figure 2: Some examples of bithermal graphs; for simplicity, each link represents two directed links. A bithermal graph is constituted of two kinds of nodes at respective temperature T A and T B , indicated here by their colors. Each A node is connected to a subset of B nodes and vice versa. (a) a (4,12,6,2) generalized star, consisting of 4 central nodes and 12 peripheral nodes, each central node connected to 6 peripheral and each peripheral connected to two central nodes ; (b) an example where the combined effect of the number of links and their weight makes the graph bithermal ; (c) A symmetric 2 dimensional, bithermal crystal with periodic boundaries;

Figure 3 :

 3 Figure 3: The connectivity matrix of a 2-thermal graph for a D-B process, where nodes of type A are numbered from 1 to M A and nodes of type B from M A + 1 to M A + M B . For a B-D process, the sum of the elements of a row is equal to unity, and the sum of the elements of a column is equal to T A or T B .

Figure 4 :

 4 Figure 4: Fixation probability π f of bithermal EGs of size M = 100 and r = 1.25. (a) π f as a function of the number M A of nodes of type A for a generalized (M A , M B ) star. circles (black) and triangles (red) correspond to numerical resolution of the fixation probabilities for B-D and D-B processes. Solid lines correspond to exact solutions given by equation (21). The dashed lines correspond to π f = 1 -1/r and π f = 1 -1/r 2 . Numerical simulations were performed by using a Gillespie Algorithm(Gillespie [4]) ; for each point, 10 5 stochastic paths were generated and used to compute the fixation probability. (b)Fixation probability π f for random bithermal EG for different values of M A (circles 5, squares 10, × 15, diamond 20, triangle up 30,triangle down 45) and D-B processes. To the left of the vertical line (#index ≤50) : Each point corresponds to a random bithermal connectivity matrix where the two blocks are related through α = T A β ⊺ . To the right of the vertical line (#index >50) : For each point, the two blocks of the bithermal connectivity matrix are random and α = T A β ⊺ . Solid lines indicate the theoretical values computed by equation (21).

Figure 5 :

 5 Figure 5: Generalization to bi-level graphs, where each node can be considered as composed of N individuals, conserving the same topology, i.e. the weight of new links is downscaled by a factor N : m ′ ij = m ij /N

Acknowledgements.

We are grateful to E. Geissler and O. Rivoire for the careful reading of the manuscript and fruitful discussions. This work was partly funded by Agence Nationale de la Recherche Française (ANR) grant "Evo-Div."

Once the fixed point is known, the fixation probability is easily determined from equation [START_REF] Shakarian | A review of evolutionary graph theory with applications to game theory[END_REF]:

For a B-D process, the computation of the fixed point follows exactly the same line of argument. The result is obtained by permuting M A and M B :

The fixation probability is given by the same expression (21). Note that D-B and B-D processes act in opposite directions. Compared to a non-structured population of the same size, a bithermal D-B process acts as a suppressor of selection, lowering the fixation probability where the B-D process is an amplifier of selection, increasing the fixation probability. The maximum π f for a D-B is obtained for T A = 1 and is equal to expression (1) ; the maximum π f for a B-D process is obtained for T A → ∞ and is equal to 1 -r -2 . A plot of both fixation probabilities and their comparison to numerical simulations is shown in figure 4 .

We stress that for bithermal graphs, the details of the connectivity are not important : the fixation probability depends only on the total number of A and B nodes. Consider for example the symmetric generalized star (M A , M B , p, q) where each A node is connected to p nodes of type B and each B nodes to q nodes of type A. For a fully connected generalized star, p = M B and q = M A ; an example of partially connected generalized star is given in figure 2a. We emphasize that for generalized stars, the fixation probability does not depend on the detail of the connectivity p, q, a result which would be hard to predict by other methods. We also note from numerical simulations that the fixation time is also only a function of M A and M B and does not depend on the detail of the connectivity.

Finally, we note that even when the constraint (20) is not respected and α = T A β ⊺ , the fixation probability of bithermal graphs is well approximated by expression (21). In this case, the point ζ computed from equations (17,18) is only a quasi-fixed point, but as we have shown earlier (Houchmandzadeh and Vallade [6]), for large communities, quasi-fixed points can be used for a good approximation of the fixation probability. It can be observed in figure (4b) that the numerical errors of fixation probabilities of connectivity matrices having exact fixed points (left of vertical line) or only quasi-fixed point (right of vertical line) are of the same magnitude, for a system as small as M = 100. We also note from numerical simulations that fixation time, is only

Discussion and Conclusion.

The approach we presented above can be generalized in various directions. We have restricted our approach to the case where each node contains only one


A Martinelli
email: agostino.martinelli@ieee.org

C Troiani
email: chiara.troiani@inria.fr

and Robotics, contract #231855 (sFLY)

Vision-Aided Inertial Navigation Using Virtual Features

Chiara Troiani and Agostino Martinelli

Abstract-In this paper we consider an aerial vehicle equipped with a monocular camera and inertial sensors. Additionally, a laser pointer is mounted on the vehicle and it produces a laser spot. The laser spot is observed by the monocular camera and it is the unique point feature used in the proposed approach. We focus our attention to the case when the vehicle moves in proximity of a planar surface and in particular when the laser spot belongs to this surface. The paper provides two main contributions. The former is the analytical derivation of all the observable modes, i.e. all the physical quantities that can be determined by only using the inertial data and the camera observations of the laser spot during a short time-interval. Specifically, it is shown that the observable modes are: the distance of the vehicle from the planar surface; the component of the vehicle speed, which is orthogonal to the planar surface; the relative orientation of the vehicle with respect to the planar surface; the orientation of the planar surface with respect to the gravity. The second contribution is the introduction of a simple recursive method to perform the estimation of all the aforementioned observable modes. This method is based on a local decomposition of the original system, which separates the observable modes from the rest of the system. The method is validated by using synthetic data. Additionally, preliminary tests with real data are provided and more complete experiments are in progress. The presented approach can be integrated in the framework of autonomous take-off and landing, safe touch-down and low altitude manoeuvres even in dark or featureless environment.

I. INTRODUCTION

In recent years, vision and inertial sensing have received great attention by the mobile robotics community. These sensors require no external infrastructure and this is a key advantage for robots operating in unknown environments where GPS signals are shadowed. Additionally, these sensors have very interesting complementarities and together provide rich information to build a system capable of vision-aided inertial navigation and mapping.

A special issue of the International Journal of Robotics Research has recently been devoted to the problem of fusing vision and inertial data [START_REF] Dias | 2nd Workshop on Integration of Vision and Inertial Sensors[END_REF]. In [START_REF] Corke | An Introduction to Inertial and Visual Sensing[END_REF], a tutorial introduction to the vision and inertial sensing is presented. This work provides a biological point of view and it illustrates how vision and inertial sensors have useful complementarities allowing them to cover the respective limitations and deficiencies. In [START_REF] Armesto | Vincze Fast Ego-motion Estimation with Multi-rate Fusion of Inertial and Vision[END_REF], these sensors are used to perform egomotion estimation. The sensor fusion is obtained by an Extended Kalman Filter (EKF ) and by an Unscented Kalman Filter (U KF ). The approach proposed in [START_REF] Gemeiner | Simultaneous Motion and Structure Estimation by Fusion of Inertial and Vision Data[END_REF] extends the previous one by also estimating the structure of the environment where the motion occurs. In particular, new landmarks are inserted on line into the estimated map. This approach has been validated by conducting experiments in a known environment where a ground truth was available. Also, in [START_REF] Veth | Fusing low-cost image and inertial sensors for passive navigation[END_REF] an EKF has been adopted. In this case, the proposed algorithm estimates a state containing the robot speed, position and attitude, together with the inertial sensor biases and the location of the features of interest. In the framework of airbone SLAM, an EKF has been adopted in [START_REF] Kim | Real-time implementation of airborne inertial-SLAM[END_REF] to perform 3D-SLAM by fusing inertial and vision measurements. It was observed that any inconsistent attitude update severely affects any SLAM solution. The authors proposed to separate attitude update from position and velocity update. Alternatively, they proposed to use additional velocity observations, such as air velocity observation. More recently, a vision based navigation approach in unknown and unstructured environments has been suggested [START_REF] Blosch | Vision Based MAV Navigation in Unknown and Unstructured Environments[END_REF]. Natraj et al. [START_REF] Natraj | Vision based attitude and altitude estimation for UAVs in dark environments[END_REF] proposed a vision based approach, close to structured light, for roll, pitch and altitude estimation of UAV. They use a fisheye camera and a laser circle projector, assuming that the projected circle belongs to a planar surface. The latter must be orthogonal to the gravity vector in order to allow the estimation of the aforementioned quantities. The attitude estimation of the planar surface becomes crucial in order to extend the operational environment of UAVs. Shipboard operations, search and rescue cooperation between ground and aerial robots, low altitude manoeuvres, require to attenuate the position error and to track the platform attitude.

Recent works investigate the observability properties of the vision-aided inertial navigation system [START_REF] Jones | Inertial Structure From Motion with Autocalibration[END_REF], [START_REF] Jones | Visual-inertial navigation, mapping and localization: A scalable real-time causal approach[END_REF], [START_REF] Kelly | Visual-inertial simultaneous localization, mapping and sensor-to-sensor self-calibration[END_REF], [START_REF] Martinelli | Closed-Form Solution for Attitude and Speed Determination by Fusing Monocular Vision and Inertial Sensor Measurements[END_REF], [START_REF] Martinelli | Vision and IMU Data Fusion: Closed-Form Solutions for Attitude, Speed, Absolute Scale and Bias Determination[END_REF] and [START_REF] Mourikis | A Multi-State Constrained Kalman filter for Vision-aided Inertial Navigation[END_REF]. In particular, in [START_REF] Martinelli | Vision and IMU Data Fusion: Closed-Form Solutions for Attitude, Speed, Absolute Scale and Bias Determination[END_REF], the observable modes are expressed in closed-form in terms of the sensor measurements acquired during a short time-interval.

In this paper we consider a vehicle which accomplishes a 3D-trajectory in the surrounding of a planar surface. The vehicle is equipped with a monocular camera and inertial sensors. This is the typical navigation problem in an indoor environment or in a city-like environment. All the approaches previously mentioned, require to extract natural features from the images provided by the camera and in particular to detect the same features in different images. The feature matching task becomes critical in outdoor environment because of possible illumination changes. In order to significantly reduce the computational burden required to perform these tasks and to make the feature matching more robust, we introduce a virtual feature by equipping our vehicle with a laser pointer. The laser beam produces a laser spot on the planar surface. This laser spot is observed by the monocular camera and it is the unique point feature used by the proposed approach.

To the best of our knowledge, this problem has never been considered so far. Compared to classical vision and IMU data fusion problems, the feature is moving in the environment but we exploit the hypothesis that it moves on a planar surface. The first question which arises is to understand which are the observable modes, i.e. the physical quantities that can be determined by only using the inertial data and the camera observation of the laser spot during a short time-interval. The results provided in section III address precisely this issue (more details about the analytical computation can be found in [START_REF] Martinelli | Vision-Aided Inertial Navigation Using Virtual Features[END_REF]). Then, the second step we consider is to analytically determine the link between the observable modes and the sensor data. This is obtained by performing a local decomposition of the original system (section IV). This decomposition separates the observable modes from the rest of the system and will allow us to introduce a simple recursive method to perform the estimation of all the observable modes (second part of section IV). The method is validated by using synthetic data (section V). Preliminary tests with real data are also provided and more complete experiments are in progress.

II. THE CONSIDERED SYSTEM

Let us consider an aerial vehicle equipped with a monocular camera and IM U sensors. The vehicle is also equipped with a laser pointer. The configuration of the laser pointer in the camera reference frame is known. The vehicle moves in the surrounding of a planar surface and we assume that the laser spot produced by the laser beam belongs to this planar surface (see fig. 1). The position and the orientation of this planar surface are unknown. The camera observations consist in the position of the laser spot in the camera frame up to a scale factor. The IMU consists of three orthogonal accelerometers and three orthogonal gyroscopes. We assume that the monocular camera is intrinsically calibrated and that the transformations among the camera frame and the IM U frames are known (we can assume that the vehicle frame coincides with the camera frame). The IM U provides the vehicle angular speed and acceleration. Actually, regarding the acceleration, the one perceived by the accelerometer (A) is not simply the vehicle acceleration (A v ). It also contains the gravity acceleration (G). In particular, we have A = A v -G since, when the camera does not accelerate (i.e. A v is zero) the accelerometer perceives an acceleration which is the same of an object accelerated upward in the absence of gravity.

We will use uppercase letters when the vectors are expressed in the local frame and lowercase letters when they are expressed in the global frame. Hence, regarding the gravity we have: g = [0, 0, -g] T , being g 9.8 ms -2 .

Finally, we will adopt a quaternion to represent the vehicle orientation. Indeed, even if this representation is redundant, it is very powerful since the dynamics can be expressed in a very easy and compact notation [START_REF]Quaternions and rotation Sequences: a Primer with Applications to Orbits, Aerospace, and Virtual Reality[END_REF].

Our system is characterized by the state [r, v, q] T where r = [r x , r y , r z ] T is the 3D vehicle position, v is its time derivative, i.e. the vehicle speed in the global frame (v ≡ dr dt ), q = q t +iq x +jq y +kq z is a unitary quaternion (i.e. satisfying q 2 t +q 2 x +q 2 y +q 2 z = 1) and characterizes the vehicle orientation. The analytical expression of the dynamics and the camera observations can be easily provided by expressing all the 3D vectors as imaginary quaternions. In practice, given a 3D vector w = [w x , w y , w z ] T we associate with it the imaginary quaternion w q ≡ 0 + iw x + jw y + kw z . The dynamics of the state [r q , v q , q] T are:

       ṙq = v q vq = qA vq q * = qA q q * + g q q = 1 2 qΩ q (1)
being q * the conjugate of q, q * = q t -iq x -jq y -kq z and Ω the angular velocity. We derive the expression of the camera observation consisting in the position of the laser spot in the camera frame up to a scale factor. The laser spot is on a planar surface whose configuration is unknown. Without loss of generality, we choose the camera frame with the z-axis parallel to the laser pointer (see figure 2). In addition, the camera frame is such that the laser beam intercept the xy-plane in [L, 0, 0]. In [START_REF] Martinelli | Vision-Aided Inertial Navigation Using Virtual Features[END_REF] we introduce a simple and efficient method in order to determine the parameter L together with the rotation to transform vectors from the original camera frame into the chosen camera frame.

Finally, we characterize the planar surface in the global frame with the equation z = ky, where k is an unknown parameter.

In these settings, by carrying out analytical computation (which uses the basic quaternion rules) we obtain the analytical expression of the position [X s , Y s , Z s ] of the laser spot in the camera reference frame. We have: The original camera frame XY Z, the chosen camera frame X Y Z and the laser module at the position [Lx, Ly, 0] and the direction (θ, φ) (in the original frame) and position [L, 0, 0] and the direction (0, 0) (in the chosen camera frame).

         X s = L Y s = 0 Z s = r z + 2q z q x L -2q y kq x L -2q y q t L -2q t Lkq z -kr y 2kq z q y -2kq t q x -q 2 z -q 2 t + q 2 y + q 2 x
(2) The camera provides the vector [X s , Y s , Z s ] up to a scale factor. This is equivalent to the two ratios Xs Zs and Ys Zs . Hence, since the latter is identically zero, the camera observation is given by h cam = Xs Zs , which is:

h cam = L(2kq z q y -2kq t q x -q 2 z -q 2 t + q 2 y + q 2 x ) r z + 2q z q x L -2q y kq x L -2q y q t L -2q t Lkq z -kr y (3)

III. OBSERVABILITY PROPERTIES

We investigate the observability properties of the system whose dynamics are given in (1) and whose observation function is given in [START_REF] Blosch | Vision Based MAV Navigation in Unknown and Unstructured Environments[END_REF]. We have also to consider the constraint q * q = 1. This can be dealt as a further observation (system output):

h const (r q , v q , q) = q * q (4)
Finally, we want to investigate whether the parameter k is identifiable or not. This is done by performing an observability analysis on the extended state S = [r q , v q , q, k] T , whose dynamics are given in [START_REF] Armesto | Vincze Fast Ego-motion Estimation with Multi-rate Fusion of Inertial and Vision[END_REF] and by the additional equation k = 0.

We apply the method introduced in [START_REF] Martinelli | State Estimation Based on the Concept of Continuous Symmetry and Observability Analysis: the Case of Calibration[END_REF]. This will allow us to detect all the observable modes, i.e. all the physical quantities that can be determined by only using the information contained in the data provided by the IMU and the camera in a given time-interval. In [START_REF] Martinelli | Vision-Aided Inertial Navigation Using Virtual Features[END_REF] we perform this analytical computation and we show that the system has six observable modes. Additionally, in [START_REF] Martinelli | Vision-Aided Inertial Navigation Using Virtual Features[END_REF], we provide the physical meaning of these observable modes. Specifically, we found convenient to express the observable modes in a new global reference frame, denoted with x, ỹ, z. In this frame the z-axis coincides with the axis normal to the planar surface. From now on, we will adopt this global frame to characterize the vehicle configuration. The x, ỹ, z-frame is obtained by rotating the x, y, z-frame about the x-axis of the angle arctan(k). Hence, it is characterized by the quaternion:

p = cos arctan(k) 2 + i sin arctan(k) 2 (5)
The state in this frame is S = [r q ṽq q, k] T , where:

q = p * q rq = p * r q p ṽq = p * v q p ( 6
)
The observable modes are (see the derivation in [START_REF] Martinelli | Vision-Aided Inertial Navigation Using Virtual Features[END_REF]):

m 1 = rz m 2 = ṽz m 3 = 2(q t qx + qy qz ) m 4 = 2(q t qy -qx qz ) m 5 = k m 6 = q * q (7)
m 1 is the z-coordinate of the vehicle, m 2 the component of the vehicle speed along the z-axis, m 3 and m 4 are related to the roll and pitch angles of the vehicle in this frame. In particular, the roll is

arctan m3 √ 1-m 2 3 -m 2 4
and the pitch is arcsin(m 4 ). m 5 is related to the orientation of the xỹ-plane with respect to the gravity. In particular, the z-axis makes an angle arctan(k) = arctan(m 5 ) with the gravity. m 6 is trivially the magnitude of the quaternion q, which is 1 since it describes a rotation. By summarizing the results of the observability analysis performed in this section (and in [START_REF] Martinelli | Vision-Aided Inertial Navigation Using Virtual Features[END_REF]) we say that the information contained in the data provided by the IM U and the camera during a given time-interval, allows us to determine the six modes m 1 , ..., m 6 . For this reason, in the rest of the paper, we will focus our attention only on these six quantities (actually, on the first five, since m 6 trivially expresses the constraint of having a unitary quaternion).

IV. LOCAL DECOMPOSITION AND RECURSIVE ESTIMATION

The goal of this section is to provide a method able to estimate the observable modes in [START_REF] Gemeiner | Simultaneous Motion and Structure Estimation by Fusion of Inertial and Vision Data[END_REF]. To achieve this goal, the first step is to determine the link between the observable modes and the sensor data. By adopting the terminology introduced in [START_REF] Isidori | Nonlinear Control Systems[END_REF], we have to perform a local decomposition of our system. We remind the reader that the local decomposition is the extension of the Kalman canonical decomposition [START_REF] Brogan | Modern Control Theory[END_REF] to the case of a non linear system. It consists in writing the equations characterizing the dynamics and the observation only in terms of the observable modes. We also remind the reader that in the non linear case it is often impossible to characterize the system with a unique decomposition. The decomposition only holds in a local region of the space of states. This is the reason why it is called local decomposition. To cover the entire space of states more than one decomposition is required (see [START_REF] Isidori | Nonlinear Control Systems[END_REF]). In the following, we will show that for our system the number of decompositions is two.

We first provide the dynamics of the state S = [r q ṽq q, k] T . We obtain [START_REF] Martinelli | Vision-Aided Inertial Navigation Using Virtual Features[END_REF]:

             ṙq = ṽq vq = qA q q * + p * g q p q = 1 2 qΩ q k = 0 (8)
A local decomposition for the dynamics is:

               ṁ1 = m 2 ṁ2 = -m 4 A x + m 3 A y + ξ(m 3 , m 4 )A z + g z ṁ3 = Ω x ξ(m 3 , m 4 ) + Ω z m 4 ṁ4 = Ω y ξ(m 3 , m 4 ) -Ω z m 3 ġz = 0 (9)
where g z is the component of the gravity along the z-axis, i.e.

g z = -g cos(arctan(k)) = -g √ 1+k 2 = -g √ 1+m 2 5
; the function ξ(m 3 , m 4 ) depends on the original state and in particular changes its sign depending on the sign of q2

x + q2 y -1 2 :

ξ(m 3 , m 4 ) ≡      1 -m 2 3 -m 2 4 if q2 x + q2 y < 1 2 -1 -m 2 3 -m 2 4 if q2 x + q2 y > 1 2 (10
) Hence, as previously said, we have two local decompositions for our original system. The validity of (9) can be checked by using [START_REF] Gemeiner | Simultaneous Motion and Structure Estimation by Fusion of Inertial and Vision Data[END_REF] and [START_REF] Isidori | Nonlinear Control Systems[END_REF]. Note that deriving (9) is troublesome. In contrast, checking its validity is very simple since it only demands to perform differentiation.

To complete the local decomposition we need to express the camera observation function in (3) in terms of the observable modes. We obtain:

h cam = Lξ(m 3 , m 4 ) m 4 L -m 1 (11)
The validity of (11) can be checked by using (3), ( 5), ( 6), [START_REF] Gemeiner | Simultaneous Motion and Structure Estimation by Fusion of Inertial and Vision Data[END_REF] and [START_REF] Jones | Visual-inertial navigation, mapping and localization: A scalable real-time causal approach[END_REF]. The equations ( 9) and ( 11) represent a local decomposition for our system. They provide the analytical expression of the link between the observable modes and the sensor data. Specifically, equation [START_REF] Jones | Inertial Structure From Motion with Autocalibration[END_REF] provides the link between the observable modes and the IM U data. Equation [START_REF] Kelly | Visual-inertial simultaneous localization, mapping and sensor-to-sensor self-calibration[END_REF] provides the link between the observable modes and the data delivered by the monocular camera. Having these equations allows us to perform the estimation of the state [m 1 , m 2 , m 3 , m 4 , g z ]. An efficient and simple approach is obtained by using an Extended Kalman Filter, EKF . To implement this filter it suffices to compute the Jacobian of the dynamics in [START_REF] Jones | Inertial Structure From Motion with Autocalibration[END_REF] and the Jacobian of the observation function in [START_REF] Kelly | Visual-inertial simultaneous localization, mapping and sensor-to-sensor self-calibration[END_REF] (see [START_REF] Bar-Shalom | Tracking and data association[END_REF]). We provide these Jacobians in [START_REF] Martinelli | Vision-Aided Inertial Navigation Using Virtual Features[END_REF].

V. PERFORMANCE EVALUATION

We evaluate the performance of the proposed strategy by using both synthetic and real data. The advantage of simulations is that the ground truth is perfectly known and this allows us a quantitative evaluation of the proposed strategy.

A. Simulations

We simulate many different trajectories in 3D and many different scenarios corresponding to different orientation of the planar surface. For all the simulations we use the proposed strategy to estimate the observable modes, i.e.:

1) the distance of the vehicle from the planar surface (d = |m 1 |); 2) the component of the vehicle speed orthogonal to the planar surface (v o = m 2 ); 3) the roll (R) and the pitch (P ) angles in the x, ỹ, zframe (i.e. the frame where the x, ỹ-plane coincides with the planar surface); 4) the orientation of the plane with respect to the gravity (α). Specifically, in all the simulations the values of the estimated d, v o , R, P and α will be compared with the ground truth values.

1) Simulated Trajectories: The trajectories are generated by randomly generating the linear and angular acceleration of the vehicle at 100 Hz. In particular, at each time step, the three components of the linear acceleration and the angular speed are generated as Gaussian independent variables whose mean values will be denoted respectively with µ a and µ ω and whose variances will be denoted respectively with σ 2 a and σ 2 ω . We set the parameters in order to be close to a real case: µ a = 0 ms -2 , σ a = 1 ms -2 , µ ω = 0 deg s -1 and σ ω = 10 deg s -1 . The initial vehicle position is at x = 0, ỹ = 0, z = 1m. The initial vehicle speed is [1, 0, 0]ms -1 in the x, ỹ, z-frame.

2) Simulated Sensors: Starting from the performed trajectory, the true angular speed and the linear acceleration are computed at each time step of 0.01s (respectively, at the time step i, we denote them with Ω true i and A true v i ). Starting from them, the IMU sensors are simulated by randomly generating the angular speed and the linear acceleration at each step according to the following:

Ω i = N Ω true i
-Ω bias , P Ωi and A i = N A true v i -A gi -A bias , P Ai where: • N indicates the Normal distribution whose first entry is the mean value and the second its covariance matrix; • P Ωi and P Ai are the covariance matrices characterizing the accuracy of the IM U ; • A gi is the gravitational acceleration in the local frame and A bias is the bias affecting the data from the accelerometer;

• Ω bias is the bias affecting the data from the gyroscope. In all the simulations we set both the matrices P Ωi and P Ai diagonal and in particular: P Ωi = σ 2 gyro I 3 and P Ai = σ 2 acc I 3 , where I 3 is the identity 3 × 3 matrix. We considered several values for σ gyro and σ acc , in particular: σ gyro = 1 deg s -1 and σ acc = 0.01 ms -2 .

Regarding the camera, the provided readings are generated in the following way. By knowing the true trajectory, and the position and the orientation of the planar surface, the true bearing angles of the laser spot in the camera frame are computed 1 . They are computed each 0.1s. The parameter L is set equal to 0.3m. Then, the camera readings are generated by adding to the true values zero-mean Gaussian errors whose variance is equal to (1 deg) 2 for all the readings.

3) Simulation Results: We also evaluated the robustness of the filter with respect to systematic errors on the imu-camera calibration and lasercamera calibration. Specifically, we performed simulations by introducing errors of one cm and one deg on the calibration parameters. These systematic errors affect the estimated α (the difference with respect to the ground truth is in the range [START_REF] Brogan | Modern Control Theory[END_REF][START_REF] Dias | 2nd Workshop on Integration of Vision and Inertial Sensors[END_REF]deg) while for all the other observable modes the effect is negligible (less than one deg for R and P and less than 1cms -1 and 1cm respectively for v o and d).

B. Preliminary experiments

In this section we provide preliminary results obtained by using a data set provided by the autonomous system laboratory at ETHZ in Zurich. The data are provided together with a reliable ground-truth, which has been obtained by performing the experiments at the ETH Zurich Flying Machine Arena [START_REF] Lupashin | A simple learning strategy for high-speed quadrocopter multi-flips[END_REF], which is equipped with a Vicon motion capture system. As previously said, the observations of the laser spot are simulated. This was possible thanks to the fact that a reliable ground truth was provided together with the inertial data. In particular, given the true trajectory, we simulated the same planar surface described in the previous section. By having the true vehicle configuration it was possible to create the observations performed by the camera on the laser spot produced by a laser pointer as in the simulations (see the last paragraph in V-A.2).

Figure 9 displays the estimated α.

C. Experiments in progress

In this last section we describe our vehicle and sensors. The robot platform is a Pelican from Ascending Technologies equipped with an Intel Atom processor board (1.6 GHz, 1 GB RAM). Our sensor suite is composed by an Inertial Measurement Unit (3-Axis Gyro, 3-Axis Accelerometer), a monocular camera (Matrix Vision mvBlueFOX, F OV : 130 deg) and a Laser Module (SparkFun TTL, wavelength: 650nm, poweroutput: 0.45-0.8mW). The Laser module and the monocular camera are mounted on a fixed baseline, and the latter is calibrated using the Camera Calibration Toolbox for Matlab [START_REF] Bouguet | Camera Calibration toolbox for Matlab[END_REF]. The calibration between IMU and camera has been performed using the Inertial Measurement Unit and Camera Calibration Toolbox by Lobo [START_REF] Lobo | Relative pose calibration between visual and inertial sensors[END_REF], whereas the Camera-Laser Module calibration technique we used is described in [START_REF] Martinelli | Vision-Aided Inertial Navigation Using Virtual Features[END_REF]. By using this latter technique we obtained the following values: The IMU provides measurements update at a rate of 100Hz, while the camera framerate is 13Hz. We are currently equipping our lab with a motion capture system to perform experiments. Fig. 12. The Pelican quadcopter equipped with a monocular camera, a laser module and passive markers which will be used to collect sensors data and ground truth.

VI. CONCLUSION

In this paper we considered an aerial vehicle equipped with a monocular camera and inertial sensors. Additionally, a laser pointer is mounted on the vehicle and it produces a laser spot.

The laser spot is observed by the monocular camera and it is the unique point feature used in the proposed approach. We focused our attention to the case when the vehicle moves in proximity of a planar surface and in particular when the laser spot belongs to this surface. The paper provided two main contributions. The former is the analytical derivation of all the observable modes, i.e. all the physical quantities that can be determined by only using the inertial data and the camera observation of the laser spot during a short timeinterval. Specifically, it is shown that the observable modes are: the distance of the vehicle from the planar surface; the component of the vehicle speed, which is orthogonal to the planar surface; the relative orientation of the vehicle with respect to the planar surface; the orientation of the planar surface with respect to the gravity. Once the observed modes have been derived, a local decomposition of the original system has also been provided. This decomposition separates the observable modes from the rest of the system and allowed us to introduce a simple recursive method to perform the estimation of all the observable modes (second contribution). The use of a virtual laser spot feature overcomes the limits of feature tracking algorithms and makes our approach suitable to work even in dark or featureless environment. The method is validated by using synthetic data. The validation with real data is in progress. We presented a low-cost low-weight sensor suite and a low computational complexity approach in the framework of aerial navigation. It can be integrated in the framework of autonomous takeoff and landing, safe touchdown and low altitude manoeuvres. However, we want to emphasize that both the paper contributions are very general and can be applied in other frameworks. In particular, in all the environment where GPS is denied and where the most of objects have planar surfaces (e.g. in an indoor or city-like environment). For instance, these paper contributions could be used in the framework of humanoid robotics (where visual and inertial sensing are often adopted and the navigation usually occurs in an indoor environment).

A very important extension of the proposed strategy is to consider laser patterns in order to improve the precision. So far, the simplest case of a point-spot was considered. The performance of the considered system depends on the extrinsic parameters of the laser module with respect to the distance from the planar surface. Therefore we want to consider the option of a laser module with a servomotor in order to improve precision according to the trajectory. We finally want to implement this approach on our quadrotor platform to perform autonomous landing and safe touchdown tasks.

Fig. 1 .

 1 Fig.1. Quadrotor equipped with a monocular camera, IM U and a laser pointer. The laser spot is on a planar surface and its position in the camera frame is obtained by the camera up to a scale factor.

Fig. 2 .

 2 Fig. 2.

Figure 3 Fig. 3 .

 33 Fig. 3. A typical vehicle trajectory in our simulations.

Figures 4 aFig. 4 .

 44 Figures 4 a and b display the estimated α respectively in the case without and with bias. The values of the biases adopted in our simulations are: Ω bias = [0.03 0.03 0.03] T (deg s -1 ) and A bias = [0.03 0.03 0.03] T (ms -2). As expected, the estimation in presence of bias becomes worse. However, the error on the estimated α in presence of bias is smaller than 1 deg.

Figures 5 ,Fig. 5 .Fig. 6 .Fig. 7 .

 5567 Figures 5, 6, 7 and 8 a and b display respectively the estimated P , R, v o and d. In each figure, both the cases of unbiased and biased inertial measurements are displayed. We initialized the filter by using a value of the initial observable state which differs from the ground truth by a relative error in the range [10, 20]%.

Fig. 8 .

 8 Fig. 8. Estimated d in absence (a) and in presence of bias (b) on the inertial data. Blue dots indicate ground true values while red discs indicate the estimated values.

  Figures 10 a and b display the estimated P and R and figures 11 a and b display the estimated v o and d. All the observable modes are estimated with very good accuracy. Additionally, we remark that the convergence of the filter occurs in less than half second for all the observable modes.

Fig. 9 .Fig. 10 .Fig. 11 .

 91011 Fig. 9. Estimated α in the experiment. Blue dots indicate ground true values while red discs indicate the estimated values.

  θ = 47.1 deg, φ = -3.1 deg and L x = -0.146 m, L y = -0.005 m. The resulting L is 0.100 m.

This is obtained also by knowing that the laser pointer has the same orientation as the camera and that it is located at the position [L, 0, 0]

This work was supported by the European Project FP7-ICT-2007-3.2.


Hugo Talbot

Christian Duriez

Hadrien Courtecuisse

Jatin Relan

Maxime Sermesant

Stéphane Cotin

Hervé Delingette

Towards Real-Time Computation of Cardiac Electrophysiology for Training Simulator

This work aims at developing a training simulator for interventional radiology and thermo-ablation of cardiac arrhythmias. To achieve this, a real-time model of the cardiac electrophysiology is needed, which is very challenging due to the stiff equations involved. In this paper, we detail our contributions in order to obtain efficient cardiac electrophysiology simulations. First, an adaptive parametrisation of the Mitchell-Schaeffer model as well as numerical optimizations are proposed. An accurate computation of both conduction velocity and action potential is ensured, even with relatively coarse meshes. Second, a GPU implementation of the electrophysiology was realised in order to decrease the computation time. We evaluate our results by comparison with an accurate reference simulation using model parameters, personalized on patient data. We demonstrate that a fast simulation (close to real-time) can be obtained while keeping a precise description of the phenomena.

Introduction

Cardiac arrhythmias are characterised by a pathological electrical activity in the myocardium (heart muscle) which can be lethal. Catheter thermo-ablation is a minimally invasive technique that can prevent fibrillation by removing the substrate or the trigger of such a phenomenon through the ablation of the responsible cardiac cells. This procedure is performed by highly skilled and experienced interventional cardiologists, yet the success rate remains limited due to the complexity of such pathologies. Moreover, due to population ageing, an important increase is predicted in such arrhythmias, without an associated increase in interventional cardiologists. Therefore, there will be an important need for novel tools in learning, training and planning, such as simulators.

One important missing step towards the development of such a simulator is the ability to simulate the electrophysiology of the heart in real-time. The electrical wave propagating inside the cardiac walls corresponds to ion exchanges between the cells and through the cell membrane. The myocardial cells are polarised, i.e. there is a potential difference between the inside and the outside of the cells. This potential difference is called action potential (AP) or transmembrane potential. Before a thermo-ablation, cardiologists first carry out a diagnostic electrophysiological study (programmed electrical stimulation) allowing a definitive diagnosis. From the acquired data, the ablation intervention can be planned. We focus here on the training of cardiologists for these catheter-based ablation of cardiac tissues in the context of cardiac arrhythmias.

Research on cardiac electrophysiology models is very active. The proposed models can be divided into three different classes: (i) biophysical models, which are complex models including the different ionic concentrations and channels, involving many parameters and simulating the electrophysiology close to the cellular scale [START_REF] Tusscher | A model for human ventricular tissue[END_REF]; (ii) phenomenological models, which are simplified models [START_REF] Richard | Impulses and physiological states in theoretical models of nerve membrane[END_REF], [START_REF] Aliev | A simple two-variable model of cardiac excitation[END_REF], [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with fiber rotation[END_REF], [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF] derived from the biophysical models, involving less parameters and capturing the AP shape and its propagation at the organ scale; (iii) Eikonal models [START_REF] Keener | An eikonal-curvature equation for action potential propagation in myocardium[END_REF], which correspond to static non-linear partial differential equations for the depolarization time derived from the previous models. These models enable to simulate wave propagation but cannot accurately account for complex physiological states and the parameters have no direct physiological meaning.

This paper addresses the problem of real-time simulation of electrophysiology thanks to a new modelling approach adapted to efficient computational strategies. The approach relies on one of the phenomenological models from Mitchell-Schaeffer (MS) [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF]. The following sections presents the MS model, the model parameters and the numerical settings. Then we detail how we adapt the model parameters and numerical settings in order to recover patient-specific features, namely the conduction velocity (CV) and the action potential duration (APD, denoting APD90). To further improve our computation times, we propose an efficient GPU implementation with CUDA. Finally, we present our first results on patient-specific data and we conclude with some ideas for future work.

2 Real-Time Cardiac Electrophysiology Simulation

Discretisation of the Cardiac Electrophysiology Model

Mitchell-Schaeffer Model. For the last fifty years, numerous models were proposed for the cardiac cell AP. Regarding the features of the model categories presented previously, the phenomenological models meet best our requirements in terms of computational efficiency at the organ scale. The model that we chose is the MS model [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF] because of the following reasons:(i) it has only 5 parameters that we detail below, (ii) each parameter has a physiological meaning and (iii) it provides a better estimation of the AP compared to other phenomenological models (as the Aliev-Panfilov model [START_REF] Aliev | A simple two-variable model of cardiac excitation[END_REF]).

The MS model is a two-variable model derived from the Fenton Karma model [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with fiber rotation[END_REF]. The equations describing the model are written in 1:

               ∂ t u = div(D∇u) + zu 2 (1 -u) τ in - u τ out + J stim (t) ∂ t z =        (1 -z) τ open if u < u gate -z τ close if u > u gate (1)
where u is a normalized transmembrane potential1 and z is the gating variable associated to the sodium ion influx, thus depicting the repolarization phase. The extra-cellular potential can be recovered using the method described in [START_REF] Chhay | How to compute the extracellular potential in electrocardiology from an extended monodomain model[END_REF]. The diffusion term is defined by an 3x3 anisotropic diffusion tensor D = d•diag(1, r, r) so that the planar conduction velocity in the fiber direction is 2.5 times greater than in the transverse plane (r = 1

(2.5) 2 ). d is the diffusion coefficient. The parameters τ in and τ out define the repolarization phase whereas the constants τ open and τ close manage the gate opening or closing depending on the changeover voltage u gate . The term J stim (t) is the stimulation current applied in the pacing area. The default values (describing the common action potential) of these parameters are given in [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF]. The initial conditions of our simulation were set as Dirichlet conditions. The pacing area (u 0 = 1 and z 0 = 1) has been extracted from the depolarization time map.

Patient-Specific Spatial Discretisation. The cardiac geometry that we used in this study is an anatomy segmented from a patient MRI. This patient suffers from chronic ischemia and a left bundle branch block (noted LBBB).

The method used for the spatial discretization of the model is the finite element method which requires a volumetric mesh of the myocardium. The volumetric bi-ventricular geometry has been meshed with linear tetrahedra using CGAL (www.cgal.org) algorithms. The edge size used is about dx = 4 mm and implies 65547 tetrahedra which is already substantial for real-time simulation. Usually the usual edge size required for cardiac electrophysiology is dx < 0.1 mm. However the element size has to be fine enough to correctly model the wave propagation, especially in the grey zones around scars. Scar regions are known to be very unstable electrical zones and can be responsible for cardiac reentries. These areas therefore need to be discretized with accuracy.

A LBBB abnormal conduction implies a late activation of the left ventricle. Regarding this pathology, the edge size of 4 mm should be sufficient for our simulation since only the times of depolarization and repolarization have to be computed accurately. However to model a possible ventricular tachycardia induced by an isthmus the edge size should be finer. As demonstrated in [START_REF] Lamecker | Adaptive tetrahedral meshing for personalized cardiac simulations[END_REF], a 0.8 mm edge size must be used for accurate simulation. Time Integration Scheme. For all our simulations, we use the second-order semi-implicit solver called Modified Crank Nicholson/Adams Bashforth (MC-NAB), detailed in [START_REF] Ethier | Semi-implicit time-discretization schemes for the bidomain model[END_REF]. This solver defines an implicit integration of the diffusion term and an explicit integration of the ionic current. The construction of the matrix system leads to a symmetric definite positive case that we can solve using many existing algorithms.

In our GPU simulations, we work with a conjugate gradient (CG) to solve the resulting linear system (written Ax = b). To improve the convergence rate, we use a preconditioner Jacobi. The Jacobi method computes the inverse matrix of the diagonal of the system matrix A. Since our matrix A is diagonal dominant, the use of this preconditioner is straightforward. The factorization provided by the Jacobi method can be updated when needed, thus allowing to interact with the model (e.g. thermo-ablation).

Numerical Study

To achieve our performance goal, the numerical settings of our simulation need to be optimized. In the literature about cardiac electrophysiology simulation, it is shown that both element size and time step have to be small enough to capture the cardiac phenomenon. Usual values for the edge size are defined such that dx ≤ 0.1 mm and the time step is lower than dt ≤ 0.01 ms (see [START_REF] Pathmanathan | The significant effect of the choice of ionic current integration method in cardiac electro-physiological simulations[END_REF]). However meshing a heart geometry with such small elements would imply millions of tetrahedra. Such small elements as well as very small time steps prevent realtime simulations. To achieve real-time performances, we have to use larger time steps and larger finite-element. In this part, we study the impact of using large elements, the limitation in terms of time step as well as solutions to reach realtime performances.

Locally-Adapted Model Parameters

To ensure the reliability of the training simulation, we need a relatively accurate computation of the CV. Nevertheless using large elements will affect the diffusion effect, i.e. the wave propagation and its CV. We studied the influence of the element size on the CV with the lumping using the MS model on a 3D bar (15cm × 2.5cm × 2.5cm, regularly meshed) on which the wave is propagating along the bar axis with a planar front. A similar study was conducted in [START_REF] Pathmanathan | The significant effect of the choice of ionic current integration method in cardiac electro-physiological simulations[END_REF] where they also considered other integration method for the ionic term than the lumping method. The results are compared to a very accurate computation of the MS model in 1D using Matlab with a time step dt = 10 -5 s and a spatial step dx = 5 • 10 -6 m that provides a CV reference value: CV ref = 0.5124 m/s. In the Fig. 1, the results show that the CV decreases when the mesh becomes coarse.

A way to compensate this integration error is to artificially increase the diffusion coefficient d. Based on this idea, a personalization of the diffusion coefficient d and the opening and closing time constants τ open and τ close has been computed to fit measured patient data, using the personalization method of [START_REF] Relan | Quantitative comparison of two cardiac electrophysiology models using personalisation to optical and mr data[END_REF]. This personalization step gave us values of the diffusion coefficient for each tetrahedron, as well as values of opening and closing time constants on each vertex.

Time-Step Limitations

Even if the time step does not affect the electrophysiology features, the explicit formulation of the ionic term implies a stability condition on the time step. This limitation of the time step has already been pointed out in [START_REF] Ethier | Semi-implicit time-discretization schemes for the bidomain model[END_REF] and [START_REF] Coudière | Stability and convergence of a finite volume method for two systems of reaction-diffusion equations in electro-cardiology[END_REF]. Since the diffusion is computed implicitly the limiting term is the ionic term. Following [START_REF] Coudière | Stability and convergence of a finite volume method for two systems of reaction-diffusion equations in electro-cardiology[END_REF], we should define dt ≤ 1 |inf (∂f )| = τinτout τin+τout = 0.286 ms, if we would use a semi-implicit Euler solver (f denoting the ionic term).

Using the MCNAB, we measured dt ≤ 0.59 ms. For more stability, we will use dt = 0.4 ms.

Optimization on CG tolerance As detailed in 2.1, we use a CG to solve our system. Our solver stops when the difference between the results given by two successive iterations is below the tolerance, noted tol. The tolerance tol ref = 10 -10 ensures us to reach the "exact" solution. However to improve the computation performances, we computed the optimal tolerance value. We found the tolerance tol opt = 10 -6 giving a L 2 -error about err = 3.79 • 10 -5 on the action potential field, which is fully acceptable. In the results part, we will show performances using both tolerances tol ref = 10 -10 and tol opt = 10 -6 in section 3.

GPU Implementation

The GPU architecture consists in several multi-processors able to carry out highly parallel tasks independently. The complexity of GPU programming results in defining an optimal distribution of the threads and minimizing the memory access latency. In our approach, we relied on the CUDA toolkit (dedicated to NVidia's GPUs) to develop the GPU version of the method. Similar results could be obtained using OpenCL and other GPU models.

The implementation of the ionic term of MS is based on classical parallelization methods. Using the lumping integration method, the ionic term is computed on each vertex. Therefore each thread is dedicated to one vertex and computes the contribution of the MS term for this vertex. However our computation strategy ensures a tiled access in memory in our GPU code.

The diffusion term div(D • grad(u)) is more complicated to implement in parallel. The contribution of the diffusion term is computed from the edges and summed on each vertex. In a parallel computation, this algorithm can lead to writing conflicts: two threads solving two adjacent edges could write on the same point simultaneously. New GPUs supporting CUDA 2.0 handle now these atomic operations. Nevertheless a solution to this problem is detailed in [START_REF] Allard | Implicit fem solver on gpu for interactive deformation simulation[END_REF] which is twice faster than atomic functions from CUDA. We decided to adapt this algorithm, originally designed for deformable finite element equations, to our diffusion model. This technique consists in first computing the neighborhood of each node. Then the contribution of the diffusion is computed and saved on all the edges. Finally, the contribution on each vertex is obtained by accumulating the values computed on the edges with the neighborhood information. In this last accumulation step, a parallel reduction is carried out using several threads by vertex as explained in [START_REF] Allard | Implicit fem solver on gpu for interactive deformation simulation[END_REF].

Our GPU implementation also includes a conjugate gradient and a Jacobi preconditioner that are used to solve very quickly the electrophysiology model on the whole mesh. Only single float precision has been considered since it is sufficient regarding the electrophysiology accuracy. Moreover simulations using single float precision are about two times faster than simulations using double precision as detailed in [START_REF] Bartocci | Toward real-time simulation of cardiac dynamics[END_REF].

Finally, it must be stressed that this GPU implementation is a novel application of an efficient method originally proposed by [START_REF] Allard | Implicit fem solver on gpu for interactive deformation simulation[END_REF] for deformable finite element computations. This parallel implementation will be crucial in our performance results.

Results and Discussion

The local adaptations on the MS model that we detailed in 2.2 will still generate some errors. We now want to evaluate this error in comparison with our ground-truth data. Moreover a GPU implementation has been done to improve the computation time of our simulation. The performances using the optimized model parameters will also be assessed with this implementation.

Depolarization Times and APD Error Maps

We compared our results with the reference simulation using very small time step (dt = 10 -5 s) and model parameters fitting patient data. We computed an average error on depolarization time of 4.05 ms (∼ 5.48%) and an average error on APD of 19.1 ms (∼ 5.62%). These results shows that the errors are bounded. As future work, the restitution analysis of this patient-specific case seems to be compulsory to fully validate our simulation.

GPU Performances

Using the optimization presented previously would not be sufficient for fast computations. However, when combined to a GPU implementation, we obtain a significant speedup. Table 3.2 sums up computation times and other metrics of performance on the CPU and on the GPU (GeForce GTX580 with 512 cores). The processor used for CPU computations is a Intel Xeon W3550. From this table, we clearly notice the performance gain offered by the parallel implementation. The GPU simulation is more than 10 times faster than the CPU one. The computation time spent in the computation of diffusion and MS becomes negligible. Now the limitation of our simulation is now the solver, since it is the most time-consuming part of our computation. We notice that reasonably decreasing the tolerance of the CG significantly improves the performances.

To discuss these results, it can be stressed that 65547 tetrahedra represent a large amount of computation in a real-time context. A first solution would be to use on a coarser mesh. It must be stressed that the electrical wave of this patient was especially fast implying a small time step. Next generation GPU cards could also be most helpful.

User Interaction

A simulation dedicated to a training system assumes that cardiologists will interact with it. Our simulations already allow to stimulate any area of the heart in real-time (constraint with u = z = 1). This allows to simulate the stimulation done by the surgeon before the ablation in order to set the definitive diagnosis. The procedure of thermo-ablation is also handled by our computations. The cardiac cells treated by thermo-ablation can not conduct the current anymore. These cells can be seen as zero conductivity area. To take this change into account, we update the factorization computed by the Jacobi method. This thermo-ablation step can therefore be done without affecting the computation times.

Conclusion

In the coming years, the growing number of cardiology students implies to improve the training methods. Medical training simulators are an efficient solution to answer this need. Their main advantage is that simulators enable to train on virtual patient, thus avoiding to early operate on patients. Moreover such simulators could be a way of improving the course of medical studies.

In this paper, we presented some contributions to speed up simulations of cardiac electrophysiology. First, we locally modified time constant parameters from the Mitchell-Schaeffer model in order to use larger time steps and a coarser mesh. Then, we calculated the optimal time step and tolerance value for our iterative solver. Finally, we implemented a GPU version of the code modeling the diffusion and the Mitchell-Schaeffer model. With these optimizations, we reached an interactive simulation close to real-time based on 3D patient-specific data. This work is a very encouraging step toward the development of a training simulator.

In this work, only monocycle simulations have been considered, i.e. simulating only one cardiac cycle. A multicycle study needs to be carried out to analyze the restitution curves of our model and to validate our simulations. This paper presents a personalized simulation using data from a patient suffering from LBBB.

In the future, we want to address the challenging issue of simulating the electrophysiology of patients suffering from a ventricular tachycardia with an isthmus. This would lead to compute models with finer elements around the regions of interest (scars and isthmus) in order to accurately capture the electrical activity around these regions. To do so, we plan to implement innovative numerical method to use finer meshes without altering our performances.

Fig. 1 .

 1 Fig. 1. (Left) initial conditions on the 3D bar (Right) Evolution of the CV (m/s) against element size (mm).

Fig. 2 .

 2 Fig. 2. Two different views of an electrical stimulation (with a catheter) applied by the user on a two-chamber geometry (ventricles).

Table 1 .

 1 Performance comparison between CPU and GPU.

		CPU	GPU	GPU (using tolopt)
	Mean computation time for 25.4 one time step (ms)	3.89	2.17
	real-time ratio (slower than real-time)	75.8	7.04	2.95

The MS is a monodomain model since it is expressed according to the transmembrane potential whereas bidomain models depends on both intra-and extra-cellular potential.