Abimael Gonzalez-Hernandez
email: abimaelgh@comunidad.unam.mx

PhD Alexandre Charlet
email: acharlet@unistra.fr

PhD Abimael González-Hernández

Oxytocin, GABA and TRPV1, the analgesic triad?

Keywords: Oxytocin, GABA, and TRPV1, the Analgesic Triad?

come

ABSTRACT

Oxytocin (OT) has emerged as an interesting endogenous molecule with potential analgesics effects modulating nociception at peripheral, spinal and supraspinal sites in the nervous system during physiological and pathological pain conditions. Experimental data suggest that physiological OT could modulate the nociceptive input activating descending endogenous mechanisms. Furthermore, pharmacological data strongly support the notion that this neuropeptide could be used as analgesic during pathological pain (i.e. neuropathic pain). Indeed, experimental pain models that resembles some pathophysiological aspects of chronic pain, showed that intrathecal OT relieves this type of pain, pointing out the relevance of oxytocinergic mechanisms modulating the nociceptive during pathological conditions. In this regard, although much effort has been made to elucidate the mechanisms/pathways which OT engage to induce its effects, the key downstream mechanisms remain unsolved. Canonically OT activates its cognate receptor (OTR), a seven-transmembrane receptor preferentially coupled to Gq proteins which in turn mobilizes the intracellular Ca 2+ leading to an array of possible metabotropic effects. Currently, we know that at spinal level the GABAergic recruitment seems to be essential in the spinal OT effect, whereas at peripheral level direct OTR activation diminishes the nociceptive input. Nevertheless, a pair of basic studies show us that OT could affect the pain transmission by novel unexplored pathways, first could be modulating directly (desensitizing) the activity of TRPV1 and second activating indirectly GABAB receptors which in turn restore the TRPV1 function in pathological pain. Together, these studies remark the pleiotropic nature of the mechanisms activated by the OT to modulate pain transmission, and in a world where the abuse of opioids drugs to treat pain becomes as a public health issue, the discovery of new molecules with novel analgesic mechanisms emerges as an interesting area that needs to be further explored.

MAIN TEXT

The hypothalamic nonapeptide oxytocin (OT) has several physiological functions, ranging from lactation to social attachment [START_REF] Gimpl | The oxytocin receptor system: structure, function, and regulation[END_REF]. More recently, a cumulus of evidences revealed that OT inhibits pain transmission at peripheral [START_REF] Juif | Long-lasting spinal oxytocin analgesia is ensured by the stimulation of allopregnanolone synthesis which potentiates GABAA receptor-mediated synaptic inhibition[END_REF][START_REF] De Araujo | Selenoether oxytocin analogues have analgesic properties in a mouse model of chronic abdominal pain[END_REF][START_REF] Tzabazis | Oxytocin receptor: Expression in the trigeminal nociceptive system and potential role in the treatment of headache disorders[END_REF][START_REF] González-Hernández | Peripheral oxytocin receptors inhibit the nociceptive input signal to spinal dorsal horn wide-dynamic-range neurons[END_REF], spinal and supraspinal levels [START_REF] Eliava | A new population of parvocellular oxytocin neurons controlling magnocellular neuron activity and inflammatory pain processing[END_REF][START_REF] Poisbeau | Oxytocin signaling in pain: cellular, circuit, system, and behavioral levels[END_REF]. These data support the prospect that spinal OT could be translated to the clinical pain practice [START_REF] Rash | Oxytocin and pain: a systematic review and synthesis of findings[END_REF][START_REF] Eisenach | Phase 1 safety assessment of intrathecal oxytocin[END_REF][START_REF] Condés-Lara | Successful pain management with epidural oxytocin[END_REF]. While the receptor(s) and intracellular mechanisms responsible for the OT-induced analgesia are under scrutiny, the main dogma rely on the activation of OT receptor (OTR). Indeed, ex vivo studies using capsaicin-sensitive dorsal root ganglia (DRG) neurons showed that OT induce membrane hyperpolarization and consequently inhibits the nociceptive transmission by Ca 2+ -dependant mechanisms [START_REF] Hobo | Oxytocin inhibits the membrane depolarization-induced increase in intracellular calcium in capsaicin sensitive sensory neurons: a peripheral mechanism of analgesic action[END_REF]. In addition, in vivo electrophysiological recordings of the second order wide dynamic range (WDR) neurons, ex vivo patch-clamp and behavioural analysis suggest that OT inhibits the nociceptive input either directly [START_REF] Eliava | A new population of parvocellular oxytocin neurons controlling magnocellular neuron activity and inflammatory pain processing[END_REF] or via the enhancement of the GABAergic transmission [START_REF] Rojas-Piloni | GABA-mediated oxytocinergic inhibition in dorsal horn neurons by hypothalamic paraventricular nucleus stimulation[END_REF][START_REF] Breton | Oxytocin-induced antinociception in the spinal cord is mediated by a subpopulation of glutamatergic neurons in lamina I-II which amplify GABAergic inhibition[END_REF], probably by inducing neurosteroids synthesis [START_REF] Juif | Long-lasting spinal oxytocin analgesia is ensured by the stimulation of allopregnanolone synthesis which potentiates GABAA receptor-mediated synaptic inhibition[END_REF]. In all cases, OTR were involved. However, the contribution of vasopressin V1A receptor (V1AR), for which OT has a physiologically-relevant affinity, cannot be rule out [START_REF] Schorscher-Petcu | Oxytocininduced analgesia and scratching are mediated by the vasopressin-1A receptor in the mouse[END_REF][START_REF] González-Hernández | Oxytocin and analgesia: future trends[END_REF][START_REF] Manzano-García | The role of peripheral vasopressin 1A and oxytocin receptors on the subcutaneous vasopressin antinociceptive effects[END_REF]. In any case, these findings point out the relevance of GABAergic transmission in the OT-induced analgesia (Figure 1).

In this context, two recent interesting studies [START_REF] Nersesyan | Oxytocin modulates nociception as an agonist of painsensing TRPV1[END_REF][START_REF] Sun | Oxytocin relieves neuropathic pain through GABA release and presynaptic TRPV1 inhibition in spinal cord[END_REF] add a new layer on the analgesic mechanisms of the OT, suggesting than apart from the role of OTR and its downstream mechanism, this peptide also interact (directly or indirectly) with the Transient Receptor Potential Vanilloid type-1 (TRPV1) at peripheral and spinal level. Certainly, since this channel has a physiological function as a detector of noxious stimuli [START_REF] Caterina | Impaired nociception and pain sensation in mice lacking the capsaicin receptor[END_REF], and seems to be crucial for the transition from acute to chronic pain [START_REF] Choi | Emerging role of spinal cord TRPV1 in pain exacerbation[END_REF], there is a strong interest to develop molecules modulating TRPV1 to treat pain. Nersesyan et al. first showed that aversive behaviours induced by intracutaneous (i.c.) injection of capsaicin into the hindpaw was reduced in animals concomitantly treated with i.c. OT.

While one can propose that the OT effect is due to OTR activation at peripheral nerve endings [START_REF] González-Hernández | Peripheral oxytocin receptors inhibit the nociceptive input signal to spinal dorsal horn wide-dynamic-range neurons[END_REF]; Figure 1), the authors challenged an audacious hypothesis involving direct OT interaction with the TRPV1 channel. Indeed, since OT could modulate the influx of extracellular Ca 2+ , the role of voltage-gated or ligand-gated channels may play a role, as previously suggested [START_REF] Sanborn | Molecular mechanisms regulating the effects of oxytocin on myometrial intracellular calcium[END_REF]. Under this line of reasoning, they showed that OT at μM concentration, elicits Ca 2+ influx through TRPV1 expressed in HEK293 cells. Furthermore, the hypothesis proposed was strengthened considering that both the TRPV1 antagonist capsazepine and a TRPV1-specific small interfering RNA (siRNA) abolished the OT-induced Ca 2+ influx. These results were then confirmed by their replication in nociceptive DRG neurons isolated from wild-type animals and the modest OT-induced analgesia observed in TRPV1 knock-out mice. Further electrophysiological, planar lipid bilayer and in silico experiments showed that OT act as a partial agonist of TRPV1 channels and induce a strong desensitization of these channels under nociceptive stimulation. This pioneer study strongly suggests that at peripheral level OT-induced analgesia pass, at least in part, through direct TRPV1 interaction.

To further explore if OT-induced analgesia in neuropathic pain may occur via TRPV1, Sun et al. first confirmed that intrathecal OT diminishes the aversive behaviour in neuropathic animals [START_REF] Miranda-Cárdenas | Oxytocin and electrical stimulation of the paraventricular hypothalamic nucleus produce antinociceptive effects that are reversed by an oxytocin antagonist[END_REF]. Most interestingly, the cerebrospinal fluid OT concentration was increased, in accordance with the spinal OT increase measured after a peripheral inflammatory sensitization [START_REF] Juif | Long-lasting spinal oxytocin analgesia is ensured by the stimulation of allopregnanolone synthesis which potentiates GABAA receptor-mediated synaptic inhibition[END_REF]. Furthermore, while the spontaneous excitatory post-synaptic currents (sEPSC) recorded in spinal cord superficial layers were not altered by OT, it prevented the capsaicin-induced increases of sEPSC, an effect blocked by the GABAB receptor antagonist saclofen, suggesting a presynaptic effect. Hence, OT enhance the spontaneous inhibitory postsynaptic currents (sIPSC) mediated by GABA. These data suggest that OT recruit GABAergic interneurons which in turn activates GABAB receptors and consequently inhibits the capsaicinsensitive neurons. It is interesting to note that TRPV1 and GABAB receptors were found to be coexpressed in primary nociceptive neurons. In addition, the TRPV1 spinal cord expression was enhanced in neuropathic rats, an up-expression corrected by intrathecal OT treatment. Accordingly, the OT-induced analgesia was partially reduced in TRPV1 knock-out mice. All in all, this study goes in the same direction than Nersesyan et al. and point that OT-induced analgesia may rely on both GABA release and modulation of TRPV1 activity.

During the past decades, OT slowly emerged as an important mediator of endogenous analgesia [START_REF] Poisbeau | Oxytocin signaling in pain: cellular, circuit, system, and behavioral levels[END_REF]. Thus, the relevance of oxytocinergic mechanisms implicated in the modulation of nociceptive transmission are under intense scrutiny. These two studies highlight a novel and interesting mechanistic explanation about how OT modulates the nociceptive transmission at peripheral and spinal level under physiological and pathological pain conditions. Although at spinal level the GABAergic enhancement induced by OT has been reported by several research groups, Sun et al., demonstrates, in line with the observation that GABAB receptors revert the TRPV1 sensitization by non-canonical signalling under pathological pain [START_REF] Hanack | GABA blocks pathological but not acute TRPV1 pain signals[END_REF], that this GABA facilitation engage the function of post-synaptic GABAB receptors which in turn modulate the TRPV1 function during neuropathic pain. In addition, the discovery by Nersesyan et al. that OT could act as an allosteric modulator of TRPV1 highlight the multifaceted actions of OT as neuromodulator. Certainly, the finding that OT partially blocked the peripheral capsaicin-induced nocifensive behaviour could explain in part why blocking OTR at peripheral level only partially reverted the OT antinociceptive effect [START_REF] González-Hernández | Peripheral oxytocin receptors inhibit the nociceptive input signal to spinal dorsal horn wide-dynamic-range neurons[END_REF]). relevant, OT-mediated TRPV1 modulation. Since OTR was recently involved in functional heterodimers, for example with vasopressin or dopamine receptors [START_REF] Terrillon | Oxytocin and vasopressin V1a and V2 receptors form constitutive homo-and heterodimers during biosynthesis[END_REF][START_REF] De La Mora | Signaling in dopamine D2 receptor-oxytocin receptor heterocomplexes and its relevance for the anxiolytic effects of dopamine and oxytocin interactions in the amygdala of the rat[END_REF], the hypothesis of a functional heterodimer formation between OTR and TRPV might be interesting to further explore.

An intriguing point in

Finally, although the source of OT for its spinal effects are neuronal and seems to play a role in endogenous analgesia [START_REF] Eliava | A new population of parvocellular oxytocin neurons controlling magnocellular neuron activity and inflammatory pain processing[END_REF], it is more unclear at peripheral level. Nersesyan et al. suggest that the OT to stimulate TRPV1 is neuronal, but we need to keep in mind that other sources could explain the findings reported, such as keratinocytes [START_REF] Denda | Oxytocin is expressed in epidermal keratinocytes and released upon stimulation with adenosine 5′-[γ-thio] triphosphate in vitro[END_REF], which may release OT upon nociceptive stimulation or self-soothing and consequently induce pain relief [START_REF] Grinevich | Oxytocin: pain relief in skin[END_REF][START_REF] Walker | C-tactile afferents: Cutaneous mediators of oxytocin release during affiliative tactile interactions?[END_REF]; Figure 1). Certainly, the pleiotropic nature of this peptide to modulate pain transmission indicates that differential routes are followed.

In conclusion, a cumulus of evidence briefly touched here suggests that a triad between OT, GABA and TRPV1 exist. They shed light on new mechanisms involved in OT-induced analgesia, which finally seems to be much more complex than initially thought.

  Sun et al. study is how exactly, at spinal level, OT engages the intracellular machinery to modulate the TRPV1 expression. At least two mechanisms can be proposed. The first might be a direct action of OT on TRPV1, as suggested by Nersesyan et al. The second might be that the OT-induced activation of OTR lead to intracellular cascade able to modulate the expression or epigenetic modifications of TRPV1, as suggested by Sun et al. In the same line, it remains unclear how exactly OT modulates the TRPV1 functions and if it remains true at physiological concentrations. Indeed, while Nersesyan et al. clearly show a direct interaction between OT and TRPV1, neither them nor Sun et al. were able to reveal a strong, behaviorally

Figure legend Figure 1 .

 legend1 Figure legend

Acknowledgements

This work was supported by the Programa de Apoyo a Proyectos de Investigaión e Innovación Tecnológica PAPIIT-UNAM Mexico IA203117 (to AGH); NARSAD Young Investigator Grant 24821, ANR-DFG GR 3619/70 and PICS 07882 grants (to AC).

Conflicts of interest

The authors declare no competing interest.


F Cazals

C Mueller

C Robert

A Roth

Towards Morse Theory for Point Cloud Data

Keywords: Morse theory, Morse homology, point cloud data, multi-scale analysis, energy landscapes Théorie de Morse, Homologie de Morse, nuages de points, analyse multi-echelle, surfaces d'énergie

Morse theory provides a powerful framework to study the topology of a manifold from a function dened on it, but discrete constructions have remained elusive due to the diculty of translating smooth concepts to the discrete setting. Consider the problem of approximating the Morse-Smale (MS) complex of a Morse function from a point cloud and an associated nearest neighbor graph (NNG). While following the constructive proof of the Morse homology theorem, we present novel concepts for critical points of any index, and the associated Morse-Smale diagram. Our framework has three key advantages. First, it requires elementary data structures and operations, and is thus suitable for high-dimensional data processing. Second, it is gradient free, which makes it suitable to investigate functions whose gradient is unknown or expensive to compute. Third, in case of under-sampling and even if the exact (unknown) MS diagram is not found, the output conveys information in terms of ambiguous ow, and the Morse theoretical version of topological persistence, which consists in canceling critical points by ow reversal, applies. On the experimental side, we present a comprehensive analysis of a large panel of bi-variate and tri-variate Morse functions whose Morse-Smale diagrams are known perfectly, and show that these diagrams are recovered perfectly.

In a broader perspective, we see our framework as a rst step to study complex dynamical systems from mere samplings consisting of point clouds.

Morse theory provides a powerful framework to study the topology of a manifold from a function dened on it, but discrete constructions have remained elusive due to the diculty of translating smooth concepts to the discrete setting.

Consider the problem of approximating the Morse-Smale (MS) complex of a Morse function from a point cloud and an associated nearest neighbor graph (NNG). While following the constructive proof of the Morse homology theorem, we present novel concepts for critical points of any index, and the associated Morse-Smale diagram.

Our framework has three key advantages. First, it requires elementary data structures and operations, and is thus suitable for high-dimensional data processing. Second, it is gradient free, which makes it suitable to investigate functions whose gradient is unknown or expensive to compute. Third, in case of under-sampling and even if the exact (unknown) MS diagram is not found, the output conveys information in terms of ambiguous ow, and the Morse theoretical version of topological persistence, which consists in canceling critical points by ow reversal, applies. On the experimental side, we present a comprehensive analysis of a large panel of bi-variate and tri-variate Morse functions whose Morse-Smale diagrams are known perfectly, and show that these diagrams are recovered perfectly.

In a broader perspective, we see our framework as a rst step to study complex dynamical systems from mere samplings consisting of point clouds.

Introduction 1.Morse Theory and Morse Homology

From dierential topology to algorithms. In the introduction of his seminal memoir Analysis situs, published in 1895, Henri Poincaré was advocating the investigation of the relative position of points, lines and surfaces, regardless of their size, to understand objects beyond the third dimension. While these ideas are recognized to have contributed to the emergence of modern topology, it is probably fair to admit that the study of the variation of a function dened on a topological space, which complements purely topological considerations and is the core of Morse theory, has proved an equally fertile eld [START_REF] Milnor | Morse Theory[END_REF][START_REF] Bott | Morse theory indomitable[END_REF]. Morse theory is indeed concerned with the study of the connexion between a function dened on a manifold and the topology of that manifold, and Morse related constructions are key to countless applications, including morphological analysis of 3D shapes in general [BDFF + 08] and molecular shapes in particular [START_REF] Cazals | Molecular shape analysis based upon the Morse-Smale complex and the Connolly function[END_REF][START_REF] Thomas | Symmetry in scalar eld topology[END_REF], image segmentation via the so-called watershed transform [START_REF] Roerdink | The watershed transform: Denitions, algorithms and parallelization strategies[END_REF], study of so-called energy landscapes in biophysics [START_REF] Wales | Energy Landscapes[END_REF], high-dimensional data mining [START_REF] Gerber | Visual exploration of high dimensional scalar functions[END_REF], clustering [START_REF] Chazal | Persistence-based clustering in riemannian manifolds[END_REF], etc.

On the one hand and from a purely mathematical standpoint, Morse theory [START_REF] Milnor | Morse Theory[END_REF] and Morse homology [START_REF] Banyaga | Lectures on Morse Homology[END_REF] classically fall in the realm of dierential topology and homology theories. In particular, homology calculations can be carried out thanks to the so-called Morse-Smale-Witten chain complex which encodes incidences between critical points whose indices dier by one [START_REF] Banyaga | Lectures on Morse Homology[END_REF]. On the other hand, given the prominence of applied areas which have borrowed upon such ideas, it is especially interesting to see which routes which have been followed to implement Morse related ideas in a computer, given that the data manipulated are inherently discrete at some point. In the following, we briey list the main approaches, and discuss the key diculties faced along the way.

Morse theory based on simulation of dierentiability. Mimicking the smooth setting, e.g. by linearly interpolating a function dened on the vertices of a simplicial complex has also been Inria investigated [START_REF] Edelsbrunner | Hierarchical Morse complexes for piecewise linear 2-manifolds[END_REF]. The construction culminates with the Morse-Smale complex, namely critical points, their connections, and their (un)stable manifolds. However, the approach is bound to low dimensional data due to the exponential memory footprint of such complexes, and coping with degeneracies such as coalescing critical points remains a challenge.

Morse theory based on arbitrary-precision arithmetic. Given the dierential nature of classical Morse theory, it is tempting to design generic algorithm taking for granted arbitrary-precision arithmetic. While such algorithms can be certied under appropriate assumptions, see [START_REF] Chattopadhyay | Certied computation of planar morsesmale complexes[END_REF] for the 2D case, they seem of limited practical use due to the stringent assumptions on the numerics.

Discrete Morse theory. Motivated by the role of so-called CW complexes in Morse theory, Forman developed a purely combinatorial Morse theory, where both the critical elements and the ow operator are dened in terms of cellspractically in a cell or simplicial complex [START_REF] Forman | Morse theory for cell complexes[END_REF]. However, dening optimal discrete Morse function is in general NP-hard [START_REF] Lewiner | Optimal discrete Morse functions for 2manifolds[END_REF], and greedy algorithms, albeit conceptually simple, do not provide intrinsic constructions [START_REF] Cazals | Molecular shape analysis based upon the Morse-Smale complex and the Connolly function[END_REF]. Forman's Morse functions were recently used to propose a version of Morse homology [START_REF] Robins | Theory and algorithms for constructing discrete morse complexes from grayscale digital images[END_REF], which equips the connections of the Morse-Smale complex with the multiplicity information required to form the Morse-Smale-Witten chain complex, from which the homology is computed. Practically, the underlying cell complex being a cubical complex [START_REF] Kaczynski | Computational Homology[END_REF][START_REF] Wagner | Ecient computation of persistent homology for cubical data[END_REF]. However, for data originating from a smooth function sampled on a 3D grid, relating the critical cells to the critical points of the smooth function is an open question.

Morse theory in the context of Voronoi diagrams. Remarkably, functions used to dene ane Voronoi diagrams and α-shapes [START_REF] Edelsbrunner | Surface reconstruction by wrapping nite point sets in space[END_REF][START_REF] Buchin | Recursive geometry of the ow complex and topology of the ow complex ltration[END_REF] have proved amenable to Morse-theoretic constructions, in particular the ow complex [START_REF] Giesen | The ow complex: A data structure for geometric modeling[END_REF]. On the positive side, the ow complex enjoys the same properties as the CW complex associated with a classical Morse function, it proved important for reconstructing (non-manifold) shapes [START_REF] Dey | Critical points of the distance to an epsilon-sampling of a surface and ow based surface reconstruction[END_REF][START_REF] Cazals | Reconstructing 3d compact sets[END_REF], and it is also amenable to calculations (critical points and their incidence) in high dimension [START_REF] Giesen | A parallel algorithm for computing the ow complex[END_REF]. However, computing the (uns)stable manifolds remains challenging even in 3D [START_REF] Cazals | Robust construction of the threedimensional ow complex[END_REF], letting alone the size of the complex.

Morse theory for point cloud data. It is only recently that Morse theoretical constructions on point cloud data were investigated. For applications to distributed sensor networks analysis, a 2D Morse-Smale decomposition was proposed in [START_REF] Zhu | Topological data processing for distributed sensor networks with morse-smale decomposition[END_REF]. The algorithm handles domains of arbitrary topology, but the construction used to identify critical points (cut locus) is bound to 2D, and no comparison is reported against certied MS diagrams of smooth functions. The topography of graphs embedded in dD is studied in [START_REF] Weinan | The landscape of complex networks[END_REF]. But the notion of critical event proposed in not discussed in the context of smooth functions, and is actually suers from a contamination eect (see below, description of step 1).

Contributions. This work makes a stride towards an algorithmic version of Morse theory enjoying the following key properties:

Ease of implementation and ability to handle high-dimensional data.

Amenability to simple interpretations for data associated with smooth Morse functions.

Encapsulation of Morse homology calculations.

Guided by the rst wish, our approach is based on the analysis of point cloud data, that is, we only assume as input a nearest neighbor graph connecting the samples of a point cloud sampling a smooth Morse function, and we only resort to sorting and traversals on graphs. To match the second wish, a critical point is plainly a sample approximating the exact critical point RR n°8331 of the function studied. Similarly, the (un)stable manifold W s (p) (W u (p)) of a critical point p will consist of a sub-graph of the input graph.

Central to our constructions is the notion of bifurcation, which aims at capturing ambiguities of the ow, such as those found on the compactication of a region dened as the intersection W u (p) ∩ W s (q) of unstable and stable manifolds for a smooth Morse function. While these constructions are not backed-up by theorems at this stage, we believe that our work makes two important contributions. First, we provide insights on subtle behaviors of discrete ow lines, which are not faced in the smooth setting. Second, we report correct Morse-Smale diagrams for challenging Morse functions in 2D and 3D, the systematic comparison of a smooth Morse-Smale diagram against its approximation being undertaken for the rst time in this work, to the best of our knowledge.

Notations

Graphs and related constructions. A set denoted by a plain capital letter; a set of sets is denoted with a calligraphic font ie F = {F i }.

A graph dened by its vertex and edge sets is denoted in sans serif font, e.g.

G = (V [G], E[G]
).

An induced sub-graph G of G is a sub-graph of G such that for all pairs of vertices (u, v) of G , (u, v) is an edge of G i it is an edge of G. The sub-graph induced by a vertex set T = V [G ] ⊂ V [G] is dened similarly: (u, v) is an edge of G i it is an edge of G. A Hasse diagram is a directed acyclic graph.
In a directed graph, the nodes pointed by outgoing (resp. incoming) edges are called the successors (resp. ancestors) of that node. We consider a set of samples denoted P = {p i } i=1,...,n , each endowed with a real function value denoted by f (p i ). We assume that these samples are connected by a symmetric nearest neighbors graph (NNG), denoted G. The lower link of a sample consists of its neighbors in the NNG having a lower function value.

Discrete ow operator. We use the NNG connecting the samples to dene a descending pseudo-gradient graph G -by linking each sample p i to a sample q i , thus approximating the negative gradient of the function f . (Several options are explored in section 3.2.) A sample p i is called regular in G -if it has an outgoing edge, and singular otherwise. Singular vertices are local minima. The discrete ow operator is the operator iteratively following the neighbor (below in G -, above in G + ) of a sample. Using this ow operator, the origin and destination of a sample The ow lines (i.e.integral curves) associated with the negative gradient of the Morse function -∇f , and the associated orbits, that is the space of ow lines quotiented out by additive reparameterizations.

p i are denoted α -(p i ), ω -(p i ) if
The stable and unstable manifolds of the critical points, from which one denes regions of homogeneous ow W (q, p) = W u (p) ∩ W s (q). When W (q, p) = ∅, this set is used to dene a partial order q p between pairs of critical points connected by an orbit. The critical points of f and this partial order dene the phase diagram of the Morse function.

The compactication W (q, p) of of W (q, p). It can be shown that the boundary of W (q, p) consists of ow lines going through critical points intermediate between q and p. These ow lines are called boundary ow lines.

The moduli space M q,p = W (q, p)/R, namely orbits joining q to p. In particular, when λ q = k and λ p = k -1, M q,p has dimension zero, i.e. consists of a nite number of orbits.

The so-called Morse-Smale-Witten chain complex (MSW) connecting pairs of critical points q and p with λ q -λ p = 1, together with the signed multiplicities of these orbits, so as to dene the boundary operator from which the homology groups are computed. Note that when the multiplicity is omitted, we plainly talk about the Morse-Smale (MS) complex.

Overview

We aim at capturing the orbits of the moduli spaces W (q, p) with λ q -λ p = 1. As recalled above, connections between such critical points are found on the boundary of regions of homogeneous ow W (q, p). Because non empty regions W (q, p) ∩ W (q, p ) precisely delimit two regions of homogeneous ow, our strategy precisely consists of identifying such regions, thanks to bifurcations.

A sample p i bifurcates if it faces a ow ambiguity in the following sense: selected neighbors of his lower link ow to dierent index k -1 critical points. (Prosaically, two neighboring samples sitting across a ridge on a mountain dene a bifurcation provided that they ow to dierent local minima.) Albeit intuitive, the notion of bifurcation is dicult to capture: on the one hand and as just discussed, bifurcations are intimately related to ow ambiguities; on the other hand, while orbits are disjoint in the smooth setting, merges and splits are observed between their discrete analogs, due to the discrete nature of the sampling. In particular, the notions of conuence samples, forks, sticks and linkers to be introduced precisely aim at retaining the proper connections between critical points, which is non trivial due to these merges and splits.

Key steps

We now present the key steps and constructions, illustrated on Figs. 1 and 2, the details being provided in the supplemental section 6. Due to the aforementioned structure of the MSW chain complex, our construction iteratively builds slices of the MS complex, a slice precisely consisting of the index k critical points, their unstable manifolds, and the connections to index k -1 critical points. To this end, the kth iteration requires (i) the index k -1 critical points, (ii) a graph called the global thickening T[k -1], and (iii) the associated pseudo-gradient graph. For the rst iteration (k = 1), this information consists of the local minima, and the global thickening T[0] is the NNG, denoted G.

To ease the presentation, for each step, we rst present the main ideas, and proceed with subtle issues (identied by the pictogram ), whose understanding requires material presented in the supplemental section 6.

Step A bifurcating sample is such that two points from its lower link (bifurcation lower link, actually), ow to dierent index k -1 critical points. As discussed above, such a sample is located in a region that lies on the boundary of two regions of homogeneous ow. Critical points will be found amongst bifurcating samples. Through a bifurcating sample p i , there is a boundary ow line owing to a σ (k-1) by passing through a σ (k) rst. As an example on Fig. 1, p i ows to a local minimum by visiting an index one saddle rst.

A conuence sample is a sample such that its entire lower link consists of bifurcating samples. Intuitively, such samples witness the conuence of ow lines, an artifact of the discretization. These intuitive denitions actually hide one complication: when checking whether a sample p i is bifurcating, one should use the subset of its lower link consisting of non-bifurcating samples. If this precaution is not taken, because the samples are processed by increasing function value, a contamination eect is observed. That is, the number of bifurcating samples increases as a function of the distance to local minima of T[k -1] (the distance is counted in numbers of edges in the global thickening). This accounts for the notion of bifurcating lower link, which is a subset of the lower link for k ≥ 2 (supplemental section 6.3). We also dene the global stream B[k], namely the sub-graph of T[k -1] induced by the bifurcating samples. Note that the remaining samples, which are regular, are of no interest for further processing, since they belong to the stable manifolds of the index k -1 critical points.

Step 2: Computing anchors and local thickenings We wish to partition the samples of B[k] into equivalence classes with respect to tuples of index k -1 critical points (Fig. 1(C)). Consider such a subset, which we call a k-anchor, denoted A (k) . All the samples from B[k] having the same anchor dene the local thickening T(A (k) ) of this anchor. Intuitively, a local thickening is found across two (ore more) regions of homogeneous ow. Intuitively, if a local thickening has a unique local minimum, it represents the stable manifold of the critical point identied with this local minimum. If not, it represents the union of the stable manifolds associated to critical points sharing the same successors along the ow, identied by the critical points σ (k) k -1 dening its anchor.

Step 3: Computing monotonic sections and their incidence graph A local thickenings witnesses bifurcations with respect to the same index k -1 critical points, regardless of the variation of the function itself. We therefore compute a covering of T(A (k) ) with a nite set of monotonic sections M (α,ω) (T(A (k) )), each of them consisting of samples sharing the same α and ω limits within T(A (k) ) (Fig. 1(C)). Furthermore, the edges belonging to the global stream but not to any monotonic section are used to build the graph G (k) M whose nodes are monotonic sections (but so-called linkers, see below), with one edge for two incident such sections.

Inria

Conuence samples do not bifurcate by construction (Figs. 6 and7, but monotonic sections consisting of such samples (except their boundary) are important since following them downstream brings to regions where monotonic sections associated to dierent anchors split. Such monotonic sections are called linkers, and are excluded from G (k) M (supplemental section 6.5). In fact, linkers allow maintaining the connections between critical points of index k and index k -1 (supplemental section 6.6). Note that linkers do not contribute critical points. Note also that linkers account for the denition of bifurcating sets using conuence samples, as seen in step 1.

Step 4: Simplifying the incidence graph of monotonic sections We wish to transform the incidence graph G (recall that the function value decreases along a ow line). Because each monotonic section has a dening anchor A (k) (a tuple of index k -1 critical points), and since anchors enjoy a partial order for the inclusion, the edition of

G (k) M into G (k)
M is based on two partial orders: the inclusion partial order ≺ i just discussed, and a function value related partial order denoted ≺ h , which qualies the relative height of monotonic sections which are connected by a path in G (k) M . Given two elements to be compared, two partial orders yield ve dierent cases: the two elements can be compared (two cases), only one partial order applies (two cases), the elements cannot be compared by any of the two partial orders (one case). Accordingly, the edges of

G (k) M are classied as canonical, composite, H-non-I, I-non-H, or non-hierarchical. The graph G (k) M is then edited into G (k)
M by keeping canonical edges only, so that G (k) M precisely consists of the edges (and associated incident nodes i.e. monotonic sections) satisfying the following two conditions, which dene the partial order denoted ≺ ih :

M i (A (k) u ) ≺ i M j (A (k) v ) and M i (A (k) u ) ≺ h M j (A (k) v ). (1)
The diculty consists of handling the edges which are not canonical. These edges are indeed related to artifacts of the discretization in general. For example, composite edges are related to the existence of monotonic sections called forks and sticks (supplemental Fig. 6), namely incident monotonic sections such that the two partial orders ≺ h and ≺ i disagree. When detected at stage k, these sections actually yield conuence samples at the k + 1th iteration. Other situations account for the three remaining types of edges (supplemental section 6.6).

Step 5: Computing the kth slice of the Morse-Smale complex The graph G (k) M is used to build the kth slice L (k) M S of the Morse-Smale complex. First, index k critical points σ (k) are identied: each monotonic section of the graph

G (k)
M which is minimal with respect to the partial order ≺ ih contributes its local minimum as an index k critical point. Second, the connections between each index k critical point σ (k) and its successors are set. Consider the anchor associated with the monotonic section dening σ (k) : we create one connexion between σ (k) and each index k -1 critical point of this anchor. Third, stable manifolds of index σ (k-1) are computed. More precisely, W s (σ (k) ) consists of the samples of B[k] together with the samples from linkers discovered during the current iteration, which ow to σ (k) . 3 Experiments

Benchmark

Rationale. To test our construction, we selected challenging functions of two and three variables whose Morse-Smale diagrams are known exactly. (Equations, critical points, and illustrations are presented in the supplemental section 7.) We devoted a particular attention to polynomial functions, whose critical points can be certied using real solving tools. More precisely, for a multi-variate polynomial P in d variables, one can form the system equating all partial derivatives to zero, and solve it if it is zero-dimensionalthe solution set consists of isolated points. This task is performed by the Maple function RootFinding[Isolate], which combines four ingredients, namely (i) a Gröbner basis calculation with the F4 algorithm [START_REF] Faugère | A new ecient algorithm for computing gröbner bases (f 4 )[END_REF], (ii) a Rational Univariate Representation (RUR) calculation [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF], (iii) a certied isolation of the roots of the polynomial produced by the RUR [START_REF] Rouillier | Ecient isolation of polynomial real roots[END_REF], and (iv) multi-precision interval arithmetic to report boxes isolating the solutions of the system [START_REF] Revol | Motivations for an arbitrary precision interval arithmetic and the mp library[END_REF]. Upon isolating the roots, one evaluates the Hessian matrix at each critical point, and gets the index from the number of negative eigenvalues.

In discussion these examples, the signature of a function is the number of critical points of each index.

Functions of two variables. The 2D Himmelblau function is a quartic polynomial dening a simple multimodal landscape, yet challenging for optimization algorithms since the four global minima are at the same height. Its signature is (4, 4, 1), and all critical points can be certied using the aforementioned real solving tools. The 2D Rastrigin function [START_REF] Mühlenbein | The parallel genetic algorithm as function optimizer[END_REF] is a challenging case due to the cosine modulation of a parabolic function. In the domain [-5.12, 5.12] 2 , its signature is (121, 220, 100). We also created two terrains (Gauss6A, Gauss6B), each dened as a sum of six Gaussians. In particular, the former illustrates the case where two index one saddles are connected to the same two local minima.

Functions of three variables. Since we did not nd any function with certied sets of critical points the literature, we used quartic polynomial functions whose critical points were certied using real solving. The associated MS complex was established by visual inspection, by tracking connections between sub-level sets.

To probe our construction beyond smooth functions, we used the Euclidean distance to a 3D point cloud, resulting in a Morse-Smale diagram known as the ow complex [START_REF] Giesen | The ow complex: A data structure for geometric modeling[END_REF]. For that case too, certied MS diagram can be obtained using multi-precision rational numbers [START_REF] Cazals | Robust construction of the threedimensional ow complex[END_REF].

Protocol

Implementation. We implemented our algorithm, whose pseudo-code is provided in the supplemental section 6.1, in generic C++. The NNG connecting the samples being taken for granted, each iteration requires sorting steps and graph traversal algorithms. That the worst-case complexities are linear in the size of the graphs manipulated (global thickening, local thickenings, monotonic sections and their incidence graph). The algorithm works in any dimension, and can be instantiated with a variety of points and distances (points in a Euclidean space under the Euclidean distance, points representing molecular conformations under the least root mean square deviation, etc).

Inria

Nearest Neighbor Graph and pseudo-gradient graph. To focus on intrinsic properties of the construction and avoid diculties related to an uneven sampling, we used point clouds sampled on a grid superimposed to the 2D or 3D domain of interest. Thanks to this grid sampling, the neighborhood of a sample was dened as its set of l neighbors, with l =

d-1 k=0 2 d-k d k .
To dened the PSG, we plainly connected each sample p i to its neighbor q j dening the steepest slope, amidst the l candidates.

Results

Timings. All the datasets processed, up to one million of points, required less than ve minutes on a desktop computer (one core Intel Xeon(R) CPU X5482 at 3.2GHz, running Linux Fedora core 17).

Smooth functions of two variables. Exact results (critical points and their connections) are obtained for all the functions under consideration, for samplings as small as 100 2 . See Fig. 3.

Smooth functions of three variables. We obtained correct results for samplings in the range 25 3 to 50 3 . One particular diculty we encountered is linked to the presence of open local thickenings (supplemental section 6.4), which rarely produced redundant critical points, ltered out as explained below. See Fig. 4.

Non smooth functions: the 3D ow complex. Correct results were also obtained, despite the non-smoothness of the distance function. At times, several critical points were also observed within monotonic sections, and were ltered out as explained below. A simple case is illustrated on Fig. 2.

Filtering out spurious critical points. By construction, each monotonic section is dened with respect to one local minimum and maximum of a local thickening (supplemental section 6.7). In rare cases, redundant local minima (and also maxima) located near one another were observed within a local thickening, and all but one were ltered out. To see how, consider two local minima p and q having a common ancestor r in the local thickening, with p nearest to r in function value (w.l.o.g.). In all the cases observed, the distance between p and r was at most three edges. To lter out p, we plainly revert the ow from p to r, which is the usual strategy to run persistence from the Morse-Smale diagram [START_REF] Cazals | Reconstructing 3d compact sets[END_REF].

Discussion and Outlook

This paper presents a framework to compute an approximation of the Morse-Smale diagram of a function, from a point cloud sampling the domain over which the function is dened. Central is the identication of bifurcating samples, i.e. samples facing ambiguities in terms of discrete ow dened on a nearest neighbor graph connecting the input points. From an algorithmic standpoint, our framework has three key advantages. First, it requires elementary data structures and operations, and is thus suitable for high-dimensional data processing. Second, in the context of non convex optimization, it is gradient free, which makes it suitable to investigate functions whose gradient is unknown or expensive to compute. Third, in case of under-sampling and even if the exact (unknown) MS diagram is not found, the output conveys information in terms of ambiguous ow, and the Morse theoretical version of topological persistence, which consists in canceling critical points by ow reversal, applies. From an applied standpoint, we present exact MS diagram for various challenging functions of two and three variables (including smooth functions and distance functions), a miles-stone rst achieved in this paper to the best of our knowledge. Despite these achievements, we foresee both theoretical and applied extensions.

On the theoretical side, two directions seem particularly promising. The rst one is related to the investigation of sucient conditions on the sampling, so as to ensure an isomorphism between the MS diagram computed and the exact one. Under such hypothesis and upon decorating the edges of the MS diagram with their (signed) multiplicity, our construction would allow computing the (persistent) homology of the underlying manifold using the Morse-Smale-Witten boundary operator a very attractive perspective since this boundary operator handles a CW complex with one cell per critical point. The second one deals with the extension of our framework to dynamical systems beyond gradient vector elds, so as to detect invariant sets such as attractors consisting of xed orbits. We believe in particular that some of the notions introduced in our work, related to the conuence of ow lines, will prove useful in this respect.

On the applied side, given the ubiquity of dynamical systems in science and engineering, our constructions should nd various applications, and we mention four of them. First, we plan to use it in the context of Hamiltonian systems, where phase transitions are typically investigated through the variation of the Euler characteristic χ of the phase space. Our approximate MS complex, extended to the MSW chain complex as just discussed, clearly contain more information. Second, our framework actually handles points sampling a space of arbitrary topology not just a contractible domain. We thus plan to use it to investigate sampled energy landscapes with non trivial topologies, as frequently encountered for molecular system that do not visit their entire conguration space (intuitively, orphan regions correspond to locii where the potential energy is prohibitively high). Third, we plan to work on the design of collective coordinates of macro-molecular systems, namely coordinates accounting for low frequency -large amplitude of molecular motions. The design of such coordinates is typically based on the identication of (signicant) stable minima of the energy landscape of the system, and our ability to report all transitions paths connecting two local minima across index one saddles is appealing to handle cases featuring saddles of comparable heights. Fourth, given the excellent performances observed on tri-variate functions sampled on a grid, we plan to use our tool in the context of 3D image analysis, and also to analyze general functions the latter in the context of scientic software such as Maple or Mathematica.

Artwork

Figure 1 Regions of homogeneous ow, bifurcations, and monotonic sections, illustrated for d = 2 and k = 1 (A) A region of homogeneous ow W (σ (2) 0 , σ

(0) 1 ) = W u (σ (2) 0 ) ∩ W s (σ (0)
1 ) on a two dimensional manifold. The compactication of this region also involves ow lines running through the intermediate critical points σ k) , are grouped into a local thickening T(A (k) ), and each local thickening is itself partitioned into monotonic regions M (α,ω) (T(A (k) )). Displayed here are the four monotonic regions associated with the anchor (σ

(0) 1 , σ (0)
2 ) involving the two local minima. Two monotonic regions associated to the same anchors (resp. two dierent anchors) are termed interior (resp. exterior) incident. (Middle) Sub-level set before the apparition of an index 2 critical point. The stable manifolds of σ (1) (pink points) and their boundary, the 2-global thickening(yellow) are depicted. (Bottom) The stable manifolds of σ (2) (orange points) are depicted. However, as the terrain does not have any maximum, these don't intersect. The six steps of the construction are recalled on Fig. 5, and the corresponding pseudo-code is presented in Algorithm 1. Complexity-wise, the following comments are in order for the execution of each iteration:

σ (0) 1 σ (1) 1 σ (1) 2 σ (2) 2 (A)(k = 1) σ (0) 1 σ (1) 2 σ (2) 1 σ (2) 2 σ (0) 2 σ (1) 1 σ (2) 3 (C)(k = 1) σ (0) 1 p k p j σ (2) 2 σ (2) 3 σ (1) 2 p i p l σ (1) 3 σ (1) 1 σ (2) 1 σ (0) 2 (B)(k = 1)
(0) 1 M (1) 2 M (1) 3 M (1) 1 σ (1) 1 σ (0) 2 σ (3) 1 σ (0) 4 σ (0) 3 σ (2) 1 σ (0) 2 σ (0) 1 σ (0) 4 σ (0) 3 σ (2) 1 σ (3) 1 σ (1) 1 RR n°8331
σ (1) 5 σ (2) 2 σ (1) 4 σ (2) 1 W s (σ (0) 2 ) W s (σ (0) 1 ) W s (σ (0) 4 ) W s (σ (0) 3 ) W s (σ (1) 1 ) W s (σ (1) 3 ) W s (σ (1) 5 ) W s (σ (1) 2 ) W s (σ (1) 4 ) W s (σ (2) 2 ) W s (σ (
The samples of T[k] by increasing function value.

Computing the alpha and omega limits of each sample within T[k] is linear in the size of the graph, in the worst-case.

Computing the anchor of each point requires the omega limit for all samples in its bifurcation lower link.

Computing local thickenings requires sorting samples using the lexicographic ordering on the tuples dening their anchors.

Identifying the local minima and maxima in the local thickenings has again linear complexity in the size of these local thickenings, in the worst-case.

Since incidences between monotonic sections are recorded while creating the monotonic sections, creating the incidence graph G

M is linear in its size.

All subsequent operations are also linear in the size of G

M .

RR n°8331

Figure 5 The six steps of our construction, together with the main denitions.

kth iteration, step 2:Computing anchors and local thickenings Main definitions:

-Hooks and anchors

A (k) -Local thickenings T(A (k) ) -Open local thickenings T(A (k) )
kth iteration, step 3:Computing monotonic sections and their incidence graph Main definitions:

- Monotonic sections M (α,ω) (T(A (k) )) -Linkers -Incidence graph of mono- tonic sections G (k) M
kth iteration, step 5: Computing the kth slice of the MS complex Main definitions:

-Index k critical points σ (k) -Stable manifolds of W s (σ (k-1) ) -Unstable manifolds of W u (σ (k-1) ) -Layer
Initialization (k = 0) -Local minima -Global thickening T[0] k → k + 1 Algorithm 1 Algorithm (0) Pre-processing
Compute nearest neighbor graph Compute pseudo-gradient graph Identify local minima Build the kth slice of the Morse-Smale complex:

(1) Compute the bifurcating and the conuence samples w.r.t. {σ (k-1) } in the global thickening (2) Group samples into equivalence classes, called local thickenings (3a) Decompose each local thickening into monotonic sections (3b) Construct graph of incidence between all monotonic sections, excepting linkers (4) Simplify the incidence graph between monotonic sections, based on the partial order ≺ ih (5a) Identify critical points σ (k) as local min of critical monotonic sections (5b) Compute (Un-)stable manifolds and set connections σ (k) → σ (k-1) (6) Update the global thickening for the next iteration Inria 6.2 Implementation tricks Monotonic sections and their extrema In the continuous case, the stable manifold of a critical point σ (k) can be decomposed into regions of homogeneous ow, each region corresponding to the intersection of provided that all points belong to the intersection one expects a monotonic regions, that is a non-empty region that represents W s (σ (k) ) ∪ W u (σ (k+1) ).

Details for step 1

Denitions. We assume that:

The critical points σ (k-1) of index k -1 have been identied.

The global thickening T[k -1] has been computed, with n k the cardinal of its vertex set.

The samples in T[k -1] have been sorted by increasing function value, so that p 1 and p n k respectively refer to the sample with minimum and maximum function values.

Due to the recursive nature of the sets constructed, we rst dene:

Denition. 1. The bifurcating samples B (k) i and the conuence samples

C (k) i ⊂ B (k) i for i = 1, . . . , n k , are subsets of vertices of the vertex set of the global thickening T[k -1], built incrementally by processing the n k vertices of T[k -1], as indicated below. A sample of T[k -1] not in B (k) n k is termed regular. The subgraph of T[k -1] induced by the bifurcating samples is called the global stream B[k].
Initialization. One has:

B (k) i
This set is initialized with the lowest sample p i such that two points from its lower link ow to dierent index k -1 critical points in the global thickening T[k -1]. Formally:

ω - |T[k-1] (p j ) = ω - |T[k-1] (p k ), with p j , p k ∈ L - |T[k-1] (p i ). (2)
Upon initialization of

B (k) i , one set C (k) i = ∅.
Iteration. Requires three steps:

Update of the conuence set:

C (k) i = C (k)
i-1 ∪ {p i } with p i a conuence sample characterized by: its entire lower link is a subset of

B (k) i-1 .
Denition of the bifurcation lower link L - |T[k-1] (p i ) of p i , by removing from the lower link all the bifurcating samples that are not conuence samples:

L - |T[k-1] (p i ) = L - |T[k-1] (p i )\(B (k) i-1 \C (k) i-1 ), if p i ∈ C (k) i-1 L - |T[k-1] (p j )
, with p j the lowest neighbour of p i , otherwise.

(3) Classication of p i as a bifurcating sample if:

ω - |T[k-1] (p j ) = ω - |T[k-1] (p k ) for at least two neighbors p j , p k ∈ L - |T[k-1] (p i ), (4)
and update

B (k) i = B (k) i-1 ∪ {p i } accordingly.

RR n°8331

Properties. The following properties are direct consequences of the iterative constructions:

Conuence samples are only observed for k ≥ 2.

The vertex set of T[k -1] is partitioned into bifurcating samples B }, is a set of critical points of index k -1. The anchor A (k) of a sample p i is the maximal hook satisfying:

∀σ (k-1) j ∈ A (k) , ∃p j ∈ L - |T[k-1] (p i ) ∪ {p i } such that ω - |T[k-1] (p j ) = σ (k-1) j . (5)
Anchors are used to group samples as follows:

Denition. 3. [Local thickening and its lower boundary] The local thickening T(A (k) ) of a k-anchor A (k) is the induced graph G |V , with V the union of vertices whose anchor is A (k) . The lower boundary of a local thickening ∂ l T(A (k) ) is the set of vertices p i , for which there exists an edge

(p i , p j ) ∈ T[k -1], with p j ∈ L - |T[k-1] (p i ).
Consider regions of homogeneous ow W (p, q i ), in the neighborhood of the critical point p. Due to the discrete nature of the sampling, we may observe edges in T[k -1] connecting regions associated with q i and q j such that W (q i , p) ∩ W (q j , p) = {p}. Such a connexion yields a local thickening, which is spurious since the region W (q i , p) ∩ W (q j , p) being reduced to a singleton, it cannot dene the stable manifold of any critical point. We identify these spurious local thickenings thanks to their thickness, and dene: for each sample p consider its neighborhood q i ∈ B The monotonic sections of a local thickening decompose it into sets of samples having the same α and ω limits for the ow in the local thickening. More formally:

Denition. 5. [Monotonic section] Given a local minimum α and a local maximum ω of a local thickening T(A (k) ), a monotonic section M (α,ω) (T(A (k) )) is the sub-graph of T(A (k) ) induced by the set of vertices such that ω - |T(A (k) ) (p i ) = α and ω +

|T(A (k) ) (p i ) = ω.
In the sequel, for the sake of conciseness, we shall simply write M i (A (k) ) for M (αi,ωj ) i,j (T(A (k) )). To handle the case of monotonic sections consisting of conuence samples, and whose local minima are not candidate for critical points, we dene: Denition. 6. [Linker] A monotonic section M j (A (k) ) is called a linkerif:

∀p i ∈ M j (A (k) ) \ ∂ l T(A (k) ), one has p i ∈ C (k) n . (6)
We nally arrive at the graph encoding incidences between all the monotonic sections of local thickenings at level k (Fig. 1(C)): Denition. 7. [Incidence graph of monotonic sections ] Two monotonic sections are said to be incident if in the stream B[k] the union of their vertices is connected. They are interior incident if they share the same anchor and exterior incident otherwise. The incidence graph of monotonic sections G (k) M is the graph dened as follows: its contains one node per monotonic section which is not a linker; two incident monotonic sections are connected by an edge.

Properties. The following properties are immediate:

The vertex set of the set of monotonic sections associated is a covering of the local thickening.

The vertex set of G (k)

M is a subset of the set of monotonic sections found at k -st iteration, provided that linker monotonic sections are not included.

The edge set of the union of monotonic sections is included in the edge set of the local thickening.

6.6

Details for step 4

Our goal is to edit

G (k) M into simplied graph G (k)
M , retaining selected edges only. The process relies on three partial orders ≺ i , ≺ h , ≺ ih , the latter mimicking the edges found in the MSW chain complex of a Morse function.

In sections 6.6.1 and 6.6.2, we respectively dene the partial orders and perform a cases study of all congurations.

Remark 1. Because the partial order ≺ i refers to the inclusion between anchors, in the sequel, we are only concerned with the comparisons of monotonic sections associated with dierent anchors.

H-non-I edge. Such an edge occurs in the context of composite edges detailed above, as the stick can also be incident to the M j (A (k) c1 ) and M j (A (k) c2 ) and the edge in will be labeled H-non-I. I-non-H edge. This edge occurs when a cycle (Fig. 9) is found in the sub-level associated with a critical value. One way to better represent this type of incidence , which aims at 'correcting' the ow induced at k+1 is the following : consider two adjacent nodes satisfying M i (A

(k) u ) ≺ i M i (A (k) v ) and M i (A (k) u ) h M i (A (k)
v ) The operation we wish to perform is to shift M i (A (k) v ) to some region where all critical points of index k are reachable. In order to do so, we relabel the samples found at the boundary of the two monotonic sections in question.

Non-hierarchical edge. A non-hierarchical edge occurs due to spurious connections in the T[k -1], which may cause erroneous paths when considering the ow operator at the next iteration (Fig. 8).

Note that unlike the forks and the sticks, in this case only the edge in G (k) M is removed, since discretization is responsible only for the incidence relationship, not for the presence of a spurious monotonic section. However, in the k + 1th iteration, the edges in T[k -1] accounting for a non-hierarchical edge must be dismissed. (1) 6 . In particular, the function value is less M

(1) 5 than on M

(1) 6 (using the partial order ≺ h : M

(1)

5 ≺ h M (1) 6 ; but M (1)
6 bifurcates with respect to two minima while M

(1) 5 bifurcates with respect to three minima (using the partial order ≺ i : M

(1) 6 ≺ i M (1) 5
). The two partial orders conict, so that the edge connecting M Properties. The following properties hold:

The vertex set associated to the union of such fork-stick couple will be a linker at k + 1 iteration (Fig. 7).

The samples that have incident edges in the T[k -1] that link two monotonic sections sharing a non-hierarchical edge in G Remark 2. [Linking σ (k) to σ (k-1) ] Each σ (k) will be linked to the σ (k-1) that form the anchor of its associated monotonic section.

Denition. 12. [Stable manifold of σ (k-1) ] The stable manifold of σ (k-1) is dened by:

W s (σ (k-1) ) = {p ∈ C(B[k]) ∪ M l (A (k) )|ω - |T[k-1] (p) = σ (k-1) }, (11)
where M l (A (1,0,0) #2 index 0 : (2,0,0) #3, index 0 : (3,0,0) #4, index 0 : (4,0,0) #5, index 1 : (3,0,0) #6, index 1 : (2,0,0) #7, index 1 : (1,0,0) #8, index 1 : (1,1,0) (creates a lled torus) #9, index 2 : (1,0,0) (lls the void in the middle of the lled torus) #10, index 1 : (1,1,0) (recreates a lled torus) #11 index 2 : (1,0,0) (rells the void of the lled torus)

  1: Identifying bifurcating and conuence samples This step aims at identifying one set B (k) i of bifurcating samples, and one set C (k) i of conuence samples, by processing the vertices of the global thickening T[k -1] by increasing function value (Fig. 1(A,B)):

M

  into a simplied graph G (k) M compatible with the MSW chain complex, by retaining edges (M (k) i , M (k) j ) between monotonic sections which potentially contribute critical points whose indices dier by one, and such that the value of the function on M (k) i less than the value of the function on M (k) j

Step 6 :

 6 Computing the global thickening for the next iteration The nal step consists in preparing the global thickening T[k] to be used during the next iteration. The vertex set of this graph consists of the vertices of the monotonic sections in G (k) M . RR n°8331 The graph T[k] is not an induced sub-graph of T[k -1], due to the removal of edges in T[k -1] responsible for the non-hierarchical incidences in G (k) M .

  circled blue region marks a neighborhood of one of the boundary ow lines, analyzed in B. The orange area marks a region where a conuence of two ow curves respectively descending from σ the smooth setting, these curves would only meet at the local minimum; in the discrete case, they coalesce before. (B) Points in the neighborhood of the boundary ow line of (A) may ow to σ ambiguity is used to dene the the bifurcating samples (red samples), which dene the set B (0) n ). Each such sample is witnessed by its bifurcation lower link L - |G (p i ) depicted in blue. Note that p l does not belong to the bifurcation lower link of p i since it bifurcates itself. (C) Samples bifurcating with respect to the same tuples of index k -1 critical points, called an anchor A (

Figure 2

 2 Figure 2 Concepts illustrated for the Euclidean distance function to the vertices of a tetrahedron (Main) Seven critical points σ (k) i are represented: index 0 in purple, index 1 in blue, index 2 in orange, index 3 in red. The three monotonic sections M (1) 1 , M (1) 2 , M (1) 3 partitioning the local thickening of the pair of local minima (σ () , σ () ) are also tagged. (Inset) Overview of the 7 critical points of the main gure.

  σ

Figure 3

 3 Figure 3 The critical points and manifolds of 2D functions (Top) The Himmelblau function has signature (4,4,1). All minima are at the same critical value.(Middle) The Rastrigin function has (121,220,100) as signature. (Bottom) The Gauss6A function has (3,5,3). One pair of minima shares three σ(1) . In the inset we can notice the two σ (1) (blue points) that are linked to the same anchor.

Figure 4

 4 Figure 4 Stable manifolds of the critical points of poly1 (Top) Stable manifolds of the four local minima up to the critical value of σ (1) 4 . Note that the W s (σ (0) i ) are separated by the global thickening T[1].

  of the Morse-Smale diagram L (k) HD kth iteration, step 4: Simplifying the incidence graph of monotonic sections Main definitions: -Partial orders ≺ ih -Reduced incidence graph G (k) M of monotonic sections kth iteration, step 6: Computing the global thickening -Global thickening T[k] kth iteration, step 1: Identifying bifurcating and confluence samples Main definitions: -Bifurcating samples -Confluence samples -Global stream B[k]

  stream B[k] into equivalence classes of samples facing identical bifurcations, we dene: Denition. 2. [Hook and Anchor] A hook of index k, also called a k-hook, denoted A (k) = {σ (k-1) j

  Denition. 4. [Open local thickening] An interior node of a local thickening is a sample p i for which all neighboring nodes that are also in B[k] share the same anchor as p i . A local thickening is open if has no interior node. Note that vertices in open local thickenings are reassigned as described below. Consider the set V [T(A (k) )] , the reassignment process ows as follows: iterate over samples in V [T(A (k) )] in increasing function value.

  k \T(A (k) ) Properties. The following properties holds: The vertices of local thickenings partition the vertex set of the steam B[k].

Figure 6

 6 Figure 6 Incidence between monotonic sections: composite edge. The split of ow lines after yields forks and sticks. (A) At k = 1, in the discrete setting, pairs and triples of local minima dene the monotonic sections M

Figure 7 Figure 8 Figure 9

 789 Figure 7 Composite edges and linkers. (A) In the smooth setting, four ow lines (in red) originating at a σ (1) converge to the index two critical point (red circle). Each black contoured region M (1) i

M

  would determine the presence of an open local thickening or a linker monotonic section at k + 1 iteration, as illustrated in 8(Left). This is what justies the removal of edges that link the two monotonic sections in question. The vertex set and the edge set of G (k) M are subsets of the vertex set, respectively the edge set of G . [Critical point] A critical point is the local minimum of a minimal monotonic section of the graph G (k) M .

  (k) ) are linkers andC(B[k]) is the complementary of V [T[k -1]].RR n°8331 7 Benchmark-→Frederic says: I suggest to organize the dataset as follows, since we may not include the Hasse diagram in the certif paper← function used in this study Himmelblau, see also http://en.wikipedia.org/wiki/Himmelblau%27s_function. f (x, y) = (x 2 + y -11) 2 + (x + y 2also also http://en.wikipedia.org/wiki/Rastrigin_function.f (x) = An + d i=1 x 2 i + A cos(2πx i ) , (x = (x 1 , . . . , x d ))(18)

Figure 10

 10 Figure 10 Bi-variate functions of the benchmark: illustrations (A) Himmelblau (picture from wikipedia) (B) 2D-Rastrigin (picture from wikipedia) (C) Gaussian-6A (D) Gaussian-6B

Figure 11 Figure 12

 1112 Figure 11 poly1(Eq. (19): overview of critical points and Morse-Smale diagram

  we follow the vectors of G -. Note in particular that a local minimum is characterized by ω - |G (p i ) = p i . The previous notions generalize for any induced sub-graph of G. Given a set of samples H ⊂ V [G], the induced descending pseudo-gradient graph G - |H is the pseudo-gradient graph dened on the sub-graph of G induced by the set V [H]. Note that this graph is dierent from the sub-graph of G -induced by the set V [H]. (If a node from H ows to a node of V [G]\H, it is orphan in the latter graph.) Likewise, the destination of a sample under the descending ow is denoted ω - |H (p i ), etc.

	2 Construction: Overview
	2.1	Ingredients of Morse homology
	Our framework to study a sampled Morse function is guided by the smooth setting, and by
	the so-called Morse homology theorem, which allows computing the homology of the underlying
	manifold from the incidences between pairs of critical points whose indices dier by one. The
	key ingredients of Morse homology are [BH04]:
		Morse theory. An index k critical point is denoted σ (k) , and its (un)stable manifolds are respectively W s (σ (k) ) and W u (σ (k) ). The index of a critical point p is denoted λ p . The sub-level
		set of a function f up to level h is denoted M [f ] h or M h for short if the function is clear.
		Conventions for artwork. For screen-shots, the colors associated with critical points are as
		follows: index 0 i.e. local minimum: purple; index 1 saddle: blue; index 2 saddle: orange; index
		3 i.e. local maximum: red.
		Inria

Table 2

 2 3D functions used in this study poly1f (x, y, z) = x 2 + y 2 + z -1 2 + x 2 + z 2 + y -2 2 + y 2 + z 2 + x -3 = x 2 + y 2 + z -5 2 + x 2 + z 2 + y -7 2 + y 2 + z 2 + x -11 , y, z) = x 2 + y 2 + z -5 2 + x 2 + z 2 + y -2 2 + 7 + x 2 + 0.2 x (21)

	2	(19)

RR n°8331

Acknowledgments. T. Dreyfus is acknowledged for his help with the software development. This work has partially been supported by the Computational Geometric Learning STREP project of the EC 7th Framework Programme (EC contract No. 255827).

Inria

The partial orders

The inclusion partial order between two k-anchors yields a partial between monotonic sections, denoted ≺ i : Denition. 8. [Partial order ≺ i ]

The height partial order, denoted ≺ h , aims at comparing monotonic sections linked by a (sequence of) monotonic section(s) within the graph

Denition. 9. [Partial order ≺ h ] The partial order is dened by:

where

Finally the partial order ≺ ih is dened by: Denition. 10. [Partial order ≺ ih ]

and

M

Given a partial order ≺ and two elements u and v, two situations are faced: either the two elements can be compared (u ≺ v or v ≺ u), or they cannot (u v). In dealing with two partial orders denoted ≺ i and ≺ h , this yields ve situations, namely: two if the two partial orders apply, two if only one of them applies, and one if none of them applies.

In the following, let us identify edges, say

v ) wrt to these two relations:

H-non-I : M

In the following, we discuss the situations where these types of edges occur.

Composite edge. This type of incidence id illustrated in Fig. 6. We term stick and fork, the vertices associated to the composite edge

Notice that the stick might be minimal wrt to ≺ i , but in terms of ow, the omega of the stick ows into the alpha of the fork in the graph induced by the vertex sets of the two monotonic sections. So the union of vertex sets of the stick and the fork form a monotonic region whose extrema are (α u , ω v ). Note that as far as the fork is concerned, as it is not minimal wrt to ≺ i it has at least two incident edges that are canonical, say M j (A

). In an ideal setting the 'premature' merge of the two c.c. associated to M j (A (k) c1 ) and M j (A (k) c2 ) shouldn't occur and one should redistribute the points in the support of the fork and stick to the incident monotonic sections found downstream. Inria Remark 3. When linkers (and implicitly conuence samples) are present, the stable manifold of a critical point is not only described by the points belonging to C(B (n k ) k ).

6.8

Details for step 6

The global thickening T[k] to be used during the next iteration is dened as follows:

Denition. 13. [Global thickening] The vertex set of T[k] is dened by:

The edge set of T[k] is dened by: