Ke Dong

Itumeleng Moroenyane

Binu Tripathi

Dorsaf Kerfahi

Koichi Takahashi

Naomichi Yamamoto

Choa Cho

Hyunjun Cho

Jonathan Adams

Choa An

Soil nematodes show a mid-elevation diversity maximum and elevational zonation on Mt. Norikura, Japan

Ecologists have long been fascinated by how communities and individual taxa respond to elevational gradients. Elevational trends in diversity and community structure have been well documented among a wide range of forms of life, including animals, plants, insects or some other larger invertebrates [START_REF] Heaney | Frontispiece: diversity patterns of small mammals along elevational gradients[END_REF][START_REF] Jacquemyn | Elevational gradients of species diversity, breeding system and floral traits of orchid species on Reunion Island[END_REF][START_REF] Hodkinson | Terrestrial insects along elevation gradients: species and community responses to altitude[END_REF][START_REF] Powers | Distribution, community structure, and microhabitats of soil invertebrates along an elevational gradient in Taylor Valley, Antarctica[END_REF] . In recent years, elevational trends in community structure and diversity of microorganisms have also been discovered, due to advances in extraction and sequencing of environmental DNA [START_REF] Bryant | Microbes on mountainsides: contrasting elevational patterns of bacterial and plant diversity[END_REF][START_REF] Singh | Elevational Patterns in Archaeal Diversity on Mt. Fuji[END_REF][START_REF] Miyamoto | The mid-domain effect in ectomycorrhizal fungi: range overlap along an elevation gradient on Mount Fuji, Japan[END_REF] . However, there has so far been very little concerted effort to determine how diversity patterns and community structure of soil nematodes vary with elevation, despite their abundance and important roles in ecosystem processes [START_REF] Bongers | Nematode community structure as a bioindicator in environmental monitoring[END_REF] .

Nematodes are considered to be the most abundant animals on Earth [START_REF] Bongers | Nematode community structure as a bioindicator in environmental monitoring[END_REF] . They occupy almost any environment that provides an available source of organic carbon, playing a significant role in nearly all the world's ecosystems [START_REF] Yeates | How plants affect nematodes 61-113[END_REF] . Because of the key position of nematodes as primary and intermediate consumers in soil food webs, assessing the controls on diversity of nematode community structure can be seen as being of great importance for assessing future environmental change effects [START_REF] Bongers | Nematode community structure as a bioindicator in environmental monitoring[END_REF][START_REF] Ritz | Utility of nematode community analysis as an integrated measure of the functional state of soils: perspectives and challenges[END_REF] .

To date, only one study has focused on diversity patterns and community structure of soil nematodes vary with elevation, finding no evidence of changes in nematode community diversity across a 1,618 m range of elevation in the USA [START_REF] Treonis | Soil nematodes and their prokaryotic prey along an elevation gradient in the Mojave Desert (Death Valley National Park, California, USA[END_REF] . It is likely that T-RFLP, the low-resolution molecular technique used in this study, greatly limited the applicability of diversity and community structure. Given the lack of existing studies of soil nematode communities in relation to elevational gradients, it is important to add new work on mountain systems to understanding what patterns exist and what environmental factors may cause these patterns. As well as being of interest in the study of nematode ecology, this would be an important contribution to arriving at general patterns and principles in ecology.

Here, we set out to conduct systematic and representative study focusing on soil nematodes, using high resolution molecular techniques over a broad elevational range. Reviewing studies of elevational diversity trends on a broad range of other taxa, Rahbek [START_REF] Rahbek | The role of spatial scale and the perception of large-scale species-richness patterns[END_REF][START_REF] Rahbek | The elevational gradient of species richness: a uniform pattern?[END_REF] found that around half showed a mid-elevation maximum. Since this is the most common elevational pattern, it is most likely that nematode community will also follow the same pattern, showing a diversity maximum at mid-elevations. Aside from the search for underlying mechanisms, we were firstly interested in knowing whether this empirical pattern is still more widespread, occurring in yet another group of organisms. Moreover, we anticipated that the tendency for other groups of organisms to show mid elevation maxima might actually interact with nematode ecology to produce this pattern in nematodes. If the diversity of trophic sources is greater in mid elevations -due to greater taxonomic diversity in a group that is a food source -this can be expected to allow a greater number of nematode niches to coexist. No relevant studies of elevational trends in other groups of soil organisms exist for Norikura, with the exception of archaea which also show a mid-elevation diversity maximum on Norikura [START_REF] Singh | Similarities and Contrasts in the Archaeal Community of Two Japanese Mountains: Mt. Norikura Compared to Mt. Fuji[END_REF] . However, on nearby Mt. Fuji, which has a very similar climate and reaches a similar height to Mt. Norikura, there are mid elevation diversity maxima in both bacteria [START_REF] Singh | A Hump-Backed Trend in Bacterial Diversity with Elevation on Mount Fuji, Japan[END_REF] and fungi [START_REF] Miyamoto | The mid-domain effect in ectomycorrhizal fungi: range overlap along an elevation gradient on Mount Fuji, Japan[END_REF] , as well as Archaea [START_REF] Singh | Elevational Patterns in Archaeal Diversity on Mt. Fuji[END_REF] . As these are groups that many soil nematodes feed on [START_REF] Yeates | Feeding habits in soil nematode families and generaan outline for soil ecologists[END_REF] , the greater diversity of these food sources could hypothetically contribute to a mid-elevation maximum in nematode diversity.

A further reason to anticipate a mid-elevation diversity maximum for nematodes is the predominance of a mid-elevation maximum in precipitation across high mountains in central Japan [START_REF] Fujimura | Of the climate and weather of Mt. Fuji in Report of the scientific survey of Mt[END_REF] . For example, Mt. Fuji has an annual precipitation maximum at 2,000-2,500 masl [START_REF] Fujimura | Of the climate and weather of Mt. Fuji in Report of the scientific survey of Mt[END_REF] . This coincides with the elevation of the predominant cloud layer, and although the elevational trend of precipitation on Norikura has not been measured, the cloud layer likewise tends to occur at around 2,000-2,500 masl (pers. obvns. by the authors). A mid elevation precipitation maximum for Norikura is potentially significant for soil nematodes because they are strongly dependent on moving within the soil water that surrounds and fills spaces between soil particles [START_REF] Stirling | Biological control of plant-parasitic nematodes: soil ecosystem management in sustainable agriculture[END_REF] . Greater precipitation, combined with cooler temperatures in mid elevations, may be expected to give moister -but still aerated -soils that are physiologically favorable to nematodes and their feeding activity. This itself may promote greater diversity of nematodes because there are fewer restrictions of survival and feeding, allowing greater overall population densities even in specialized niches, promoting niche specialization, and thus species diversity. The precipitation/ water balance maximum in mid elevations has been suggested as a cause of diversity maxima in other parts of the world, in other groups that are particularly sensitive to water balance [START_REF] Wolf | Diversity patterns and biomass of epiphytic bryophytes and lichens along an altitudinal gradient in the northern Andes[END_REF][START_REF] Olson | The distribution of leaf litter invertebrates along a Neotropical altitudinal gradient[END_REF] .

There has been much discussion of the possible reasons why mid-elevational diversity maxima occur so commonly [START_REF] Colwell | The mid-domain effect: geometric constraints on the geography of species richness[END_REF][START_REF] Mccain | The mid-domain effect applied to elevational gradients: species richness of small mammals in Costa Rica[END_REF][START_REF] Lomolino | Elevation gradients of species-density: historical and prospective views[END_REF] . One explanation for the mid elevation diversity maximum is the 'community overlap' hypothesis, which assumes that there are often two relatively distinct environments on a mountain -an 'upper mountain' environment subject to more extreme low temperatures or a different precipitation regime (and often above the treeline), and a 'lower mountain' environment that has a distinct ecology (and is often by contrast forested) [START_REF] Lomolino | Elevation gradients of species-density: historical and prospective views[END_REF] . If each of these two environments has its own distinct set of species, at the transitional zone where they intersect both groups of species will occur together or nearby one another, giving a mid-elevation diversity maximum [START_REF] Lomolino | Elevation gradients of species-density: historical and prospective views[END_REF] .

Another widely favoured explanation for mid elevation diversity maxima is the 'mid-domain effect' (MDE), which states that if all species ranges are scattered randomly between the limits of the top and bottom of a mountain, there will be a 'bulge' of maximum numbers of overlapping species in the mid elevations [START_REF] Colwell | The mid-domain effect: geometric constraints on the geography of species richness[END_REF] . A recent advance of MDE theory has been to include a midpoint attractor -a unimodal gradient of environmental favorability, using a Bayesian simulation model to estimate the location and strength of the attractor from empirical species distribution data along the elevations, within geometric constraints [START_REF] Colwell | Midpoint attractors and species richness: Modelling the interaction between environmental drivers and geometric constraints[END_REF] . It has been suggested that gradients of environmental favorability, together with the geometric constraints imposed by the base of a mountain and its summit, will more parsimoniously explain elevational species richness patterns.

We were interested in testing whether the diversity pattern of soil nematode communities would differ with elevation and show a mid-elevation maximum, and whether this pattern would be consistent with the leading hypotheses. Besides testing these predictions for overall trends on mountains, we were also interested in a more generalized way in understanding which particular environmental variables contribute to whatever elevational diversity trend exists, since studies have shown that nematode communities are sensitive to the surrounding environment and several environmental variables. For example, salinity, temperature, pH, soil moisture and organic nutrients have been identified as important factors influencing diversity of nematode communities [START_REF] Powers | Distribution, community structure, and microhabitats of soil invertebrates along an elevational gradient in Taylor Valley, Antarctica[END_REF][START_REF] Freckman | Low-diversity Antarctic soil nematode communities: distribution and response to disturbance[END_REF] .

We were also interested in investigating the community structure of soil nematode in relation to elevation, testing the hypothesis that niche specialization will produce a succession of communities along elevation gradients. It is known that for many taxa, there are distinct and predictable community assemblages of organisms which occur at particular elevational levels on mountains [START_REF] Patterson | Contrasting patterns of elevational zonation for birds and mammals in the Andes of southeastern Peru[END_REF][START_REF] Mccoy | The distribution of insects along elevational gradients[END_REF][START_REF] Singh | Strong elevational trends in soil bacterial community composition on Mt. Halla, South Korea[END_REF] . This is a result of the fact that many types of organism have clearly defined spatial ranges that are linked to climate and the ecosystem conditions that parallel elevation. The overlap of these ranges is what produces the characteristic communities of each elevational zone. Given that this elevational turnover pattern in species and communities is so widespread, we hypothesized that nematodes, too, would have finely adapted environmental niches in relation to particular elevations, and that this would produce a series of distinctive communities along the elevational gradient.

Lastly, considering that empirical range data of OTUs would be collected and used in the simulation of null model in midpoint attractor model, we were interested in testing for the elevational pattern of species ranges in term of Rapoport's elevational hypothesis. This is a derivative of Rapoport's rule [START_REF] Stevens | The latitudinal gradient in geographical range: how so many species coexist in the tropics[END_REF] , which notes that higher latitude species tend to have larger geographical ranges. The elevational variant of this rule predicts that species found at higher elevations will tend to occur over a broader range of elevations than those from lower down the mountain [START_REF] Stevens | The elevational gradient in altitudinal range: an extension of Rapoport's latitudinal rule to altitude[END_REF] . The pattern has been documented along elevational gradients in a variety of different taxa, and it has been suggested that it occurs because climates at higher elevations are more variable, so species that occur there must be more tolerant of temperature variation -which also translates into larger elevational ranges [START_REF] Stevens | The elevational gradient in altitudinal range: an extension of Rapoport's latitudinal rule to altitude[END_REF] .

Results

Broad taxonomic features of soil nematode community on Mt. Norikura. In total, 55 samples were collected from 11 elevations on the mountain, along elevational isoclines bands separated by ~200 m of elevation. Four samples were excluded from further nematode community analysis owing to failure of amplification of nematode 18S rRNA genes (Supplementary Table 1).

A total of 70,227 high quality 18S rRNA gene sequences were assigned to 1,121 operational taxonomic units (OTUs) at ≥99% similarity level from 51 samples after the removal of low quality, chimeric, and rare sequences. On average, 98.6 OTUs (±4.4 SE) were found in each sample, which were collected along a ~2,200 meters elevational range (Supplementary Tables 2 and3). The greatest OTU richness was found at 1,105 meters above sea level (masl) with 149 OTUs, while the lowest OTU richness was found at 2,941 masl with 22 OTUs (Supplementary Table 2).

The nematode community on Mt. Norikura as a whole was dominated by Prismatolaimidae which accounted for 15.3% of the total reads, followed by Mononchidae (10.7%), Qudsianematidae (10.6%), Rhabditidae (5.1%), Criconematidae (5.0%), Plectidae (4.9%), Chromadoridae (4.5%), Cephalobidae (3.9%) and Tripylidae (3.8%). Remaining 18.1% of the total reads were shared by 24 nematode families together, and 17.9% nematode sequences remained unclassified on family level (Fig. 1).

Classifiable reads were binned into four feeding guilds (bacteria feeding guild [BF], fungi feeding guild [FF], plant feeding guild [PF] and omnivore/predator guild [OP]) according to their family taxonomy and all feeding guilds of nematode were found on Mt. Norikura (Fig. 2). BF and OP together accounted for a major portion of 72.3%, while FF were found low across all elevations of 2.4% in total. The percentage of PF was low on high elevations compared to the lower elevations below 2,500 masl (Fig. 2).

Elevational change in diversity and community zonation of soil nematodes. OTU richness,

Shannon index and Faith's PD were used to assess the ecological pattern in nematode diversity. All of the diversity indices were significantly correlated with elevation (p < 0.001), and the richness/diversity showed a mid-elevation Relative abundance (%) of nematode feeding guilds at eleven different elevations. BF, bacteria feeding; FF, fungi feeding; PF, plant feeding; OP, omnivore/predator. Note the sum of the relative abundance of the four feeding groups is less than 100 percent, due to a sizeable proportion of unclassified sequences. maximum, with maximum diversity for all three indices falling between 1,000 and 2,000 masl (Fig. 3). Quadratic models gave the best fit for all three diversity indices (data not shown).

An NMDS performed on the Bray-Curtis similarity matrix of nematode community structure showed significant variability in relation to the elevational gradient (ANOSIM: Global R = 0.761, p < 0.001) (Fig. 4). The samples also formed clear clusters according to the elevation that they came from, visually indicating that samples belonging to different elevational zones harbored distinct communities.

Nestedness analysis indicated that the communities of soil nematodes along the elevational gradient follow a nested structure (p ≪ 0.001). A packed matrix order categorizing nestedness of each sample from low to high was generated, with samples from elevations above 2,500 masl intensively clustered towards the top of the packed matrix order (Supplementary Table 3). This indicates that the communities of higher elevations (above 2,500 m) tend to be a subset of the communities of lower elevations (below 2,500 m).

Mid domain effect on soil nematode diversity. Each of the empirical ranges of the 1,121 OTUs (mean midpoint of 0.536, mean range of 0.267) was stochastically placed on the domain rescaled from a 3,026 m elevational gradient without replacement of 200 simulations. The result indicated that the simulation is driven by a midpoint attractor at 0.417, with a standard deviation of 0.299 (Fig. 5a). Regression exhibited a good fitness of modeled richness against empirical richness (R 2 = 0.746, p < 0.001, Supplementary Fig. 1), indicating that gradients in elevational favorability, subject to geometric constraints, may account for the unimodal elevational richness pattern of soil nematodes. Note that the empirical richness, which is underlain by the inferred elevational ranges of OTUs, is not the same concept as the OTU richness of each replicate in Fig. 3. Midpoint and range of OTUs were plotted in Fig. 5b.

Physicochemical factors influencing soil nematode communities. Nine physicochemical environ-

mental factors for the 55 samples were measured directly, or for elevation calculated on the basis of the moist air lapse rate [START_REF] Mark | Fundamentals of Atmospheric Modeling[END_REF] , with the parameters being statistically different between elevations (p < 0.01) (Supplementary Table 4). Mean annual temperature (MAT), which was calculated using a mean lapse rate of 0.6 °C/100 masl, showed a complete correlation with elevation. The other physicochemical factors including total carbon (TC), total nitrogen (TN), available phosphate (P 2 O 5 ), nitrogen in ammonium (NH 4 -N), nitrogen in nitrate (NO 3 -N), potassium (K), pH and soil texture all showed linear unimodal patterns against elevation (Fig. 6).

The multiple regression analysis showed the relative influence of physicochemical factors on nematode diversity indices. Seven variables were tested (see Methods and Materials; Supplementary Fig. 2) and TN and MAT were significantly correlated with all the richness/diversity indices (Table 1).

Canonical correspondence analysis (CCA) suggested that among the physicochemical factors relating to elevation, pH (pseudo-F = 3.3, p < 0.005), MAT (pseudo-F = 2.2, p < 0.005), TN (pseudo-F = 1.9, p < 0.005) and C/N ratio (pseudo-F = 1.4, p < 0.005) were the significant contributors to variability of nematode communities on Mt. Norikura (Fig. 7). A total of 66.9% of fitted variation with two axes was explained in an accumulative variance for the interaction between communities and the four significant variables.

Mean elevational range of nematode OTUs. Mean elevational range of nematode OTUs was calculated,

for each level on the mountain. The average species range extent for nematode community members became greater with increasing elevation (p < 0.01) (Fig. 8).

Discussion

We found the diversity of soil nematodes was highly correlated with elevation, and showed a mid-elevation maximum on Mt. Norikura. This finding adds a new example to the often-observed unimodal diversity pattern that is found along elevational gradients, although this in itself does not provide an explanation for its occurrence. There are various candidate hypotheses that could explain this pattern.

As explained in the introduction to this paper, mountains in central Japan tend to have rainfall maxima in their mid elevations [START_REF] Fujimura | Of the climate and weather of Mt. Fuji in Report of the scientific survey of Mt[END_REF] , and it is possible that the extra soil moisture explains the survival and coexistence of more nematode species in the mid elevations. However, the only available data on precipitation are for Fuji, where the maximum is 2,000-2,500 masl, while maximum diversity of nematodes is only found lower down the mountain at around 1000-2000 masl (Fig. 3). This weakens the case that an elevational gradient in precipitation is a key factor bringing about the observed trend in nematode diversity. The 'community overlap' hypothesis, which assumes that two distinct environments on a mountain each have their own distinct set of species, was obviously not effective in explaining the observed diversity pattern of soil nematode community, because our results indicated that the communities of the lowermost and uppermost elevations on Mt. Norikura are not distinct. Instead, nestedness analysis shows that the communities of higher elevations tend to be a subset of the communities of lower elevations.

On the other hand, the midpoint attractor model, by simulation in a null model, successfully reproduced the empirical richness pattern found on Mt. Norikura, indicating gradients of elevational favorability, subject to geometric constraints, may account for the unimodal elevational richness pattern of soil nematodes. MDE has often been put forward to explain the unimodal patterns found for various taxonomic groups along mountain elevation gradients [START_REF] Singh | Elevational Patterns in Archaeal Diversity on Mt. Fuji[END_REF][START_REF] Mccain | The mid-domain effect applied to elevational gradients: species richness of small mammals in Costa Rica[END_REF][START_REF] Brehm | The role of environment and mid-domain effect on moth species richness along a tropical elevational gradient[END_REF][START_REF] Wang | Contrasting patterns in elevational diversity between microorganisms and macroorganisms[END_REF] .

Besides the simulation in null model explaining the ecological diversity pattern, we were also interested in understanding which particular environmental factors may contribute to the deterministic process in terms of 'niche theory' and the empirical effects of legacy of past evolution and the current capacity of the system to sustain diversity. Amongst the physicochemical factors we used in analysis of the elevational pattern of nematode diversity, the inferred MAT and measured TN concentration in soil were the strongest predictors of soil nematode diversity, predicting most of the variation. On Mt. Norikura, the TN concentration is greatest in soils in the mid elevations (correlating with the diversity maximum), which -following global scale patterns -may result from slower decay of leaf litter products at moderately cool temperatures [START_REF] Sarathchandra | Effect of nitrogen and phosphate fertilisers on microbial and nematode diversity in pasture soils[END_REF] , resulting in accumulation of nitrogen in the upper parts of soils, which receive and accumulate partially broken down leaf litter [START_REF] Post | Global patterns of soil nitrogen storage[END_REF] . At lower elevations, with their increased temperatures, more rapid leaf litter decay relative to NPP tends to reduce the N concentrations in the upper soil layers [START_REF] Post | Global patterns of soil nitrogen storage[END_REF] . At higher elevations above 2,500 masl on Norikura, the tree canopy disappears and the vegetation opens up to scrub and then montane tundra [START_REF] Miyajima | Changes with altitude of the stand structure of temperate forests on Mount Norikura, central Japan[END_REF] . This may also be coupled with drier conditions above the main cloud layer [START_REF] Fujimura | Of the climate and weather of Mt. Fuji in Report of the scientific survey of Mt[END_REF] . The lower net primary productivity at the cooler and drier higher elevations, with reduced litter fall more than compensating for slower decay, would explain why TN concentration declines towards upper elevations of Norikura, and with it the diversity of soil nematodes. Although TN was significantly correlated with nematode diversity in this study, greater nitrogen concentration does not always equal greater nematode diversity. Beyond certain levels of surplus nitrogen -for example in farmland ecosystems -nematode diversity decreases [START_REF] Sarathchandra | Effect of nitrogen and phosphate fertilisers on microbial and nematode diversity in pasture soils[END_REF][START_REF] Xu | Response of soil fauna to simulated nitrogen deposition: A nursery experiment in subtropical China[END_REF] . This is thought to be because too much organic nitrogen addition results in very high concentrations of microorganisms, leading to such dense nematode communities that certain nematodes which are adapted of competition between nematode species reduce the overall nematode diversity in that system [START_REF] Sarathchandra | Effect of nitrogen and phosphate fertilisers on microbial and nematode diversity in pasture soils[END_REF] . Soil temperature (for example MAT) has a marked effect on nematode growth and can significantly influence the composition of nematode communities [START_REF] Bakonyi | Soil nematode community structure as affected by temperature and moisture in a temperate semiarid shrubland[END_REF] . Driven by adiabatic cooling, temperature on Mt Norikura is predicted to decline fairly steadily with increasing elevation, starting from the warmest MAT about 10 °C at the lowest elevation site we sampled [START_REF] Takahashi | Effects of climate on the radial growth of tree species in the upper and lower distribution limits of an altitudinal ecotone on Mount Norikura, central Japan[END_REF] . While predicted temperature is closely correlated to nematode diversity on Norikura, there is no obvious reason why it should produce a mid-elevation diversity maximum.

We found that the OTU structure of nematode communities was highly variable between different elevations, displaying a clear 'progression' of community structure towards successively higher elevations, as is often found for other groups of organisms along elevational gradients [START_REF] Heaney | Frontispiece: diversity patterns of small mammals along elevational gradients[END_REF][START_REF] Miyamoto | The mid-domain effect in ectomycorrhizal fungi: range overlap along an elevation gradient on Mount Fuji, Japan[END_REF][START_REF] Patterson | Contrasting patterns of elevational zonation for birds and mammals in the Andes of southeastern Peru[END_REF][START_REF] Singh | Strong elevational trends in soil bacterial community composition on Mt. Halla, South Korea[END_REF] . In this study, physicochemical factors of soils and climate were highly correlated with elevation and showed significant influences on community structure of soil nematodes, suggesting that nematodes on Mt. Norikura have quite finely divided niches in relation to the physicochemical factors that vary with elevation. Consistent with previous studies, we found MAT, TN, C/N ratio and pH were all important environmental variables influencing nematode community structure (Fig. 7) [START_REF] Li | Effects of organic and other management practices on soil nematode communities in tea plantation: a case study in southern China[END_REF][START_REF] Song | Responses of soil nematodes to water and nitrogen additions in an old-field grassland[END_REF][START_REF] Zhao | Soil nematode assemblages in an acid soil as affected by lime application[END_REF] . This suggests that indeed, nematode OTUs tend to have predictable distributions which are linked to identifiable environmental factors.

When we classified the nematodes on Mt. Norikura in terms of family-level feeding guilds, it was striking that the upper elevations had very few plant feeders, and that the OTUs with very broad elevational ranges belonged to families classified as bacterial feeders or omnivores. This trend could be partly a result of the lack of plant cover at the uppermost elevations on Norikura where an open tundra-like vegetation prevails [START_REF] Miyajima | Changes with altitude of the stand structure of temperate forests on Mount Norikura, central Japan[END_REF] . While the fungal-based energy channel is one of the main decomposition process in soil [START_REF] Yeates | Nematodes as soil indicators: functional and biodiversity aspects[END_REF] , fungal feeders were consistently a minor part of the community over the whole elevational gradient on Mt. Norikura. This may suggest that a faster decomposition model based on the bacterial energy channel predominates on Mt. Norikura.

We calculated the mean elevational ranges of soil nematodes occurring at each 200 m sampling level on the mountain, and found that the mean elevational range for soil nematodes was actually greater in nematode assemblages occurring at higher elevations. The general prediction of a Rapoport pattern as applied to elevational gradients is based on the idea that there will be greater environmental variability on long and short timescales on high mountains [START_REF] Stevens | The latitudinal gradient in geographical range: how so many species coexist in the tropics[END_REF][START_REF] Stevens | The elevational gradient in altitudinal range: an extension of Rapoport's latitudinal rule to altitude[END_REF] . Detailed climatic and microclimatic data for Norikura do not exist, but it is clear that direct heating of the sparsely vegetated upper slopes of mountains can produce very high soil temperatures during the day in summer, followed by drastic cooling at night and during winter [START_REF] Miyajima | Changes with altitude of the stand structure of temperate forests on Mount Norikura, central Japan[END_REF][START_REF] Takahashi | Effects of climate on the radial growth of tree species in the upper and lower distribution limits of an altitudinal ecotone on Mount Norikura, central Japan[END_REF] . This may require much more generalized physiological niches for soil nematodes. Additionally, if primary productivity decreases with elevation, only species that have very generalized niches may be able to gain enough food to survive at viable population densities at high elevations. The same nematode species may also be able to survive further down the mountain -by virtue of their generalized feeding niches -allowing them to have a wide elevational range. However, the climatic variability hypothesis in explaining elevational range size was rarely based on real climatic data and largely speculative. Though our results were consistent with the Rapoport's elevational hypothesis, its ecological implications cannot be further interpreted in this limited context of this study.

In conclusion, our study shows that nematodes exhibit various elevational patterns that have been noted from other groups of organisms, indicating how extensive these patterns are in nature. Notably, they may follow a mid elevational diversity maximum, which has been found extensively for many different groups of organisms, and this maximum may follow the deterministic influence of environmental variables as well as the predictions of 'mid domain effect' incorporating a mid-point attractor. We also found the nematode OTUs in this study have more extensive elevational ranges prevailing towards higher elevations, which may provide examples to Rapoport's elevational hypothesis. Further studies are necessary to validate how general these patterns are for nematodes in mountain environments. The modification of the 18s primer for nematodes to use with MiSeq offers a considerable improvement over 454, with much larger yields of sequence data.

Methods

Description of sampling sites. Our study site was Mt. Norikura, an extinct volcano in central Japan. Mt. Norikura reaches 3,026 masl and has a cool temperate monsoon climate at its base (MAT of 8.5 °C at 1,000 masl) with abundant year-round precipitation (2,206 mm at 1,000 masl) concentrated in summer [START_REF] Miyajima | Changes with altitude of the stand structure of temperate forests on Mount Norikura, central Japan[END_REF][START_REF] Takahashi | Effects of climate on the radial growth of tree species in the upper and lower distribution limits of an altitudinal ecotone on Mount Norikura, central Japan[END_REF] . Its surface covering of soils is uniformly derived from andesitic ash deposited over underlying lavas in a series of large eruptions around 15,000 years ago [START_REF] Miyajima | Changes with altitude of the stand structure of temperate forests on Mount Norikura, central Japan[END_REF][START_REF] Takahashi | Effects of climate on the radial growth of tree species in the upper and lower distribution limits of an altitudinal ecotone on Mount Norikura, central Japan[END_REF] . Mt. Norikura has since then undergone natural ecological succession to give a series of vegetation zones: a montane deciduous broad-leaved forest zone between 800 and 1,600 masl; a subalpine coniferous forest zone between 1,600 and 2,500 masl; an alpine dwarf pine Pinus pumila scrub zone between 2,500 and 2,700 masl: and an open herbaceous cushion plant tundra up to about 3,000 m [START_REF] Miyajima | Changes with altitude of the stand structure of temperate forests on Mount Norikura, central Japan[END_REF] .

Sampling and gathering soil nematodes. Sampling was carried out over ten days from late July to early August 2014. In total, 55 samples were collected from 11 elevations on the mountain, along elevational isoclines separated by ~200 m of elevation. At each elevational level, five separate samples spaced 100 m apart in a line were taken. Within each individual sample, five soil cores were combined to make a composite sample. The soil cores were 10 cm in diameter and 10 cm deep and were taken at four corners and central point of the quadrat of each sample (Supplementary Fig. 3). Samples were transported to the laboratory at Shinsu University within 5 hours of being gathered. At the laboratory, the contents of each bag were gently mixed, then gently passed through a 5mm sieve.

Nematodes were collected using a modified Baermann funnel method, with nematodes travelling downwards through waterlogged soil and being collected in a tube below [START_REF] Van Bezooijen | Methods and techniques for nematology[END_REF] . Immediately on arrival of the soil sample at the laboratory, 200 g of soil was wrapped loosely with 2 mm medical gauze and placed into a funnel, which was then filled with distilled water at 20 °C, with the apparatus kept in an air conditioned room at 20 °C during extraction. After 24 hours, 30 ml water in the bottom of the funnel, containing both nematodes and soil particles was collected to centrifuge and the sediment was used to extract nematode DNA.

Soil analyses. Soil analysis of each was carried out at Shinsu University, using standard SSSA protocols (Soil Science Society of America). The parameters analyzed were TC, TN, P 2 O 5 , NH 4 -N, NO 3 -N, K, pH, soil texture, MAT and C/N ratio. Total percentage of silt and clay content are used to indicate soil texture in this paper. MAT was calculated by using a mean lapse rate of 0.6 °C /100 masl, so it showed a complete correlation with elevation.

DNA extraction, PCR amplification and sequencing of 18S rRNA gene. DNA was extracted using the MoBio Power Soil DNA isolation kit (MoBio Laboratories, Carlsbad, CA, USA) following manufacturer's instructions, and random empty vials were chosen and run to serve as controls. The isolated DNA was stored at -80 °C until the PCR stage.

Miseq sequencing was used to assess the nematode community through massively parallel sequencing of the 18S rRNA gene. A 350-bp region of the 18S small-subunit rRNA gene was amplified with the primer pairs NF1 and 18Sr2b [START_REF] Porazinska | Evaluating high-throughput sequencing as a method for metagenomic analysis of nematode diversity[END_REF] with adapter sequences for the Illumina MiSeq. Polymerase chain reaction was performed in a 50-μl reaction mixture composed of 1 or 2 μl of DNA extract, 0.4 μM of each primer, 0.2 mM of each dNTP mix, 1X Taq Reaction Buffer and 1.25 U Solg TM Taq DNA Polymerase (SolGent co., Ltd., Korea). PCR conditions were 10 min at 95 °C, followed by 30 cycles of 60 s at 95 °C, 45 s at 50 °C, and 180 s at 72 °C. Final elongation was at 72 °C for 10 min using the C1000 TM thermal cycler (Bio-Rad Laboratories, Inc., Hercules, CA, USA).

Index PCR was performed for the purified PCR products using the Nextera XT Index kit (Illumina, Inc., San Diego, CA, USA). Each of the 50-μl reaction mixtures was composed of 5 μl of each index primer, 2 × PCR Solution Premix Taq TM DNA polymerase (Takara Bio Inc., Otsu, Shiga, Japan), and 5 μl of the purified DNA. PCR conditions were 3 min at 95 °C, followed by 10 cycles of 30 s at 95 °C, 30 s at 55 °C, and 30 s at 72 °C. Final elongation was at 72 °C for 5 min.

After purification by the AMPure XP beads (Beckman Coulter, Inc., Brea, CA, USA), each amplicon was normalized to 4 nM with 10 mM Tris-HCl (pH = 8.5) and pooled with an internal control PhiX (30%). Heat-denatured pooled amplicons were loaded to a v3 600 cycle-kit reagent cartridge (Illumina, Inc.), and 2 × 300 bp paired-end sequencing was performed by the Illumina MiSeq at the Graduate School of Public Health, Seoul National University.

DNA analysis.

The sequences were demultiplexed and trimmed with a read quality score above 20 by the MiSeq Reporter v2.5 (Illumina, Inc.). The sequence data obtained was processed following the Miseq SOP in Mothur [START_REF] Schloss | Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities[END_REF] . Sequences with any ambiguous bases, sequences with more than 8 homopolymers and sequences with lengths less than 200 bp were removed using the screen. Seqs command in Mothur. Putative chimeric sequences were detected and removed via the Chimera Uchime algorithm contained within Mothur in de novo mode. Rare sequences (less than 10 reads) were removed to avoid the risk of including spurious reads generated by sequencing errors [START_REF] Brown | Scraping the bottom of the barrel: are rare high throughput sequences artifacts?[END_REF] .High quality sequences were assigned to OTUs at ≥99% similarity level.

Taxonomic classification of each OTU was obtained by classifying alignments against SILVA Release 115 databases [START_REF] Quast | The SILVA ribosomal RNA gene database project: improved data processing and web-based tools[END_REF] Mid domain effect on soil nematode diversity. Colwell's midpoint attractor model [START_REF] Mark | Fundamentals of Atmospheric Modeling[END_REF] , which simulates the interaction between a unimodal gradient of elevational favorability and the geometric constraints imposed by domain limits, was applied to explain the MDE on the elevational diversity gradient of soil nematode communities. The midpoint attractor model uses two parameters A ([0 < A < 1], the mean location on the unit line domain) and B ([0 < B < 1], the standard deviation around the location) to describe the elevational favorability (attractor).

To model the expected pattern of species richness under the influence of the attractor, each of the empirical ranges in a data set is placed on the domain stochastically, without replacement of 200 simulations, with its midpoint drawn at random from a proposed attractor distribution N (A, B). Algorithm 2 [START_REF] Mark | Fundamentals of Atmospheric Modeling[END_REF] were used for placing ranges within the domain. To enforce the geometric constraint and maintain the empirical range size frequency distribution, the midpoint is sampled from this distribution only over the interval of feasible midpoints, given the size of each range, such that the range does not extend beyond either the lower or upper domain limit (Colwell, 2000; Supplementary Appendix).

Statistical analysis.

We used ANOVA for normal data and Kruskall-Wallis tests for non-normal data to test whether physicochemical factors differed among different elevations. As the midpoint attractor model does not incorporate physicochemical factors into the estimation of the parameters, to evaluate in more detail the effect of the seven remaining physicochemical factors on OTU richness, Shannon Index and Faith's PD, we performed multiple regression analyses. Non-significant predictor variables were removed sequentially until only significant variables were left in the model. Before applying multiple regression to the dataset, we looked for redundant physicochemical factors using the Varclus procedure in the Hmisc package 49 in the R platform. High correlation was found between TN and TC (Spearman's ρ 2 ≥ 0.91), between TN and P 2 O 5 (Spearman's ρ 2 ≥ 0.73) and between TN and soil texture (Spearman's ρ 2 ≥ 0.61; Supplementary Fig. 2). Therefore, we removed TC, P 2 O 5 and soil texture from the analysis, and used the remaining seven physicochemical factors, {i.e. TN, NH 4 -N, NO 3 -N, K, pH, MAT and C/N ratio} for multiple regression analyses.

Diversity indices such as Shannon Index and Faith's PD and OTU richness were calculated using the Mothur platform [START_REF] Schloss | Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities[END_REF] . To assess the best fitting model of correlations between elevation and richness/diversity and environmental variables, linear and polynomial (quadratic) models were tried out using SigmaPlot v 10.0 (Systat Software, San Jose, CA). Model selection was carried out based on adjusted R 2 and RMSE (root mean square error).

We performed a Non-metric Multidimensional Scaling (NMDS) using Primer v6 to explore patterns in community structure using a Bray Curtis similarity matrix. Abundance data of OTUs were square root transformed to calculate the distance between soil nematode communities. Community structure was compared between elevations using ANOSIM with 999 random permutations.

To individually assess the influence of physicochemical factors on community structure of soil nematode, CCA was used in CANOCO ver. 5. Forward selection was used to select significant explanatory variables with 999 permutations and only significant (p < 0.05) variables were included in the models.

Nestedness analysis was performed by BINMATNEST with default input parameters [START_REF] Rodríguez-Gironés | A new algorithm to calculate the nestedness temperature of presence-absence matrices[END_REF] to test whether the community of one set samples tending to be a subset of the community present in another set of samples. The significance value of nestedness was tested using default input parameters and null model 3 which calculates the p-value for row and column totals following Dong et al. [START_REF] Dong | Soil fungal community development in a high Arctic glacier foreland follows a directional replacement model, with a mid-successional diversity maximum[END_REF] . The samples were re-ordered following the packed matrix order which indicates a high-to-low categorized nestedness.

Furthermore, we investigated the relative abundance of the dominant nematode families across all sites. From their sequence reads, nematodes were identified to family level and binned into their respective functional feeding guilds using previous published comprehensive lists [START_REF] Yeates | Feeding habits in soil nematode families and generaan outline for soil ecologists[END_REF] . Four feeding groups were created: bacteria feeding, fungi feeding, plant feeding, and omnivore/predator.

Figure 1 .

 1 Figure 1. Relative abundance (%) of nematode families on different elevational isoclines, based on total sequence reads.

Figure 2 .

 2 Figure 2.Relative abundance (%) of nematode feeding guilds at eleven different elevations. BF, bacteria feeding; FF, fungi feeding; PF, plant feeding; OP, omnivore/predator. Note the sum of the relative abundance of the four feeding groups is less than 100 percent, due to a sizeable proportion of unclassified sequences.

Figure 3 .

 3 Figure 3. Soil nematode diversity measured along the elevational gradient. OTU richness, Shannon index, Faith's PD were used. Linear, quadratic and cubic regression models were fitted to assess the relationship of elevation with diversity indices of the nematode communities. Model selection was carried out based on adjusted R 2 and root mean square error. Only the best fit quadratic regression was shown. Significance levels were less than 0.001, shown as ***. Adj R 2 is coefficient of determination.

Figure 4 .

 4 Figure 4. NMDS of Bray Curtis similarity of overall community structure in relation to elevation.

Figure 5 .

 5 Figure 5. The Bayesian midpoint attractor model applied to the soil nematode data set rescaled to [0, 1] unit line, where 0 represents sea level. (a) Mean species richness (dark blue line) and 95% confidence interval (light blue band) for 200 simulations. Dark green vertical line, midpoint attractor; light green rectangle, standard deviation; blue line, modelled OTU richness; black dots, empirical OTU richness. (b) Midpoint-range plot for the same data. The x-axis is the location of the range midpoint for each OTU on the elevational domain, and the y-axis plots the elevational span of the range (range size). The large triangle sets the geometrically feasible midpoint limits for ranges of a give size. Black and grey points and associated horizontal line segments illustrate the empirical midpoint and range values for the 1,121 OTUs of soil nematodes. Because many OTUs have identical ranges and midpoints in the data set, the shading of each point is proportional to the number of coincident OTU midpoints. The white-to-blue color scale in the 16 small triangles is proportional to the mean number of modelled points falling in each triangle, averaged over the 200 runs of the simulation.

Figure 6 .

 6 Figure 6. Physicochemical factors shown in relation to elevation. Model selection was carried out based on adjusted R 2 and root mean square error and only the significant ones are reported. The values of TC and TN are shown in percentages. Total of percentage silt and clay content are used to indicate soil texture.

Figure 7 .

 7 Figure 7. Canonical Correspondence Analysis ordination plot of soil nematode community structure based on 18S gene OTUs and a vector overlay of the environmental variables. Environmental variables were shown in arrows and only the significant ones were presented. Different colors of symbols denote different elevations.

Figure 8 .

 8 Figure 8. The mean elevational range of high elevation nematodes have broader elevational ranges and therefore broader temperature ranges. High correlations (p < 0.01) were found at all sequence similarities. MAT, mean annual temperature.

Table 1 .

 1 Linear

model regression between environmental variables and diversity indices. Coefficients are shown for significant predictor variables. Significance level is shown at ***p < 0.001; **p < 0.01 and *p < 0.05. The analysis was performed with all measured variables and only the significant ones are shown. Abbreviations: MAT, mean annual temperature; TN, total nitrogen; TC, total carbon.

  using the classify command in Mothur at 80% cutoff with 1000 iterations. The Miseq sequence data used in this study are deposited in the MG-RAST server (Meyer et al., 2008) under project ID 17384 (http://metagenomics.anl.gov/linkin.cgi?project=17384).

Scientific RepoRts | 7: 3028 | DOI:10.1038/s41598-017-03655-3

Acknowledgements

This work was supported by the Polar Academic Program (PAP) grant funded by Korea Polar Research Institute (KOPRI) under the number 0409-20140091. Jonathan Adams was additionally supported for fieldwork by the Foreign Faculty Funding of Seoul National University (3344-20130047). Ke Dong was supported by a scholarship under the project of BK21 PLUS (Brain Korea 21 Program for Leading Universities & Students). We particularly thank Graduate School of Public Health, Seoul National University, for providing us with the sequencing facility. We would like to show our gratitude to Dr. Robert K. Colwell for testing the MDE using midpoint attractor model (24) and sharing his R code with us.

Author Contributions

Itumeleng Moroenyane, Koichi Takahashi and Jonathan Adams contributed to field sampling. Ke Dong, Koichi Takahashi, Naomichi Yamamoto, Choa An and Dorsaf Kerfahi contributed to lab experiments. Ke Dong, Binu Tripathi, Itumeleng Moroenyane and Hyunjun Cho did data analysis. Ke Dong, and Jonathan Adams wrote the main manuscript text. All of the authors reviewed and contributed to the manuscript.

Additional Information

Supplementary information accompanies this paper at doi:10.1038/s41598-017-03655-3

Competing Interests: The authors declare that they have no competing interests.

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


V I E Anireh

E N Osegi
email: nd.osegi@sure-gp.com

HTM-MAT: An online prediction software toolbox based on cortical machine learning algorithm

Keywords: Machine learning, predictions, cortical learning algorithms, hierarchical temporal memory

HTM-MAT is a MATLAB® -based toolbox for implementing cortical learning algorithms (CLA) including related cortical-like algorithms that possesses spatiotemporal properties. CLA is a suite of predictive machine learning algorithms developed by Numenta Inc. and is based on the hierarchical temporal memory (HTM). This paper presents an implementation of HTM-MAT with several illustrative examples including several toy datasets and compared with two sequence learning applications employing state-of-the-art algorithms -the recurrentjs based on the Long Short-Term Memory (LSTM) algorithm and OS-ELM which is based on an online sequential version of the

Extreme Learning Machine. The performance of HTM-MAT using two historical benchmark datasets and one real world dataset is also compared with one of the existing sequence learning applications, the OS-ELM. The results indicate that HTM-MAT predictions are indeed competitive and can outperform OS-ELM in sequential prediction tasks.

Motivation and Significance

Predictive systems are very vital in the monitoring and operation of modern day businesses, infrastructure and equipment's. For instance in the power industry, utility operators need to determine the expected load profile in advance in order to make proper load schedule for the day; bank operators need to determine the likelihood that a fraudulent bank transaction will occur; web system operators need to detect rogue like behaviour in order to avert security breaches within a given web portal. In order to gain insight into the behaviour of the aforementioned tasks, such systems require predictive coding techniques which in turn play a major role in predicting the hidden causes of incoming sensory information [START_REF] Huang | Predictive coding[END_REF]. In recent times, cortical learning algorithms (CLA) have been proposed in [START_REF] Hawkins | Hierarchical temporal memory including HTM cortical learning algorithms[END_REF] and also applied to the solution of several real challenging problems including retina analysis [START_REF] Boone | Image processing and hierarchical temporal memories for automated retina analysis[END_REF], sign language and gaze gesture recognition [START_REF] Rozado | Optimizing hierarchical temporal memory for multivariable time series[END_REF][START_REF] Rozado | Gaze gesture recognition with hierarchical temporal memory networks[END_REF] and weld flaw detection [START_REF] Rodriguez-Cobo | Feasibility study of Hierarchical Temporal Memories applied to welding diagnostics[END_REF] have been investigated with very promising results. This suite of algorithms has very interesting properties and represents a major step in the quest for a more biological plausible neural network for machine learning tasks. However, until this present moment, an open source cortical learning toolbox that can facilitate research and understanding of CLA-like algorithms in a variety of machine learning predictive tasks is lacking in the technical computing community. In this paper we present HTM-MAT, a software tool which allows such predictive tasks to be accomplished. HTM-MAT is a barebones minimalist version of the cortical learning algorithm proposed in [START_REF] Hawkins | Hierarchical temporal memory including HTM cortical learning algorithms[END_REF] that allows regression and classification functions to be implemented in a predictive way. The theory is based on the hypothesis of the use of spatiotemporal sparse memory predictive neural architecture for online machine learning tasks. The primary purpose of developing HTM-MAT is to make readily available a domain friendly, reusable and comprehensive cortical-learning toolbox for machine learning experiments to be performed on both synthetic and real world data in a technical computing environment.

It is intended to serve as a reusable, straightforward and easy to re-invent/reuse machine learning tool for academic researchers, students and machine learning experts. HTM-MAT depends on the Hierarchical Temporal Memory (HTM) which is the underlying theory that defines its core operations.

HTM theoretical background Hierarchical Temporal Memory (HTM) is originally a memory-prediction framework, and an emerging Artificial

Intelligence (AI) approach that is based on the notion that time plays a crucial role in the operation of neural networks [START_REF] Hawkins | Hierarchical temporal memory including HTM cortical learning algorithms[END_REF][START_REF] Hawkins | On intelligence[END_REF][START_REF] Hawkins | Why neurons have thousands of synapses, a theory of sequence memory in neocortex[END_REF][START_REF] Cui | Continuous online sequence learning with an unsupervised neural network model[END_REF]. In its current implementation, it uses the cortical learning algorithms -a suite of algorithms that favors the use of a more biological plausible paradigm. HTM advocates that the neocortex is the primary entity and the main seat of intelligence in the brain.

The core operations that make HTM interestingly unique are centered on the properties of sparse distributed representations (SDRs) being the fundamental data structure in HTM systems [START_REF] Ahmad | Properties of Sparse Distributed Representations and their Application to Hierarchical Temporal Memory[END_REF]. The primary operations in HTM include [START_REF] Cui | Continuous online sequence learning with an unsupervised neural network model[END_REF]:

 The activation and prediction of cell states/columns based on the notion of permanence  Computation of a similarity metric called the "overlap"  A learning rule for updating permanence's

The prediction and activation states at time step, t, are typically computed using Eq. ( 1) and Eq. ( 2) respectively as:

         otherwise D if t d j i d t ij 0 || || 1 1    (1)                  otherwise i t j i and t W j if t j i and t W j if t ij a 0 0 1 1 1 1 1   (2)
where, , is fundamentally referred to as the "overlap".

In order to improve on the HTM predictions, a Hebbian-like learning rule is applied to the dendrite segment:

  d j i t d j i d j i D p A D p D       1  (3)
where, p + and p -are permanence increments and decrements respectively.

The percentage of active columns typically follows a 2% rule for maximum fixed sparse distribution. The assumption that the 2% rule holds and can be determined apriori remains a fundamental issue and yet to be resolved by HTM researchers. More so, it is still not clear if the brain uses this rule for learning sparsely. Nonetheless, the HTM model still works reasonably well in certain sequence learning tasks such as in time series related problems.

Concept of an online sequential cortical learning algorithm Online sequential learning emphasizes the predictions of patterns that show some interesting regularities or properties in data e.g. the prediction of the likelihood of an anomaly in equipment given a history of equipment conditions. HTM-MAT is a sequence learning framework based on cortical learning principles and is generically designed to operate in an online manner.

A typical online learning process for cortical-like algorithms consists of four key stages. First, a dataset containing the source file (sensory input data) is encoded into a suitable (mixed-integer) representation. This is achieved by using a scalar encoder. Second, using the concepts introduced earlier in Section 1.1, a sparse distributed representation (SDR) of the encoder representations is generated using a spatial pooler algorithm. Third, the SDR formed in the spatial pooler stage is further processed to generate a hierarchy of sparse representations using a temporal pooler algorithm. Fourth, a greedy prediction search is performed to extract the possible sparse pattern meeting a pre-specified overlap threshold. Full details of the spatiotemporal processing inherent in cortical-like applications including a comprehensive mathematical treatment of the formulation of sparse distributed structures can be found in [2, 9 and 10]. Note, that in HTM-MAT, we do not consider the expensive method of combinatorics proposed in [START_REF] Ahmad | Properties of Sparse Distributed Representations and their Application to Hierarchical Temporal Memory[END_REF], but rather use a Monte Carlo simulation to evolve a sparse distributed representation. An implementation of HTM-MAT is provided in [START_REF] Osegi | [END_REF].

Related sequential learning algorithms

The online sequential extreme learning machine (OS-ELM) is also an emerging sequence learning framework that is based on the very popular Extreme Learning Machines (ELM).OS-ELM has been shown to be very fast and highly accurate [START_REF] Liang | A fast and accurate online sequential learning algorithm for feedforward networks[END_REF]. However, this model does not fully account for the notion of time and hierarchy in space. It also does not account for sparse distributed representations which are core features of truly biological-plausible machine learning applications. In addition OS-ELM does not have any inherent mechanism for handling categorical inputs and requires supervised operation. Implementations of OS-ELM can be found in [13]. Another sequence learning algorithm, Long Short-Term Memory (LSTM) neural network originally developed in [START_REF] Hochreiter | Long short-term memory[END_REF] represent the state-of-the-art in predictive mining tasks specifically in the area of language modelling and grammar recognition/prediction. It can operate both in supervised, semi-supervised and unsupervised modes. In particular, LSTM networks can be used as a reference benchmark for evaluating the predictive power of other machine learning algorithms or neural networks where the task is to predict in advance, the possible output response state given a history of previous states. LSTM model neural networks have been used in a variety of tasks including language modelling [START_REF] Palangi | Deep sentence embedding using long short-term memory networks: Analysis and application to information retrieval[END_REF][START_REF] Sundermeyer | LSTM Neural Networks for Language Modeling[END_REF][START_REF] Kim | Character-Aware Neural Language Models[END_REF], grammar/word learning tasks [START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF] and speech recognition [START_REF] Graves | Hybrid speech recognition with deep bidirectional LSTM[END_REF][START_REF] Graves | Speech recognition with deep recurrent neural networks[END_REF]. An implementation of LSTM can be found in [START_REF] Karpathy | [END_REF].

HTM-MAT HTM-MAT is a MATLAB ® toolbox for performing machine learning experiments on streaming or sequential data. The primary feature of HTM-MAT is that it uses a sparse mixed-integer representation for encoding data during training and learning. Using a sparse mixed integer representation gives HTM-MAT the advantage of a character-level or symbolic explanation facility i.e. HTM-MAT can generate more meaningful and expressive

results that still follow a sparse distributed representation using the core HTM principles.

Since HTM-MAT uses an unsupervised machine learning algorithm, it can also derive a solution to the problem in a more compact and adaptive manner. Though, HTM-MAT is specifically intended for researchers and students in academic environments, it may be possible to adapt it to industry or commercial applications.HTM-MAT can also be helpful to users new to the use of cortical-like machine learning algorithms for predictive classification tasks.

Software description HTM-MAT exploits the functional programming capabilities of MATLAB ®. HTM-MAT interacts with the user through a set of functional routines typically instantiated in a main functional class. The architectural concept is as shown in Fig.1. Basically, the user interacts with HTM-MAT (top-layer) through a set of functions (middle-layer). These functions include routines for handling spatial and temporal pooling operations as well as for classifying the recognition units. The lower layer (DATASET FILE

) is used to define the dataset needed for the prediction or classification task.

Fig1.Architecture of HTM-MAT application

Full details of the core functions used in HTM-MAT can be found in the user guide that comes with the software package [START_REF] Osegi | [END_REF].

Illustrative Examples

A sample implementation of a prediction task using a synthetic or benchmark dataset is as shown in Table 1. In this task, the goal is to implement a code routine to predict the target value of the last exemplar given a history of previous inputs. This routine consists of the following blocks of code:  HTM-System Parameters: First the data source (dataset) name is defined followed by a definition of several core system parameters such as the number of Monte Carlo (MC) iterations (iters), minimum overlap (min_overlap), permanence threshold (perms_th), the number of considered sequences (seq_size) and the percent overlap adjustment for greedy-max predictions.

 Spatiotemporal processing: This requires a function call to ST_HTM.m functional class using most of the parameters defined in the HTM-System Parameters block. This operation returns the full mean absolute percentage error (MAPE), a mixed-integer representation of the original input data sequence (original_data) and sparse generated data (sparsed_data), the overlap (overlap), temporal sparse data (sparsed_data_T), minimum values obtained meeting the overlap threshold (min_val) and the permanence.

 Error-plots: Plots the MAPE errors generated by the system after the pre-specified number of MC iterations.

 Char-Integer transformation: This transformation converts the mixed-integer character representations of original and sparse generated data into an equivalent numeric (integer-based) representation.This is very important for extracting the predictions (using a prefix-search) and performing the necessary

HTM-MAT

HTM-MAT FUNCTIONS

DATASET FILE 6 classifications.

 Argmax-predictions: This block allows greedy character-wise predictions to be performed on the recognition units. The predictions are stored in a variable called predicted_out.

 Performance metrics: This block computes the performance metrics such as the percentage classification accuracy, mean absolute percentage error (MAPE), and the root-mean squared error (RMSE) given the original (in_p) and sparse recognition (out_s) integer units.

Table1. A sample HTM-MAT implementation data_name = 'heart_data'; % dataset %% HTM system parameters: iters = 50; min_overlap = 2; perms_th = 0.21; desired_localActivity = 2; seq_size = 700; per_adjust = 99; %% Spatiotemporal processing: [MAPE_n,original_data,sparsed_data,overlap,sparsed_data_T,min_val,permanence] = ST_HTM(iters,... min_overlap,perms_th,desired_localActivity,data_name,seq_size); %% Error Plots: plot(MAPE_n); %% Char-Integer Transformation: in_p = original_data.*1; out_s = sparsed_data.*1; %% Argmax-predictions: predicted_out(pred_idx_n,:) = char(out_s(pred_idx_n,:)); %% Performance metrics: %Compute the errors: RMSE = RMSE_func(in_p,out_s)

In Table 2 is shown a comparative prediction of the HTM-MAT with two applications employing state-of-the-art algorithms using four toy datasets; these applications, recurrentjs developed in [START_REF] Karpathy | [END_REF] and OS-ELM applications using three different activation functions [13] have been earlier introduced in sub-section 1.2.1. The toy datasets include a multiplication operation obtained from the very popular times-table used in high-school mathematics, and three fictitious datasets of word representations (applicable to only HTM-MAT and recurrentjs since they both can handle character-level input data). These datasets are included in the HTM-MAT toolbox [START_REF] Osegi | [END_REF]. In the first toy dataset (times_trainv1), the target value of the last exemplar was intentionally set to zero in order to find out if the algorithm can accurately compute the correct answer. Table 3 also shows comparative experimental results of HTM-MAT with the OS-ELM applications using two historical benchmark datasets and one field data for the predictive classification task. The historical datasets can be obtained from [START_REF][END_REF] and consists of a heart dataset (heart_data) for categorizing a heart condition as either present [START_REF] Huang | Predictive coding[END_REF] or absent (0) and a fraud dataset (australian_data) for detecting the presence [START_REF] Huang | Predictive coding[END_REF] or absence of a fraud (0) . The field data represent two-hourly pressure readings obtained from a sensor attached to an oil pipeline where a 1 indicates a likely threat on the facility. These results are interpreted in terms of the root-mean-squared error (RMSE). Simulations were performed for a maximum of 50 MC iterations for HTM-MAT. The results in Table 2 indicate that HTM-MAT is indeed capable of numeric and word predictions that are indeed a reflection of the pattern in the dataset. The results in Table 3 show that HTM-MAT outperforms the OS-ELM on all tasks presented and achieves 100% accuracy on the real world dataset. It is important to note that the OS-ELM parameters was especially hand-tuned to obtain the best possible set of predictions which was found to be achievable when the number of hidden neurons is greater than 100 and the number of initial data and block size is set to 3 and 0 respectively.

Impact

HTM-MAT is an open source MATLAB toolbox that allows cortical learning algorithms to be explored in a technical computing environment. HTM-MAT to the best of our knowledge presents the first attempt to implement a comprehensive toolbox of HTM cortical like algorithms in a technical computing environment such as MATLAB ®.

HTM-MAT is a reusable framework that can easily be adapted to include novel neural-like algorithms that are biologically plausible. This implies that HTM-MAT is indeed a framework for exploring any kind of neural network or technology that is deeply rooted in cortical-like principles. However, this remains an open challenge for neural network experts but a world of possibilities for software developers. It is expected that HTM-MAT will facilitate research in this direction and help advance the state-of-the-art in predictive applications.

Conclusions

This paper presents HTM-MAT, an open source toolbox implementation of the cortical learning algorithm (CLA) which is based on the hierarchical temporal memory (HTM) theory. HTM-MAT is inspired by how the neocortex operates and proposes an online spatiotemporal sparse-distributed representational neural learning structure. This paper also presents illustrative examples using simple tasks on how to use HTM-MAT. It also compares its prediction performance with results reported using similar applications employing state-of-the-art algorithms on three historical benchmark datasets. The results of simulations show that HTM-MAT is indeed capable of precise word/numeric predictions and can give better results than a popular sequence learning application.

  M×N binary matrix representing the permanence of a connected synapse belonging to a dendrite segment d = a cortical segment j i, = cell and column states t  = an M×N binary matrix denoting the activation state of the cortical network  = the cortical segment activation threshold or permanence N = number of cortical columns M = number of neurons (cortical cells) per column t W = percentage of columns with most synaptic inputs The operation

Table2.

  Comparative word/numeric predictions of HTM-MAT, recurrentjs and OS-ELM using toy datasets ELMrbf -OS-ELM using radial basis activation function b OS-ELMsig -OS-ELM using sigmoid activation function c OS-ELMsin -OS-ELM using sine activation function

Table 3

 3 Root-mean squared errors of HTM-MAT and OS-ELM variants using two historical and one real dataset

	Dataset	OS-ELMrbf	OS-ELMsig	OS-ELMsin	HTM-MAT
		(RMSE)	(RMSE)	(RMSE)	(RMSE)
	heart_data	1.5275	0.6734	1.5816	0.2582
	australian_data	0.8236	0.7164	0.8249	0.0381
	pressure_data	0.5075	0.5075	0.5075	0.0000


Longbo Zhang

Chen Yu
email: yuchen@hust.edu.cn

Hai Jin
email: hjin@hust.edu.cn

Dynamic Spray and Wait Routing Protocol for Delay Tolerant Networks

Keywords: Delay Tolerant Network (DTN), mobile networks, epidemic routing, Spray & Wait routing, dynamic control

Delay Tolerant Networks (DTN) is one of the mobile wireless networks that cannot set up the end-to-end communication path between the source and destination nodes pair in most of time. In this paper, we propose an improved routing protocol named as Dynamic Spray & Wait (DS&W), which is based on Spray and Wait routing strategy attempting to gain the higher delivery ratio benefits of replication-based routing as well as the lower resource utilization benefits of forwarding-based routing. But different from S&W, DS&W routing strategy focus on dynamic control to reduce the scale of the messages' flooding by calculating the delivery success probabilities to destinations. DS&W routing protocol can dramatically reduce the overhead ratio in DTNs. Simulation results also evident that DS&W protocol outperforms other existing DTN routing protocols.

Introduction

Delay Tolerant Networks (DTN) [START_REF]Delay tolerant networking research group[END_REF] often refer to sparse mobile ad hoc network, where a routing path does not necessarily exist. In DTNs, both nodes and links may be unreliable. Many emerging communication networks fall into this paradigm including wildlife tracking and habitat monitoring sensor networks (IPN) [START_REF] Juang | Energy-efficient Computing for Wildlife Tracking: Design Tradeoffs and Early Experiences with Zebranet[END_REF], interplanetary networks [START_REF] Burleigh | Delaytolerant Networking: An Approach to Interplanetary Internet[END_REF], nomadic communities networks [START_REF] Udn | Providing Connectivity to the Saami Nomadic Community[END_REF], etc. Many challenges affect the routing in DTNs such as the changing network topology, low delivery ratio and high delay. The problems can be mitigated by using flooding-based scheme or by predicting the network topology [START_REF] Singh | Temporal Behavior Analysis of Mobile ad hoc Network with Different Mobility Patterns[END_REF] to ensure the delivery ratio and delay. Many works shows that flooding-based scheme has a better performance on delivery ratio and delay than predicting scheme, but it will generate numbers of redundant copies and increase the overhead ratio which would have negative impact on the efficiency, if the flooding scale can not be controlled in a suitable scale, it would lead to network congestion. Thus, the efficiency of a flooding-based scheme DTN routing protocol relies essentially on the appropriate scale of message flooding.

In this paper, we propose DS&W, a routing protocol for DTNs, which dynamically control the scale of message flooding by predicting the delivery probabilities. The contributions presented in this paper are: we propose a suitable method to control the scale of message flooding; based on the dynamic control, we propose DS&W, a novel routing protocol which relies on the control of the scale of flooding. With DS&W protocol, the scale of message flooding is dynamically changeable, and different from other flooding-based scheme routing protocols, it could reduce the overhead ratio with little effect on delivery ratio and delay.

The paper is organized as follows. Section 2 presents the state of the art for DTN routing and the traditional method of the message flooding. Section 3 presents the dynamic control and DS&W protocol. Section 4 provides the simulation results of DS&W routing and related discussion. Finally, section 5 concludes the paper.

Related Works

Several routing protocols have been proposed for DTNs. These protocols can be classified into two categories: single-copy routing protocol [START_REF] Spyropoulos | Single-copy Routing in Intermittently Connected Mobile Networks[END_REF] and multi-copy routing protocol [START_REF] Spyropoulos | Multiple-copy Routing in Intermittently Connected Mobile Networks[END_REF]. With the single-copy routing protocols, only one copy of the message is transported by the topology knowledge of the network. In order to increase the delivery success probabilities, the topology knowledge should be predicted by historical encounters between nodes, location information or others [START_REF] Chen | Enabling Disconnected Transitive Communication in Mobile ad hoc Networks[END_REF]. The Direct Delivery protocol and the First Contact protocol are the use of this strategy. With the multi copy routing, the node carrying the message sends a copy to each encountered node. This is repeated until the destination receives the message. In this case, the contacts are assumed to be totally opportunistic. PRoPHET protocol [START_REF] Lindgren | Probabilistic Routing in Intermittently Connected Networks[END_REF], Epidemic Routing protocol [START_REF] Vahdat | Epidemic Routing for Partially Connected ad hoc Networks[END_REF], and Spray & Wait routing protocol [START_REF] Spyropoulos | Spray and Wait: An Efficient Routing Scheme for Intermittently Connected Mobile Networks[END_REF] are the use of this strategy. Due to the low delivery ratio and long delay of the single-copy routing, the multi-copy routing is the mainstream routing protocol.

PRoPHET protocol determines the delivery probability in each node, which is related to the history of the encounters. Transfer is done only for an encountered node, which has a good probability to deliver the message to its final destination. The delivery probability used by each node is recalculated according to three rules:

(a) When the node M encounters another node E, the probability for E is in-

creased: ( ) ( ) ( ) ( ) encounter old old new L E M P E M P E M P × - + = , 1 , , (1)
where L encounter is an initialized constant. (b) The predictabilities for all destinations D other than E are aged:

( ) ( ) k old new D M P D M P γ × = , , (2)
where γ is the aging constant and k is the number of time units that has elapsed since the last aging. (c) Predictabilities are exchanged between M and E and the transitive property of predictability is used to update the predictability of destinations D for which E has a value on the assumption that M is likely to meet E again:

( ) ( ) ( ) ( ) ( ) ( ) β × × × - + = D E P E M P D M P D M P D M P old old new , , , 1 , , (3)
where β is a scaling constant.

In PRoPHET protocol, if a node i encounter another node j, node i will decide whether to send the message to node j based on the delivery probability of them.

Epidemic routing can be described as following: the source node, which has N neighbor nodes, sends N copies of the message to its N neighbor nodes. Afterwards, each neighbor node judges whether the received message copy has been received before or not. If it has, the neighbor node discards the copy; otherwise the neighbor node will keep this message copy and send the same copy to all its neighbor nodes until the message reaches the destination node. In the ideal situation, epidemic routing is considered to be the most effective strategy, but in reality, the overhead ratio is high and it is faced with serious network congestion.

In order to solve these problems, fixing the number of copies is the most common method. Spray & Wait routing protocol uses this solution. In Spray & Wait routing, there are two phases: spray phase and wait phase. In spray phase, the source node of a message initially starts with L copies of this message, which is related with the number of nodes and the requirements for delay. Any node A that has n>1 message copies (source or relayed), when it encounters another node B with no copies, it hands over to

⎣ ⎦ 2 n B
copies and keeps

⎡ ⎤ 2 n B
copies for itself; when A is left with only one copy, it switches to direct transmission. L can be calculated from:

( ) ( ) 1 1 1 2 6 2 . 1 2 2 2 3 3 - = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - + + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + - M M L M M M a L H L H M M π (4)
where M is the number of nodes; a is multiples of the time units;

∑ = = n i n i H 1 1 τ τ
is the n th Harmonic number of order of τ.

The multi-copy routing protocol can achieve higher delivery ratio and lower delay than the single-copy routing protocols, but it will increase the overhead ratio and wastes energy. Although Spray & Wait routing to a certain extent solves the above problems, these problems still exist. Considered Spray & Wait routing, the number of copies L is not appropriate in all cases such as the following situation: in spray phase, a message sends copies to other nodes after a period of time, set the number of nodes L 1 , then it will continue to flooding message to the left node until the number of nodes reaches L; then the routing protocol goes in the wait phase. If the message is successfully forwarded by a node in L 1 , we can consider the spray phase after L 1 is useless.

Dynamic Control and DS&W Routing Protocol

We can find the limitation of Spray & Wait routing is that it has to fix the copies of message to control the overhead ratio; the spray phase will not stop until the message get to the destination or wait phase. The flooding scale L related with the network size can guarantee the maximum delivery delay and stationary message flooding, but it will increase some messages' flooding. For some messages, a smaller L can not only guarantee the delivery ratio and delay, but also decrease the overhead ratio. Since the redundant message flooding in spray phase, the routing efficiency drops. So, if the spray phase can be controlled dynamically, it will have a better performance than Spray & Wait protocol. In order to ensure the delivery ratio and delay, the scale of message flooding should also not be too small. So how to ensure the changeable flooding scale appropriate is very important. Next, we will give our method using the probability of nodes to realize the dynamic control.

Considering PRoPHET routing protocol, it gives a method to evaluate the probability to deliver the message. In the spray phase, we can calculate the total probability of the nodes. When the probability is high enough to ensure the delivery ratio and delay while the spray phase is not end, we can force it end, so that the spray phase can be controlled dynamically. The following section will prove that our assumption is correct, the method of dynamic flooding is effectively. Based on the dynamic flooding, we will give the pseudo code of DS&W routing protocol.

In order to prove our assumption, we do the following works. For the delay, we divided it into two parts: T s and T w , means the spray phase and wait phase. It is certain that if we reduce the number of copies T s will decrease and T w will increase. Assume that the node can meet another node after t unit , and the node's average delivery probability as p . In spray phase, expect of T s is:

unit L s t T 2 log =
. In wait phase, all nodes carrying packets is connected with the destination with the probability of

( ) L p - - 1 1
, so expect of T w is:

( ) L unit w p t T - - = 1 1
. The delivery delay is:

( ) L unit unit L p t t T - - + = 1 1 log 2 (5)
If we decrease the number of nodes to λL (λ∈(0, 1)), then in spray phase, expect of T s ' is:

unit L s t T λ 2 log = ′
. In wait phase, when the number of node is L if the transmission node is in the node of λL, the wait phase is

( ) L unit w p t T - - = ′ 1 1
, if the node is not in the nodes of λL, the message will be transported by the λL nodes, then

( ) L unit w p t T λ - - = ′ 1 1
, so the expect of T w is:

( ) ( ) ( ) L unit L unit w p t p t T λ λ λ - - - + - - = ′ 1 1 1 1 1
. The delivery delay is:

( ) ( ) ( ) L unit L unit unit L p t p t t T λ λ λ λ - - - + - - + = ′ 1 1 1 1 1 log 2 (6) ( ) ( ) unit L L t p p T T T ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ - - - - - - - + = - = ∆ λ λ λ λ 1 1 1 1 1 1 log ' 1 2 (7)
Derivative of ∆T, we find that 0

1 < ′ ∆ = λ T and 0 0 > ′ ∆ = λ T
, the function graph is similar to a parabola.

0 1 = ′ ∆ = λ T
, which means there is a λ that makes 0 > ∆ λ T , and T' is smaller than T, and limit of 0 = ∆ λ T is negative infinity. So when λ is close to 1, the deliver delay is shorter than spray & wait routing, and there should be a minimal value of the delay. After that when λ gets smaller, the delay gets larger. By calculating the total probability, we can make the λ in a suitable scope. By the research, we reach the following conclusions: with the dynamic control, the routing reduces the overhead while has little influence on the delivery ratio and delay.

By the calculation and analysis, the routing forwarding algorithm can be described as: (1) maintain the delivery probability of each node which is related to the history of the encounters; (2) send messages to its neighbor node, calculate the total delivery probability when nodes are connection, and update the record of the messages and the number of copies; (3) if the total delivery probability is higher than p max or L is equal to 1, switches to wait phase. Next is the pseudo code for DS&W forwarding algorithm:

Performance Evaluation

To demonstrate and evaluate the performance of DS&W, we use ONE1.4.1 [START_REF]Opportunistic Network Environment[END_REF], an Opportunistic Network Environment simulator which provides a powerful tool for generating mobility traces, running DTN messaging simulations with different routing protocols, and visualizing both simulations interactively in real-time and results after their completion. Table 1 shows the delivery ratio of S&W routing and DS&W routing with different p max . We find that when the total probability is from 0.9 to 1, the delivery ratio is similar to S&W routing, this corresponds with our conclusion: if the total probability is high enough, there is negligible affection on delivery ratio. If the total probability gets lower, the delivery ratio would have slightly bad performance.

Figure 1 and 2 show the performance comparison between DS&W routing and S&W routing. We can find that when p=0.95, the overhead ratio reduces about 30.07%, and the delay reduces about 5%; when p=0.93, the overhead ratio reduces about 33.24%, and the delay is basically unchanged; when p=0.9, the overhead ratio reduces about 38.22%, and the delay increases about 13%. It shows a good performance in reducing the burden of network and the delay corresponds with our conclusion: when λ is close to 1, the deliver delay is shorter than Spray & Wait routing, and there should be a minimal value of the delay, after that when λ gets smaller, the delay gets larger. So, when p max is higher than 0.93, the routing can reduce overhead ratio obviously without influencing the delay, when it is from 0.93-0.9, it will increase the delay, but it is still acceptable. Further, we consider the situation when p max =0.85, overhead ratio reduces about 49.5%, and the delay increases about 42.4%. If the network is troubled by the network congestion and is not stringent on delay, this performance still can be acceptable; but when p get lower to 0.8, the overhead ratio reduces about 58.75%, the delay increases about 61.04%, this is unacceptable. So, if p gets smaller, the routing should be regarded as failed. Through our analysis, we can summarize that when p∈(0.9, 1), the DS&W routing shows a good performance in reducing the burden of network while with little influence to delay; when p∈(0.85, 0.9), the DS&W routing protocol can be suitable for the network that is troubled by the burden of network and not particularly stringent requirements on delay.

Conclusions

In this paper, we propose a new routing protocol, named DS&W, which is suitable for mobile delay tolerant networks. DS&W routing protocol is based on multi-copy routing protocol and predicts probabilistic using history of encounters and transitivity. DS&W uses the total probability of the nodes carrying message to dynamically control the number of copies. In this way, the overhead-ratio can be reduced, and deliv-ery-ratio and delay will not change much. Simulation results show that DS&W routing protocol outperforms Spray & Wait routing protocol.

Fig. 1 .

 1 Fig. 1. Overhead ratio comparison between DS&W and S&W

Fig. 2 .

 2 Fig. 2. Delay comparison between DS&W and S&W

Table 1 .

 1 125 nodes are distributed in Helsinki, capital of Finland, about , the nodes are divided into four groups. Group 1 and 4 simulate the pedestrians' nodes, they all have 10MB memory, 1Km/h move speed, 100m transmit range, but move in different area. Group 2 and 3 simulate the car nodes, they have 15MB memories, 150m transmit range, move only on the rode of city, but have different speeds: 2Km/h and 3Km/h. The flowing figures show the performance comparison of three routing protocols. Delivery-ratio of DS&W routing and S&W routing

	Algorithm: Core of DS&W routing
	Spray phase
	Connect with neighbor
	if (neighbor=destination) transport message; jump end;
	else transport message;
	total probability=1-(1-total probability)*(1-neighbor. total probability);
	L=L/2;
	if (total probability> max p or L=1) jump wait phase;
	else jump spray phase;
	Wait phase
	Connect with neighbor
	if (neighbor=destination) transport message; jump end;
	else jump wait phase;
	End;

Table 2 .

 2 Comparison of low p-DS&W and S&W

	Routing	DS&W-0.85 DS&W-0.8	S&W
	Overhead-ratio	10.3094	8.4273	20.4340
	Delay (ms)	4267.6299	4809.0593 2996.4775

Acknowledgement The work is supported by National Natural Science Foundation of China (No.61003220), Research Fund for the Doctoral Program of Higher Education of China (No.20090142120025), Fundamental Research Funds for the Central Universities (HUST:2010QN051) and Natural Science Foundation of Hubei Province of China (No.2010CDB02302)