Ruben Van De Walle
email: ruben.vandewalle@hotmail.com

Garben Logghe

Nina Haas

François Massol

Martijn L Vandegehuchte

Dries Bonte

Van De

Orc Walle Ruben

Arthropod food webs predicted from body size ratios are improved by incorporating prey defensive properties

Keywords: Feeding trials, hunting strategy, insects, invertebrate community, marram grass dunes, trophic interactions

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

Introduction

Biodiversity is in a worldwide crisis, with unprecedented declines in the species richness of many taxa [START_REF] Ceballos | Accelerated modern human-induced species losses: Entering the sixth mass extinction[END_REF]. An understudied aspect of this biodiversity crisis is the loss of species interactions [START_REF] Fortuna | Habitat loss and the disassembly of mutalistic networks[END_REF]. These interactions even disappear at faster rates compared to species loss [START_REF] Valiente-Banuet | Beyond species loss: the extinction of ecological interactions in a changing world[END_REF] and therefore contribute substantially to the degradation of ecosystem functions and services [START_REF] Griffiths | The value of trophic interactions for ecosystem function: dung beetle communities influence seed burial and seedling recruitment in tropical forests[END_REF][START_REF] Keyes | An ecological network approach to predict ecosystem service vulnerability to species losses[END_REF]. Not only mutualistic interactions, such as pollination and facilitation [START_REF] Montoya | Emerging perspectives in the restoration of biodiversity-based ecosystem services[END_REF][START_REF] Traveset | Effects of habitat loss on the plantflower visitor network structure of a dune community[END_REF], but also trophic interactions, such as herbivory and predation, can have important effects on ecosystem functioning [START_REF] Lavorel | A novel framework for linking functional diversity of plants with other trophic levels for the quantification of ecosystem services[END_REF][START_REF] Schmitz | Effects of Predator Hunting Mode on Grassland Ecosystem Function[END_REF]. To realistically gauge the effect of disappearing trophic interactions, a precise quantification and understanding of food web structure is necessary [START_REF] Novak | Predicting community responses to perturbations in the face of imperfect knowledge and network complexity[END_REF]. This is, however, challenging for smaller and less studied organisms such as invertebrates, even though they often attain high numbers and diversity in natural systems [START_REF] Gongalsky | Soil macrofauna: Study problems and perspectives[END_REF][START_REF] Traugott | Empirically Characterising Trophic Networks[END_REF].

Body size is generally used to identify potential trophic interactions, with the assumption that predators usually feed on prey that are (slightly) smaller than themselves [START_REF] Gravel | Inferring food web structure from predator-prey body size relationships[END_REF]. Body size further constrains the potential diet of a species between a lower and upper limit [START_REF] Cohen | Body Sizes of Animal Predators and Animal Prey in Food Webs[END_REF][START_REF] Hirt | Rethinking trophic niches: Speed and body mass colimit prey space of mammalian predators[END_REF]. On the one hand, prey that are much larger than the predator will be very difficult to handle and subdue, with pack-hunting and host-parasite systems as notable exceptions [START_REF] Hirt | Rethinking trophic niches: Speed and body mass colimit prey space of mammalian predators[END_REF][START_REF] Kalinkat | Body masses, functional responses and predator-prey stability[END_REF]. Very small prey items, on the other hand, would not yield sufficient energy to compensate for the invested effort, except if harvesting can be done efficiently [START_REF] Naisbit | Phylogeny versus body size as determinants of food web structure[END_REF]. However, ratios between predator and prey body sizes can vary greatly between and within taxonomic entities [START_REF] Brose | CONSUMER-RESOURCE BODY-SIZE RELATIONSHIPS IN NATURAL FOOD WEBS[END_REF][START_REF] Naisbit | Phylogenetic signal in predator-prey bodysize relationships[END_REF].

Empirical evidence on predator-prey size ratios is growing fast. Nevertheless, data is still scarce for terrestrial invertebrates [START_REF] Gongalsky | Soil macrofauna: Study problems and perspectives[END_REF][START_REF] Traugott | Empirically Characterising Trophic Networks[END_REF] despite mentionable efforts [START_REF] Brose | Foraging theory predicts predator-prey energy fluxes[END_REF][START_REF] Eitzinger | Testing the validity of functional response models using molecular gut content analysis for prey choice in soil predators[END_REF][START_REF] Miller-Ter Kuile | Predator-prey interactions of terrestrial invertebrates are determined by predator body size and species identity[END_REF]. Furthermore, a positive relationship between body size and trophic level is usually very clear in both marine and freshwater environments 13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License [START_REF] Boukal | Trait-and size-based descriptions of trophic links in freshwater food webs: current status and perspectives[END_REF][START_REF] Sholto-Douglas | 13C/12C and 15N/14N isotope ratios in the Southern Benguela Ecosystem: indicators of food web relationships among different size-classes of plankton and pelagic fish; differences between fish muscle and bone collagen tissues[END_REF], but less pronounced in terrestrial ecosystems [START_REF] Brose | CONSUMER-RESOURCE BODY-SIZE RELATIONSHIPS IN NATURAL FOOD WEBS[END_REF][START_REF] Potapov | Trophic position of consumers and size structure of food webs across aquatic and terrestrial ecosystems[END_REF]. The difference is due to inherent dissimilarities between aquatic and terrestrial ecosystems: first, the range of body sizes of the main terrestrial and aquatic primary producers (vascular plants vs. unicellular algae respectively) is different [START_REF] Brose | CONSUMER-RESOURCE BODY-SIZE RELATIONSHIPS IN NATURAL FOOD WEBS[END_REF][START_REF] Perkins | Energetic equivalence underpins the size structure of tree and phytoplankton communities[END_REF][START_REF] Potapov | Trophic position of consumers and size structure of food webs across aquatic and terrestrial ecosystems[END_REF][START_REF] Shurin | All wet or dried up? Real differences between aquatic and terrestrial food webs[END_REF]; second, many terrestrial invertebrate predators have developed strategies to handle larger prey species by use of toxins or suctorial mouth parts [START_REF] Brose | CONSUMER-RESOURCE BODY-SIZE RELATIONSHIPS IN NATURAL FOOD WEBS[END_REF]. Thus, whereas some small terrestrial herbivores consume large plants (e.g. most aphids and caterpillars) and some small terrestrial predators can attack relatively large prey (e.g. small spiders killing large insects), most aquatic prey are fed upon by larger predators, in turn fed on by even larger predators and so on [START_REF] Potapov | Trophic position of consumers and size structure of food webs across aquatic and terrestrial ecosystems[END_REF].

Apart from body size differences, hunting strategy is another factor influencing the likelihood of predator-prey interactions [START_REF] Brose | Predator traits determine food-web architecture across ecosystems[END_REF][START_REF] Miller | Fear on the move: Predator hunting mode predicts variation in prey mortality and plasticity in prey spatial response[END_REF][START_REF] Potapov | Feeding habits and multifunctional classification of soil-associated consumers from protists to vertebrates[END_REF]. Ambush predators use the element of surprise to hunt prey that happen to pass by [START_REF] Hirt | Rethinking trophic niches: Speed and body mass colimit prey space of mammalian predators[END_REF][START_REF] Kersch-Becker | The role of spider hunting mode on the strength of spider-plant mutualisms[END_REF][START_REF] Nentwig | A comparison of prey lengths among spiders[END_REF]. Ambush predators are incentivized to attack a larger range of prey species since their overall encounter rates are lower than those of active predators, as the former stay stationary for longer periods of time (invertebrates and ectotherm vertebrates) or occupy a smaller territory (endotherm vertebrates) [START_REF] Hirt | Rethinking trophic niches: Speed and body mass colimit prey space of mammalian predators[END_REF][START_REF] Strand | Host Specificity of Insect Parasitoids and Predators[END_REF]. Active predators, however, confront their prey head-on, a strategy that is more successful with a larger body size [START_REF] Hirt | Rethinking trophic niches: Speed and body mass colimit prey space of mammalian predators[END_REF][START_REF] Nentwig | A comparison of prey lengths among spiders[END_REF].

The outcome of a predation event is not only influenced by the predator, but also by the prey species.

Prey use different mechanisms to try to escape predation, ranging from evasive behaviour (e.g. leafand grasshoppers, gazelles) to the development of defensive properties such as camouflage (e.g. owls, stick insects), hard body armour (e.g. beetles, turtles) or chemical defences (e.g. dart frogs, shield bugs) [START_REF] Jeschke | Animal Prey Defenses[END_REF][START_REF] Peschel | Do oribatid mites live in enemy-free space? Evidence from feeding experiments with the predatory mite Pergamasus septentrionalis[END_REF]. There is considerably more selective pressure to specialize towards particular prey species than there is to develop specific defensive traits against particular predator species, because it is more necessary to avoid overlap of resources than it is to avoid sharing natural enemies [START_REF] Rossberg | Food webs: Experts consuming families of experts[END_REF]. As a consequence, defence mechanisms seem to be largely uniform across each arthropod taxonomic group [START_REF] Witz | Antipredator Mechanisms in Arthropods: A Twenty Year Literature Survey[END_REF]. Some studies even suggest that prey phylogeny is more important than body size when determining predator-prey interactions [START_REF] Naisbit | Phylogeny versus body size as determinants of food web structure[END_REF].

Food webs are important tools to study ecosystem functioning [START_REF] Montoya | Ecological networks and their fragility[END_REF][START_REF] Thébault | Cascading extinctions and ecosystem functioning: contrasting effects of diversity depending on food web structure[END_REF]. Despite some notable exceptions (e.g. [START_REF] Mclaughlin | Temporal Variability in Predator-Prey Relationships of a Forest Floor Food Web[END_REF], empirical studies on terrestrial arthropod food webs are mostly restricted to one-on-one interactions (e.g. among specific pest species and their main predators) or highly simplified food webs from agricultural systems or extreme environments (e.g. [START_REF] Curtsdotter | Ecosystem function in predator-prey food webs-confronting dynamic models with empirical data[END_REF][START_REF] Digel | Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License of forest soil food webs: higher omnivory and more trophic levels[END_REF][START_REF] Sint | Resolving the predator first paradox: Arthropod predator food webs in pioneer sites of glacier forelands[END_REF]. Size-based predictions of predator-prey interactions from soil food webs were shown to provide useful insights into the realised food web structure, but their accuracy can be optimized by including species traits [START_REF] Potapov | Multifunctionality of belowground food webs: resource, size and spatial energy channels[END_REF]. Furthermore, validation of the assumed interactions between species is still needed for such soil food webs [START_REF] Potapov | Multifunctionality of belowground food webs: resource, size and spatial energy channels[END_REF][START_REF] Potapov | Feeding habits and multifunctional classification of soil-associated consumers from protists to vertebrates[END_REF].

In general, food webs can be constructed by compiling fundamental trophic niches, which include all pairwise trophic interactions that a species can potentially establish with other species [START_REF] Torres-Campos | The return of the trophic chain: Fundamental vs. realized interactions in a simple arthropod food web[END_REF]. These niches can be inferred from literature or ratios between putative prey and predator body sizes (e.g. [START_REF] Hines | A meta food web for invertebrate species collected in a European grassland[END_REF]. Not all potential interactions are, however, necessarily realized in a given food web [START_REF] Torres-Campos | The return of the trophic chain: Fundamental vs. realized interactions in a simple arthropod food web[END_REF]. Predator and prey species can for instance be separated in space or time (e.g. through differences in seasonal activity, microhabitat choice or behaviour) or potential prey can have defensive properties that deter predators from attacking [START_REF] Ruxton | Avoiding Attack: The Evolutionary Ecology of Crypsis, Warning Signals and Mimicry[END_REF][START_REF] Torres-Campos | The return of the trophic chain: Fundamental vs. realized interactions in a simple arthropod food web[END_REF]. Additionally, smaller juveniles might eat smaller individuals and even species that are no prey for larger adults [START_REF] Cuff | Density-independent prey choice, taxonomy, life history, and web characteristics determine the diet and biocontrol potential of spiders (Linyphiidae and Lycosidae) in cereal crops[END_REF]. Since average adult body size is often used for the construction of food webs, interactions by juveniles or exceptionally large or small individuals are usually ignored [START_REF] Gongalsky | Soil macrofauna: Study problems and perspectives[END_REF]. Realised food webs can thus be more complex than those based on adult size, which begs the question to what extent a purely size-based approach is adequate to construct food webs in terrestrial ecosystems [START_REF] Petchey | Size, foraging, and food web structure[END_REF]. Food webs can be based on realized interactions through isotope studies, field observations or gut content analyses (mostly through DNA metabarcoding) (e.g. [START_REF] Jacob | The Role of Body Size in Complex Food Webs[END_REF][START_REF] Miller-Ter Kuile | Predator-prey interactions of terrestrial invertebrates are determined by predator body size and species identity[END_REF]. However, these methods have some shortcomings, the most important one being the inability to distinguish predating and scavenging when gut content is used as a proxy for realized interactions [START_REF] Kamenova | Invasions Toolkit: Current Methods for Tracking the Spread and Impact of Invasive Species[END_REF][START_REF] Miller-Ter Kuile | Predator-prey interactions of terrestrial invertebrates are determined by predator body size and species identity[END_REF].

In order to reconstruct a natural but practically tractable food web, we focused on the arthropod community associated with a keystone plant species of coastal sand dunes: marram grass (Calamagrostis arenaria). This grass is an engineering species of coastal foredunes where harsh conditions, such as salt spray, strong winds and extreme absolute temperatures, restrict the number of species that are capable of surviving in this unique ecosystem [START_REF] Mclachlan | Ecology of coastal dune fauna[END_REF]. Nevertheless, species communities in marram grass dunes are much more diverse and complex than those in most agricultural systems, making this ecosystem appropriate for studying the reliability of size-based food web construction [START_REF] Moore | Impact of agricultural practices on soil food web structure: Theory and application[END_REF].

The aim of this study is to assess whether the arthropod food web in marram grass dunes can be accurately predicted through the assessment of body size ratios. This was done by constructing two food webs representing the trophic interactions in marram grass dunes. A theoretical food web was constructed based on body size ratios, expert knowledge and literature data on microhabitat, phenology and trophic interactions. The second food web was derived from experimental feeding trials with arthropod pairs in small vials. We hypothesise that the theoretical food web based on body size ratios is able to predict trophic interactions, but would be improved by incorporating species traits like prey defences and predator hunting strategy.

Material & Methods

Construction of the theoretical food web

To construct the arthropod food web of marram grass dunes, we used data from coastal dune areas along the Channel and the North Sea, covering the North of France, Belgium, the Netherlands and the United Kingdom (Supp. Mat. 1). Arthropods were sampled by sweep-netting and manually searching the base of 638 marram grass tussocks during the summers of 2017-2019. In total more than 15 000 specimens were caught and identified, using a stereomicroscope. The overarching food web was constructed by pooling the data of all sites and countries. Interactions between two species were based on the species' feeding ecology, body size differences between predator and prey, microhabitat (i.e. whether the species reside primarily among the marram grass leaves or on the ground at the base of the tussock), seasonal activity periods (all collected from literature; see Supp. Mat. 2 for a full overview of the used literature) and expert knowledge. This method has been used before (e.g. [START_REF] Hines | A meta food web for invertebrate species collected in a European grassland[END_REF]. We refer to the resulting network as the theoretical food web.

Construction of the empirical food web

This study did not require ethical approval or a permission for fieldwork. (Auchenorrhyncha) and true bugs (Heteroptera), which were provided with small pieces of marram grass to avoid mortality. Water was presented to all individuals in the form of a small ball of wet paper.

Collection of arthropods

During

Pairwise feeding trials

Feeding trials were used to test pairwise interactions between and within species. These trials were conducted in 60 mL vials that were closed off with mesh to allow free air flow. The bottom of each arena was covered in sand originating from coastal dunes to simulate a slightly more natural environment. Neither marram grass, nor any other form of structure was added to exclude any influence these could have on the interaction. Water was provided on a daily basis. Each trial consisted of two randomly chosen individuals interacting with each other. These could also be individuals from the same species but with different sizes, in order to test for cannibalism. Additionally, the pairs of interacting individuals could be any combination of herbivores, omnivores and carnivores. However, two strict herbivores were never paired. Testing of two predators (either omnivores or carnivores) was included to study intra-guild predation. At later stages of the experiment, some specific interactions were tested in order to maximize the number of tested potential links.

After introducing the arthropods to the vial, it was stored in a shady area at room temperature and left alone. The outcome of the trials was checked and recorded daily with minimal interference. Trials were terminated when at least one of the two interacting individuals had died, i.e. there was no fixed time limit for the trials. In case two known carnivores (e.g. two spiders, either belonging to the same or two different species) did not consume each other, the same trial was counted twice: once with individual one as predator and individual two as prey and vice versa. The dead arthropod was stored in 70% ethanol for later identification, whereas the surviving individual, after being deprived of food for at least one day, could be re-used for subsequent trials. At the end of the series of experiments, all tested arthropods were identified with a stereomicroscope to the most detailed taxonomic resolution according to morphology and/or life stage. This resulted in 617 specimens identified to species level, 63 to genus level and 38 to family level. The food web that resulted from pooling the pairwise trophic interactions will be referred to as the experimental food web.

Determining body size

Since it was not possible to determine the dry weight of individuals, we used body length as an alternative measure for body size [START_REF] Ganihar | Biomass estimates of terrestrial arthropods based on body length[END_REF][START_REF] Martin | The biogeography of insects' length-dry mass relationships[END_REF][START_REF] Moretti | Handbook of protocols for standardized measurement of terrestrial invertebrate functional traits[END_REF]. Before the start of each trial, a digital photograph of every individual was taken with a camera mounted on a stereomicroscope and used to measure its body length using ImageJ version 1.52a [START_REF] Schneider | NIH Image to ImageJ: 25 years of Image Analysis[END_REF]. Body length was measured from cephalon to abdomen, excluding appendages such as antennae, chelicerae, legs or wings. The body length across all tested species and individuals ranged from 1.6 to 12.3 mm (Supp. Mat. 4).

Statistical analyses

Data preparation

Predators were categorized as active or ambush hunter. A separate third category for web-building spiders was included because the use of a web possibly enables spiders to catch larger prey [START_REF] Enders | The Influence of Hunting Manner on Prey Size, Particularly in Spiders with Long Attack Distances (Araneidae, Linyphiidae, and Salticidae)[END_REF][START_REF] Nentwig | A comparison of prey lengths among spiders[END_REF] and spiders appeared as outliers in size-based automated food web construction attempts [START_REF] Bohan | Automated discovery of food webs from ecological data using logic-based machine learning[END_REF]. However, since almost all web-building spiders did not construct a web during the feeding trials, this category served more as a correction for the absence of a web, i.e. we expect lower prey size for web-building spiders in the context of standardized experiments in which web-building was not possible. Since closely related species are known to have similar defence mechanisms, order-level taxonomy of the prey was used as a proxy for defensive traits.

Because of large differences in defence mechanism within the order of the Hemiptera, suborders Auchenorrhyncha (escape predation by rapid movement) were used instead. 42 feeding trials were omitted from the data, because predation could not be reliably determined, e.g. it seemed likely that the "predator" was scavenging or an individual escaped. After removing these data points, the data set consisted of 586 trials, of which 170 ended in predation and 416 did not.

Testing pairwise interactions in the lab

To test the link between predation success, predator/prey body size ratio, hunting strategy of the predator (ambush, active, web-building) and prey taxonomy, we used generalized linear mixed models with a logit link function and binomial distribution. The outcome of the feeding trial, i.e. success (1) or failure (0), was used as response variable, while the log-transformed length of the prey (included for the purpose of correction) and the log-transformed predator/prey size ratio together with categorical variables for hunting strategy and prey taxonomy were used as explanatory variables. Orders that included fewer than 15 tested individuals (Blattodea, Collembola, Dermaptera, Diptera, Hymenoptera, Isopoda, Myriapoda, Neuroptera, Opiliones, Orthoptera, Trombidiformes and Sternorrhyncha), were grouped together into a single category due to convergence issues. This artificial group was then used as a baseline for the comparison with the other included groups (Araneae, Auchenorrhyncha, Coleoptera, Heteroptera, Pseudoscorpionida and Psocoptera). The importance of hunting strategy and prey taxonomy for predator-prey interactions was determined by modelling them as fixed effects, with or without an interaction with the predator-prey body size ratio. According to model selection based on the corrected Akaike Information Criterion (AICc), several models could be considered to best describe the data (Supp. Mat. 3). The model with the lowest AICc was used because model averaging does not allow post-hoc comparison for prey taxonomy and hunting strategy.

The model used to analyse the data from the feeding trials was: Individual 'prey ID' was used as a random variable to correct for re-using some prey individuals in different trials. 'Predator ID' was included as a random variable in the initial models but excluded from final models because it explained almost no variation.

Comparing the empirical and theoretical food webs Some species were not represented in both food webs, since the theoretical food web was based on biodiversity data and the experimental food web depended on the sampled species. The exclusion of these species resulted in the two final food webs containing the same 85 species (Fig. 1; see Supp. Mat. 4 for the full species list). The number of interactions in the subset theoretical food web was much higher than in the empirical food web (726 vs. 143 interactions) because not all possible interactions could be tested in the experimental trials. Therefore, only interactions that were tested in the feeding trials were selected and compared between both food webs (Fig. 1). The interactions were classified into four categories of agreement, with the first and second digit indicating whether or not the interaction was observed in the theoretical or the empirical food web respectively: (1-1) interactions that did or (0-0) did not occur in both food webs, (0-1) interactions that were only observed during the feeding trials or (1-0) interactions only described in the theoretical food web. These categories were linked to the taxonomy of the prey and predator species separately via the construction of mosaic plots. The mosaic plot helps visualise the observed and expected interactions more clearly. The relative frequency of the prey taxonomic groups within each category of agreement between the two food webs are represented by the size of the boxes. We use 0-0 and 1-1 to represent the theoretical and experimental food web agreeing on the absence or presence of an interaction, respectively. 0-1 represents a link that is not present in the theoretical food web but happened during the experimental Since we were interested in the presence or absence of interactions between different species, the length of the arthropods was not needed in this analysis, which allowed us to include 12 extra lab trials (for which we could not reliably determine the size of one of the tested individuals due to wrong scale bars on the photographs), resulting in a data set with 598 tested interactions. R version 4.0.4 (R Core Team, 2021) was used to perform all data analyses. Generalized linear mixed models (glmm) were fitted with 'lme4' [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF] and 'lmerTest' [START_REF] Kuznetsova | lmerTest Package: Tests in Linear Mixed Effects Models[END_REF]. The 'anova' function was used to perform likelihood ratio tests between nested models. The 'emmeans' function, with Tukey adjustment of p-values, from 'emmeans' package [START_REF] Lenth | emmeans: Estimated Marginal Means, aka Least-Squares Means[END_REF] was used to assess and test parameter differences between categories of hunting strategies and taxonomic groups.

The 'emtrends' function from the same package was used to test differences in the effect of predatorprey body size ratio between different hunting strategies. Network metrics were calculated using the 'igraph' package [START_REF] Csardi | The igraph software package for complex network research[END_REF]. Mosaic plots were made with the 'mosaicplot' function of the R base package 'graphics' [START_REF] Core | R: A language and environment for statistical computing[END_REF]. Gephi software [START_REF] Bastian | Gephi: An open source software for exploring and manipulating networks[END_REF]

Results

Pairwise interactions

On average, the body size ratio between (potential) predator and prey was larger for successful predation events compared to trials where predation did not happen (Fig. 3). The probability of successful predation depended on the ratio of predator to prey size (LRT, χ 2 1 = 10.27, p = 0.001), but was generally lower when the prey was larger (LRT, χ 2 1 = 9.43, p = 0.002). Overall, the results indicated that individual predators prefer to consume prey that are ca. 20% shorter than themselves (Fig. 2). The probability of predation was further influenced by predator hunting strategy (LRT, χ 2 2 = 13.43, p < 0.001) and prey taxonomy (LRT, χ 2 6 = 110.3, p < 0.001), while the interaction between predator hunting strategy and body size ratio did not significantly increase the goodness-of-fit (LRT, χ 2 2 = 4.37, p = 0.113).

Active (z ratio = 3.182, p = 0.004) and ambush hunters (z ratio = 2.932, p = 0.009) engaged more into predation compared to web-building spiders; however, no difference was found between the first two hunting strategies (z ratio = 0.420, p = 0.907). Nonetheless, differences between active and ambush hunting strategy were found when only spider data was used for the analysis (Supp. Mat. 6). Post-hoc comparisons among taxonomic groups showed that barklice (Psocoptera), leafhoppers (Auchenorrhyncha) and spiders (Araneae) are most likely to be preyed upon when tested during the feeding trials, while beetles (Coleoptera) and pseudoscorpions (Pseudoscorpionida) are least likely to be consumed (Fig. 4, Supp. Mat. 5).

Food web comparison

The theoretical (based on body size ratios and expert knowledge) and experimental food webs (based on the feeding trials) agreed on most interactions (436/598 or 72.9%; 370 interactions were absent and 66 were present in both food webs), although there were still some discrepancies between the two (81 interactions that were present only in the theoretical and 81 interactions only in the empirical food web).

Spiders were captured by more species during the feeding trials than was expected from the theoretical food web. In contrast, beetles were consumed by fewer species than expected (Fig. 5). The vulnerability of leafhoppers, aphids, springtails (Collembola) and barklice were well estimated (Fig. 5).

Finally, Fig. 5 suggests that Hymenoptera and Neuroptera were, like spiders, eaten by more species than expected. However, these results are based on 2 and 1 observation(s) respectively, so they should be interpreted with caution. As predators (Fig. 6), spiders and beetles again show the most remarkable results. For spiders, the theoretical and experimental food webs often disagreed (Fig. 6). This means that spiders were more often able to subdue prey in the feeding trials than expected based on expert knowledge and body size ratios and conversely, some of the prey indicated for spiders in the theoretical food web were not attacked. Unlike the mismatch between the two food webs regarding beetles as prey, expectations for beetles as predators based on the theoretical food web corresponded rather well to the feeding trials (Fig. 6). Finally, pseudoscorpions were less frequently able to kill their prey during the feeding trials than expected from the theoretical food web (Fig. 6).

Discussion

Overall, the theoretical food web based on body size, activity period, microhabitat and expert knowledge converged quite well with the food web based on experimental feeding trials for both predator and prey species. Additionally, the pairwise feeding trials confirmed the general consensus that predators are typically larger than their prey [START_REF] Brose | CONSUMER-RESOURCE BODY-SIZE RELATIONSHIPS IN NATURAL FOOD WEBS[END_REF][START_REF] Cohen | Body Sizes of Animal Predators and Animal Prey in Food Webs[END_REF]. Potential prey with a body length larger than 80% of that of the predator were in most cases not attacked.

Despite these results, several mismatches between the two food webs for some taxonomic groups were highlighted by the mosaic plots.

First, the vulnerability of certain prey taxa such as beetles (Coleoptera) and pseudoscorpions is consistently overestimated in the theoretical food web. This is likely due to their defensive traits, which deter large predators from attacking [START_REF] Ruxton | Avoiding Attack: The Evolutionary Ecology of Crypsis, Warning Signals and Mimicry[END_REF][START_REF] Shinohara | Functional diversity and trade-offs in divergent antipredator morphologies in herbivorous insects[END_REF][START_REF] Wang | Biomechanical Strategies Underlying the Robust Body Armour of an Aposematic Weevil[END_REF]. Our feeding trials provide some indirect evidence for this hypothesis by indicating prey taxonomy as an important predictor of their outcome. While prey taxonomy is linked with a spectrum of traits, we postulate that the most relevant trait for our study is defence strategy. For instance, hard exoskeletons like those of beetles effectively lower the predation probability, even when attacked by a much larger predator [START_REF] Peschel | Do oribatid mites live in enemy-free space? Evidence from feeding experiments with the predatory mite Pergamasus septentrionalis[END_REF][START_REF] Wang | Biomechanical Strategies Underlying the Robust Body Armour of an Aposematic Weevil[END_REF]. Taxonomic groups without any particular defensive structures, such as aphids and springtails, were frequently consumed, further supporting this hypothesis. Apart from their strong exoskeleton, Coleoptera also have low protein content [START_REF] Reeves | Arthropod prey vary among orders in their nutrient and exoskeleton content[END_REF], which could further explain why they were not predated very often. The low predation probability of pseudoscorpion prey might illustrate that, sometimes, offence is the best defence. Pseudoscorpions do not have a hard exoskeleton, but are venomous and can rapidly paralyze or even kill a potential predator [START_REF] Von Reumont | Quo Vadis venomics? A roadmap to neglected venomous invertebrates[END_REF]. Because predators without alternative prey in our experimental setup still refused to consume these well-defended prey, their defences will likely be even more effective in nature. In case of prey species relying more on evasive behaviour, probabilities of predation were likely inflated in our trials since consistent escape was futile. About 60% of the spider trials and 50% of the pseudoscorpion trials that were differently classified by the theoretical food web, had a beetle as prey species, which resulted in lower-than-expected predation rates (Fig. 6). Since spiders are considered generalist predators, predictions for interactions between spiders and beetles were mainly based on size differences between predator and prey, but in reality, it appears difficult for spiders to kill even smaller beetles. For pseudoscorpions, however, interactions with beetles are mentioned in some studies [START_REF] Del-Claro | Ecological and evolutionary pathways of social behavior in Pseudoscorpions (Arachnida: Pseudoscorpiones)[END_REF][START_REF] Harvey | The systematics and biology of pseudoscorpions[END_REF]. However, [START_REF] Harvey | The systematics and biology of pseudoscorpions[END_REF] mentioned beetle larvae rather than adults as prey for pseudoscorpions. In nature, beetle adults are probably mainly predated by larger vertebrate predators [START_REF] O'connor | From carrion to Christmas beetles: the broad dietary niche of the red fox in a hybrid coastal ecosystem in south-eastern Queensland[END_REF][START_REF] Oosten | Comparative Breeding Biology of Three Insectivorous Songbirds in Dutch Dune Grasslands[END_REF][START_REF] Petracci | Diet of Sanderling in Buenos Aires Province , Argentina[END_REF].

Second, spiders both ate more prey and were preyed upon more frequently -i.e. intraguild predation was more prevalent-during the feeding trials than expected based on the theoretical food web. This may be explained by the inclusion of microhabitat, a factor known to affect spider diet [START_REF] Potapov | Feeding habits and multifunctional classification of soil-associated consumers from protists to vertebrates[END_REF], in the theoretical food web. For the feeding trials, arthropod species were matched randomly.

Consequently, two (spider) species living in different microhabitats could be placed together, resulting in more frequent spider-on-spider predation. Given their high protein content and relatively soft exoskeleton [START_REF] Reeves | Arthropod prey vary among orders in their nutrient and exoskeleton content[END_REF], spiders make high risk -high reward prey for other predators.

Targeting smaller spider individuals might help mitigate this risk. Our data support this, since we noticed unpredicted links in the feeding trials between predators that were substantially larger than their prey, contributing to the overall size-conditionality of predator-prey relationships. Since it is known that young spiders disperse more than adults [START_REF] Humphrey | Fluid mechanic constraints on spider ballooning[END_REF], it can be expected that they encounter adult spiders inhabiting other microhabitats, i.e. that these interactions, though rare, do occur in nature.

Our feeding trials limited the natural behaviour of the species and consequently their predation success. For instance, ambush predator species are known to have more successful predation events when facing larger-bodied prey than actively hunting species [spiders [START_REF] Nentwig | A comparison of prey lengths among spiders[END_REF]; Verdeny [START_REF] Verdeny-Vilalta | Foraging mode affects the evolution of egg size in generalist predators embedded in complex food webs[END_REF], mammals [START_REF] Hirt | Rethinking trophic niches: Speed and body mass colimit prey space of mammalian predators[END_REF]]. This hypothesis was supported for spiders, but not across all considered taxa. We speculate that this may be due to the absence of structural complexity within the vials resulting in decreased ambush predation success, since the predator could not really hide. The positive effect of structural complexity on predation success is known from terrestrial systems [START_REF] Donadio | Linking predation risk, ungulate antipredator responses, and patterns of vegetation in the high Andes[END_REF][START_REF] Mullin | The foraging ecology of the gray rat snake (Elaphe obsoleta spiloides). III. Searching for different prey types in structurally varied habitats[END_REF], but it is especially well established in aquatic systems [START_REF] Flynn | Effect of habitat complexity and predatory style on the capture success of fish feeding on aggregated prey[END_REF][START_REF] Schultz | Influence of seagrass meadows on predatorprey habitat segregation in an Adriatic lagoon[END_REF], where predators are even reported to switch from ambush to active hunting according to the structural complexity (Říha et al., 2021). Webbuilding spiders can be considered a specific type of ambush predator, using webs to tackle larger prey [START_REF] Enders | The Influence of Hunting Manner on Prey Size, Particularly in Spiders with Long Attack Distances (Araneidae, Linyphiidae, and Salticidae)[END_REF][START_REF] Nentwig | A comparison of prey lengths among spiders[END_REF]. However, since almost all web-building spiders did not build a web during our experiment, a significantly lower chance of predating any prey species was observed for these spiders. Direct observations of insects caught in the web, or gut content DNA metabarcoding;

Schrojenstein [START_REF] Schrojenstein Lantman | Body size and tree species composition determine variation in prey consumption in a forest-inhabiting generalist predator[END_REF] offer better possibilities to study the prey species diversity of webbuilding spiders, but then have as disadvantage that sizes are difficult to assess after consumption.

Trophic interactions involving prey without structural defence traits were well predicted because, as discussed above, these species could not evade predation. Consuming such prey species is easy if predators are able to find the prey in a natural environment. Our results therefor imply that existing information about size, feeding behaviour, life cycle and niche differentiation can be used to generate a coarse yet useful estimate of the trophic interactions among arthropods. We demonstrated this for terrestrial arthropods in marram grass dunes, but we support the wider applicability of such theorybased food web estimates in similar habitats. Size-structured predation remains an important factor when considering intra-guild predation, and can only be quantified by direct methods such as feeding trials, because indirect methods based on molecular markers cannot correct for intraspecific variation in body size. We therefore encourage the use of this under-rated method to study intraspecific variation within traits influencing trophic interactions; because "Humans are different; why do we think arthropods are the same?" (adapted from [START_REF] Gongalsky | Soil macrofauna: Study problems and perspectives[END_REF].

  the summer of 2019, 718 arthropod individuals belonging to 155 different species were collected in several marram grass dune areas along the Belgian coast: 93 individuals were caught in De Panne (51.0949°N, 2.5635°E; 27/07), 107 in Wenduine (51.2972°N, 3.0717°E; 03/08), 270 in Blankenberge (51.3205°N, 3.1446°E; 06/07, 15/07, 10/08 and 24/08) and 248 in Oostduinkerke (51.1368°N, 2.6997°E; 03/07, 17/08 and 31/08). All arthropods were caught in marram grass tussocks (Calamagrostis arenaria) by hand or aspirator. After sampling, all arthropods were deprived of food for at least 24h to create a baseline, except for aphids (Sternorrhyncha), leafhoppers 13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  Sternorrhyncha (reproduce rapidly to survive), Heteroptera (chemical defence mechanisms) and

  success ~ log(Length prey) + log(Length predator/Length prey) + Hunting strategy + Prey taxonomy + log(Length predator/Length prey)*Hunting strategy Post hoc tests with Tukey adjustments were used to compare different levels of the fixed effects.

  trials, while 1-0 means the opposite. Over-or under-representation of a taxon in one of the categories is determined by the size of standardized residuals(Meyer et al., 2006) and is indicated by the colouring of the boxes: blue indicates over-while red indicates under-representation. The statistical significance of the residuals is eventually indicated by the shading (legend of Fig.5).

Figure 1 :

 1 Figure 1: a) The theoretical food web and b) the empirical food web. Node colours indicate the species' taxonomic group: beetles (purple), spiders (green), true bugs (blue), leaf hoppers (black), flies (orange), ants (pink), cockroaches (beige), other (grey). Interaction arrows point from predator to prey and the colour indicates predator taxonomic group.

Figure 2 :

 2 Figure 2: Overview of (potential) predator and (potential) prey size for all tested pairwise interactions. Colours indicate whether or not predation occurred.

Figure 3 :

 3 Figure 3: Box-and violin plots represent predator-prey body size ratio in relation to the outcome of the feeding trials. Boxes represent median, 1 st and 3 rd quartiles. Whiskers represent 1.5 times the interquartile range.

Figure 4 :

 4 Figure 4: Effect of prey taxonomic group on the relationship between predation probability and predator-prey body size ratio.

Figure 5 :

 5 Figure 5: Mosaic plot representing the relative frequency of the prey taxonomic groups within each category of agreement between the two food webs. 0-0 and 1-1 represent the theoretical and experimental food web agreeing on respectively the absence or presence of an interaction respectively. 0-1 represents a link that is not present in the theoretical food web but happened during the experimental trials, while 1-0 means the opposite. The relative frequency is indicated by the area of the boxes. Blue colouration indicates over-while red indicates under-representation. Example: a blue 0-1 box indicates that the taxon is consistently more eaten during the feeding trials compared to what was expected from the theoretical food web. Shading represents the statistical significance (light colours for α = 0.05 and dark colour for α = 0.0001).

Figure 6 :

 6 Figure 6: Mosaic plot representing the relative frequency of the predator taxonomic groups within each category of agreement between the two food webs. Coding and colours are as in Figure 5.

  13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  was used for network visualization. 13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

13652656, ja, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2656.13905 by Portail BibCNRS INSB, Wiley Online Library on [25/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Acknowledgements

This research was conducted with financial support by the INTERREG 2 seas project Endure. RVDW was funded by the University of Lille. GL was funded by FWO (grantnr: 1130223N). DB was funded by BOF (grantnr: BOF/24J/2021/066). FM would like to thank CNRS for general support.

We thank everyone who has assisted with the fieldwork and/or identification of the arthropods: Paulien Vanhauwere, Laurian Van Maldegem, Gillis Sanctobin, Noëmie Van den Bon, Pieter Vantieghem, Margaux Pottiez, Arthur Van Roey, Elias Van Den Broeck, Tine Cooreman, Jurjen Logghe and the Endure consortium.

Finally, we would like to thank two anonymous reviewers for their constructive and valuable comments on the manuscript.

Data availability

The raw data are available via Zenodo at https://doi.org/10.5281/zenodo.7546960 [START_REF] Van De Walle | Arthropod food webs predicted from body size ratios are improved by incorporating prey defensive properties[END_REF].

Competing interests

The authors have declared that no competing interests exist.

Author contributions

Ruben


G L J A Rikken

Dielectric magneto-chiral anisotropy in triglycine sulfate

Recently we have predicted, solely on the basis of symmetry arguments, the existence of strong magneto-chiral anisotropy in the displacement current in chiral dielectrics and reported its experimental observation in chiral ferroelectrics near their ferroelectric-paraelectric phase transitions (G.L.J.A. Rikken and N. Avarvari, Nat. Comm. 13, 3564 (2022)). Here we present a microscopic model that describes both direct and inverse dielectric magneto-chiral anisotropy (dMChA) in the chiral ferroelectric triglycine sulfate (TGS), based on the specifics of the charge movement in a time-varying electric field. We predict a strong inverse dMChA in TGS that should be within reach of experimental observation and we find good agreement between our model and the experimental observation of the direct dMChA in TGS.

Chirality is vital in many areas of physics, chemistry and biology, where entities exist in two nonsuperimposable forms (enantiomers), one being the mirror image of the other. Chirality corresponds to an absence of inversion symmetry and if time-reversal symmetry is also absent because of a magnetization or an external magnetic field, an entire class of effects called magneto-chiral anisotropy (MChA) becomes allowed. Optical MChA corresponds to a difference in the absorption and refraction of unpolarized light propagating through the chiral medium parallel or anti-parallel to the field [1], [2]. Initially observed in the visible wavelength range [3], [4], [5], its existence was later confirmed across the entire electromagnetic spectrum, from microwaves [6] to X-rays [7]. MChA was further generalized to other transport phenomena [8]. It was experimentally observed in the electrical transport in bismuth helices [8], in carbon nano tubes [9], in bulk organic conductors [10], in metals [11], [12], in superconductors close to the transition temperature [13] and in semiconductors [14] as an electrical resistance  that depends on the handedness of the conductor and on the relative orientation of electrical current I and magnetic field B:

  (B I) =  0 (1 +  B • I) (1)
with  = -  referring to the right-and left-handed enantiomer of the conductor. Electrical MChA should be considered a generalization of chirality-induced spin selectivity (CISS) which involves strong non-reciprocities in spin polarized electronic transport through very thin layers or across interfaces of chiral materials [15]. For a recent review of this rapidly expanding field, see [16]. Another manifestation of MChA was recently observed in the propagation of ultrasound in a chiral crystal [17], further illustrating its universal character. Inverse MChA, a longitudinal magnetization induced in a chiral medium by an unpolarized flux, first proposed in the optical domain [18], and later for an electrical current [19], has also been observed [20]. MChA has become a prominent representative of the wider class of non-reciprocal transport phenomena in broken-symmetry systems, that play an import role in topological quantum systems and in Berry phase physics [21]. For a recent review of MChA, see [22].

Recently we have claimed, using solely symmetry arguments, that MChA should also exist in the displacement current in chiral dielectrics, and have conjectured in particular that this type of MChA should be relatively strong near the ferroelectric-paraelectric phase transition of chiral ferroelectrics [23]. We have named this effect dielectric MChA (dMChA) and have experimentally observed it in the chiral ferroelectrics triglycine sulfate (TGS) and Rochelle salt (RS). Quite surprisingly, dMChA was was found to be orders of magnitude stronger than the any resistive MChA ever reported before [23]. No microscopic mechanism nor theoretical estimate for its magnitude were given in Ref. [23]. In this Letter we will present a microscopic model that describes both direct and inverse dMChA in TGS. This model is found to be in good agreement with the reported experimental values for direct dMChA in TGS, predicts accessible values for inverse dMChA in TGS and provides the basis to understand dMChA in any chiral dielectric.

Insulating dielectrics when submitted to a time-varying electric field E() will still carry a displacement current density J() because of the movement of bound charges. This displacement current density is given by J = Ṗ =  0 χ Ė where P is the polarization density and χ the permittivity. Similar to Eq. 1 we can write a general form for the symmetry-allowed displacement current density in a chiral dielectric subject to a magnetic field up to first order as

J =  0 χ Ė(1 +  0   B • χ Ė) (2)
where   quantifies the dielectric MChA, now expressed in terms of the current density. The main difference with Eq. 1 is that the latter also applies to constant currents whereas the parameter   in Eq. 2 is only meaningful in the presence of time-varying electric fields.

In the case of a periodic driving electric field,   is independent of the field's frequency, apart from any nontrivial time-dependence of the charge dynamics, as is also the case for  in Eq. 1. Eq. 2 shows that the dM-ChA appears at the second harmonic frequency of this driving field. dMChA can therefore be quantified by the current anisotropy ratio

 ≡ ( 2 () - 2 (-))  .
Very large dMChA, of opposite signs for the two handednesses of the crystals, was observed in TGS with a magnetic field parallel to the polar -axis, with values for  up to 10 -3  -1 at a frequency of 50 kHz, which corresponds to values of   up to 3 • 10 -5  2  , at temperatures just below the ferroelectric-paraelectric phase transition at 322 K [23]. Quite surprisingly, this value for   in TGS is several orders of magnitude larger than any value reported before for MChA in (semi)conducting materials [22] which may open a road to practical applications of this novel effect. Note that its existence is solely derived from symmetry arguments and that it should therefore exist in all chiral dielectrics, which represents an enormous materials class, of particular interest in biological and pharmaceutical contexts.

When neglecting spin effects, to first order, the only magnetic field effect on charge movement is through the Lorentz force, wich is only effective for charge movements perpendicular to the magnetic field, the latter being parallel to the driving electric field for dMChA. A chiral medium however can support charge movements perpendicular to these fields, as symmetry arguments tell us that a longitudinal magnetization density will be induced in any chiral medium by a time varying electric field:

M =   E  =   J   0 χ (3)
where   = -  . Eq. 3 is the generalization of inverse electrical MChA [19] to chiral dielectrics. The simplest physical picture behind this formal symmetry argument is to consider the chiral medium as consisting of helicoidal units with their axis parallel to E In this picture, charge movement along E will be intrinsically accompanied by an azimuthal charge movement around E. This azimuthal component in the charge movement around the electric field then generates a longitudinal magnetization. This electrically induced magnetization couples to the external magnetic field, leading to a magneto-chiral energy density difference

∆  = -M • B = -    0 χ J • B (4)
depending on the relative orientation of current and magnetic field and on the handedness of the medium, exactly the characteristics of dMChA. Eq. 4 is general and will apply to all chiral dielectrics. The modern theory for electrical polarization in general, and ferroelectricity in particular, is based on a quantum-mechanical Berry phase formalism [25], and a general theory for dMChA will probably also have to be based on that. However, in molecular dielectrics, charges are quite localized and a classical approach should be adequate. In the following we will first present a specific classical microscopic model for charge movement in TGS that allows to calculate ∆  . Subsequently, we will present a general thermodynamic method based on the Landau-Ginzburg-Devonshire model to translate ∆  into the current anisotropy factor .

The combination of chirality and ferroelectricity is by no means exotic, as 5 of the 10 polar point groups that support ferroelectricity are also chiral. The first ferroelectric ever discovered, Rochelle salt, is chiral, and introducing chirality is one of the appproaches in the search for new molecular ferroelectrics [26]. The most studied and widely used chiral ferroelectric, TGS, is a so-called gyroelectric crystal [27], a chiral ferroelectric where the handedness of its chirality is coupled to the sign of the remanent polarization along the polar axis, in this case the crystalline  axis [24] [28]. This remanent polarization, and thereby the handedness of the crystal, can be reversed by an electric field parallel to this axis, with a strength above the coercive field   . TGS is most accurately described as (NH

+ 3 CH 2 COOH) 2 •NH + 3 CH 2 COO -• SO 2-
4 and is an order-disorder type ferroelectric. It is the NH + 3 -C dipole of one of the two glycine cations that is principally responsible for the ferroelectricity, the projection of this dipole on the crystal -axis having two bistable opposite values, as illustrated in Fig. 1. Switching between the two polarities happens through a rotation of this dipole around the -axis. This charge movement corresponds essentially to a displacement of charge along the -axis, which is to first order not influenced by a magnetic field B along the -axis, and indeed no effect of such a magnetic field on the phase transition was detected [29]. The acompanying chirality related azimuthal charge movement around the polar axis, the existence of which was inferred above from symmetry arguments, can be recognized in the transfer of the two H26 protons between neighboring glycine molecules [30] (See Fig. 1). Such proton tautomerism is quite common in molecular crystals, and has also been identified as the underlying mechanism for ferroelectricity in some crystals [31]. This movement can be approximated as a rotation of the H26 protons in the - plane, over an angle of ap-

   =  2 ∆     =  2 ∆     0   (5)
This leads to the MChA energy density

∆  = -    = - 2 2 ∆  0      ≡      (6)
where  is the H26 proton density (3 • 10 27 m -3 ) and.

  = -  .
In the above, we have assumed that the response to an applied electric field of the NH + 3 -C dipole-H26 protons in a given unit cell is independent of what happens elsewhere in the crystal. This is clearly not the case close to the ferroelectric-paraelectric phase transition, where the charge dynamics becomes correlated over a large number of unit cells. This is most evident from the value of the permittivity close to the transition, which is orders of magnitude larger than what would be expected for an ensemble of independent, thermally disordered electric dipoles. Clearly, close to the transition temperature, a correlation across many unit cells exists in the charge movements and we can quantify this through a correlation number , given by the ratio between the experimentally observed AC polarization and the theoretical polarization of a dipole gas

 ≡  0        (7)
where   = 1 4 Debye is the NH + 3 -C dipole strength. In Ref. [23] typical values of  ≈ 1500 were reported close to T  , from which follows that  ' 10 6 . Much higher values of  have been reported in TGS crystals with lower defect concentrations [29], implying even larger correlation numbers in those crystals. This correlation will also manifest itself in the MChA energy density which becomes ∆  =      ≡   . For the TGS crystals used in [23] we find a value for the inverse dM-ChA parameter  ≈ 35 • 10 -2  . For comparison, the inverse MChA parameter in the chiral semiconductor tellurium is experimentally found to be 6 • 10 -5   [20] and in the chiral metal CrNb 3 S 6 the maximum reported value is 3 • 10 -4   [32]. So the predicted inverse dMChA in TGS is several orders of magnitude bigger than any reported experimental value for inverse MChA in (semi)conductors. This stems mainly from the strong amplifying effect of the charge movement correlation, expressed by the large correlation factor . The experimental observation of inverse dMChA in TGS therefore seems feasible.

The next step will be to calculate the effect of the MChA energy density on the macroscopic polarization response of the medium. For this we will consider the thermodynamics of the ferroelectric-paraelectric phase transition which can be described by the phenomenological Landau-Ginzburg-Devonshire (LGD) model, which defines a Gibbs' free energy density of the form

  (   ) =  0 + 1 2 ( ) 2 + 1 4  4 + 1 6  6 - (8) with ( ) = (  - )
where  is the so-called Curie constant and for a second order transition, like the one in TGS,   0 [START_REF]Physics of Ferroelectrics: A Modern Perspective[END_REF]. We extend this model heuristically through the addition of the correlated MChA energy density defined above:

  =   +   .
Although such an addition is not a priori justified, it is considered valid as long as it respects the symmetry of the problem [34], which is the case here. The inclusion of this time-dependent additional term will be valid as long as the frequency is much lower than the inverse intrinsic relaxation time scale of the material's response, which for TGS is in the sub-microsecond range [35]. We can therefore write the free energy for dielectric MChA as

  =   +     (9)
From the steady state condition  = 0 we obtain

 -   =  +  3 +  5 (10)
We look for an approximate solution of this equation for  through the Ansatz  =   +  () where   is the static remanent polarization, found from   ( = 0) = 0, assuming that  ¿   and retaining only terms up to  2 . A series expansion of the solution of the resulting quadratic equation for  gives

 =  -    + 3 16 ³  -   ´2  2   + + 9 16 ³  -   ´3  3  2  +  (11)
As () =  cos ,  will have components at , 2 3 etc. By collecting the second harmonic terms in Eq. 11 up to the second order, we arrive at

 = (  2 () - 2 (-))    ' 3 2   2  2  ( 12
)
In the LGD model  = ( 0 ) -1 and with the experimental values of   = 10 -2  2 ,  ≈ 1500,  = 3 • 10 5  -1 and  = 10   and the result obtained above from our microscopic model for the inverse dMChA parameter  ≈ 35 • 10 -2  , Eq. 12 predicts  = 4 • 10 -4  -1 which is to be compared with the typical experimental values reported in [23] of  ≈ 5 • 10 -4  -1 . From the different approximations in the derivation of Eq. 12, an uncertainty in the calculated result of 50% can be estimated, and the agreement between calculation and experiment is therefore satisfactory. From this agreement we can therefore conclude a posteriori that our heuristic extension of the LGD model is valid. The good agreement also automatically confirms our calculation of the inverse dMChA parameter of TGS. Note that the term    in Eq. 10 plays the role of a dynamic "magneto-chiral" electric field. By defining the corresponding magneto-chiral polarization   ≡  0  we can rewrite Eq. 12 in a more easily interpretable form as

 •  = 3 2         (13)
Our model therefore identifies the ratios between three material factors that determine the direct dMChA strength: the ratio between the electrically induced polarization and the remanent polarization and the ratio between the magneto-chiral polarization and the remanent polarization. Although   ¿   , the collective charge dynamics close to the transition results in a very large value of , and thereby of   , which means that Eq. 13 can still result in significant values of . Values for  (1) of 10 5 , up to frequencies of 100 MHz have been reported for TGS [36], suggesting that even larger  values than those reported in [23] are possible in this material. This would open the door to practical applications of this effect, for instance the possibility of non-destructive readout of a ferroelectric memory. Whereas the microscopic part of our model is specific for TGS, the basic constituents, expressed by Eqs.4, 10 and 13 will apply to all chiral ferroelectrics, with material specific values for   . This is supported by the observation in Ref. [23] of dMChA in another chiral ferroelectric, Rochelle salt, albeit one order of magnitude weaker than in TGS. TGS is a ferroelectric of the order-disorder type, in which one or more electric dipoles fluctuate in the paraelectric phase, and become collectively ordered in the ferroelectric phase, generating a macroscopic polarization. As we outlined above, the electrically induced magnetization is not directly linked to the alignement of these dipoles, but results from other, secondary, charge rearrangements in the ferroelectric phase. Although the approach of our calculation is of general validity, it is difficult to quantitatively generalize it to an arbitrary orderdisorder ferroelectric as the secondary charge rearrangement will subtly depend on the details of the crystal structure. RS belongs to another classe of ferroelectrics called displacive, where, at least partially, the macroscopic polarization is induced by the linear displacement of charged entities [37]. As for the case of order-disorder ferroelectrics, it is not this primary source of macroscopic polarization that is responsible for dMChA, but other, more subtle, secondary charge rearrangements will be at its origin. Calculations of the strength of dMChA in Rochelle salt and in other chiral displacive ferroelectrics will therefore also have to be made on a case-by-case basis, depending on the details of the charge movement in those materials. Such moving charges may be electrons, protons or even larger charged entities as in the case of RS, where the tartrate ions are displaced. As large dielectric constants, corresponding to large correlation factors  are ubiquitous in molecular ferroelectrics, we believe that many other chiral ferroelectrics with similarly large dMChA will exist. Although dMChA is not restricted to ferroelectrics, our analysis shows that its strength scales strongly with the value of the permittivity , in particular for our model   ∝  3 , which suggests that experimental observation of dMChA in 'normal' dielectrics with typical dielectric constants in the range of 10-30 will be an experimental challenge. Taking advantage of the linear frequency dependence of the anisotropy factor may be a way to (partially) meet this challenge.

We have limited our model to diamagnetic chiral dielectrics, as to our knowledge, no chiral ferromagnetic ferroelectrics have been reported so far, wich the exception of helimagnetic ferroelectrics, in which the chirality resides in a helicoidal spin structure [38]. In such media even larger dMChA could be expected, the remanent magnetization interacting strongly with the electrically induced magnetization resulting in a large ∆  . More subtle magneto-electric couplings will also contribute to the dMChA. Such materials represent therefore interesting topics for further experimental and theoreti-cal studies.

In summary, we have presented a microscopic model for direct and inverse dielectric MChA in TGS, and a general thermodynamical approach to calculate direct dMChA in any chiral ferroelectric. The predictions of our model for direct dMChA in TGS are in good agreement with recent experiments and its prediction for inverse dMChA is within reach of experimental observation.

FIG. 1 :

 1 FIG. 1: View of the unit cells of the two enantiomers of the TGS crystal, with the corresponding macroscopic polarizations Ps

Acknowledgements This work was supported by the French National Agency for Research project SECRETS (ANR PRC 20-CE06-0023-01).


Kevin Feghoul

Mehdi Deise Santana

El Amrani

Mohamed Daoudi

Ali Amad

Mehdi El Amrani

Spatial-Temporal Graph Transformer for Surgical Skill Assessment in Simulation Sessions

Keywords: Graph Convolutional Networks, Transformer, Surgical Skill Assessment, Hand Skeleton, Simulation, Education

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

Introduction

Surgical skill assessment refers to the process of evaluating and measuring surgeon's technical proficiency and competence in executing surgical procedures. It delivers targeted feedback that enables efficient skill development through the provision of guidance, ultimately resulting in better patient treatment. Traditionally, evaluation has been performed by senior surgeons using both global and task-specific checklists [START_REF] Goh | Global evaluative assessment of robotic skills: validation of a clinical assessment tool to measure robotic surgical skills[END_REF][START_REF] Martin | Objective structured assessment of technical skill (osats) for surgical residents[END_REF]. However, classical surgical skill assessment checklists have several limitations, such as having a restricted scope, being prone to evaluator bias, lacking standardization, being a time-intensive and expensive process. Therefore, the development of automated tools to evaluate surgical skills is of significant interest. Collection and analysis of tool motion or video data can lead to an accurate assessment of the trainee's surgical proficiency. The proficiency can be quantified numerically through metrics such as the average OSATS score [START_REF] Martin | Objective structured assessment of technical skill (osats) for surgical residents[END_REF], or categorized into novice or expert levels, providing a clear and objective evaluation.

The conventional approach for automatically evaluating surgical proficiency relies on analyzing instrument motion, which can be obtained from various data sources such as video object tracking [START_REF] Pérez-Escamirosa | Objective classification of psychomotor laparoscopic skills of surgeons based on three different approaches[END_REF], video spatial-temporal features [START_REF] Zia | Video and accelerometer-based motion analysis for automated surgical skills assessment[END_REF], and robotic kinematics [START_REF] Ismail Fawaz | Evaluating surgical skills from kinematic data using convolutional neural networks[END_REF][START_REF] Wang | Deep learning with convolutional neural network for objective skill evaluation in robot-assisted surgery[END_REF]. Other techniques focus solely on utilizing video data. For instance, Funke et al. [START_REF] Funke | Video-based surgical skill assessment using 3d convolutional neural networks[END_REF] proposed to use a Temporal Segment Network [START_REF] Wang | Temporal segment networks for action recognition in videos[END_REF] by fine-tuning a pre-trained 3D Convolutional Neural Network on a stack of video frames. In [START_REF] Liu | Towards unified surgical skill assessment[END_REF], the authors proposed a unified multi-path framework for automatic video-based surgical skill assessment, taking into account various aspects of surgical skills, such as surgical tool usage, intraoperative event patterns, and other skill proxies. To capture the relationships between these factors, a path dependency module has been specially designed.

In recent years, Graph Convolutional Networks (GCNs) have become the de facto choice for modeling relational data due to their ability to capture both the local and global structure of graphs. This has resulted in GCNs achieving stateof-the-art performance in various tasks related to spatial-temporal data [START_REF] Guo | Attention based spatial-temporal graph convolutional networks for traffic flow forecasting[END_REF][START_REF] Slama | Str-gcn: Dual spatial graph convolutional network and transformer graph encoder for 3d hand gesture recognition[END_REF][START_REF] Yu | Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting[END_REF]. Similarly, Transformers [START_REF] Vaswani | Attention is all you need[END_REF] have revolutionized the field of natural language processing and have become the go-to method for various natural language processing (NLP) tasks. In addition to language-related applications, the Transformer architecture has also been applied to tasks beyond NLP, such as skeleton-based action recognition, and has produced outstanding results, as demonstrated in studies such as [START_REF] Mazzia | Action transformer: A self-attention model for short-time pose-based human action recognition[END_REF][START_REF] Plizzari | Spatial temporal transformer network for skeleton-based action recognition[END_REF][START_REF] Zhang | Stst: Spatial-temporal specialized transformer for skeleton-based action recognition[END_REF].

In this study, we explored the potential of using hand skeleton sequences for surgical skill assessment. Our framework offers several advantages, including (1) being lighter and easier to train than models that process entire video sequences, and (2) providing an affordable alternative to expensive robotic surgical systems that can provide kinematics data, since the hand skeleton can be extracted from affordable mobile phones. Additionally, hand skeleton detection is performed in real-time, which ensures its practicality and suitability for use in real-world scenarios. As far as our knowledge extends, this is the first attempt to use hand skeleton sequences for evaluating surgical proficiency. Considering the graph structure of the hand skeleton and the dynamic spatial-temporal patterns in sequences of hand movement, we propose the STGFormer framework that combines the strengths of spectral GCNs for learning spatial-temporal represen-

Spatial-Temporal Graph Transformer -STGFormer

Hand Skeleton Sequence

Transformer Encoder FC

Surgical Skill Classifier

Output probabilities Spatial-Temporal GCN The contributions of this work are twofold and can be summarized as follows:

(1) we propose to use sequences of hand skeleton for the task of surgical skill assessment. This approach offers several advantages, such as being non-invasive, objective, and extensible to operating rooms. Moreover, hand skeletons can be extracted from inexpensive devices, such as a smartphone. By analyzing hand dynamics, practitioners can gain valuable insights into their performance, which can be used for improvement and ultimately lead to better patient treatment;

(2) we developed a new spatial-temporal model that learns the dynamic spatialtemporal correlations of hand skeletons. It consists of a spectral GCN for spatialtemporal feature learning follows by a Transformer encoder for capturing global temporal dependencies. This combination of proven techniques leads to the best prediction performances compared to existing state-of-the-art models.

Proposed Approach

This section introduces our STGFormer framework, which is illustred in Figure 1. The framework consists of two essential components: (1) a spectral GCN responsible for learning spatial-temporal representation from hand skeleton sequences, and (2) a Transformer encoder designed to capture global temporal patterns.

Spectral Graph Convolutional Networks

In order to learn higher-level feature representations, we constructed a spatialtemporal graph and employed a spectral domain GCN.

Graph Construction In this study, we constructed an undirected spatialtemporal graph G = (V, E) to obtain high-level representations of a hand skeleton sequence consisting of N joints over T frames. The set of nodes in the graph is represented by V , while E denotes the set of edges. The construction process is outlined as follows: Nodes: the nodes in the graph consist of all joints in the sequence, expressed as V = {v ti | t = 1, .., T, i = 1, .., N }. Each node v ti is initialized with its 3D coordinate information. In this study, N is equal to 21. Edges: the set of edges E is defined as the union of intra-skeleton connections, E intra , and inter-frame connections, E inter , in the graph, defined as follows:

E intra = {v ti v tj | (i, j) ∈ H, t ∈ {1, .., T }} (1)
E inter = {v ti v (t+1)i | i ∈ {1, .., N }, t ∈ {1, .., T -1}} (2)
In Eq. 2., H represents the set of naturally connected hand joints.

Graph Learning We trained a spectral deep GCNs based on the previously constructed graph G. We define the graph convolution operator as in [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF]:

H (l+1) = σ( D -1 2 A D -1 2 H (l) W (l) ) (3)
where A = A + I n denotes the adjacency matrix of the undirected graph G with inserted self-connections, I n represents the identity matrix, D ii = j A ij is the diagonal degree matrix, W (l) is a learnable weight matrix, and σ(.) an activation function. H l represents the matrix of activations in the l th layer; H 0 = X, where X is the matrix of input node feature.

Transformer Encoder

In order to capture complex temporal patterns in the hand skeleton sequences, we feed the final high-level representation, previously extracted from the GCN, to a Transformer encoder. To be more specific, we concatenate the newly obtained representations of every joint j ti into a vector h 0 t for each frame. This concatenation process is illustrated in Eq. 4. Next, we combine the initial representations h 0 t from all frames into a vector h 0 as depicted in Eq. 5. This vector h 0 serves as the input for a Transformer encoder.

h 0 t = [j t1 , j t2 , .., j tN ] ( 4
)
h 0 = [h 0 1 , h 0 2 , .., h 0 T ] ( 5
)
The Transformer is an advanced neural network architecture that relies on the self-attention mechanism, enabling the model to effectively process input sequences and generate predictions. Unlike traditional recurrent neural networks, which are limited by sequential processing, the Transformer can simultaneously attend to different parts of the input sequence, making it highly efficient at capturing long-term dependencies.

The self-attention mechanism in the Transformer calculates a weighted sum of the input sequence, with the weights being learned during the training process. This allows the model to assign importance to different positions in the sequence, focusing on the most relevant information for prediction. By considering the entire input sequence rather than just past representations, the Transformer can effectively capture contextual information and make accurate predictions.

A crucial component of the Transformer is the Multi-Head Attention (MHA) module. It enhances the model's ability to capture long-range dependencies and enables simultaneous attention across multiple representation subspaces at different positions. The MHA achieves this by utilizing Query-Key-Value (QKV) pairs. Each QKV triple is transformed into separate linear projections, and the scaled dot-product attention mechanism is applied. The scaled dot-product attention can be defined as follows:

Attention(Q, K, V ) = sof tmax( QK T √ d k )V (6)
where 1

√ d k
is used to counteract the vanishing gradient problem cause by the softmax function.

Each head of the MHA module is computed in parallel. The MHA module can be mathematically represented by the following equations:

M ultiHead(Q, K, V ) = Concat(head 1 , ..., head h )W O (7)
with

head i = Attention(QW Q i , KW K i , V W V i ) (8)
where

W Q i ∈ R dm×d k , W K i ∈ R dm×d k , W V i ∈ R dm×dv , W O ∈ R hdv×dm rep
resent the query, key, value, and output projection learnable weight matrices, respectively. h and d m correspond to the number of heads and the output dimension of the encoder block. In this study, we choosed

d k = d v = d m /h.

Surgical Skill Classifier

After the forward pass through k-th Transformer encoder layer, the learned representation h k , as shown in Eq. 9, is utilized as input to a fully connected neural network. This network is responsible for making predictions about whether the hand skeleton sequence is related to a senior surgeon or a surgical resident.

h k = [h k 1 , h k 2 , .., h k T ] (9)
3 Experimental Results

This section presents the dataset collected for the surgical skill assessment task, as well as the results obtained using our proposed approach and several stateof-the-art deep learning-based models.

Dataset

Fig. 2: Illustration of the circular cutting exercise using the VirtaMed simulator.

Data Collection

We gathered data from a total of 16 participants, consisting of 4 experienced surgeons and 12 surgical residents. The participants executed a circular cutting exercise using the VirtaMed medical simulator, as depicted in Figure 2. The first step of the cutting exercise was to use a laparoscope, as illustrated in Figure 3a, to enter the virtual environment and position the view at the correct location. Following that, the participants utilized an atraumatic grasper tool (Figure 3b) to apply tension to the tissue and execute a precise cut along a circular incision between two lines using a pair of scissors (Figure 3b).

For each participant, we recorded their hand movements while they performed the exercise using a smartphone equipped with 4K recording capability. The circular cutting exercise was conducted in a simulated environment at the PRESAGE medical simulation center (Plateforme de Recherche et d'Enseignement par la Simulation pour l'Apprentissage des Attitudes et des GEstes), which is a department affiliated with the Faculty of Medicine at the University of Lille. This simulation center accurately replicates surgical training scenarios, making it an ideal setting for developing surgical skills. It is common for surgical novices to practice on medical simulator, where they perform tasks from curricula such as the Fundamentals of Laparoscopic Surgery (FLS) [START_REF] Peters | Development and validation of a comprehensive program of education and assessment of the basic fundamentals of laparoscopic surgery[END_REF].

FLS Program

The circular cutting exercise is an important part of the training of residents and is included in the FLS program. Initiated in 2004, the FLS program was designed to deliver standardized training for laparoscopic procedures and encompasses theoretical knowledge as well as practical skills. Surgeons often need to complete the FLS program to obtain certification in laparoscopic surgery.

The circular cutting exercise is a crucial component of the FLS program, along with few other simulation exercises, and holds significant importance within the training curriculum. Despite seeming straightforward, it remains an important aspect of the training curriculum. This exercise helps residents develop precise control over laparoscopic instruments, particularly scissors, and improves their hand-eye coordination. It also enables them to understand the tactile feedback and resistance encountered when cutting tissue using these instruments, and enhances their depth perception skills by accurately assessing the distance and thickness of simulated tissue.

Data Preprocessing

We used the method from [START_REF] Zhang | Mediapipe hands: On-device real-time hand tracking[END_REF] to extract the hand skeleton of both hands from the recorded videos of each individual. The hand landmark model outputs a set of 21 3D coordinates for each frame, based on the hand intra connectivity structure, as illustrated in Figure 1. We opted to rely solely on right hand landmarks, as left hand detection was unreliable and played a minor role in this task. Indeed, the right hand was primarily responsible for the cutting, while the left hand primarily held the tissue with limited movement. Afterward, we normalized each hand skeleton sequence by subtracting the coordinate of the first wrist joint (v 00 ) from each joint. Finally, we generated non-overlapping sliding windows of 20 seconds, which correspond approximately to 600 data points. As a result, we have a varying number of data sequences for each subject, which are directly dependent on the time taken to complete the exercise, with a duration of the recordings ranging from 1 minute and 33 seconds to 6 minutes and 17 seconds, with an average duration of 3 minutes and 6 seconds.

Results

Evaluation framework

In line with the JIGSAWS [START_REF] Gao | Jhu-isi gesture and skill assessment working set (jigsaws): A surgical activity dataset for human motion modeling[END_REF] dataset, which is widely used as a benchmark for evaluating surgical skill assessment, our study also takes into account the surgeon's experience as a valuable indicator of surgical proficiency. In our case, given the existence of two distinct groups of practitioners, namely senior surgeons and surgical residents, we formulate the surgical skill assessment as a binary classification task. Our evaluation strategy involved utilizing a subject-independent 6-fold crossvalidation to enhance the robustness of our evaluation. This approach was necessary because the data sequences of hand movements from the same subjects are likely to exhibit correlations. To ensure fairness in distributing the limited number of surgeons in our dataset across each fold, we generated all possible combinations of two surgeons, resulting in a total of six combinations. This ensured that each surgeon had an equal presence in both the training and test sets.

Additionally, we ensured that surgical residents were evenly distributed across the six folds to maintain homogeneity. To assess the performance of our model, we employed accuracy as well as the weighted average F1-score. The inclusion of the F1-score allowed us to account for imbalanced class distributions within our dataset.

Surgical Skill Classification We compared our approach with eight state-ofthe-art models that we re-implemented. Our approach was compared to classical deep learning-based methods such as TCN [START_REF] Bai | An empirical evaluation of generic convolutional and recurrent networks for sequence modeling[END_REF], LSTM [START_REF] Hochreiter | Long short-term memory[END_REF], DeepGRU [START_REF] Maghoumi | Deepgru: Deep gesture recognition utility[END_REF], and Transformer [START_REF] Vaswani | Attention is all you need[END_REF], trained directly on sequences of raw hand landmarks. In addition, we compared our approach with state-of-the-art spatial-temporal graphbased models, including GCN [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF], ST-GCN [START_REF] Yu | Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting[END_REF], and ASTGCN [START_REF] Guo | Attention based spatial-temporal graph convolutional networks for traffic flow forecasting[END_REF]. The ST-GCN consists of multiple spatial-temporal convolutional blocks, each of which includes two temporal gated convolution layers and one spatial graph convolution layer in the center. The ASTGCN consists of multiple blocks, each composed of a spatial-temporal attention mechanism and a spatial-temporal convolution that utilizes graph convolutions to capture spatial patterns and standard convolutions to describe temporal features simultaneously. We also compared our framework with a model trained on handcrafted features, namely the SoCJ descriptor [START_REF] De Smedt | Skeleton-based dynamic hand gesture recognition[END_REF], which extracts a descriptor from the hand skeleton based on its geometric shape. These features are then input into a LSTM model.

In Table 1, we presented the results of our STGFormer model and above mentioned state-of-the-art baselines. Our STGFormer achieves the best performance in terms of both evaluation metrics, achieving 83.29%, and 81.41% in terms of accuracy, and F1-score respectively, as shown in Table 1, which represent an improvement of 1.37% and 1.28% when compared to the best state-of-the-art model.

The SoCJ approach, which involves extracting spatial descriptors, exhibits the lowest F1-score among the evaluated methods. In addition, even when compared to a LSTM model trained directly on raw data, the SoCJ descriptor proves to be inefficient, highlighting the limitations of the descriptor extraction process for our particular task.

As part of our ablation study, we observed that STGFormer outperformed both the GCN and Transformer models by a significant margin, achieving an accuracy and F1-score improvement of at least 1.37% and 1.28% respectively. This outcome clearly demonstrates the effectiveness of combining graph-based and transformer-based approaches in the context of learning surgical skill evaluation. These results highlight the importance of incorporating spatial and temporal information for accurate and robust assessment of surgical skills.

Therefore, based on these findings, we can draw several conclusions regarding the effectiveness of using temporal data either individually or in combination with spatial data. Firstly, the superiority of STGFormer over the GCN and Transformer models suggests that leveraging both spatial and temporal information provides a more comprehensive understanding of surgical skill performance. By capturing the interplay between spatial relationships and temporal dynamics, STGFormer can extract more informative features, leading to improved accuracy and F1-score.

Secondly, the performance gap between STGFormer and the other models implies that solely relying on either spatial or temporal data may not be sufficient for accurate surgical skill assessment. Spatial information alone might not capture the dynamic nature of the surgical procedure, while temporal information alone might lack the contextual understanding provided by spatial relationships. Therefore, combining both spatial and temporal data, as done in STGFormer, proves to be crucial for achieving superior performance for the particular task of surgical skill assessment.

Conclusion

This study demonstrates the feasibility of utilizing hand skeleton sequences for accurate surgical skill assessment. The successful development of automated surgical skill assessment holds significant importance in training aspiring surgeons and enhancing their proficiency in performing safe interventions. In order to achieve this goal, we proposed a novel approach called STGFormer, which effectively captures spatial-temporal correlations and long-range dependencies in the hand skeleton sequences of practitioners as they perform tasks within a simulated environment. Extensive experiments were conducted on a dataset comprising both senior surgeons and surgical residents, and our STGFormer framework achieved an accuracy of 83.29% and a weighted average F1-score of 81.41%. These results strongly support the efficiency of our approach to accurately distinguish between senior surgeons and surgical residents, highlighting its potential as a valuable tool for evaluating surgical skills.

In a future study, we plan to extend our research by investigating a multimodal approach that combines hand skeleton sequences with RGB data to improve the accuracy of surgical skill assessment. This integration aims to leverage the complementary information provided by both modalities, further enhancing the robustness and effectiveness of our assessment framework.

Fig. 1 :

 1 Fig. 1: Illustration of our STGFormer based surgical skill assessment framework, which is composed of two key components: Spatial-Temporal Graph Transformer and Surgical Skill Classifier.

Fig. 3 :

 3 Fig. 3: (a) Laparoscope; (b) Atraumatic Grasper / Scissors.

Table 1 :

 1 Surgical skill assessment: comparison with state-of-the-art methods.

	Method	Acc	F1-score
	SoCJ [2]	80.39	77.55
	TCN [1]	80.08	78.25
	LSTM [7]	81.21	79.36
	DeepGRU [11]	81.42	79.48
	Transformer [18]	80.53	78.19
	GCN [9]	81.92	80.13
	ST-GCN [21]	79.14	79.54
	ASTGCN [6]	79.30	79.49
	STGFormer (ours)	83.29	81.41