Didier Le Ray
email: leray@lnf.cnrs-mrs.fr

François Clarac

Daniel Cattaert

FUNCTIONAL ANALYSIS OF THE SENSORY MOTOR PATHWAY OF RESISTANCE REFLEX IN CRAYFISH. I-MULTI-SENSORY CODING AND MOTOR NEURON MONOSYNAPTIC RESPONSES

The in vitro preparation of the fifth thoracic ganglion of the crayfish was used to detail the negative feedback loop involved in the control of passive movements of the leg. Releasesensitive primary afferents issued from the coxo-basipodite chordotonal organ (CBCO), a proprioceptor which strand is released by upward movement of the leg, monosynaptically connect the depressor motor neurons (Dep MNs). Extracellular identification of sensory units from the CBCO neurogram allowed us to determine the global coding of a sine wave movement, imposed from the most released position of the CBCO strand. Intracellular recordings from sensory terminals (CBTs) and ramp movement stimulations applied to the CBCO strand allowed us to characterize two groups of release-sensitive CBCO fibers. The first group, divided into two subgroups (phasic and phasico-tonic), is characterized by discontinuous firing patterns: phasic CBTs fired exclusively during release movements; phasico-tonic CBTs displayed both a phasic firing and a tonic discharge during the more released plateaus. The second group was continuously firing whatever the movement, with higher frequencies during the release phase of the movement stimulation. All CBTs displayed a marked sensitivity for release movements while only the phasico-tonic ones showed a clear sensitivity to maintained positions. A previous study described different subtypes among the 12 Dep MNs of the hemi-ganglion: only nine Dep MNs (one assistant and eight resistant) are monosynaptically connected by CBCO afferents. Systematic intracellular recordings from all resistant Dep MNs, performed in high divalent cation saline, allowed us to describe two shapes of monosynaptic resistance reflex responses. A phasic response was characterized by bursts of EPSPs occurring exclusively during CBCO strand release movements. A phasicotonic response was characterized by a progressive depolarization occurring all along the release phase of the stimulation: during maintained released positions, the amplitude of the sustained depolarization was position-dependent; in addition, each release movement produced a phasic burst of EPSPs in the MN. The parallel study of the Dep MN properties failed to point out any correlation between the type of reflex response recorded from the MN and the MN intrinsic properties, which would indicate that the type of MN response is entirely determined by the afferent messages it receives. This question is addressed in a companion paper.

INTRODUCTION

In our knowledge on motor control regulation, the stretch reflex is a very widespread phenomenon described in all groups of vertebrates and invertebrates [START_REF] Granit | Receptors and sensory perception[END_REF][START_REF] Pringle | Proprioception in arthropods[END_REF]. Studied in details in the past, this reflex was first considered as rigid and stereotyped, but the concept of reflex has evolved recently in parallel with the data accumulated on the central pattern generator (Barnes and Gladden 1985;[START_REF] Clarac | How do sensory and motor signals interact during locomotion[END_REF][START_REF] Pearson | Common principles of motor control in vertebrates and invertebrates[END_REF]).

In the crustacean walking legs, each joint possesses one or two mechano-receptors, the chordotonal organs, able to record the different parameters of position and movements (Bush and Laverack 1982;[START_REF] Whitear | The fine structure of crustacean proprioceptor. I. The chordotonal organs in the legs of the shore crab, Carcinus maenas[END_REF][START_REF] Wiersma | Movement receptors in decapod crustacea[END_REF]. They are composed of an elastic strand in which ten to hundred bipolar cells are inserted. [START_REF] Mill | Chordotonal organs of Crustacean appendages[END_REF] defined four different classes of sensory afferents in the dactyl chordotonal organ. Two are sensitive to changes of the chordotonal organ strand length; they operate in opposite direction, one being stretchsensitive, the other release-sensitive. The two others code the strand length (stretched or released states of the strand). Considering the overall chordotonal discharge, these receptors are mainly activated during extreme angular movements and positions, and minimally activated for mid-ramp movements and positions. Bush (1962) demonstrated that chordotonal organs are able to induce a stretch reflex, called a resistance reflex, which consists in the activation of the passively stretched muscle. The stretch reflex is not limited to its proper joint but can spread out on other joints or even on other legs [START_REF] Clarac | Motor coordination in crustacean limbs[END_REF][START_REF] Vedel | Combined reflex action by several proprioceptive inputs in the rock lobster walking legs[END_REF].

Intensity and gain of this reflex appear to be very variable and may lead to complete reversal.

In the stick insect, [START_REF] Bässler | Reversal of a reflex to a single motoneuron in the stick insect Carausius morosus[END_REF], in the crab, Di Caprio and [START_REF] Di Caprio | Reversal of a walking leg reflex elicited by a muscle receptor[END_REF], and in the crayfish, [START_REF] Skorupski | Phase-dependent reversal of reflexes mediated by the thoracico-coxal muscle receptor organ in the crayfish, Pacifastacus leniusculus[END_REF] demonstrated that the same mechano-receptor is able depending upon "the central state" of the CNS to change a resistance reflex into an assistance reflex. Such variability has been similarly found in mammals [START_REF] Pearson | Common principles of motor control in vertebrates and invertebrates[END_REF] and in human [START_REF] Duysens | Interlimb reflexes during gait in cat and human[END_REF].

In the arthropods, intracellular studies have demonstrated that resistance reflex involves monosynaptic and polysynaptic connections [START_REF] Bässler | The femur tibia control system of stick insects. A model system for the study of the neural basis of joint control[END_REF][START_REF] Burrows | Monosynaptic connections between wing stretch receptors and flight motoneurones of the locust[END_REF]. The variability of the reflex has often been associated with interneurons which are inserted within the reflex pathway, between sensory terminals and MNs [START_REF] Burrows | Local circuits for the control of leg movements in an insect[END_REF][START_REF] Büschges | Non-spiking pathways antagonize the resistance reflex in the thoraco-coxal joint of stick insect[END_REF]Büschges and Wolf 1995;[START_REF] Le Ray | Neural Mechanisms of reflex reversal in coxo-basipodite depressor motor neurons of the crayfish[END_REF]. However in the crayfish, monosynaptic connections between the coxo-basal chordotonal organ (CBCO) and the levator (Lev) and depressor (Dep) MNs that command the coxo-basal joint, have been demonstrated to be responsible for efficient resistance reflexes [START_REF] El Manira | Monosynaptic connections mediate resistance reflex in crayfish (Procambarus clarkii) walking legs[END_REF].

The CBCO contains 40 afferents, half are sensitive to CBCO stretch (corresponding to downward movement of the leg in intact animal) and half to CBCO release (corresponding to upward movement of the leg in intact animal); the levator pool contains 19 MNs and the depressor pool 12 MNs. To study the determinants of the monosynaptic sensory-motor relationship, we have used a saline enriched in calcium and magnesium (x 2.5) in order to suppress the majority of the polysynaptic connections [START_REF] Berry | Criteria for distinguishing between monosynaptic and polysynaptic transmission[END_REF] and keep intact the direct CBCO afferent-motoneuronal pathways.

Up to now, sensory afferents have been considered depending upon the temporal parameters of their firing as phasic and tonic units, corresponding roughly to movement and position coding fibers. In the past, the classification of MNs was mainly established upon the types of muscle fiber they innervate: a MN is considered as phasic when its electrical stimulation elicits a twitch contraction, and as tonic when it induces a slow postural response.

In this paper we analyze in detail the coding of movement parameters by sensory afferents and the output response of the motor neurons (MNs). Using new computer data processing, it is possible to identify individually each sensory afferent in a sensory nerve discharge, based upon the shape of each extracellular action potential. To simplify this study we focused on the release-sensitive CBCO afferent fibers and the 9 Dep MNs that received monosynaptic inputs from the CBCO fibers (three of the 12 Dep MNs do not receive any direct CBCO inputs; Le Ray and Cattaert 1997). Therefore, we restrained mechanical movements imposed to the CBCO strand to its most released range, where this reflex is the most efficient. In this paper, we classified the different release-sensitive sensory afferents, and the different postsynaptic Dep MNs according to a qualitative and quantitative analysis of their responses and properties. In the following paper, using paired intracellular recordings, the wiring of the sensory-motor connections will be established and the role of specific synaptic parameters controlling MN response will be presented.

METHODS

Preparation.

Results are based on more than 130 intracellular recordings from CBTs and Dep MNs performed on adult male and female crayfish, Pacifastacus leniusculus and Procambarus clarkii. Animals were maintained in aquarium at 18°C and fed once a week.

The in vitro preparation (Fig. 1A) consisted of the last three thoracic ganglia along with the two couples of antagonistic motor nerves innervating the two proximal joints of the 5 th leg (Promotor/Remotor and Depressor/Levators). The CBCO, which encodes the vertical movements of the leg, was dissected out together with its sensory nerve (CB nerve). The preparation was pinned down dorsal side up in a Sylgard-lined Petri dish and superfused with oxygenated crayfish saline.

Stimulations / Recordings.

Extracellular recordings and nerve stimulations were performed using platinum electrodes contacting the nerves, isolated from bath with Vaseline, and directed to a 4-channel differential AC amplifier (A-M Systems). Single and paired intracellular recordings from CBTs and MNs (Fig. 1B) were realized with thin-walled glass micro-electrodes filled with a potassium chloride solution (3M) and having a 25-30 M resistance.

Salines.

A normal saline (NaCl: 195 mM, KCl: 5 mM, CaCl 2 : 13 mM, MgCl 2 : 2 mM), was used to perform the sensory coding analysis. Because strong central synaptic inputs masked the monosynaptic sensory inputs from the CBCO, the Dep MN responses were analyzed with saline in which divalent cation concentrations were increased (CaCl 2 : 34 mM, MgCl 2 : 6.4 mM, with the sodium concentration reduced accordingly); it was perfused exclusively on thoracic ganglia, in order to raise the activation threshold of all central neurons without affecting the CBCO sensory inputs (Fig. 8). Saline solutions were buffered with 3 mM N-[2-Hydroxyethyl] piperazine-N'-[2-ethanesulfonic acid] (Hepes) and pH adjusted at 7.7 at 15°C.

Analysis.

Signals were analyzed using the CED programs SPIKE2 and SIGAVG. The "Wave Marker" SPIKE2 program was used to classify the different profiles of extracellular spikes recorded from the CBCO nerve. Circular statistics were used to study the sensory coding of the cyclic imposed movement. Each sensory unit was identified over several movement cycles. Each cycle is defined from the starting of the release phase (chosen as zero) to the next release phase (chosen as 360 degrees). Each unit is figured by a unit vector  v i , with a phase angle according to its occurrence in the stimulus cycle. The resulting vector  R is the sum of all unit vectors, divided by the number (n) of spikes (Fig. 2A):

  R v n i n   ,  R R R x y   2 2 ,         ArcTan R R y x ,
with R x and R y representing the x and y coordinates of the resulting vector  R .

The resulting vector length (  R ) has a value between 0 and 1, and is a measure for the dispersion of the spikes within a stimulus cycle. A value of zero means that the spikes are evenly distributed; a value of one would indicate that all spikes occurred at the same phase.

The frequency distribution of  R over all experiments (Fig. 2B) is bimodal, the separation between the two peaks being 0.3. Therefore, units with  R value below 0.3 were considered as non-specific units (Fig. 2C), whereas the others were classified as release or stretchsensitive units.

SPIKE2 scripts were used to analyze, for each part of the mechanical stimulation, the instantaneous and mean firing frequencies from CBTs, and the Dep MN responses. SIGAVG program was used to record MN response to various intracellular current pulses, in order to realize the V/I curves from Dep MNs. Statistical analysis and regressions were performed with the GraphPad Prism statistic programs.

RESULTS

As shown in figure 1, C and D, the mechanical stimulation of the CBCO strand is able to produce a reflex activity in the in vitro preparation. This is characterized by the discharge of the sensory units that can be recorded extracellularly in the CBCO neurogram and intracellularly in the CBTs (Fig. 1C). This sensory firing induces, in the postsynaptic Dep MN, a depolarization of its membrane potential (Fig. 1D) because of the summation of each compound EPSP.

I -Analysis of the CBCO sensory coding of the leg movements.

Two complementary studies have been performed to describe the coding of the CBCO sensory informations. First, from the neurogram, an extracellular approach allowed to characterize in a general manner the different CBCO units. Because of the limitations of this method we realized a second analysis, the intracellular characterization of the CBCO fibers.

To facilitate our work, we focused the intracellular study on the release-sensitive afferents, impaled in their terminals.

I A -Extracellular study from the CBCO neurogram.

The extracellular study has been performed from a sinusoidal movement stimulation of the CBCO strand. In the angular range of the imposed sinusoidal movement, and depending on the preparation, about the half of the CBCO units were activated as shown on figure 2C and figure 3. In our conditions, release-sensitive afferents were predominantly recruited (11/20  2 release-sensitive and 7/20  3 stretch-sensitive CBCO fibers). Table of figure 2C gives the number of the sensory afferents activated by sinusoidal movements imposed to 10

CBCOs. A total number of 195 extracellular spikes has been studied to determine their sensitivity to CBCO strand movement. 55.5% of the recorded units were release-sensitive fibers (n = 108) against 28% of stretch-sensitive units (n = 55). In both cases, most of the CBCO units coded mainly one direction of the imposed movement, their discharge lasting during a part of the opposite movement. In all experiments, some CBCO fibers (16.5%, n = 32) were classified in the "Non-Specific" group because their response was not clearly related to one phase of the imposed movement (R  0.3; see Methods). Figure 3 presents a typical experiment in which different types of CBCO units were characterized from the CBCO neurogram, using the « Wave Marker » Spike2 program. The different identified extracellular spike templates are presented with the corresponding event histograms, which represent the occurrence of each extracellular spike in the stimulation cycle (bin size: 3 ms). In this experiment, 14 distinct release-coding and five stretch-sensitive CBCO units could be differentiated. Nevertheless, three sensory units did not code preferentially one direction of the imposed movement (R  0.3) All the release-sensitive fibers responded during the whole range of the release movement. Therefore, it appears that the CBCO coding is distributed, i.e.

sensory fibers do not code an extremely specific range of the movement. The phase () of the mean of each release-sensitive unit varied between 92.3 and 152.5 degrees, which indicates some differences in the coding of the various sensory units. Moreover, although some interunit differences exist (event histogram peaks occurred between 65 and 125 degrees), most of the units fired when the velocity of the imposed movement was maximum (90 degrees).

Thus, sinusoidal movement stimulation and extracellular studies of sensory neurograms allowed us to determine how many CBCO afferents were excited in the angular range of the stimulation applied to the CBCO strand. Ramp movement stimulation during intracellular recordings from identified CBTs were used to separate position and movement sensitivities of the sensory neurons.

I B -Classification of the release-coding sensory afferents -Intracellular study.

Two groups of release-sensitive fibers were defined on the basis of their response to ramp stimulations. Some were active only during one direction of movement; this group was subdivided into two subgroups (phasic and phasico-tonic), depending upon their position sensitivity. Others displayed a continuous although modulated discharge whatever the direction of the movement. In this study, we analyzed 19 phasic and 23 phasico-tonic, from the first group, and 22 continuous firing fibers from the second group.

Phasic CBCO units.

CBCO units that responded exclusively during the release movements of the CBCO strand (see example in Fig. 4A), and were silent during plateaus (see detail in insert), were classified as phasic units. Instantaneous frequency of this kind of sensory terminals varied from 5  3 to 75  8 Hz during release movements, and mean frequency varied from 10  5 to 50  10 Hz. As shown by the mean frequency histogram of figure 4A example (bottom), there was no clear relationship between the mean firing frequency and the angular sector in which movement was applied. Moreover, the pattern of discharge was very variable from one cycle to the next. This explains the apparent discrepancy between the raw data (Fig. 4A, top) and the mean frequency (histogram in Fig. 4A). When low velocity (0.25 mm.s -1 ) ramp stimulation was applied (Fig. 4B), the mean frequency of every phasic CBTs was dramatically reduced to about one third of that obtained with high velocity (1.25 mm.s -1 ) ramp stimulation.

Phasico-tonic CBCO units.

Phasico-tonic fibers (Fig. 5) were characterized by both a phasic response during release movements and a response to position. In general, the tonic response developed from a particular release position (different for each unit) and increased with the degree of release up to the maximum release position. As shown on figure 5, A and B, phasico-tonic types of CBCO units can be separated into two groups depending upon the hysteresis of position coding: The first one (CBT 1, Fig. 5A) was characterized by spike emission for the most released positions, whatever the direction of the movement applied between maintained positions. However for a given position, the discharge was larger after a release movement than after a stretch movement. The second group (CBT 2, Fig. 5B) stopped firing during the stretching phase; even for most release position, a stretch movement induced an inhibition of tonic firing that persisted during the two seconds of maintained position. Nevertheless, firing patterns of both types (CBT 1 and 2) were similar during the CBCO release (see detail in insert): the high frequency phasic firing was followed by a lower frequency firing during the plateaus. Histograms of both mean frequencies, and linear regressions, showed that phasic (white histogram) and tonic (dark histogram) firings were increasing with the release of the CBCO strand. It was noticeable that in both cases the phasico-tonic units fired with a higher frequency during movements than during plateaus (up to 48  7 Hz during ramps against up to 12  5 Hz during plateaus). The velocity sensitivity of the phasico-tonic CBTs was exclusively limited to the phasic part of their response as shown by frequency histograms on figure 5C. As is the case for phasic CBCO units, when the velocity of the ramp was increased from 0.25 to 1.25 mm.s -1 (i.e. multiplied by 5) the firing frequency of the movement-sensitive component of the phasico-tonic fiber was multiplied by 3 (white histogram). By contrast, the firing frequency of the tonic component appeared to be non significantly affected by the change in velocity (dark histogram).

As figure 6 shows, the phasico-tonic CBCO units presented an adaptation of the tonic firing when the mechanical stimulation was stopped in a released position for a long time (> 10s, Fig. 6A). At the end of the last ramp, the tonic firing began with a frequency of 18  3

Hz that was reduced to a basal value of 4  2 Hz within 5 seconds. Figure 6B shows the mean firing frequency calculated over 5 s intervals (histogram) and instantaneous frequency (filled circles) plotted as a function of time, after the mechanical stimulation was stopped in the most released position (time 0), for one phasico-tonic fiber. The basal instantaneous frequency value was obtained only after six seconds in released position.

Continuous discharge CBCO units.

Never described in detail in previous reports, this kind of fibers was characterized by a continuous firing (at least 10  4 Hz) that was modulated by movements imposed to the CBCO strand. Nevertheless, release sensitive units presented a specificity for the release phase of the mechanical stimulation (where firing frequency could raise to a mean value of 65

 15 Hz during ramps). As shown on figure 7A, continuous discharge fibers fired with a higher frequency during the CBCO release but were still spiking during the stretch of the sensory organ. In this type of units, a movement sensitivity was also evident. The mean frequency histograms (Fig. 7B) show a marked difference between the movement response during the releasing ramps (Fig. 9A). These bursts of EPSPs were almost constant from the first to the last releasing movement of a stimulation cycle. This kind of response, where EPSP bursts stopped since the beginning of the plateaus (see detail in insert of Fig. 9A), was defined as a « phasic response ». The second type of reflex responses (Fig. 9B) was characterized by the same phasic bursts of EPSPs, which summated on a general depolarization of the MN that occurred during the whole release of the CBCO. In this case, the amplitude of the EPSP burst weakly increased from the first to the last releasing ramp (1  0.4 mV during the first ramp, and 3  0.3 mV during the last). In parallel, the Dep MN depolarization was increasing from the beginning until the end of the CBCO release (from 0 to 1.7  1 mV). This kind of response (see detail in Fig. 9B insert) defined a « phasico-tonic » Dep MN response. As shown in figure 9C, although the phasic bursts of EPSPs swiftly disappeared, the tonic depolarization of the phasico-tonic Dep MNs persisted (with an almost constant value) when the stimulation was stopped in a released position. The Dep MN completely repolarized only when the CBCO stretch started.

Tonic / phasic properties of MN and monosynaptic response to movement.

Since it was often claimed that phasic MN had larger diameter and higher conduction velocity, we used these parameters to classify the Dep MNs. Figure 10A In order to classify MNs on the basis of their firing properties, we used injection of long lasting depolarizing current pulses (2s-duration, 0.3 Hz). Some MNs fired continuously all along the depolarizing pulse, with a discharge frequency depending upon the amount of injected current; they are defined as tonic MNs. Other MNs fired only few spikes at the beginning of the current pulse and rapidly decreased their firing frequency. Their discharge stopped before the end of the pulse; they are defined as phasic MNs. It appeared that the firing pattern of the Dep MN and its reflex response were not correlated since both patterns could be recorded in Dep MNs displaying different monosynaptic reflex response (data not shown).

Active membrane properties of MN and monosynaptic response to movement.

The Dep MN V/I curves were performed in a 5 mV range around the membrane potential in order to estimate the role of the membrane properties in the reflex response, since CBT-related EPSPs generally may reach up to 5 mV in MNs. V/I curves, issued from injection of depolarizing and hyperpolarizing current pulses, revealed four types of electrical behavior in high divalent cation saline (Fig. 10B). Some Dep MNs (35%) displayed almost a linear response to current injection. The majority (51%) presented a negative rectification for negative current pulses (V1) and a weak positive rectification for positive current pulses

DISCUSSION

In this paper, we have characterized the properties of both elements involved in the "stretch reflex" in crustacean walking legs: the CBCO sensory afferents and the connected MNs. We will discuss successively the functional aspects of the present results, the types of sensory coding and MN responses.

The crayfish thoracic ganglion in vitro preparation.

The crayfish thoracic walking system is dissected out in such a way to eliminate most of the sensory inputs that can excite the MNs. Only the CBCO is kept intact in the Petri dish with its sensory nerve that innervates the fifth thoracic ganglion, only leaving functional the monosynaptic reflex pathway. Perfusing the high Ca 2+ and high Mg 2+ saline upon the ganglion, by rising the threshold for spiking, we got sure to strongly limit the polysynaptic pathways which allowed us to study much clearly the monosynaptic connections between the CBCO sensory afferents and the MNs. Consequently high Ca 2+high Mg 2+ saline emphasized the monosynaptic MN responses (Fig. 8). In all experiments, we were very careful to have very stable conditions and reproducible responses before performing the recordings (about two hours after the onset of the dissection). The CBCO was dissected out and disposed in the Petri dish in such a position that the puller exactly imposed movements to the strand in the same way as if it was in the intact leg. The CBCO strand length was measured in the fully released position, in order to apply stretches that were one third of the original length. All stimulations were applied in this angular sector that corresponds to the angular range covered during the locomotor movements of the leg: The CB joint allows movements within the angular sector from the maximum of levation when the leg contacts the thorax, to the maximum of depression when the leg contacts the ventral carapace. Accidentally the CB joint can reach these extreme positions, for example when the animal is upside down. Nevertheless, movements involved in locomotion or in posture are restricted to about 60 degrees, between the horizontal plane and the complete levation. During levation the CBCO strand is completely relaxed, and during depression it is slightly stretched; we focused our study on this functional part of CB joint movements.

Sensory unit coding.

By the past, the crustacean walking leg chordotonal organs have been classified in two types: the first type present a symmetric discharge, as the pro-dactylopodite (PD), the carpopropodite (CP 1 ), the mero-carpopodite (MC 1 ) and the CB chordotonal organs, i.e. they have fairly the same number of stretch-sensitive and release-sensitive sensory units; the second type only possess release-sensitive sensory fibers, as the CP 2 , the MC 2 and the basi-ischiopodite chordotonal organs [START_REF] Mill | Chordotonal organs of Crustacean appendages[END_REF][START_REF] Mill | The fine structure of the PD proprioceptor of Cancer pagarus. I. The receptor stand and the movement sensitive cells[END_REF]. Other chordotonal organs are not disposed at a joint: the cuticular stress detectors (CSD 1 and CSD 2 ), located near the CB joint, are linked to a small window of soft cuticle and respond to mechanical stretching and releasing of their strand when the cuticle is stimulated [START_REF] Marchand | Physiology and in vivo activity of cuticular stress detector afferents in crayfish[END_REF]. In 1991, El Manira et al. described the CBCO: it is composed of 40 cells, half of them are sensitive to the stretch of the strand and the others are sensitive to its release. In this study, the analysis of the sensory coding was performed using two types of mechanical stimulations imposed to the CBCO strand.

The sinusoidal movements were used because it resembles the real movement performed by the CB joint during walking activity. Ramp movements were used to finely distinguish between position and velocity-sensitive components of the responses. However, ramp movements did not allow analysis from the neurogram because of the large number of units synchronized during each ramp. Therefore, ramp stimulations were exclusively performed during intracellular recordings from CBTs. By contrast, sine wave movements allowed the analysis of the neurogram because no synchronization occurred. From one experiment to the other, we were very careful to always perform ramps which had the same speed and amplitude characteristics. However, because the released CBCO length varied from 1 to 1.8 mm, the number of ramps imposed to the strand was modified accordingly (Figs. 4 to 7). This procedure results in ramp amplitudes being different relatively to the released CBCO length. However, the results were very comparable whatever the imposed ramp amplitudes.

Our results confirm the previous observations that described the unidirectional responses of the sensory fibers (Bush 1965a, b;[START_REF] Mendelson | Some factors in the activation of crab movement receptors[END_REF][START_REF] Mill | Chordotonal organs of Crustacean appendages[END_REF]) and, from more sophisticated tools, give new insights concerning the chordotonal organ fine sensory coding.

By the past, most of the studies were performed from extracellular recordings (Bush 1965a, b;Clarac 1970), the analysis being based only upon the amplitude of the action potentials, and remained sometimes quite doubtful concerning the discrimination of very close amplitude units. Our extracellular analysis, based upon the extracellular action potential recognition, thoroughly studied the discharge of almost all the sensory units excited during movements imposed in a physiological range. Sensory terminal intracellular recordings allowed us to identify specifically different types of sensory units. They can be divided into three subtypes which are the phasic, the phasico-tonic and the continuously firing sensory fibers. This classification is mainly based upon the position sensitivity because all of them are strongly sensitive to movement. The phasic units only responded to movement, their response being unchanged whatever the angular sector where it was measured. By contrast, in both phasicotonic and continuously firing units, the response to movement was considerably affected by the angular sector where it was measured. Therefore, among movement-sensitive responses (dynamic responses) we can distinguish the position-unaffected dynamic (phasic units) from the position-modulated dynamic responses (phasico-tonic and continuously firing units). In addition, the dynamic response may be velocity-sensitive (phasic and phasico-tonic units) or not (continuously firing units). Surprisingly, we did not find purely position-coding units as was usually described (Bush 1965a;Clarac 1970;[START_REF] Matheson | Responses and locations of neurones in the locust metathoracic femoral organ[END_REF]), but all units coding positions were also activated during movements (phasico-tonic CBTs). Whatever the imposed movement (ramp or sine wave), it appeared that each sensory fiber is activated in a wide angular sector. We never recorded any fiber coding specifically for a small joint angular sector.

As a conclusion of this study, it seems that the coding of movement is shared by all the sensory fibers, i.e. on the whole angular sector of the CB joint, upward movements should activate about 20 release-sensitive fibers and downward movements should activate about 20 stretch-sensitive fibers. Moreover, during movements, the firing frequency of any CBT is always three to ten times higher than during maintained positions. Consequently, the major information conveyed by this sensory system is dynamic, and whatever the postsynaptic neurons they should be more excited during movements than during maintained positions. The CBCO movement-sensitivity may therefore be responsible for the resistance reflex being a dynamic reflex, activated by any velocity variation rather than by static positions. Concerning the coding of position, CBTs seemed to present a low selectivity because during sine wave movements all release-sensitive sensory afferents were activated on the whole range of the releasing movement. However, their maximum of discharge occurred at angular positions that were different from one CBT to the other (Fig. 3).

In the angular sector where movements were imposed to the CBCO strand, the releasesensitive sensory units were activated in a greater number than stretch-sensitive ones (almost twice the number of stretch-sensitive units). Then, in the angular sector involved in locomotion and posture, the monosynaptic pathway excites much more the Dep MNs than the Lev ones. Thus the resulting monosynaptic reflex effect is the body upholding.

Monosynaptic motoneuronal responses.

Recently, different studies upon the crustacean MNs have demonstrated that, in the thorax, the abdomen or the stomatogastric ganglia, motoneuronal pools are quite complex and not homogenous. Their intrinsic properties and their connections have been described to be responsible for this heterogeneity [START_REF] Harris-Warrick | Dynamic biological networks: The stomatogastric nervous system[END_REF][START_REF] Le Ray | Neural Mechanisms of reflex reversal in coxo-basipodite depressor motor neurons of the crayfish[END_REF][START_REF] Skorupski | Heterogeneity and central modulation of feedback reflexes in crayfish motor pool[END_REF][START_REF] Wine | The structural basis of an innate behavioural pattern[END_REF]). In the same MN pool, it seems that not all the MNs have the same role: [START_REF] Skorupski | Heterogeneity and central modulation of feedback reflexes in crayfish motor pool[END_REF] demonstrated that two groups must be considered in the promotor MN pool, one involved in the resistance reflex and the second in the assistant pathway. More recently, Le [START_REF] Le Ray | Neural Mechanisms of reflex reversal in coxo-basipodite depressor motor neurons of the crayfish[END_REF] showed in the Dep MN pool the existence of one assistance-specialized Dep MN. In this study, we classified the Dep MNs according to their monosynaptic reflex response and their intrinsic properties. We did not take into account the electrotonic coupling that have been described to be very weak between Dep MNs (Chrachri and Clarac 1989). Classically MNs have been classified into tonic and phasic depending upon the type of muscle fibers they innervate and their activity (Kennedy and Takeda 1965a, b). Tonic MNs display continuous discharge whereas phasic MNs are activated only during rapid movements [START_REF] Kennedy | The organization of invertebrate nervous systems[END_REF]. Generally, phasic MNs have cell bodies and axon diameters larger than tonic ones, and they convey spikes more rapidly.

Functionally the resistance reflex is acting during posture and should involve preferentially postural MNs, i.e. tonic MNs. In the present work, Dep MN classification was not related neither to the muscle fibers innervated, nor to their activity in motor behavior. However, we attempted to classify Dep MNs according to their conduction velocity (Fig. 10A) and their firing pattern. Surprisingly this study demonstrates that this classification did not match the different types of reflex responses because resistance reflex Dep MNs were found in the whole range of the conduction delays (Fig. 10A). However this lack of correlation could be explained by the conduction delays not being representative of tonic or phasic behavior of a given MN. Concerning the three MNs which are not directly connected by the CBCO sensory afferents, they were also found in the whole range of Dep MN conduction delays (Fig. 10A).

Nevertheless, this group may also be heterogeneous, some of them being excited through one or several interneurons and other not. Moreover the properties of interneurons involved in polysynaptic pathways may be very different from one MN to the other. Such an organization has already been demonstrated in the case of the insect locomotor system [START_REF] Burrows | Positive feedback loops from proprioceptors involved in leg movements of the locust[END_REF][START_REF] Sauer | Information processing in the femur-tibia control loop of stick insects. 1. The response characteristics of two nonspiking interneurons result from parallel excitatory and inhibitory inputs[END_REF] and in the case of the crayfish reflex pathways [START_REF] Clarac | Interneuronal control of a walking leg reflex in an in vitro crayfish preparation[END_REF], and might be very important in the reflex modulation.

In A: Procedure used to analyze the different sensory units extracellularly recorded from the CBCO nerve. In this illustration, three units were identified by their shape and amplitude. The activity of each unit was then represented using circular statistics. Within the stimulation cycle, each occurrence of the unit (unit 1 in the example) was represented by an unit vector

(  v i ).
The vectorial sum of these unit vectors was calculated over several cycles and divided by the number n of occurrences; the resulting vector (  R ) indicated the coding characteristics (angle  and vector length R) of the unit (see text for details). B: Bimodal distribution of resulting vector lengths, calculated from all analyzed CBCO units. This distribution was used to distinguish between specific (R > 0.3) and non-specific (R  0.3) CBCO units. C: Using the calculated  and R values, CBCO units were classified in 10 different experiments: in the angular sector we studied, most of the release-sensitive fibers were activated. Stretch-sensitive units were generally, and non-specific units were always activated in a lower ratio.

  The signals were amplified by an Axoclamp 2B (Axon Instruments). Intracellular current pulses and nerve stimulations were controlled by an 8-channel digital stimulator (A.M.P.I.). All signals were monitored on a eight channel oscilloscope, a four channel digital oscilloscope (Yokogawa DL 1200) and stored on D.A.T. tapes (BioLogics digital tape recorder) and digitized on a PCbased computer through an A/D interface (from Cambridge Electronic Device, CED 1401PLUS). Intracellular and extracellular recordings were digitized at 5-10 kHz and written to disk. An electromechanical puller (Ling Dynamic systems Ltd.) commanded by a homemade controller was used to stretch and release of the CBCO strand, according to sinusoïdal (Fig. 1, A, C and D) or ramp (Fig. 4) protocols. Ramp stimulations are very suitable for separating position from movement sensitivities in intracellular recordings from one CBT. However, this procedure tends to synchronize the firing at the onset of each ramp movement which makes it very difficult to separate the different units within the neurogram. Therefore sine wave movements were used during experiments in which CBCO neurograms were analyzed. Stretching stimulations were performed from the most released position of the CBCO strand and total movement amplitude was one third of the released CBCO strand length (1 mm -1.8 mm). It almost corresponds to the angular range in which the leg moves during locomotion. Imposed movements elicit stretch and release of the CBCO strand, as real leg movements do. The movement control voltage traces were visualized on the oscilloscopes and stored on both tape and computer.

(

  figure 8 (left), in a normal saline the Dep MN monosynaptic reflex response induced by

  presents the conduction delays of 59 Dep MNs recorded from 12 experiments. The different MNs are grouped according to the monosynaptic reflex response they displayed, i.e. assistance, phasic resistance, phasico-tonic resistance or no response. In the four cases, no correlation appeared between conduction delays and the kind of monosynaptic response. The dispersion observed in the delays in each group is quite large: from 2.35 to 10 ms for the no responding Dep MNs (circles); from 2.5 to 8.75 ms for the phasico-tonic resistance MNs (filled squares); from 2.19 to 8 ms for the phasic resistance MNs (white squares); and from 3.44 to 6.88 ms for the assistance Dep MN (stars). This dispersion of conduction delay values for the unique assistance Dep MN would indicate some variability from one animal to the other. Nevertheless, in a single experiment where all Dep MNs were recorded, the distribution of conduction delays was similar and did not point out any relation with the reflex response displayed by the different Dep MNs.

(

  V2). A little part of Dep MNs (8.5%) showed a double negative rectification (V1 and V2). The last type of electrical behavior, a positive rectification V1 and a negative rectification V2, was detected in few cases (only 5.5% of the impaled Dep MNs). V1 and V2 values were calculated as the differences between the linear regression (around 0) and the real values for -5 and +5 nA injected current amplitudes.In order to represent the rectification properties of all recorded Dep MNs in a single diagram, we plotted V2 as a function of V1 for each MN (Fig.10C). This diagram revealed clusters of Dep MNs that would indicate distinct behavioral groups of MNs. However, each cluster contained phasic and phasico-tonic (filled squares) resistance Dep MNs, except for cluster 5 (positive V1 and negative V2), that contained exclusively the assistance Dep MN (stars). Non responding Dep MNs (open circles) were also found in clusters 1 and 2 (moderate negative V1), but not in clusters 3 and 4 (larger negative V1).

Fig. 1 :

 1 Fig. 1: Diagram of the experimental arrangement used in the analysis of sensory-motor

Fig. 2 :

 2 Fig. 2: Extracellular study -Methodology.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Characterization of the CBCO units recruited by sinusoidal stimulation.

Fig. 5 :Fig. 8 :

 58 Fig. 5: The phasico-tonic release-sensitive CBTs.

Fig. 9 :

 9 Fig. 9: Monosynaptic resistance Dep MN responses.

Fig. 10 :

 10 Fig. 10: Dep MN properties and reflex responses.

Only the assistance responding Dep MN presented particular active electrical properties (cluster 5).


Hang Zhou

Rémi Douvenot

Alexandre Chabory

Modeling the Long-Range Wave Propagation by a Split-Step Wavelet Method

A split-step wavelet method for simulating the long-range wave propagation is introduced. It is based on the fast wavelet transform. Compared to the split-step Fourier method, this method improves the computation efficiency while keeping a good accuracy. The propagation is performed iteratively by means of a pre-computed matrix containing the individual propagations of the wavelets. A fast computation method of this matrix is also presented. For the radiowave propagation in the low troposphere, a local image method is proposed to account for an impedance ground. Inhomogeneous atmospheres and irregular grounds are also considered. Finally, numerical tests of long-range propagations are performed to show the accuracy and time efficiency of this method.

Introduction

Modeling of the long-range wave propagation is a major issue for numerous applications. The parabolic equation method (PE), introduced by Feit and Fock [START_REF] Feit | Light propagation in graded-index optical fibers[END_REF], is considered as one of the most reliable numerical method. It is based on an approximation of the wave equation valid along a paraxial direction and neglecting backward propagation. This method has been widely studied to model the propagation of acoustic waves under the sea, electromagnetic waves in the low atmosphere, or to solve the Schrödinger equation, e.g., [START_REF] St Mary | A modified wide angle parabolic wave equation[END_REF][START_REF] Kuttler | Theoretical description of the parabolic approximation/Fourier split-step method of representing electromagnetic propagation in the troposphere[END_REF][START_REF] Lee | Parabolic equation development in the twentieth century[END_REF][START_REF] Taha | Analytical and numerical aspects of certain nonlinear evolution equations. II. Numerical, nonlinear Schrödinger equation[END_REF][START_REF] Muslu | Higher-order split-step Fourier schemes for the generalized nonlinear Schrödinger equation[END_REF].

The two mostly-used methods for applying PE are based on either finite-differences (FD) [START_REF] Levy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF] or the split-step Fourier method (SSF) [START_REF] Hardin | Applications of the split-step Fourier method to the numerical solution of nonlinear and variable coefficient wave equations[END_REF]. The FD has the advantage of a straightforward implementation of complex boundary conditions. SSF is more numerically efficient because it permits larger grid increments. Consequently, SSF is widely used for long-range wave propagation. The computation is done marching on in distances. At each step, the wave is transformed from the spatial domain to the spectral domain by means of a fast Fourier transform (FFT) along the vertical coordinate.

For electromagnetic propagation in the low troposphere, refractivity, ground boundary condition, and irregular relief can be considered in SSF. The phase-screens method is applied to take into account refractivity in the spatial domain [START_REF] Levy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF]. The discrete mixed Fourier transform [START_REF] Dockery | An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation[END_REF] includes the reflection over an impedance ground in the spectral transform. Finally, the relief can be considered by different algorithms [START_REF] Levy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF][START_REF] Barrios | A terrain parabolic equation model for propagation in the troposphere[END_REF][START_REF] Donohue | Propagation modeling over terrain using the parabolic wave equation[END_REF][START_REF] Janaswamy | A curvilinear coordinate-based split-step parabolic equation method for propagation predictions over terrain[END_REF].

Propagation methods in 3D have also been developed with FD and SSF for underwater acoustics [START_REF] Siegmann | A wide-angle three-dimensional parabolic wave equation[END_REF][START_REF] Sturm | On the use of higher-order azimuthal schemes in 3-D PE modeling[END_REF][START_REF] Sturm | Numerical study of broadband sound pulse propagation in three-dimensional oceanic waveguides[END_REF][START_REF] Belonosov | An iterative solver for the 3D Helmholtz equation[END_REF]. These 3D methods model the lateral couplings. For electromagnetic propagation, 3D extensions of SSF have also been developed [START_REF] Zaporozhets | Bistatic RCS calculations with the vector parabolic equation method[END_REF][START_REF] Janaswamy | Path loss predictions in the presence of buildings on flat terrain: A 3-D vector parabolic equation approach[END_REF][START_REF] Ginestet | Modélisation de la propagation d'une onde electromagnétique sur des scènes de grande taille par résolution de l'equation parabolique 3D vectorielle[END_REF]. The 3D method has a high accuracy and the capability to model lateral effects. However, if a large computation domain is considered, the computation burden in time and memory becomes the main concern.

Lately, Zhou et al. have developed the discrete SSF method (DSSF) in both 2D [START_REF] Zhou | Comparisons of discrete and continuous propagators for the modelling of low tropospheric propagation[END_REF][START_REF] Zhou | Modeling the atmospheric propagation of electromagnetic waves in 2D and 3D using Fourier and wavelet transforms[END_REF] and 3D [START_REF] Zhou | A 3-D split-step Fourier algorithm based on a discrete spectral representation of the propagation equation[END_REF][START_REF] Zhou | Modeling the atmospheric propagation of electromagnetic waves in 2D and 3D using Fourier and wavelet transforms[END_REF]. This improvement of SSF brings self-consistency as defined by [START_REF] Chew | Electromagnetic theory on a lattrice[END_REF][START_REF] Teixeira | Lattice electromagnetic theory from a topological viewpoint[END_REF]. Therefore numerical spurious solutions due to the a-posteriori discretization are avoided. The 2D DSSF is used as the reference in this paper.

Gabor bases/frames [START_REF] Einziger | Gabor representation and aperture theory[END_REF][START_REF] Lugara | Frame-based Gaussian beam summation method: Theory and applications[END_REF] decompositions have also been proposed for propagation modeling. They have been combined with analytic formulations of the beam propagation. Gabor-based beam algorithms and frame-based beam summation methods are utilized in various applications involving radiation and scattering in complex environments [START_REF] Lugara | Frame-based Gaussian beam summation method: Theory and applications[END_REF][START_REF] Maciel | Systematic study of fields due to extended apertures by Gaussian beam discretization[END_REF][START_REF] Maciel | Discretized Gabor-based beam algorithm for time-harmonic radiation from twodimensional truncated planar aperture distributions .I. Formulation and solution[END_REF][START_REF] Chabory | Novel Gabor-based Gaussian beam expansion for curved aperture radiation in dimension two[END_REF].

Wavelets are widely used for signal processing, numerical analysis and data compression [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. As a Gabor atom, a wavelet is a short length oscillation function that has a space-frequency localization property. In the past twenty years, several applications of wavelets in the wave propagation domain have been studied. Wavelets applied to integral equations [START_REF] Steinberg | On the use of wavelet expansions in the method of moments (EM scattering)[END_REF][START_REF] Sarkar | Wavelet Application in Engineering Electromagnetics[END_REF] and time-domain methods [START_REF] Fröhlich | An adaptive wavelet-vaguelette algorithm for the solution of PDEs[END_REF][START_REF] Holmström | Solving hyperbolic PDEs using interpolating wavelets[END_REF][START_REF] Hong | On a wavelet-based method for the numerical simulation of wave propagation[END_REF][START_REF] Sarkar | Wavelet Application in Engineering Electromagnetics[END_REF] have been proposed to improve their computational efficiency. In integral equations in electromagnetics, wavelets are notably used to transform a full-matrix linear system into a sparse one. This latter can be solved efficiently by means of Krylov iterative techniques. However, for simulating the long-range propagation, large computer resources and computation time are still required. Efforts to used wavelets for solving PE [START_REF] Iqbal | A split step wavelet method for radiowave propagation modelling in tropospheric ducts[END_REF][START_REF] Iqbal | Numerical modeling of radio wave propagation in horizontally inhomogeneous environment using split-step wavelet method[END_REF] have yielded a method as accurate as SSF. However, the computation complexity is the same as SSF since Fourier transforms are used.

In this paper, the split-step wavelet (SSW) method is proposed and applied to electromagnetic waves propagation in the low troposphere. Compared to DSSF, the major difference is that the fast wavelet transform (FWT) [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF][START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] is used in place of the FFT, our main purpose being to reduce the computation time. This method is based on a wavelet transform with respect to the vertical coordinate, which has a lower complexity than FFT. Moreover, after a compression step, the wavelet decomposition efficiently yields a sparse representation of a signal. The propagation is performed by a linear combination of the wavelets individually propagated, these individual propagations being stored in a pre-computed matrix. An efficient method for filling the propagation matrix is presented. The atmosphere and relief are accounted in the space domain. In addition, for considering an impedance ground, a local image method is proposed. This method reduces the computational complexity compared to the classical image method.

The article is organized as follows. In section 2, the configuration of the problem is presented. In section 3, the discrete wavelet transform (DWT) and FWT methods for decomposition and compression are briefly introduced. The choice of the wavelet basis is discussed and an example of wavelet transform and compression is presented. In section 4, the proposed SSW method is comprehensively introduced. Domain truncation, ground condition, irregular atmosphere, and relief are also considered. Then, the computational complexities of SSW and DSSF are compared. Section 5 is devoted to numerical tests demonstrating the accuracy and computational efficiency of SSW.

Configuration and discretization

In this work, we assume an e jωt time-dependence of the fields, where ω is the angular frequency. The aim is to simulate the propagation in a vertical plane (x, z) with y the direction of invariance. The source is located at x < 0 and the field u(x, z) is known at x = 0. The propagation is computed in the region x > 0, z ≥ 0. Any field can be decomposed in a transverse electric (TE) and a transverse magnetic (TM) components with respect to z. For the sake of simplicity, only the TE component is presented. The study of the TM component would be similar.

For numerical reasons, the computation domain is discretized and of finite size. The vertical domain is limited to z ∈ [0, z max [ and the range is limited to x ∈ [0, x max [. The discretization steps are ∆x and ∆z. The number of points along the 2 directions are N z = z max /∆z and N x = x max /∆x. The following uniform grid is used:

x = p x ∆x for p x = {0, . . . , N x -1}, z = p z ∆z for p z = {0, . . . , N z -1}. (1)
Finally, u[p z ] denotes u(p z ∆ z ). In the following parts, calculations are performed in this discrete domain.

Introduction to the discrete wavelet transform

In this section, DWT and FWT are briefly introduced. These methods are thoroughly described by Mallat [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. The choice of the wavelet basis for electromagnetic propagation is justified. Furthermore, a data compression criterion is introduced. Finally, decomposition and compression are applied on an example.

Discrete wavelet transform

The DWT leads to a representation of a discrete function as a linear combination of elementary functions. They are obtained by dilations and translations of a scaling function φ and a wavelet function ψ [START_REF] Mallat | Multiresolution approximations and wavelet orthonormal bases of L 2 (R)[END_REF]. For a scale l ∈ Z and a spatial index p ∈ Z, the scaling and wavelet functions are

φ l,p [p z ] = 2 -l/2 φ[2 -l p z -p] and ψ l,p [p z ] = 2 -l/2 ψ[2 -l p z -p], (2)
respectively.

As proven by Mallat [START_REF] Mallat | Multiresolution approximations and wavelet orthonormal bases of L 2 (R)[END_REF], for any L ∈ Z, (φ L,p ) p∈Z and (ψ l,p ) p∈Z,l∈ [1,L] form an orthonormal basis of the space of the finite-energy discrete functions. A discrete function u can thus be represented as the sum of an approximation component associated with the scaling function and multiresolution detail components associated with the dilated wavelets functions. The L-level multiresolution representation is given by

u[p z ] = p∈Z a L [p]φ L,p [p z ] + L l=1 p∈Z d l [p]ψ l,p [p z ], (3)
where a L [p] and d l [p] are the approximation and detail coefficients, respectively. Due to the finite size of the domain along z, for each level l, the p-indices are limited to [0, N p (l) -1] with N p (l) = N z /2 l . The total number of coefficients is N z . For the sake of clarity, in the following, the coefficients are represented by a vector U of dimension N z . The double index (l, p) is used to represent the elements of U , such that

U (l,p) = a L [p] for l = 0, p ∈ [0, N p (L) -1], d l [p] for l ∈ [1, L], p ∈ [0, N p (l) -1]. (4)
Besides, for the sake of conciseness, scaling functions and wavelets are both designated as wavelets hereafter. Thus, the representation is expressed by

u[p z ] = L l=0 p U (l,p) χ l,p [p z ], (5)
with

χ l,p [p z ] = φ L,p [p z ] for l = 0, ψ l,p [p z ] for l ∈ [1, L]. (6)
Because of the orthogonality of the basis, these coefficients can be calculated as the inner products of u and the basis elements. However, to calculate these coefficients in practice, the FWT has been developed by Mallat [39]. This operation has a O(N z ) complexity [START_REF] Mallat | Multiresolution approximations and wavelet orthonormal bases of L 2 (R)[END_REF], to be compared to the complexity of FFT in O(N z log 2 (N z )) [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. The inverse FWT, which is based on the same principle as the FWT, presents the same complexity.

Finally, a data compression is applied on the wavelet representation of the signal, taking advantage of the orthonormality of the basis. By forcing to 0 the coefficients that are weaker than a threshold V s , a regular signal can be approximated by a sparse set of coefficients within a chosen accuracy. After compression, the number of non-zero coefficients is usually much smaller than the original size of the signal N z . The vector U obtained after compression, denoted as Ũ , is a sparse vector with a maximal root mean square error (RMSE) of 20 log V s .

Note that the choice of the wavelet basis is important for the efficiency of the method. We choose the wavelet family "symlets" of order 6. These wavelets can be used for FWT. Furthermore, they have a compact support, a good regularity, and are almost symmetric [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. Therefore, they are good candidates for the SSW algorithm.

Wavelet decomposition and compression of an electromagnetic field

An example of FWT and compression applied to a field is shown in this section. A 2D complex source point (CSP) [START_REF] Deschamps | Gaussian beam as a bundle of complex rays[END_REF] with a frequency f 0 = 300 MHz is considered at x s = (x w0 +jkW 2 0 /2), y s = 0 m, z s = 1000 m, with x w0 = -50 m and W 0 = 5 m. This field is sampled on N z = 2048 points with ∆z = 1 m. Its amplitude is plotted in Figure 1 in blue.

The compression rate (CR) is defined as

CR = Number of zeros Total number of coefficients . (7)
For L = 3, V s = 10 -3 , the comparison of the coefficients before and after compression is shown in Figure 2. The number of non-zero compressed coefficients is 22, which is much smaller than N z . Furthermore, we see that the detail coefficients at level 1 are all set to zero. The signal is only represented by the coefficients of level 2 and 3. The signal is automatically described by the wavelets at appropriate dilations.

The signal reconstructed from the wavelet coefficients is given in Figure 1. The RMSE between the two signals (original and reconstructed from the wavelets after compression) is -61.6 dB.

The next section presents how to efficiently propagate this sparse field representation in complex environments. 4 Split-step wavelet method

Based on the wavelet transform and compression introduced in the previous sections, the SSW method is derived in the present section. The overview of the method is exposed. Then, specific points are addressed: the construction of the propagation matrix, the apodisation, the local image method, the accounting of the relief and atmosphere, and the computational complexity.

Overview of the split-step wavelet method

SSW is performed going back and forth from a spatial to a wavelet representation of the wave so as to model the propagation iteratively at increasing distances. The field is supposed to be known at x = 0. Then, the propagation from x to x + ∆x, with x = p x ∆x, is simulated step by step as follows:

1. The field u(x, z) is represented as a sparse wavelet vector Ũ (x) after applying FWT (denoted as W ) and a compression with threshold V s (denoted as C) as introduced in Section III.D, such that

Ũ (x) = CW u(x, z). (8)
2. The free-space propagation step is considered in the wavelet domain. The coefficients vector U fs (x + ∆x) after free-space propagation on a distance ∆x is given by

U fs (x + ∆x) = M Ũ (x), ( 9
)
where M is the pre-computed free-space propagation matrix, detailed in section 4.2, that models the waveletto-wavelet propagations. This matrix is compressed with a threshold V M . The elements of M , i.e. M (l,p),(l ,p ) , satisfy

U fs(l,p) (x + ∆x) = l ,p M (l,p),(l ,p ) Ũ(l ,p ) (x). (10)
3. The field u fs (x + ∆x, z) propagated in free space is recomposed by an inverse FWT (denoted as W -1 ) from the propagated wavelets.

u fs (x + ∆x, z) = W -1 U fs (x + ∆x). ( 11
)
4. Apodization, atmosphere, and relief can be applied in the spatial domain as in the classical SSF algorithm [START_REF] Dockery | An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation[END_REF], here represented by the operator D.

As a conclusion, the propagation from x to x + ∆x is simulated step by step as

u(x + ∆x, z) = DW -1 M CW u(x, z). (12)

Pre-computation of the free-space propagation matrix

In this section, we present the pre-computation of the matrix M . We take into account the translation and dilation properties of the wavelet basis so as to reduce its computation cost. The aim is to obtain M (l,p)(l ,p ) for all l, l , p, and p . The indices l and l correspond to the wavelet levels before and after propagation, respectively. The indices p and p correspond to the positions along the vertical axis before and after propagation, respectively. The computation of the matrix comprises two steps:

• Wavelet propagation on a distance ∆x: For the first step, only one wavelet χ l ,0 on level l is considered. It is then propagated on ∆x using DSSF. (Any other propagation method can be used here.) As illustrated in Figure 3, for any p , the propagated wavelets χ l ,p (∆x, z) can be deduced from χ l ,0 (∆x, z) by means of translations of p .

• Calculation of the matrix coefficients: The coefficients M (l,p)(l ,0) are obtained by inverse FWT and compression.

The other ones exploit the translation properties of M given here.

For l = l , the number of wavelets on the levels l and l are the same, so a translation of 1 on p results in a translation of 1 on p. Thus, M satisfies

M (l,p+1)(l ,p +1) = M (l,p)(l ,p ) for l = l . ( 13
)
Figure 3: Wavelets at level l propagated on ∆x on positions 0 and p .

Secondly, for l = l + 1, they are twice as many wavelets on level l than on level l. Thus, a translation of 2 on p results in a translation of 1 on p. Thus, we have

M (l,p+1)(l ,p +2) = M (l,p)(l ,p ) for l = l + 1. ( 14
)
Finally, for l = l -1, they are half as many wavelets on level l than on level l. In other words, a translation of 1 on p results in a translation of 2 on p. All the even coefficients p are obtained this way. The odd coefficients are obtained from translations of the inverse FWT of χ l ,1 . (An additional inverse FWT is required.)

The rule can be generalised for any l and l by

if l ≤ l , M (l,p+2 l -l )(l ,p +1) = M (l,p)(l ,p ) , ( 15
) if l > l , M (l,p+1)(l ,p +2 l-l ) = M (l,p)(l ,p ) . (16)
These properties are used to efficiently fill the matrix M . The complete procedure is given in algorithm 1.

Algorithm 1 Creation of the propagation matrix M .

Input: Maximum level of decomposition L, wavelet type for l ∈ [1, L] do χ l ,0 (0, z) ← inverse FWT on the wavelet coefficient (l , 0) set to 1, others set to 0. χ l ,0 (∆x, z) ← propagation of χ l ,0 (0, z) on ∆x by DSSF. M (l,p)(l ,0) ← FWT and compression applied to χ l ,0 (∆x, z). for l ∈ [1, L] do if l ≤ l then for q l ∈ [0, N p (l)] do M (l,p)(l ,q l ) ← duplication and translation of M (l,p)(l ,0) following [START_REF] Sturm | Numerical study of broadband sound pulse propagation in three-dimensional oceanic waveguides[END_REF].

if l > l then for q ∈ [0, 2 l-l [ do χ l ,q l (∆x, z) ← translation of χ l ,0 (∆x, 0) on q coefficients. M (l,p)(l ,q) ← FWT and compression applied to χ l ,q (∆x, z). for q l ∈ [0, N p (l)] and q l = q mod (2 l-l ) do M (l,p)(l ,q l ) ← duplication and translation of M (l,p)(l ,q) following ( 16).

return M

The cases l = 0 and l = 0 corresponding to the scaling functions have not been included for the sake on conciseness. They must be treated as additional levels L in both loops. The threshold used for the compressions during the creation of the matrix M is denoted as V M .

Finally, only L + 1 DSSF propagations and 2 L FWT are necessary to calculate the entire matrix M .

Domain truncation

In order to remove spurious reflections over the top boundary and reduce the amount of computation, many methods could be used such as apodization layers, perfectly matched layers, or non-local boundary conditions. In this paper, we use an apodization layer. A Hanning window [START_REF] Levy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF] is applied to the field at the upper half part.

Ground boundary condition: the local image method

In DSSF, if an impedance condition is applied at the ground boundary, the Fourier transform is replaced by the discrete mixed Fourier transform [START_REF] Dockery | An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation[END_REF]. We cannot reproduce this strategy here because the wavelet transform is used.

To consider a perfectly conducting or an impedance ground in SSW, we introduce a local image algorithm. The key idea is to add the symmetric field affected by its reflection coefficient in a thin image layer under the ground level. This field is reset at each propagation step. The local image algorithm is illustrated in Figure 4. The method follows these three steps:

1. At range index p x , we know the field u described on N z points. Then it is extended to u t by adding a local image layer under u on N i points. The choice of N i is discussed later. In this layer, the field is the image of the upper part. The whole field u t is defined by

u t [p x , p z ] =      u[p x , p z ] for p z ∈ [1, N z -1], (1 + Γ)u[p x , 0] for p z = 0, Γu[p x , -p z ] for p z ∈ [-N i , -1], (17)
with Γ the Fresnel coefficient of the impedance ground, considered as constant at each step. Note here, the ground wave is not considered in this method. Besides, in the case of a PEC ground, Γ = -1.

2. u t is propagated in free space from x to x + ∆x by using the SSW propagation strategy.

3. The propagated field u p [p x + 1, p z ] is obtained as

u p [p x + 1, p z ] = u t [p x + 1, p z ], for p z ∈ [0, N z -1]. ( 18
)
The previous algorithm is repeated at each step. As only the image of the field in the close vicinity of the ground is used, we call this algorithm the local image method.

At the bottom of the image layer, no apodization is applied. Therefore, an error appears at the bottom of the domain due to the boundary. However, N i is chosen so that this error does not reach the domain of interest (i.e., p z ≥ 0). At each step, the field in the image layer is updated using [START_REF] Zaporozhets | Bistatic RCS calculations with the vector parabolic equation method[END_REF]. Thus, the error in the image layer is cleared up at each step and never reaches the domain p z ≥ 0.

N i is chosen large enough so that the spurious reflections at the bottom of the image domain do not reach the domain p z ≥ 0 after one computation step. If only propagated waves are accounted (∆z > λ/2), the chosen criterion is the maximum width of the wavelets after propagation on ∆x. As the widening of one wavelet on one step ∆x is much smaller than the total computation domain, N i N z . This makes the local image method computationally efficient.

Atmosphere and relief

To consider a slowly varying atmosphere, the phase screen method [START_REF] Dockery | An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation[END_REF] is applied to the field. Irregular relief is considered by using a staircase model [START_REF] Levy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF] exactly as for SSF. More precise algorithms [START_REF] Barrios | A terrain parabolic equation model for propagation in the troposphere[END_REF][START_REF] Donohue | Propagation modeling over terrain using the parabolic wave equation[END_REF][START_REF] Janaswamy | A curvilinear coordinate-based split-step parabolic equation method for propagation predictions over terrain[END_REF] could also be applied.

Complexity comparison

Compared to the classical DSSF method, the Fourier transform is replaced by the wavelet decomposition and recomposition. At each range step, the complexity of FFT is O(N z log N z ), and the complexity of FWT is

O(N z + N i ) ≈ O(N z ) because in practice N i N z .
The computation complexity of the propagation step in DSSF is O(N z ). The propagation in SSW follows [START_REF] Dockery | An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation[END_REF]. Since M and Ũ are sparse. The effective multiplication cost N e is

N e = (l,p)∈Inz( Ũ ) N nz M (l,p),(l ,p ) ≤ max (l,p) N nz M (l,p),(l ,p ) N M N nz ( Ũ ) Ns , ( 19
)
where I nz and N nz give the indices and the number of non-zero coefficients in a vector, respectively. N M is the maximum number of non-zero coefficients in one column of M . N s is the number of non-zero coefficients in Ũ . The complexity of this step in SSW is O(N e ) with N e ≤ N M N s . In practice, N M and N s are much smaller than N z due to the high compression rate of the wavelet decomposition. The larger the signal and matrix thresholds V s and V M , the smaller N M and N s . The comparison of the computational complexities of SSW and DSSF is given in Table 1. With a high CR, SSW is faster than DSSF. As the complexities of the algorithms are linear with N x , the total complexities are obtained by multiplying by N x . Note here, for the creation of the propagation matrix M , the complexity of the DSSF propagations is O(LN z log N z ), and the complexity of the wavelet decompositions is O(2 L N z ). They are also negligible for large values of N x .

Numerical tests

In this section, several numerical tests are performed to validate the SSW method. A propagation without reflection validates the matrix propagation strategy. Then a simulation with a ground reflection validates the local image method. Finally, a complex case with varying atmosphere and relief is presented. A parametric study to discuss the computation time on this scenario concludes this section. All the simulations are performed on a 1.8 GHz core i7 CPU.

Propagation without reflection

The propagation of a CSP without gound reflection is studied. The aim is to test the accuracy of SSW with different threshold values. The accuracy of SSW is given by comparison to DSSF.

The parameters of the CSP are: frequency f = 300 MHz, x s = x w0 + jkW 2 0 /2, y s = 0 m, z s = 2000 m, with x w0 = -50 m and W 0 = 5 m. The computation grid is x max = 1 km, ∆x = 10 m, z max = 4096 m, and ∆z = 1 m. Thus, N x = 100 and N z = 4096. The chosen wavelets are symlets 6, with the maximum level L = 3.

In order to test the compression thresholds on the matrix M and on the signal, three cases are considered: no compression (V M = V s = 0), only matrix compression (V M = 2 × 10 -3 and V s = 0), and only signal compression (V M = 0 and V s = 2 × 10 -3 ). For all cases, the fields at the maximum range and the differences to DSSF are plotted with respect to altitude in Figure 5a and Figure 5b, respectively. The fields at the maximum range show a very good agreement. The differences with DSSF are only due to the compression. Without compression, the RMS difference of SSW to DSSF (not visible on Figure 5b) is -165.4 dB.

The evaluation of the RMS difference normalised at each step in range is plotted in Figure 5c in logarithmic scale. For V M = 0 and V s = 0 (green line), the error accumulates in range. The error yielded by the compression on M is of order of V M N x (red dashed line). In the cases V M = 0 and V s = 0 (black line), the error increases in range. However, it is much lower than V s N x (about -23 dB smaller at the maximum range). A least squares fit yields to an error proportional to V s N α

x , with α ≈ 0.3. Consequently, V s can be chosen much larger than V M in practice. Note that this specific value for α only stands for simple scenarios, when no ground is accounted. Many simulations with different configurations have led to a coefficient α always lower than 0.5. Therefore, the error due to the signal compression is considered of the order of V s N 0.5

x . A demonstration of this empirical result has not been obtained yet.

Propagation over a planar impedance ground in a homogeneous atmosphere

The aim of this section is to test SSW with the local image method described in Section 4.4 and to confirm the range dependency of the compression error. The same source as the previous test is chosen with a height of z s = 30 m. The simulation parameters are SSW with compressions on both M and the signal is tested. The expected error V E on the field at the final range corresponds to the sum of both the compressions, such that

V E = V M N x + V s N α
x . The chosen thresholds are V M = 2 × 10 -5 and V s = 4.47 × 10 -4 for the errors due to both compressions to be the same, equal to -40 dB. Thus, the total expected error is of the order of -34.0 dB.

The propagated field along range and altitude is shown in Figure 6a, where the interference pattern due to the ground reflection is visible. During the simulation, the CR of the propagation matrix is 86.4 % and the average CR of the signal is 74.2 %. The RMS difference to DSSF on the last vertical is -47.3 dB. This RMS difference is -13.3 dB better than the expected value.

To explain this difference, the error to DSSF with respect to distance is plotted in Figure 6b with compression on the matrix only (in green) and with compression on the signal only (in black). For both cases, the expected final error is -40 dB. When some energy goes out of the computation domain, after 100 iterations corresponding to a distance of 20 km, the error increase is slower. This may be explained by the fact that some error is removed in the apodization layer. Thus, the actual error due to compression on the last vertical is usually less than the expected error V E .

Finally, the computation time of SSW is 3.2 s whereas the DSSF simulation lasts 8.6 s. Therefore, SSW is faster than DSSF. Besides, the computation time for creating M is 1.2 s. As a conclusion, SSW with an impedance ground is successfully tested in terms of accuracy and shows its advantage concerning the computation time.

Long-range propagation over an impedance irregular relief in an inhomogeneous atmosphere

In this section, we test a long-range propagation over an irregular relief in an inhomogeneous atmosphere.

In the vertical direction, a surface-based duct is considered. It is modelled by a trilinear modified refractivity profile, as illustrated in Figure 7a In this test, the propagation range is x max = 100 km. The range step is 200 m, giving N x = 500. We choose the same compression thresholds V M = 2 × 10 -5 and V s = 4.47 × 10 -4 . That yields an expected error of the order of -34.0 dB.

The normalised field propagated using SSW in the vertical plane is shown in Figure 7c. The refractive effects of the surface-based duct and the reflection over the irregular ground are properly simulated. The fields at the maximum range and the difference to DSSF are plotted in Figure 7d. The RMSE of the final field of SSW to DSSF is -42.0 dB, 8 dB below the expected error. Therefore, SSW works well for a long range simulation over an irregular relief in an inhomogeneous atmosphere. As for the previous vase, the computation time is 3.2 s for SSW and is 8.6 s for DSSF. SSW has a better computation efficiency than DSSF, even on complex scenarios.

Discussion about the computation time

As stated in Section 4.6, the complexity of SSW depends on the number of vertical points and on the compression rates applied to both the propagation matrix and the signal at each distance step. In this section, this complexity is illustrated. We use the simulation setup of Section 5.3. All else being equal, the vertical number of points N z or the compression rates are varied. The chosen values and the corresponding computation times obtained with DSSF and SSW are summarized in Table 2. If no compression is applied (V M = 0 and V s = 0), SSW is much slower than DSSF. Indeed, even if the wavelet transform is faster than the Fourier transform, the propagation step requires a full matrix mutliplication in the absence of compression.

If smaller compression thresholds are chosen with SSW, some computation time can be saved. Following Section 5.1, the choices V M = 2 × 10 -4 and V s = 4.47 × 10 -3 lead to an acceptable expected error of -16 dB. The actual simulation gives a computation time shortened by a factor of 2 compared with the previous section and a final RMSE with respect to DSSF of -22.2 dB.

When N z is greater, the computation time increases less with SSW than with DSSF, which is consistent with the complexities of Table 1. More optimised codes would certainly match better the theoretical complexities. Finally, note that the multiplication cost O(N e ) in SSW could still be reduced because the sparse matrix-vector multiplication we use does not take advantage of the sparsity of the vector yet.

Conclusion

The split-step wavelet method (SSW) for the simulation of long-range electromagnetic wave propagation has been introduced. This method is based on the discrete wavelet transform and the discrete electromagnetic theory. First, the discrete wavelet transform (DWT) and the fast wavelet transform (FWT) methods have been briefly introduced, and the choice of the wavelet family has been discussed.

Then, the proposed SSW method has been comprehensively introduced. The field is represented as a sparse set of coefficients after applying FWT and thresholding. It is then propagated in free-space by means of a pre-computed propagation matrix. The propagated field is recomposed by inverse FWT from the propagated wavelets. Relief, atmosphere, and apodisation are treated in the same way as in split-step Fourier (SSF) techniques. The total field is obtained iteratively.

The strategy to efficiently create the pre-computed propagation matrix has been introduced. Finally, to consider a perfectly conducting or an impedance ground condition, a local image method has been proposed. It is based on the localisation property of the wavelets and does not significantly increases the computation burden.

The computation complexity of SSW has been compared to the discrete SSF (DSSF). SSW is shown to be more efficient than DSSF due to the high compression rate of the wavelet decomposition and the low complexity of the fast wavelet transform.

Finally, numerical tests of wave propagation have been presented to show the accuracy and efficiency of this method. First, propagation in free-space with different thresholds have been tested and compared. Tests of propagation over an impedance ground with both matrix and signal compressions have been performed. We have shown that the error due to compression can be anticipated. Finally, a long-range propagation over an irregular relief in an inhomogeneous atmosphere has been performed. The result of SSW shows a very good match with DSSF, and the simulation parameters are varied to illustrate the complexity of SSW. In conclusion, SSW works well for a long-range simulation with reduced computation time.

Even if the choice of the symlets 6 wavelets has been carefully thought, an exhaustive comparison of the wavelets could lead to a more efficient choice. Moreover, the computation time could be shortened by reducing the size of the pre-computed data and improving the top boundary condition. A theoretical study to obtain an upper bound of the compression error is also necessary. Finally, the method will be extended to 3D configurations in future works.

Figure 1 :

 1 Figure 1: Amplitude of the signal radiated by the CSP.

Figure 2 :

 2 Figure 2: Coefficients of the wavelet decomposition: non-zero coefficients are in red, null coefficients in blue.

Figure 4 :

 4 Figure 4: Local perfect electric conductor ground condition for SSW.

  (a) Normalised electric fields (dB) at the maximum range obtained by DSSF (doted line) and SSW (full lines) with different compressions. (b) Differences of SSW to DSSF for matrix and signal compressions. (c) RMS difference (dB) of SSW fields with DSSF at each step in range.

Figure 5 :

 5 Figure 5: Propagation of a CSP in free-space.

( a )

 a Normalised electric field (dB) in the vertical plane obtained with SSW. (b) RMS difference (dB) of SSW fields with DSSF at each step in range.

Figure 6 :

 6 Figure 6: Propagation of a CSP with reflections over an impedance ground.

  . The duct parameters are: M 0 = 330 M-units, z b = 100 m, z t = 200 m, z max = 4096 m, with gradients c 0 = 0.118 M-units/m, c 2 = -0.1 M-units/m. The relief is chosen as 2 small triangular hills of height 100 m and 200 m. The relief is displayed in Figure 7b. The characteristics of the impedance ground are chosen as r = 20 and σ = 0.02 S/m.

  A relief with 2 small triangular hills of heights 100 m and 200 m. (c) Normalised electric field (dB) obtained by SSW. (d) Normalised electric fields (dB) at the maximum range obtained by DSSF and SSW and the difference SSW-DSSF.

Figure 7 :

 7 Figure 7: Propagation over a planar impedance ground in an inhomogeneous atmosphere and over an irregular relief.

Table 1 :

 1 Computational complexity of SSW and DSSF on one step ∆x. SSW DSSF wavelet dec. / DFT O (N z ) O(N z log N z ) propagation O(N e ) O(N z ) wavelet rec. / IDFT O (N z ) O(N z log N z )

Table 2 :

 2 Computation times for varying N z and compression thresholds.

				comp. time	comp. time
	V M	V s	N z	DSSF (s)	SSW (s)
	0	0	4096	8.6	21.1
	2 × 10 -4 4.47 × 10 -3 s 4096	8.6	2.0
	2 × 10 -5 4.47 × 10 -4 s 1024	2.2	0.9
	2 × 10 -5 4.47 × 10 -4 s 8192	34.6	8.0

Acknowledgements

This work was supported by the Fundamental Research Funds for the Central Universities, project number 3122016U001.