C De Bacco

S Franz

D Saad

C H Yeung

Shortest node-disjoint paths on random graphs

A localized method to distribute paths on random graphs is devised, aimed at finding the shortest paths between given source/destination pairs while avoiding path overlaps at nodes. We propose a method based on message-passing techniques to process global information and distribute paths optimally. Statistical properties such as scaling with system size and number of paths, average path-length and the transition to the frustrated regime are analyzed. The performance of the suggested algorithm is evaluated through a comparison against a greedy algorithm.

Introduction

Among the various computationally-hard constraint satisfaction problems, routing and path optimization have attracted particular attention in recent years due to their non-localized nature and interdisciplinary relevance. The node-disjoint path (NDP) problem on graphs studied here, aims at finding a set of paths linking specified pairs of nodes (communications) such that no two paths share a node; the problem is classified among the NP-complete class [START_REF] Karp | Reducibility among combinatorial problems[END_REF] of hard combinatorial problems. This has not only been studied as a purely theoretical problem by mathematicians in the series of graph minors [START_REF] Robertson | Graph minors. xiii. the disjoint paths problem[END_REF] under the name of subgraph homeomorphism problem, but also by practitioners due to its wide applicability to various fields. For instance, in communication systems where the network performance is often strictly related to capacity limits, traffic congestion and the rate of information flow; and in problems of virtual circuit routing where switches located at nodes may become bottlenecks. Moreover, due to their distributive nature NDP is more resilient to failure and represents one aspect of optimal routing where network robustness is the main objective.

One specific communication application where efficient and effective NDP algorithms are essential is in the area of optical networks where transmissions using the same wavelength cannot share the same edge or vertex, hence all communications of the same wavelength must be non-overlapping (disjoint). Consequently, such an algorithm impacts on the achievable network capacity and transmission rate. In this field of routing and wavelength assignment [START_REF] Chatterjee | Review and performance analysis on routing and wavelength assignment approaches for optical networks[END_REF], the objective is to find a routing assignment that minimizes the number of wavelengths used. Different techniques that exploit disjoint paths heuristic algorithms have been proposed to tackle this problem; for instance, greedy algorithms [START_REF] Chen | A new model for optimal routing and wavelength assignment in wavelength division multiplexed optical networks[END_REF][START_REF] Manohar | Routing and wavelength assignment in optical networks from edge disjoint path algorithms[END_REF], approximations based on rounding integer linear programming formulations [START_REF] Banerjee | A practical approach for routing and wavelength assignment in large wavelength-routed optical networks[END_REF][START_REF] Kolliopoulos | Approximating disjoint-path problems using greedy algorithms and packing integer programs[END_REF], post-optimization methods [START_REF] Belgacem | A post-optimization method for the routing and wavelength assignment problem applied to scheduled lightpath demands[END_REF], bin packing algorithms [START_REF] Skorin-Kapov | and wavelength assignment in optical networks using bin based algorithms[END_REF], various heuristic genetic algorithms such as ant colony optimization [START_REF] Blesa | Ant colony optimization for the maximum edge-disjoint paths problem[END_REF] and differential evolution [START_REF] Storn | Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces[END_REF].

Another important application of NDP is in the design of very large system integrated circuits (VLSI), where one searches for non-overlapping wired paths to connect different integrated hardware components, to avoid cross-path interference. Similarly, in wireless ad-hoc communication networks [START_REF] Sumpter | Maximizing number of satisfiable routing requests in static ad hoc networks[END_REF][START_REF] Modiano | Minimum energy disjoint path routing in wireless ad-hoc networks[END_REF]14], where each node can act as a router, path overlaps imply signal interference and low transmission quality, whereas longer paths imply poor signal to noise ratio due to multiple relays and higher transmission power; hence the need to consider both path length and transmission overlaps to be minimized is essential for routing problems. Solutions to the NDP problem also provide fault tolerant routes due to the optimal separation of communication paths all over the network, so that if a node (router) fails, as frequently happens in wireless networks due to the mobility of hosts, only few communications will be affected [START_REF] Li | Stable node-disjoint multipath routing with low overhead in mobile ad hoc networks[END_REF][START_REF] Jain | Exploiting path diversity in the link layer in wireless ad hoc networks[END_REF]. This feature is particularly important when quality of service (QoS) is one of the main requirements in the set up of a communication network, along with the load-balancing feature of NDP that prevents network congestion by establishing non-overlapping routes. This is especially relevant to connection-oriented networks [START_REF] Tanenbaum | Computer Networks 4th Edition[END_REF] that are strongly affected by node failures and congestion [START_REF]Research challenges in qos routing[END_REF].

Practical algorithms for various applications often depend on the specific network topologies considered [START_REF] Aggarwal | Node-disjoint paths on the mesh and a new trade-off in vlsi layout[END_REF] and mostly focus on the optimization version of the problem, i.e. maximizing the number of paths routed [START_REF] Chekuri | Poly-logarithmic approximation for maximum node disjoint paths with constant congestion[END_REF]. The satisfiability version of the problem, i.e. whether all paths can be routed successfully without overlap, is not considered; theoretical studies often give bounds to the achievable approximation instead of providing a practical algorithm for individual instances and fail to calculate path lengths and possible overlaps at the same time as part of the optimization process observables. Given that paths are constrained to be contiguous and interaction between paths is non-localized, a local protocol is insufficient and global optimization is required. The computational complexity is determined by the fact that such a global optimization problem has to consider all variables simultaneously in order to minimize a cost function with non local interactions between variables.

Unlike other constrained satisfaction problems on networks, NDP has received little attention within the statistical physics community. In this paper we consider a random version of NDP on regular graphs (Reg), Erdős Rényi (ER) [START_REF] Erdős | On the evolution of random graphs[END_REF] and a dedicated type of random graph (RER) described in Section 4, with the aim of testing the efficacy of statistical physicsbased methods derived in the context of spin glass theory [START_REF] Mézard | Spin glass theory and beyond[END_REF] such as belief propagation or message-passing (MP) cavity method [START_REF] Mézard | Information, physics, and computation[END_REF][START_REF] Mézard | The cavity method at zero temperature[END_REF] as viable alternatives to greedy algorithms; we also study statistical and scaling properties of quantities of interest as a function of network size and number of paths. We study sparse regular, ER and RER random graphs as they are the most interesting for the problem at hand, but the methodology can be easily extended to accommodate other sparsely connected architectures. Clearly, due to the hard constraint of node disjoint paths, typically no solutions would be found in graphs having a non negligible number of nodes with degree k = 1, 2. Moreover, graphs with a small number of high degree nodes (hubs) or with high modularity measure, such as scale-free or planar graphs, are not interesting for the node-disjoint routing problem since when a paths passes through one of these special nodes it leads directly to graph fragmentation, hence frustration. The situation would be very different for constraints on edges instead, but this variant of the problem is left for future work. Finally, the requirement for the graph to be sparse is suggested by restrictions on the validity of the cavity method which is based on fast decaying correlation functions, i.e. a negligible number of loops in the graph.

Numerical simulations indicate that MP outperforms greedy breadth-first search algorithms not only in finding better solution but also in reaching a higher frustration threshold. Moreover, we find scaling of the expected total length of the NDP as a function of the system size and graph connectivity that goes as M log V V log γ (k-1) with exponent γ that depends on the type of graph, where V is the number of nodes and M the number of paths. We find good agreements between theory and simulation data for graphs of average degrees k = 3, 5, 7 and sizes V = 1000, 2000, 4000, 5000, 10000. Finally, we study statistical properties of physical quantities observed a posteriori, i.e. when a solution is found, such as path length distribution, degree distribution and maximum cluster size for the case of regular graphs.

The reminder of the paper is organized as follows: in Section 2 we will introduce the model used followed by the algorithmic solution in Section 3. Results obtained from numerical studies will be presented in Section 4 followed by conclusions and future research directions in Section 5.

Model

Given an undirected graph (or network) G = (V, E) characterized by V = |V| nodes and E = |E| edges we define a set of M communications C as paths on edges of the graph, each of which originates from a source node S and terminates in a receiver node R. We introduce a variable Λ µ i to characterize each node i ∈ V:

Λ µ i =            +1 if i is a sender for communication µ -1 if i is a receiver for communication µ 0 if i is neither a sender nor a receiver for communication µ (1)
The full node characterization is specified by a vector Λi := (Λ For a given set of M sender-receiver pairs (S µ , R µ ) with µ = 1, . . . , M we address the problem of finding a set of communications that optimize a cost function which penalizes path length and prevents communications overlap (traffic). The state of the network can be specified by introducing a variable I µ i j for each edge (i j) ∈ E and for each communication µ, which specifies whether communication µ passes through edge (i j) and in which direction:

I µ i j =           
+1 if µ passes through (i j) from i to j -1 if µ passes through (i j) from j to i 0 if µ does not pass through (i j)

(2)
Notice that in this formalism I µ i j = -I µ ji . We term these variables currents and define for each edge (i j) a vector Īij := (I 1 i j , . . . , I M i j ) that collects information on all currents involved in that edge. Currents are subject to Kirchhoff law:

∑ j ∈∂ i I µ i j -Λ µ i = 0 ∀µ = 1, . . . , M . (3)
For a given path optimization problem we seek the communication configuration C * that minimizes a cost function c({ Īij }), which penalizes path length and traffic congestion:

c({ Īij }) = ∑ (i j)∈E f (|| Īij ||) (4)
where

f (|| Īi j ||) is a monotonically increasing function of || Īij || := ∑ M µ=1 |I µ i j |
, that penalizes both congestion and path length; where |I µ i j | denotes the absolute value of I µ i j . We would like now to search for approximate solutions to this problem by message-passing equations [START_REF] Mézard | Information, physics, and computation[END_REF]. To derive a distributed algorithm it is useful first to consider tree-like graphs T , for which one can derive exact recursive equations, and later on use these equations as an approximation for arbitrary graphs G. If T is a tree, the removal of any edge (i j) ∈ E divides T in two disjoint subtrees T i and T j (see figure 1). We can define Êij ( Ī) as the optimized cost on T i when the current Ī flows through on the edge (i j); in this way we can write the message sent from node i to his neighbor j when the current Īij flows through the edge (i j) as

E i j ( Īij ) := Êij ( Īij ) + f (|| Īij ||). i j T i T j
Figure 1: Subtree division. T i is the subtree rooted in i when the edge (i j) is removed. Conversely, T j is the subtree rooted in j after the same edge removal.

The messages E i j ( Īij ) admit the min-sum [START_REF] Mézard | Information, physics, and computation[END_REF] recursion relation:

E i j ( Īij ) = min Īki |constraint          ∑ k∈∂i\ j E ki ( Īki )          + f (|| Īij ||) , (5)
where the symbol ∂i stands for the set of neighbors of node i and the constraint is the Kirchoff law (3).

In the following we will use the recursion equation ( 5) on arbitrary random graphs to approximate the constrained minimum of c({ Īij }), the cost defined in equation [START_REF] Chen | A new model for optimal routing and wavelength assignment in wavelength division multiplexed optical networks[END_REF]. Namely:

c * := min Ī 1 E ∑ (i j)∈E { E i j ( Ī) + E ji (-Ī) -f (|| Ī||) } (6)
where the last subtracted term is introduced to avoid double counting the cost of edge (i j).

Unfortunately, the computational complexity of this algorithm is exponential in the number of communications M. In fact, messages E i j ( Īij ) can a priori take 3 M values corresponding to all possible currents passing through a single edge (i j). Therefore, we cannot generally treat even moderately large values of M [START_REF] Yeung | From the physics of interacting polymers to optimizing routes on the london underground[END_REF]. The problem can be simplified if we introduce the hard constraint that paths cannot overlap on nodes (and thus neither on edges). This has the important consequence of reducing the configuration space from 3 M to 2M + 1 and the computational complexity becomes linear in M. This restricted version of the path optimization problem is called the node-disjoint path problem (NDP), as we already mentioned in the introduction, and is the problem we address here. Notice that since we impose the node-disjoint constraint for the communications, then one communication at most flows through the edges, so that

|| Īi j || = ∑ M µ=1 |I µ i j | ∈ {0, 1}
. This corresponds to taking:

f (|| Ī||) =            ∞ if || Ī|| ≥ 2 1 if || Ī|| = 1 0 if || Ī|| = 0 (7)
so that the cost function ( 4) represents indeed the total path length.

In order to solve equation ( 5) iteratively we define a protocol for taking into account only the allowed configurations at each edge given the current value Ī passing through it and Λi at vertex i.

If | Λi | = 0 then: E il ( Īil = 0) = min          ∑ j∈∂i\l E ji ( Ī ji = 0), (8)
min

j 1 , j 2 ∈∂i\l;µ∈M         E j 1 i (I µ j 1 i = +1) + E j 2 i (I µ j 2 i = -1) + ∑ k∈∂i\l, j 1 , j 2 E ki ( Ī ji = 0)                  E il (I µ il = ±1) = min j∈∂i\l          E ji (I µ ji = ±1) + ∑ k∈∂i\l, j E ki ( Īki = 0)          + 1 (9)
If Λ µ i = ±1 then:

E il ( Īil = 0) = min j∈∂i\l          E ji (I µ ji = ∓1) + ∑ k∈∂i\l, j E ki ( Īki = 0)          (10) E ji (I ν ji = ±1) = + ∞ (ν µ) (11)
E ji (I µ ji = ∓1) = + ∞ (12) E ji (I µ ji = ±1) = ∑ j∈∂i\l E ji ( 0) + 1 (13)
The constant +1 that appears equations ( 9) and ( 13) are the costs assigned for a unit of current passing through the considered edge (i.e. f (1) = 1). This cost is the one required for the shortest paths but can be generalized to other arbitrary types of costs. Equation (8) represents the case where i is a transit node and no current passes through edge (i j), then the allowed configurations are that either no currents pass through the remaining neighboring edges (first term inside curly brackets) or one current enters and then exits i through a pair of neighboring edges, all others edges being unused (second term inside brackets). In figure 2 you can see a diagram representing the different allowed configurations for a transit node. Equation (9) represents the case where i is a transit node and the communication µ passes through edge (i j); in this case the only allowed configuration is that where the same communication µ enters/exits from one of the other neighboring edges, all others being unused. Similar considerations are used to formulate the equations [START_REF] Blesa | Ant colony optimization for the maximum edge-disjoint paths problem[END_REF][START_REF] Storn | Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces[END_REF][START_REF] Sumpter | Maximizing number of satisfiable routing requests in static ad hoc networks[END_REF][START_REF] Modiano | Minimum energy disjoint path routing in wireless ad-hoc networks[END_REF] for senders and receivers.

The procedures of applying the algorithm can be summarized as follows: • Initialize messages at random.

"!! i i i i j 2 j 1 j l l l k j
• Pick in random order all i ∈ V and update messages using ( 8), ( 9) and (10-13) until convergence is reached (i.e., message changes are below a given threshold).

• Use the converged messages to calculate physical observables. We can see that the green path has to be re-routed to avoid both blue and red communications. Also the red one cannot take its shortest path because the sender of the green communication is on that path.

Obtaining a solution

Once the iterative equations ( 8), ( 9) and (10-13) have converged, the resulting messages can be used to calculate the solution. We define the energy per link [START_REF] Mézard | The cavity method at zero temperature[END_REF]:

E Link i j ( Ī) := { E i j ( Ī) + E ji (-Ī) -|| Ī|| } ( 14
)
where the last term on the right || Ī|| is subtracted to avoid double counting as it appears in both of the previous two terms. To find a solution we calculate:

E * Link i j := min Ī E Link i j ( Ī) (15)
for each link in the graph and store the current values that minimize the energy per link for each edge:

Ī * i j := arg minĪ E Link i j ( Ī) . (16)
Eventually, we sum over all (i j) ∈ E to find the different paths and total length:

L tot := ∑ (i j)∈E || Ī * i j || . (17)
In the cavity formalism [START_REF] Mézard | The cavity method at zero temperature[END_REF] this is equivalent to calculating the quantity:

E i := min Īki |constraint ∑ k∈∂i E ki ( Īki ) , (18)
which represents the energy per node. Finally, the total energy (or path length) is the combination of the two, which in the case of a k-regular graph establishes the formal relation:

E t = ∑ i∈V E i - k 2 ∑ (i j)∈E E Link i j . ( 19
)
Notice that the calculation of E Link * i j is carried out link by link as if the energies per link were statistically independent. It is not intuitively clear that doing this will result in the optimal paths which do not overlap and are also fully connected from the source to the receiver. This is a consequence of having used messages which implicitly contain global information on the constraints and path lengths, so that the energies per link are indeed globally interdependent albeit in a non-obvious manner.

To fully characterize the solutions statistically we calculate also other observables as explained below. Finally, we calculate the paths and corresponding lengths in a sequential order using a greedy breadth first-search local algorithm (BFS) to compare the results obtained against our MP-based algorithm.

Algorithmic complexity

The node-disjoint constraint is very restrictive and algorithmically helpful in comparison to other routing models where overlaps are allowed but minimized [START_REF] Yeung | From the physics of interacting polymers to optimizing routes on the london underground[END_REF][START_REF] Yeung | Competition for shortest paths on sparse graphs[END_REF]. This hard constraint is indeed paramount in reducing considerably the algorithm's computational complexity. If we allow for overlaps we need to span a configuration space of the order of 3 M at each cavity iteration, leading to a complexity of O(N (3 M ) k-1 ), where the exponent k -1 explores the different flux combinations for each of the k -1 independent neighboring sites of the considered message; there are order of N such messages. In order to tackle this issue proper approximations have to be introduced as in [START_REF] Yeung | From the physics of interacting polymers to optimizing routes on the london underground[END_REF][START_REF] Yeung | Competition for shortest paths on sparse graphs[END_REF] where they use techniques from polymer physics [START_REF] Yeung | From the physics of interacting polymers to optimizing routes on the london underground[END_REF] or convexity properties of the cost function [START_REF] Yeung | Competition for shortest paths on sparse graphs[END_REF]. On the contrary when the overlap is prohibited we reduce the configuration space from 3 M to (2M + 1), as this is the number of allowed configurations (the term 2M is derived from the number of possible currents I µ = ±1 and the additional +1 is due to the configuration of all 0), hence there is no need for approximations because the entire configuration space can be efficiently calculated by the cavity equation. Actually, the use of cavity MP implicitly requires one important approximation as it assumes that when node i ∈ V is removed, all its neighbors are statistically independent. This is equivalent of having fast decaying correlation function between these neighboring nodes. This hypothesis is verified in trees and in locally tree-like sparse graphs.

For the same reason is important to distinguish between edge and node overlaps. In this work we chose to consider constraints on nodes motivated by the reduced complexity as explained before; in case of edge constraints one has to consider a much bigger configuration space where all configurations with different communications entering and exiting the same transit node must be considered in the optimization routine. For this reason approximations should be introduced as in the case of the models which minimize overlap. The edge-disjoint variant of the problem will be left for future work.

We performed single instance simulations to find optimal microscopic solutions; to obtain macroscopic averages one would usually use population dynamics, one of the most commonly used numerical tools in statistical mechanics literatures [START_REF] Mézard | The bethe lattice spin glass revisited[END_REF][START_REF] Mézard | Information, physics, and computation[END_REF] for studying similar models. Population dynamics is considered when the thermodynamic limit V → ∞ is taken and the system size is not fixed a priori as in the single instance algorithm. In our model the use of population dynamics does not make much sense since the parameter M enters explicitly in the expressions of the messages because it represents the domain of the fluxes, which is of size 2M + 1. But when we fix M at the same time we are fixing a system size V, because we extract random pairs (S , R) with density M/V. Hence, it is impossible to decouple the messages domain from the system size, preventing us to properly employ the thermodynamic limit through population dynamics. There is also another problem, that such a macroscopic oriented approach would introduce averages over all possible configurations (S , R), including both frustrated and unfrustrated configurations with much higher energies. Thus the macroscopic averages are highly biased by the fewer frustrated configurations and more complex algorithm should be designed to discard such cases. For these two reasons we did not consider in the following the population dynamic counterpart of the algorithm but focused only on averages over single instances.

Greedy algorithm

To test the performance of the algorithm we compared the results obtained with those given by a greedy algorithm (or its variant) that is often used in literatures to solve the NDP problem in different contexts [START_REF] Sumpter | Maximizing number of satisfiable routing requests in static ad hoc networks[END_REF][START_REF] Modiano | Minimum energy disjoint path routing in wireless ad-hoc networks[END_REF][START_REF] Chen | A new model for optimal routing and wavelength assignment in wavelength division multiplexed optical networks[END_REF][START_REF] Manohar | Routing and wavelength assignment in optical networks from edge disjoint path algorithms[END_REF]. The greedy protocol considers only local information around the sources and then builds up a solution step by step recursively, hence reducing considerably the complexity but at the same time completely ignoring other communication positions in the network.

A typical greedy algorithm works in the following way: start by choosing an arbitrary pair (S , R), find the shortest path linking the two nodes and then remove nodes belonging to this path from the available network nodes. Choose a second pair and repeat the procedures until either all the M paths from sources to destinations have been established or no solution can be found due to frustration. Clearly, the performance of this algorithm is strictly dependent on the order in which we choose the pairs. For instance, in the extreme case the first pair selected is the one with the longest shortest path among all the M communications; this implies that we have effectively a more restricted graph and choice of paths, leading for a long second path and even more restricted choice of paths later on.

The results

We performed numerical simulations on three types of random graphs. Standard regular random graphs (Reg): each node has fixed degree k; Erdős Rényi random graphs (ER) [START_REF] Erdős | On the evolution of random graphs[END_REF]: edges are drawn at random between each pair of nodes with probability p = ⟨k⟩/V; a decorated random graph (RER): starting from a regular random graph of degree k 1 (which is the minimum degree of this graphs), we then randomly add new edges as in the ER model until the final average connectivity is ⟨k⟩ = k > k 1 . Notice that the degree distribution in this case can be obtained by writing k = k 1 + d where d is Poisson distributed with ⟨d⟩ = kk 1 . The parameters used were average degrees ⟨k⟩ = 3, 5, 7 and system sizes of V = 1000, 2000, 4000, 5000, 10000. We calculated averages over [50 -500] realizations for both the MP and the greedy algorithms; we used a smaller number of realizations for cases of higher complexity (as for V = 10 4 and k = 7). Nevertheless, results in all cases are stable and with small error bars with respect to the symbols used. We omitted the error bars from the figures for clarity. We found a system size scaling that is a cubic function of the variable x := M log V V log γ (k-1) . A qualitative explanation of the scaling is as follows. The average path length in random graphs goes as ⟨l⟩ ∼ log V/ log⟨k⟩ (see [START_REF] Albert | Statistical mechanics of complex networks[END_REF] for an extensive review of graphs properties) and in our case we have M paths to consider.We can refine the dependence on k using instead ⟨l⟩ ∼ log V/ log (k -1). Now, suppose all communications take their shortest path, the quantity x = M log V/ log γ (k -1) would be a good estimate of graph occupancy for the NDP, where the exponent γ has been introduced as a free parameter to account for the approximation in the expression for ⟨l⟩ as a function of k for different types of graphs. Furthermore, if we divide by the number of available nodes V we can define the occupancy ratio as M log V V log γ (k-1) = x. Therefore in this simple case we would expect L tot /V increasing linearly in x. If overlaps are prohibited, for a sufficiently high value of M the communications are increasingly forced to take longer routes, leading to a faster than linear increase in the scaling variable x. From numerical simulations we found for the NDP a cubic increase L tot V = ax + cx 3 in the scaling variable x = M log V V log γ (k-1) . For small M this function agrees well with the linear shortest path Figure 4: Expected total normalized path length -regular graphs (Reg). We obtained results for regular graphs of fixed degree k = 3, 5, 7 and for system size V = 2000, 4000, 5000, 10000. We found a global scaling rule with respect to the variable x := M log V l /V log γ (k -1), with l = 1.00 and γ = 0.87. In black (solid line) we draw the cubic fit ax + cx 3 , with a = 0.910 ± 0.001 and c = 0.660 ± 0.009. The dotted-dashed line represent the shortest path (not considering overlap) trivial solution. Error are not reported because smaller or comparable to point sizes.

behavior but for values of x > 0.2 the steeper increase of L tot becomes predominant. Figures 4 and5 show a good data collapse of the normalized expected total length per node L tot /V as a function of the scaling variable x for different graph connectivities for Reg and ER graphs respectively. We notice a first regime where the curves follow the linear behavior of the dashed line representing the shortest paths. The term "sparse regime" is used since paths are sufficiently far apart, M is small, and then no re-routing is needed as each communication will simply take its shortest path. For x > 0.2 the curves show the cubic steeper behavior that represents the increase in path lengths to avoid overlaps. The term "dense phase" reflects the increase in path density; M is sufficiently high so that shortest-path choices induce conflicting demands and communications are rerouted, taking longer paths to avoid overlaps. Finally, for large M we can identify different frustration points, represented by vertical lines in figure 4, that connect the largest M for which solutions have been found with the points where frustration is reached and the length is set to 0 by convention. We see that the frustration points do not collapse and that the bigger the graph size V and the higher the Figure 5: Expected total normalized path length -Erdős Rényi (ER). We obtained results for ER graphs of expected degree ⟨k⟩ = 3, 5, 7 and for system size V = 1000, 2000, 4000, 5000. We found a global scaling rule with respect to the variable x := M log V l /V log γ (k -1), with l = 1.00 and γ = 0.69. In black (solid line) we draw the cubic fit ax + cx 3 , with a = 0.77 ± 0.01 and c = 0.96 ± 0.03. The dotted-dashed line represent the shortest path (not considering overlap) trivial solution.

connectivity k, the earlier frustration sets in (as a functions of x). Arguably, this is due to algorithmic convergence rather than theoretical arguments. In fact the higher V and k are, the higher the corresponding algorithmic complexity, and hence the larger the number of iterations required to reach convergence. Due to the prohibitive computation cost we ran a smaller number of instances for higher values of V and k, and without increasing the preset maximum convergence time. We suspect that convergence can be reached in these cases albeit in a much longer times, and hence a solution could be found as well in theory, but has not been found due to the computational limits imposed. Hence we can not provide a precise measure for the frustration transitions nor make further statements regarding their collapses for different systems sizes and connectivities.

In figure 6 we can see the scaling behavior for Reg, ER and RER of given average connectivity and different system sizes; we fixed γ = 0 arbitrarily to highlight the dependence on V. We can notice how different types of graph, although having different average lengths, follow the same cubic scaling in x. The steeper slope of the ER graphs shows the smaller number of paths choices in this type of graphs that forces the path to rewire in increasingly more convoluted patterns and hence also reach frustration 6: Scaling in total length. We plotted the scaling of the total length per node as a function of M log V/V for Reg, ER and RER graphs of different system sizes; ⟨k⟩ = 3 for Reg and ER and ⟨k⟩ = 4 for RER. We can see different slopes in the cubic fits of the various curves, for instance ER graphs achieve shorter lengths but with higher cubic slope, meaning that their length increases faster with traffic due to the smaller number of path choices.

We found that our MP algorithm outperforms the greedy BFS both in finding a better solution (smaller L tot ) and in reaching higher values of the frustration transitions. In figures 7, 8 and 9 we plotted the expected normalized total length for both greedy and MP algorithms, for the three graph types, fixed connectivity and different system sizes. We focused on the case ⟨k⟩ = 3 for Reg and ER and ⟨k⟩ = 4 for RER because of its lower complexity compared to higher k; nevertheless, simulations for different k values indeed agree with the suggested scaling and exhibit the same behavior. Initially, in the x range where solutions exist the greedy algorithm gives the same total length as the global algorithm up to a certain value of M (and x). The explanation is that in this interval the graph is sparse, communications typically do not interact and shortest paths can be selected. This also shows that for a small number of paths the global procedure reduces to act similarly to the greedy algorithm does, e.g. when rerouting is required it involves only two paths, the optimal solution will adopt the shortest path for one and will reroute the second. When M increases, we see that the global optimization algorithm outperforms the greedy approach in both finding the optimal solution regime where the global algorithm gives better solution (i.e. shorter total length) we see that it is more efficient to globally reroute paths rather than taking the shortest path of selected paths and adapt the other paths. This means that the optimal solution is not a simple superposition of the first n-shortest path of the M communications, but is a more complex solution. Figure 7: Expected total normalized path length -greedy vs global optimization algorithms -Reg graphs. We compared results obtained by a greedy BFS algorithm and our global optimization for system size V = 2000, 4000, 5000, 10000,γ = 0.87 and degree k = 3. We can identify the sparse interval where the algorithms give similar results because communications are far apart and take the shortest paths. As the number of communications grow we observe an intermediate regime where global optimization performs better than BFS; and finally the dense regime where the greedy BFS algorithm fails to find a solution whereas the global optimization algorithm succeeds up to a critical M value. Cubic fits are also plotted (solid black line for global optimization, dashed line for the greedy BFS algorithm) whereas the dotted line represents the shortest path (not considering overlap) trivial solution, i.e. the sum of the M shortest path lengths, which is linear in x. Vertical lines show the frustration points where no solution is found and the total path length is set to zero. Inset: Ratio L greedy /L MP -1 is plotted as a function of x. Notice the worse performance of the greedy algorithm.

Figure 10 shows the failure ratio defined as the number of unsuccessful instances (for which a solution is not found) over the total number of realizations as a function of the scaling variable x. We notice that the greedy algorithm reaches the frustration point (as a function of x) earlier than the corresponding global MP algorithm, regardless the system sizes or graph type. This shows that, if a solution exists, a global management of the entire set of communications is required in order to find an optimal solution. Whereas if each communication acts selfishly, seeking the corresponding shortest path, unsolvable overlaps between communications emerge at lower x values. Both algorithms show an increased failure rate as the system size increases, presumably due to the unscaled limit on the number of iterations allowed and possibly inherent finite-size effects.

0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 0.5 L tot /V MlogV/ Vlog γ (k-1)

A posteriori statistics: maximum cluster size and degree distribution (regular graph case).

To better understand the optimization process and characterize the solutions obtained we carried out a statistical analysis of the solution a posteriori. Given the clearer statistical interpretation of the results obtained (due to the limited number of possible connectivity values and their evolution, and the higher frustration threshold for a given connectivity), we chose to study regular graphs. In this case one can gain more insights into the type of routes formed and the reduced effective graphs that emerge for any number of communications. By a posteriori we mean that once a solution was found, by an MP or greedy algorithm, we removed from the graph G all nodes and edges taking part in the paths and then calculated statistical properties of the remaining graph G ′ . In particular, we calculated the maximum cluster size and the degree distribution.

The existence of a solution to a given set of communications is strictly related to the connectivity of the graph. Each time a solution for a subset of communications is found, edges and nodes involved in the solution paths are effectively removed; and properties of the reduced graph provide information on its ability to accommodate more source-destination pairs and the efficiency of the obtained solution in making use of the topology. Figure 11 shows the max cluster sizes ratio of G ′ as a function of the scaling variable x. This quantity is defined as the ratio between the number of nodes in the maximum connected cluster over the number of nodes in the same graph G ′ , the graph obtained after edge and node removal of the obtained solution paths. For both the greedy and the global MP algorithm we see an abrupt step change at some x value, between a graph that has a giant connected component and a situation where no solution exists, such that we set the ratio to zero by convention. Moreover, this drop is more abrupt and occurs for smaller x values in the case of the greedy algorithm. This means that the greedy procedure does not distribute paths evenly on the graph and creates small disconnected clusters; the greedy algorithm is therefore more sensitive to V . We can notice that the greedy algorithm fails to find solution earlier than MP for all types of graph considered. ER reaches frustration sooner due to less path choices, whereas RER has higher frustration threshold because of the higher connectivity. The MP data shown are results of averages calculated over a smaller number of instances than for the greedy algorithm, hence the lines connecting them are less smooth.

small changes in connectivity compared to the global MP algorithm, for which the drop is more gradual at first and occurs at higher x values. This reconfirms the previous results that the greedy behavior is fragile and sensitive to the position of the communication pairs and the order in which they are selected.

We evaluate the a posteriori degree distribution P(k) by calculating the connectivities k i ∀i ∈ G ′ for the different k values, and from these derive the average degree ⟨k⟩ as a function of the scaling variable x. Results shown in Figure 12 for different system size and k = 3 show consistent trends; starting from a 3-regular graph we end up, close to the frustration transition point and after the node and edge removal, with about 20% of the nodes with k = 3, whereas ∼ 40% have degree k = 2 and ∼ 30% have k = 1. The decay of ⟨k⟩/k is also plotted for the same process. Also here we see a good data collapse (the different curves can only be distinguished close to frustration).

From graph theory [START_REF] Albert | Statistical mechanics of complex networks[END_REF] we know that when ⟨k⟩ ∼ 1 the graph is likely disconnected, at least to two giant components; the numerical results show that frustration is reached when has value ∼ 50-60%, which corresponds to an average degree ⟨k⟩ ∼ 1.5-1.8, still higher than the connectivity threshold 1. This can be explained by the fact that tighter constraints on edge availability are imposed in the case of the NDP problem, resulting in frustration even before the graph disconnects (i.e., disconnection is sufficient but not necessary for frustration). Indeed in our model it is insufficient to have just a good number of available links, but they should also constitute clusters of connected links to accommodate new communication paths. Hence the average connectivity value observed at the frustration point of ⟨k⟩ > 1.

Path length and stretch distribution.

Another interesting quantity to consider is the path length distribution close to the critical threshold, and its comparison with the shortest path distribution. Using the rescaled variable

L log V/ log γ (k-1)
, where L is the length per communication, we present in figure 13 the distribution obtained for different system sizes. We see a good data collapse for graphs of different system sizes and connectivities to a Gaussian-like distribution with left fat tails, as confirmed by the log-plot on the right panel. This can be explained by the fact that the shorter of the M shortest paths are less likely to be rerouted. A graph with a high number of communications exhibits a path length distribution with higher length averages (with respect to the shortest Figure 12: Degree distribution and ⟨k⟩. We plotted P(k) for k = 0, . . . , 3 as a function of x = M log V/V log γ (k -1) for a 3-regular graph and system size V = 2000, 4000, 5000, 10000. From top to bottom: the curve at the top is P(3) while the one at the bottom represents P(0). The black line represents ⟨k⟩/3 for the different system sizes.

path distribution) as well as higher variances because solution path lengths are more broadly spread. We notice that the left tails are similar for all connectivity values whereas the right tails are broader for lower connectivities close to the frustration point. This can be explained by the fact that short paths are less likely to be rerouted and occur in roughly the same proportion in graphs of different connectivities; hence the similarity in the fat left tails. Regarding the right tails -many paths are rerouted through longer routes by the MP-algorithm, but graphs with higher degree allow for more communications with shorter routes due to the higher routing flexibility they offer.

Figure 14 shows the stretch, defined as the difference between the shortest path and the path length obtained through MP optimization, for M close to frustration ( M log V V log γ (k-1) ∼ 0.5) for different system sizes. We can see that for graphs of degree k = 3 only 34 -38 % of the communications follow the shortest path, all other communications are routed through longer paths. A higher fraction of shortest-path communications is found for higher connectivity graphs, presumably due to the higher routing flexibility they offer. Looking at the tails we can see that there is a non-negligible fraction of paths that stretch considerably compared to the average shortest path length.

Conclusion

We studied the shortest node-disjoint path problem on regular, ER and RER random graphs using message-passing cavity equations. We found that the suggested MP algorithm shortest is normalized and plotted for each communication in the dense regime for degrees k = 3, 5, 7. For k = 3 we have 40% of communications correspond to the shortest paths while others take routes which are long compared to the average path length log V/ log γ (k -1). In the case of higher connectivities k = 5, 7 there is a higher fraction of paths with the same length as the shortest path due to the routing flexibility offered. Right -The same data is plotted on a log scale to highlight the tails behavior.

outperforms the greedy breadth-first search approach both in finding better solutions (shorter total path length) and in finding solutions for higher values of M. This shows that a global strategy is needed to optimally route paths which do not overlap at nodes but also have minimal path lengths. We found a scaling rule for the total length that goes as a cubic function of the occupancy ratio M log V V log γ (k-1) , with γ varying with the graph topology. This behavior resembles the shortest path length for small M but increases faster than linearly for a higher number of paths.

We also studied statistical properties of physical observables a posteriori in the case of graphs. We found good regular graphs of different system sizes and connectivities for quantities such as maximum cluster size, degree distribution and, length and stretch distributions.

We believe this approach is theoretically interesting due to its relevance to hard combinatorial complexity problems but also offers a new direction for solving important practical routing problem in communication, in particular in optical and wireless ad-hoc networks and VLSI design. This study offers the first step for realizing the potential in this new direction.

Figure 2 :

 2 Figure 2: Transit node cavity diagram. Left and center represent the two terms inside the min brackets in equation (8). Right represents equation (9) .

Figure 3 :

 3 Figure3: Edge-disjoint routes. For a single instance of a graph of size V = 20 and M = 3 communications we have on the left the shortest paths and on the right the optimal non-overlapping solutions. We can see that the green path has to be re-routed to avoid both blue and red communications. Also the red one cannot take its shortest path because the sender of the green communication is on that path.

Figure 8 :

 8 Figure 8: Expected total normalized path length -greedy vs global optimization algorithms ER. System sizes V = 1000, 2000, 4000, 5000, γ = 0.69 and degree ⟨k⟩ = 3. Inset: the ratio L greedy /L MP -1 is plotted as a function of x. Notice the worse performance of the greedy.

Figure 9 :

 9 Figure 9: Expected total normalized path length -greedy vs global optimization algorithms RER. sizes V = 1000, 2000, 4000 and degree ⟨k⟩ = 4. Inset: the ratio L greedy /L MP -1 is plotted as a function of x. Notice the worse performance of the greedy.

Figure 10 :

 10 Figure10: Failure rate for greedy and MP algorithms. We plotted for V = 2000 and for Reg and ER of degree ⟨k⟩ = 3 and RER of ⟨k⟩ = 4, the failure rate as a function of M log V V . We can notice that the greedy algorithm fails to find solution earlier than MP for all types of graph considered. ER reaches frustration sooner due to less path choices, whereas RER has higher frustration threshold because of the higher connectivity. The MP data shown are results of averages calculated over a smaller number of instances than for the greedy algorithm, hence the lines connecting them are less smooth.

Figure 13 :Figure 14 :

 1314 Figure13: Path length distribution. Left -Path lengths L ν in the dense interval M log V/ log γ (k -1) ∼ 0.5 are normalized and plotted for k = 3, 5, 7 graphs against the corresponding shortest paths (violet curves). We see a more broadly spread distributions with higher averages for all graphs with respect to the shortest paths, signaling the path rerouting due to the MP optimization. Right -path lengths are plotted in log scale to highlight the left fat tails where the shortest paths are similar irrespective of connectivity, whereas other paths are rerouted by the MP algorithm and are considerably longer, almost by a factor of two (right tails).

Acknowledgement

This work is supported by the Marie Curie Training Network NETADIS (FP7, grant 290038), the EU FET FP7 project STAMINA (FP7-265496) and the Royal Society Exchange Grant IE110151. This work is partially supported by Research Grants Council of Hong Kong (605010 and 604512).


Keywords: male size, competition, agility hypothesis, assortative mating, sexual selection opportunity, Venturia canescens

It is now acknowledged that mating system is a dynamic characteristic within species or populations. Indeed, behavioural characteristics shaping the mating system can vary at the interindividual level, providing raw material for sexual and natural selections, and at the intraindividual level, because of phenotypic plasticity and ontogeny of individual behaviours. Social environment is an additional source of variation in sexual behaviours. By focusing on males of the parasitoid wasp Venturia canescens, we assessed how sexual behaviours (courtship and mating) vary between and within individuals, considering males' size as inter-individual source of variation and age and mating experience as intra-individual sources of variation.

Consequences of the variations in sexual behaviours at the individual level are appraised by quantifying the costs in longevity and offspring production incurred by (repeatedly) mated males. At the population level, we tested "agility hypothesis" predicting that smaller males win the male-male competition to access females, resulting in a non-random assortative mating by size and a negative sexual selection on body size. We showed that V. canescens males are polygynous and that mating incurred longevity cost. Old males had a higher probability to court females, but no mating advantage. Smaller males tended to have a higher mating probability, but size was not under sexual selection and there was a random mating pattern according to individual size. Male-male competition intensity did not influence pairing pattern but tended to influence mating dynamic. Overall, our results suggest that sexual behaviours tended to vary at the inter-and intra-individual levels and the effect of social environment remains to be more precisely described. These variations may impact individuals' reproductive success with potential ecological and evolutionary consequences that are discussed.

Introduction

Mating system defines the number of sexual mates individuals of both sexes have during a breeding season [START_REF] Shuker | The evolution of insect mating systems[END_REF]). The mating system of a species or a population results from individuals' ability to access mates and reproduce, and thus from individuals' sexual behaviours [START_REF] Emlen | Ecology, sexual selection, and the evolution of mating systems[END_REF]. The mating system, often considered as a fixed characteristic of a species or a population, is dynamic because individuals' sexual behaviours and their ability to access mates vary in space and time according different sources of variation [START_REF] Shuster | Mating systems and strategies[END_REF]. First, inter-individual variation in sexual behaviours and characteristics represents a raw material on which natural and sexual selections operate. Second, there is an intra-individual variation in sexual behaviour and characteristics linked to both phenotypic plasticity and ontogeny. Finally, the environment (both social and physical) that directly affects individuals' sexual behaviour and varies in space and time. These three sources of variation (inter-individuals, intra-individuals, and environmental) shape the evolution of mating systems [START_REF] Thornhill | The evolution of insect mating systems[END_REF]. These sources of variation are well illustrated in insects [START_REF] Godfray | Parasitoids: Behavioral and evolutionary ecology[END_REF][START_REF] Shuker | The evolution of insect mating systems[END_REF] showing for instance their impact on variability of sexual behaviours in males and in their ability to secure reproduction.

Inter-individual variations in reproductive opportunities are clearly linked to mate choice and intra-sexual competition for access to a partner. Under the classical sexual selection theory, males compete with each other for females, and females use signals from males to choose their partner [START_REF] Thornhill | Sexual Selection and Nuptial Feeding Behavior in Bittacus apicalis (Insecta: Mecoptera)[END_REF]Metzger, Bernstein, et al., 2010a). Variation among males in their access to females thus depends on their competitiveness and attractiveness. Both may be determined by their body size or the size of their weapons [START_REF] Beukeboom | Size matters in insects-An introduction[END_REF], their ability to find females (Metzger, Fischbein, et al., 2010b), their coloration [START_REF] Hawkes | Sexual selection drives the evolution of male wing interference patterns[END_REF], the vigour and duration of their parade [START_REF] Balaban-Feld | Changes in courtship behaviour following rejection: The influence of female phenotype in Drosophila melanogaster[END_REF].

Intra-individual variations in reproductive opportunities are related to the change in individual's ability to access a partner during its reproductive lifetime. This change is primarily affected by ageing meaning that male's body condition and vigour decrease with age [START_REF] Rockstein | Aging in insects[END_REF]. Ageing can also generate variation in production of male pheromones that are used as signals by females to choose their partner [START_REF] Ingleby | Insect Cuticular Hydrocarbons as Dynamic Traits in Sexual Communication[END_REF]. In the butterfly Bicyclus anynana male sex pheromone composition varies with age, and females prefer mid-aged over young and old males [START_REF] Nieberding | Cracking the olfactory code of a butterfly: The scent of ageing: The scent of ageing in a butterfly[END_REF]. Performance of males, and thus their investment into ejaculates, are often predicted to decrease with age (senescence). Ejaculate quality diminishes with age in the tropical butterfly Bicyclus anynana [START_REF] Kehl | Old male sex: Large ejaculate, many sperm, but few offspring[END_REF]. Intraindividual variations in reproductive opportunities (e.g. number of courtships and matings) are also associated to changes in reproductive cost (e.g. in terms of energy expenditure) during reproductive lifetime (likely not independent of age and experience). Numerous examples support the hypothesis that increasing the number of reproductive opportunities affects residual reproductive success of males. In the ladybird, Menochilus sexmacutatus, frequent mating attempts made by males result in reduction in their body weight, and males engaging in consecutive matings suffer a longer latency before a new mating [START_REF] Shandilya | Cost of mating in male Menochilus sexmaculatus (Fabricius) (Coleoptera: Coccinellidae)[END_REF].

Environmental variations, through social environment that corresponds to change in operational sex ratio (OSR), may influence the intensity in male-male competition [START_REF] Emlen | Ecology, sexual selection, and the evolution of mating systems[END_REF]. These variations in social environment also interact with individual variations (both inter-and intra-individual) to indirectly affect males' ability to access a partner [START_REF] Kelly | The causes and evolutionary consequences of variation in female mate choice in insects: The effects of individual state, genotypes and environments[END_REF]. In the bushcricket, Poecilimon ampliatus, male body condition and male social environment affect male acoustic signals used to attract females. Heavy males modify the characteristics of their songs in response to different social conditions, when competing against light rivals. At the opposite, light males show less plasticity in their acoustic signals in response to social conditions when they compete with heavy males [START_REF] Anichini | To compete or not to compete: Bushcricket song plasticity reveals male body condition and rival distance[END_REF].

At the population level, variations in individual sexual behaviours frequently result, for instance through mate choice, in non-random pairing [START_REF] Shuster | Mating systems and strategies[END_REF] such as assortative mating patterns [START_REF] Arnqvist | Assortative mating by size: A meta-analysis of mating patterns in water striders[END_REF]. These latter may be based on several characteristics : individual colour [START_REF] Mason | Natural Insect Populations and Assortative Mating[END_REF], environmental rearing conditions [START_REF] Geiselhardt | Looking for a similar partner: Host plants shape mating preferences of herbivorous insects by altering their contact pheromones[END_REF] for an example in the mustard leaf beetle), or size which remains the most documented [START_REF] Otronen | Size Assortative Mating in the Yellow Dung Fly Scatophaga Stercoraria[END_REF], for an example in the yellow dung fly). Moreover, assortative mating patterns may be affected by the environmental conditions encountered. Indeed, social environment, by conditioning intensity of intra-sexual competition, sexual behaviours, and availability of the sexual partners, may thus influence the assortative mating patterns. In the ant Tranopelta gilva assortative mating occurs only when mating competition is weak (low ant density and a slightly male-biased sex ratio) [START_REF] Moura | Sex ratio and density may affect temporal variation in sizeassortative mating in a Neotropical ant species[END_REF].

In this study, we assessed how sexual behaviours (mating and courtship characteristics) vary between and within males to better understand mating systems in the parasitoid wasp, Venturia canescens (Hymenoptera: ichneumonidae). Our first aim was to describe courtship and mating characteristics (latency, number and duration) in males. We focused on precopulatory behaviours because V. canescens females are monandrous [START_REF] Collet | Insects and incest: Sib-mating tolerance in natural populations of a parasitoid wasp[END_REF]Metzger et al., 2010b), so that no sperm competition nor cryptic female choice occur. We expected that V. canescens males are polygynous since guarding behaviours have never been reported in these wasps. More precisely, we predicted a scramble competition polygyny (according to Thornhill & Alcock, 1983 classification) since emerging females are widely dispersed [START_REF] Van Santen | Courtship description in the solitary parasitoid wasp Venturia canescens (Hymenoptera)[END_REF], making it unlikely that males monopolise oviposition sites, i.e. "rendez-vous" sites where sexual partners encounter (Metzger, et al., 2010b). Moreover, polygynous systems should lead to differential costs in terms of post-mating residual longevity and offspring sired; the more males have access to reproduction, the shorter they live and the more daughter production may be impaired. As an inter-individual source of variation, we focused on male body size because it classically determines competitiveness (and possibly mate attractiveness), longevity and hence potential number of matings during lifespan.

In V. canescens, male size is a relevant trait; it is variable [START_REF] Metzger | Rôle de l'information dans l'acquisition des partenaires sexuels et le choix du sexe de la descendance chez les hyménoptères parasitoïdes -Exemple de Venturia Canescens Gravenhorst (Hymenoptera: Ichneumonidae)[END_REF] and may be a criterion of male quality used by females for mate choice when ejecting males at the beginning of mount (E.D. and B.C. pers. obs.). We considered potential sources of intra-individual variations in sexual behaviours by investigating effects of age and mating experience of males.

We also experimentally studied the impact of social environment variation on sexual behaviours through variation in operational sex ratio (OSR). This question makes sense because males and females encounter, and mate, in "rendez-vous" sites (e.g. fruits containing hosts), (Metzger, et al., 2010b). As a consequence, adult density and thus individual interactions are expected to vary, leading to different operational sex-ratio (OSR) in these sites. Because of the inverse sexual size dimorphism in V. canescens, this experiment allowed us to test the 'agility hypothesis' predicting that smaller males are superior in scramble competition and mate more frequently or rapidly than larger males [START_REF] Andersson | Evolution of reversed sexual size dimorphism and role partitioning among predatory birds, with a size scaling of flight performance[END_REF][START_REF] Blanckenhorn | Behavioral Causes and Consequences of Sexual Size Dimorphism[END_REF].

This advantage for smaller males should be more acute under stronger male competition (i.e. more male-biased OSR).

At last, we investigated the effect of OSR on the intensity of sexual selection on male size, predicting a negative selection on this trait and an assortative pattern where smaller males mate with larger females (i.e. with higher egg load), and larger males with females independently of their size.

Materials and methods

Insect and rearing facilities

Venturia canescens is a koinobiont (i.e. hosts continue to develop after parasitization) and solitary (i.e. a single adult at most emerges from a given parasitized host) parasitoid wasp of a broad range of moth species [START_REF] Salt | The hosts of Nemeritis canescens a problem in the host specificity of insect parasitoids[END_REF]. To access mating, males court with stereotyped behaviours. Courtships are composed of five steps: localisation, orientation, positioning, mounting and mating. During orientation and positioning, the male whirrs its wings and swings sideways. Females may reject a male when it tries to mount or mate them, leading to a failure from the male point of view [START_REF] Van Santen | Courtship description in the solitary parasitoid wasp Venturia canescens (Hymenoptera)[END_REF]. Experiments and observations were conducted with sexual (parthenogenetic arrhenotokous) individuals, which have a haplodiploid sex determination, allowing us to quantify males' offspring by counting the number of daughters that come from fertilized eggs. This measure is a reliable proxy of male reproductive success as females do not control offspring sex-ratio [START_REF] Metzger | Does constrained oviposition influence offspring sex ratio in the solitary parasitoid wasp Venturia canescens?[END_REF]. In the following experiments, we quantified sexual behaviours by calculating: male latency to court a female, average courtship duration, latency to mate and to re-mate, and duration of each mating.

Wasps used in this study came from a lab rearing (24 ± 1°C, 50 ± 10% relative humidity, DL 12:12). The strain was established from numerous individuals caught in the field during the whole summers 2017, 2018 and 2020 near Valence in an organic orchard (Drôme, France). As a host, we used the pyralid Ephestia kuehniella. All experiments were conducted under the same conditions between 10 am and 4 pm with individuals coming from different rearing boxes (to prevent from testing inbred individuals) and kept individually after emergence in plastic tubes, with food ad libitum (honey). Two hours before starting the observation, the food was removed to standardize the feeding status of individuals. In Experiments 1 and 2, individuals were anaesthetised with CO2 and immediately placed in an arena (see Figure 1). Each arena was filmed using a Raspberry Camera HQ with Raspberry Pi 4, and videos were analysed with the BORIS software [START_REF] Friard | BORIS: A free, versatile open-source event-logging software for video/audio coding and live observations[END_REF].

Experiment 1: Mating system and male sexual behaviours This experiment had two goals. First, we described male mating system. For that, we recorded sexual behaviours of single males in presence of three females. The presence of three females gave each male the opportunity to copulate several times during the observation period, and thus to test whether they are polygynous. Second, we estimated the costs of mating in terms of longevity and number of offspring produced by males.

To be able to identify the females a male mated with, each female was marked by a dot of waterpaint (blue, green or orange) on the thorax. This procedure does not affect mating behaviour in V. canescens [START_REF] Fauvergue | Sterile males in a parasitoid wasp with complementary sex determination: From fitness costs to population extinction[END_REF]. To avoid the establishment of owner versus intruder status, each anaesthetised wasp was confined (thus isolated from the other females) in a corner of an arena by a removable wall (Figure 1A). When all individuals were awake, we opened the wall first for the 3 females (Figure 1B), allowing them to be acclimated to the experimental conditions. Thirty minutes later, we released the male, and recorded male sexual behaviours for 90 minutes (Figure 1C).

To appraise male reproductive success, at the end of the observation, we gave each mated female the opportunity to oviposit in non-limited healthy hosts. We then calculated the proportion of daughters among the emerging parasitoids. We also quantified male longevity by keeping them individually in plastic tubes without food, and checking four times a day (at 8 am, 12 am, 4 pm and 8 pm) whether they were still alive. Longevity was estimated from the median between the time found dead and the last time it was observed alive. Body size was considered as an inter-individual source of variation in sexual behaviours (courtship and mating). The size of individuals was approximated by the measure of their hind left tibia length [START_REF] Pelosse | Does synovigeny confer reproductive plasticity upon a parasitoid wasp that is faced with variability in habitat richness?[END_REF]. We used a dissecting microscope with a camera (moticam 1000, software Motic Image Plus 2.0). A total of 72 males were tested; due to technical issues, videos were incomplete for 13 males, explaining the sample size differences according to measured trait.

Experiment 2: Intra-individual variations in sexual behaviours

We investigated whether age affects males' sexual behaviours (courtship and mating characteristics). For that, one virgin male and one virgin female were placed in the arena. Once awake, they were released from their closed space by opening the removable walls at the same time and sexual behaviours were observed and recorded for 45 minutes. A previous experiment showed that males live an average of 16 days under the same laboratory conditions (data not shown), so we considered ages: 1 day old ("Young") and 7 days old ("Old"). The virgin females were 1 day old. Size of males was measured as in Experiment 1. A total of 208 males were tested.

Experiment 3: Test of "agility hypothesis", effects of OSR on sexual selection and assortative mating patterns This experiment aimed at testing i) the 'agility hypothesis', ii) the intensity and direction of sexual selection on size, expected to be negative under the 'agility hypothesis', iii) the existence of an assortative mating where small males mated with larger females, and if this mating pattern varies with the intensity of male-male competition. For this purpose, we created population cages with a fixed wasp density but different operational sex ratio (OSR), assuming that the more male-biased the OSR, the more intense the male-male competition. In each cage, we observed matings and measured the size of mated partners and of all unmated individuals (see below) using two proxies of size: hind left tibia length (as in Experiments 1 and 2) and the length of the first segment of the abdomen.

Fifty randomly chosen individuals, one or two days old, were placed in a population cage (30*30*30 cm). In each cage, the OSR was 0.3, 0.5 or 0.7 corresponding to 35 females/15 males, 25/25 or 15/35, respectively. We observed each cage for a maximum of 2 hours. We recorded the time between the start of the observation and the formation of the pairs, and removed them immediately from the cage. This procedure of removing the copulating pairs was repeated until 2 hours had elapsed or half of the possible matings had taken place (8 matings for OSR = 0.7 and 0.3, 13 matings for OSR = 0.5; this allowed the individuals still in the cage to have several choices for mating). Three replicates of each OSR condition were conducted.

Data analysis

Experiment 1 showed that males are polygynous (see Results). Since only 4 males mated three times during the observation period, we grouped them with those that mated 2 times, and thus considered 3 categories of males for the statistical analyses: virgin, mated once and re-mated.

To appraise body size as an inter-individual source of variation in mating characteristics, we fitted different GLM. First, we tested the effect of male size i) on the probability to be mated at the end of the observation (GLM, binomial distribution, logit link), ii) on latency before first mating, and on first mating duration (GLM, Gamma distribution and inverse link). For these last two models, only mated individuals were analysed.

Then, we tested the effect of inter-individual variability of males on courtship characteristics (latency and average duration). As for mating, this variability was represented by male size differences and analysed by the means of GLM (Gamma distribution and inverse link).

We investigated the intra-individual variability in sexual behaviour by comparing males' characteristics of first and second matings (mating duration and latency before mating of remated males). We used a Wilcoxon test on paired data. Intra-individual variability of the courtship characteristics was analysed by comparing average duration of courtships and latency before the first courtship displayed by a given male when virgin and then mated.

Cost associated with mating in terms of longevity was analysed by the means of a GLM (Gamma distribution and inverse link) with the male reproductive status (3 modalities: virgin, mated (one time) and re-mated), size and their interaction as explanatory variables. Cost in terms of sex ratio was analysed with a mixed GLM (binomial distribution, logit link) with mating rank (i.e. if the mating was the first or subsequent (as second and third matings were pooled)) as explanatory variable and male identity as random factor. For testing the intraindividual variability, we compared sex-ratio among re-mated males with a Wilcoxon test on paired data.

In Experiment 2, because of the variability between males in courtship occurrence, we first analysed the probability to court as a function of male age, size and their interaction (GLM binomial distribution, logit link). Among males that displayed courtships, we conducted the same analyses for the sexual behaviours (mating and courtship characteristics) as in Experiment 1, but we added age as an explanatory variable together with its interaction with size. Male age was a qualitative variable, with two modalities ("Young" and "Old").

For Experiment 3, we tested the 'agility hypothesis' by testing whether small males mated faster than large males, and whether this relationship between latency and male relative size depended on OSR (GLM Gamma distribution, inverse link; latency before mating as response variable, relative size, OSR (as factor) and their interaction as explanatory variables). Male relative size (MRS) was calculated as follows: MRS = 𝑥-𝑥̅ 𝑥̅ , where x is the size of the male and 𝑥̅ is the mean male size in the cage.

We estimated the intensity of sexual selection on female and male size for each OSR condition following the method presented in Arnold & Wade [START_REF] Arnold | On the measurement of Natural and Sexual Selection: Applications[END_REF] (see Appendix 1).

To assess the assortative mating pattern, we calculated the correlation between the two proxies of male size (tibia and abdomen sizes) to decide whether we should conduct the analysis on both or a single body measure (Spearman's rank correlation). We then modelled assortative mating pattern by using mixed linear models (lmer function from lme4 package, [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF] with male size as response variable, and female size and OSR (as fixed factors) as explanatory variables; cage (i.e. population cage) was included as a random factor. As male and female sizes could be both response variables, we swapped the role of the variables male size and female size in our previous model. Conclusions are qualitatively the same, we present only results with male size as response variable.

We presented the results from the simplest model selected after a backward procedure [START_REF] Crawley | The R book[END_REF]. We accompanied results of the models with effect sizes (effectsize package, [START_REF] Ben-Shachar | Effectsize: Estimation of Effect Size Indices and Standardized Parameters[END_REF]. Effect sizes are presented with their 95% confidence interval.

Analyses were performed with R (R Core Team, 2020).

Results

Experiment 1: Mating system and male sexual behaviours Our results show that Venturia canescens males are polygynous. Fifty percent of the males mated, and 52.78% of these mated males re-mated at least once during the observation period of 1h30 (Figure 2). Smaller males tended to have a higher probability to be mated at the end of the observation (χ² = 2.987, df = 1, p = 0.08; medium effect size Hedges' g = 0.40, IC = [-0.06, 0.86]), suggesting that body size is likely an inter-individual source of variation in mating probability. Nevertheless, among mated males, size did not influence latency before mating and mating duration (latency: F = 0.02, 1 and 32 df, p = 0.88; duration: F = 0.81, 1 and 32 df, p = 0.37).

There is no evidence that size is a source of inter-individual variation in courtship characteristics. We did not detect any effect of male size on the latency before the first courtship (F = 0.426, df = 1 and 63, p = 0.517; 8 over 72 males displaying no courtship were excluded from the analysis). The average duration of courtship of virgin males was independent of their size (F = 1.75, 1 and 54 df, p = 0.19).

Analyses performed among the mated males suggested intra-individual variations in the average courtship duration according to mating experience. The average courtship duration decreased after mating (Figure 3A, 6.0 (±3.0 SE) sec when virgin versus 4.0 (±1.0) sec when mated) (N = 29, V = 373, p= 0.0004; very large effect size Glass' delta = 1.78, IC = [0.63, 2.90]). On the contrary, latency before courtship did not change with mating experience (N = 32, V = 246, p = 0.747).

Among re-mated males (N=19), we detected intra-individual variation in mating duration. The second mating was longer than the first one (Figure 3B, V = 0, p = 0.0007, N = 15; very large effect size Hedges 'g = -1.64, IC = [-2.45, -0.81]). The conclusion was similar when we included two extremes values of mating duration (Appendix 2). There was no difference in mating latency (N = 19, V = 60, p = 0.28).

There was a tendency that matings were costly in terms of longevity. Virgin males lived 21% longer than mated males in average (Figure 4, median: 3.71 days for mated versus 4.46 days for virgin males; F = 3.23, df = 1 and 65, p = 0.077). This was confirmed by large effect size (Hedges' g = -0.65, IC = [-1.13, -0.16]). Nevertheless, there was no difference in longevity within mated males (median = 3.71 days is identical for mated and re-mated). Irrespective of the mating status, longevity was also affected by the male size; larger males live longer (F = 30.81, df = 1 and 65, p <0.0001).

The reproductive success (offspring sex-ratio, SR) was not affected by the rank of mating when considering all mated males (first mating: mean SR = 0.52±0.34; re-mating: mean SR = 0.75±0.39; χ² = 0.784, df = 1, p = 0.376). We did not detect intra-individual variation of SR, when comparing SR from the first and second matings within re-mated males (first mating: mean SR = 0.55±0.32, second mating: mean SR = 0.59±0.32; Wilcoxon paired test, V = 62.5, p = 0.796, N = 17).

Experiment 2: Intra-individual variations in sexual behaviours Old males had a higher probability to court than young males (59.22% of old males versus 31.43% of young males courted the females, χ² = 16.44, df = 1, p < 0.0001; medium effect size:

Odd Ratio = 3.15,IC = [1.73,5.83]). Nevertheless, male age influenced neither courtship characteristics: average courtship duration (F = 1.73, df = 1 and 91, p = 0.192) and latency before the first courtship (F = 0.068, df = 1 and 92, p = 0.79), nor mating: the mating probability (χ² = 0.262, df = 1, p = 0.61), mating latency (F = 1.33, df = 1 and 29, p = 0.26), mating duration (F = 1.73, df = 1 and 29, p = 0.199). Male size did not affect any of the response variables in this experiment. Experiment 3: Test of "agility hypothesis", effects of OSR on sexual selection and assortative mating patterns Venturia canescens presented an inverse sexual dimorphism in size: males are smaller than females irrespective of the proxy used to estimate the size (mean tibia length for males = 1.58 ± 0.16 mm; for females = 1.73 ± 0.13 mm; t = -10.02, df = 354, p<0.0001; mean length of the first segment abdomen for males = 1.06 ± 0.10 mm, for females = 1.29 ± 0.10 mm, t = -21.06, df = 354, p<0.0001). The two proxies of body size were positively correlated both in males and in females (in males: ρ = 0.89, in females: ρ = 0.86). We therefore present the results concerning the length of the tibia (results with the first segment of the abdomen are in Appendix 3 and 4).

The latency before mating was influenced neither by relative size of males (MRS) in a cage (F = 0.719, df = 1 and 83, p = 0.399) nor by the male-male competition (i.e. OSR) (F = 2.546, df = 1 and 82, p = 0.085, Figure 6). This indicated that our data did not support the agility hypothesis. This conclusion is confirmed by the fact that the intensity of sexual selection, I, on size remained weak (and non-significantly different from zero) whatever the sex and the OSR conditions (Table 1). When conducting analyses for each cage separately, we detected an intensity of sexual selection different from 0 only in the cages 1, 2 and 8 for males, and in cage 8 for females. The direction of sexual selection was either weakly negative or positive (Table 1). These values of I also indicate that there is no evidence of any relationship between male body size and the mating probability.

We did not detect any clear pattern of assortative mating. Indeed, there was no correlation between male and female size within pairs formed in the cages (Figure 5, χ² = 0.03, df = 1, p = 0.86) whatever the OSR (interaction between female size and OSR: χ² = 0.933, df = 2, p = 0.627). The conclusions were the same when we analysed the female size as response variable and male size as explanatory variable (male size effect: χ² = 0.142, df = 1, p = 0.707, interaction between male size and OSR: χ² = 0.234, df = 2, p = 0.89).

Discussion

Our results provide a more detailed picture of the mating system in Venturia canescens males, of the variations in sexual behaviours at inter-and intra-individual levels, and their consequences in terms of reproductive costs. At last, we assess the intensity and direction of sexual selection on body size as a potential trait involved in mating.

Males are polygynous (this study) and females are monandrous [START_REF] Collet | Insects and incest: Sib-mating tolerance in natural populations of a parasitoid wasp[END_REF]. Such a polygynous mating system implies that some males have no access to reproduction while others are more successful, opening opportunities for assessing reproductive costs. A cost in terms of realised fecundity (i.e. a decrease of the offspring number per reproductive event during lifetime) is often described in non-social insects [START_REF] Walker | Sperm Utilization Strategies in Nonsocial Insects[END_REF]), but we did not detect such a cost in male V. canescens, at least in the case of two matings close together in time.

Indeed, the offspring sex-ratio did not change with the rank of successive matings in V. canescens, suggesting that sperm depletion did not occur, while it is observed in males of several parasitoid species [START_REF] Damiens | Why do sperm-depleted parasitoid males continue to mate?[END_REF]. However, we cannot rule out that sperm depletion could have occurred after a greater number of matings over a short period of time, preventing V. canescens males from replenishing their sperm stock between matings. This is the case in Trichogramma evanescens in which males transfer 50% of their sperm after 8 matings [START_REF] Damiens | Male reproductive strategy in Trichogramma evanescens: Sperm production and allocation to females[END_REF]. Alternatively, in our experiments a weak sperm depletion could have occurred, compensated for by a longer mating duration (potentially associated with a greater transfer of spermatozoids) as observed in the second mating (Experiment 1). This positive correlation between copulation duration and quantity of transferred sperm, is far from being ubiquitous in insects [START_REF] Weggelaar | Increased copulation duration does not necessarily reflect a proportional increase in the number of transferred spermatozoa[END_REF] but was described, for instance, in the bushcricket Requena verticalis [START_REF] Simmons | Microsatellite analysis of sperm-use patterns in the bushcricket requena verticalis[END_REF]. Regarding longevity, we observed a cost incurred by mated males that was not accentuated with multiple matings. This cost is classically observed in male arthropods [START_REF] Scharf | Reproductive Trade-Offs and Direct Costs for Males in Arthropods[END_REF]. In the Australian tettigoniid Kawanaphila nartee, breeding males live 6% shorter than non-breeding males, and this difference increases when the resources are low, leading to a 48% loss in longevity linked to breeding [START_REF] Simmons | Costs of breeding and their effects on the direction of sexual selection[END_REF]. Proximate factors responsible for the longevity decrease are diverse. For instance, in our laboratory experiments, we could not assess the costs of mate-finding and male-male competition (at least on the two first experiments). Therefore, we probably under-estimated the cost of reproduction, since flight is known to be highly energy demanding in V. canescens [START_REF] Amat | Fuelling flight in a parasitic wasp: Which energetic substrate to use?[END_REF] and that male-male interactions could reduce the lifespan (B.C. data not shown). A complete picture of reproductive costs in males would require further investigations.

Age is often reported as an intra-individual source of variation in sexual behaviours. For instance, females of the butterfly Bicyclus anynana prefer to mate with mid-aged males, so young and old males have a lower mating success [START_REF] Nieberding | Cracking the olfactory code of a butterfly: The scent of ageing: The scent of ageing in a butterfly[END_REF]. Our experiments suggested that courting probability increased with male age while other courtship characteristics were not affected by age. These results could be explained by a physiological change with age (e.g. hormonal modification or motivation). Facing a reduced residual lifespan, males would be more active in courtship. Another explanation for the higher courting probability in older males would be that their courtships are of lower quality. They therefore would need more courtships to seduce a female and achieve mating. At last, this behavioural change could be due to a greater experience in old males as shown in Drosophila melanogaster male [START_REF] Balaban-Feld | Changes in courtship behaviour following rejection: The influence of female phenotype in Drosophila melanogaster[END_REF]. However, this last explanation is unlikely in our study since old and young males had the same mating experience when tested in Experiment 2.

Nevertheless, we can imagine that male experience during its life impacts its motivation to mate, as we found that courtship duration decreased after the first mating (Experiment 1).

Motivation thus may vary with age according to males' reproductive history. Surprisingly, this effect of male age on courtship behaviours did not result in difference of mating probability between young and old males. So, age was not a relevant cue for the choice of the sexual partner by the females or did not confer any advantage for the males. Unfortunately, we did not compare the offspring production of young and old males. This precludes to firmly conclude about an effect of ageing in males.

Because of the inverse sexual dimorphism in size, we predicted that male size was under negative selection and expected that smaller males were superior in scramble competition for access to females ('agility hypothesis', [START_REF] Andersson | Evolution of reversed sexual size dimorphism and role partitioning among predatory birds, with a size scaling of flight performance[END_REF]. In experiment 1, in absence of competition, we detected a weak evidence of greater mating probability in small males suggesting that agility hypothesis is potentially supported in V. canescens. However, our results in population cages clearly invalidated this hypothesis: smaller males did not have shorter latency before mating even under strong intra-sexual competition (i.e. when OSR is male-biased). This conclusion is also corroborated by the fact that our measures of sexual selection intensity and direction indicated that size was not the target of sexual selection. We also expected that operational sex-ratio (OSR) influenced assortative mating patterns. However, we only observed random pairing patterns according to size whatever the intensity of male-male competition. This absence of assortative mating by size is coherent with our results concerning the weak sexual selection on size. The removal of mating pairs from the cages could have prevented us from observing skewed mating success according to body size, for instance if some males mated several times. To reinforce our conclusions, an experiment in which individuals are individually marked and thus identified in each pair would be relevant to follow potential multiple matings.

Our experimental results confirmed that males are polygynous but did not allow us to validate the hypothesis that a competitive advantage of the smaller males (agility hypothesis) is a determinant of scramble competition polygyny (as described in [START_REF] Thornhill | The evolution of insect mating systems[END_REF] in population of V. canescens. It would be relevant to conduct complementary field experiment -with the difficulties associated with monitoring minute insects -to test whether the rank of males' arrival at the 'rendez-vous' sites (e.g. fruits) could be determinant of mating opportunity.

The males that are the fastest to locate odour combination of females and their hosts, and the fastest to reach these sites should have an advantage in accessing to virgin females. B: then the 3 awake females are released in the arena by opening the walls (red arrows) but the male is still confined to his corner during 30 minutes; C: finally, the male is released by opening the wall (red arrow) and the observation begins for 90 minutes. Tables Table 1: Sexual selection opportunity, I, on size for males and females in each OSR and Cage. p-values come from permutation tests. OSR: Operational sex-ratio in population cage (Experiment 3)

Males

Figures

  Figures

FIGURE 1 :

 1 FIGURE 1: Experimental arena seen from above for Experiment 1, with removable walls represented in green. A: each anaesthetised individual is placed in a closed corner of the arena;

FIGURE 2 :

 2 FIGURE 2: Distribution of the number of matings per male (Experiment 1).

FIGURE 3 :

 3 FIGURE 3: Difference in sexual behaviours when a male is virgin versus mated. A: difference in mean courtship duration before and after the first mating. B: difference in mating duration between first and second matings for re-mated males. Two extremes values (mating durations = 241.5 and 495 sec.) were removed from the dataset (Experiment 1).

FIGURE 4 :

 4 FIGURE 4: Relationship between longevity, male status and size (Experiment 1); lines represent the predictions of the selected statistical model; shaded areas represent the confidence interval of prediction.

FIGURE 5 :

 5 FIGURE 5: Relationships between size of sexual partners when mated under different conditions of operational sex ratio (OSR) (Experiment 3).

FIGURE 6 :

 6 FIGURE 6: Relationship between latencies before mating and OSR. Kaplan-Meier plot showing the probability to mate according to time (Experiment 3). Shaded areas represent confidence interval of prediction.

Figure 1 Figure 6

 16 Figure 1

This work was partly funded by Scientific Breakthrough Program (IDEX Lyon) "Micro-be-have".

data and script are available here: https://lbbe-box.univ-lyon1.fr/d/845958507a4c4732a1cc/ and will be deposited on an official platform if the

Appendices

Appendix 1 -Estimation of the intensity of sexual selection on female and male size for each

OSR condition

The intensity of sexual selection on size was calculated as: Ii = Si √𝑣𝑎𝑟 (𝑧) where Si is the selection differential on size and var(z) is the overall variance in size. Selection differential is defined as: Si = 𝑧̅ i -𝑧̅ , where 𝑧̅ i is the mean size of paired individuals, i.e. the mean size after selection, and 𝑧̅ is the mean size of all the individuals, i.e. the mean size before selection. We tested the null hypothesis, I = 0 through a randomization test, as described in Sokal & Rolf [START_REF] Sokal | Biometry[END_REF] with 10,000 permutations. The results presented below were obtained with the measure of the first segment of the abdomen as size proxy. There was no evidence of relation between male and female sizes (χ² = 0.347, df = 1, p = 0.56), confirming that there was no assortative mating by size in V. canescens. The effect of interaction between female size and OSR on male size was not significant (χ² = 0.536, df = 2, p = 0.77) irrespective of the OSR. The same conclusions were obtained when we analysed the female size as response variable and male size as explanatory variable (male size effect: χ² = 0.27, df = 1, p = 0.60, interaction between male size and OSR: χ² = 0.57, df = 2, p = 0.75).

Figure Appendix 3: Assortative mating pattern between male and female size (length of the first segment abdomen used as size proxy) and effect of OSR on this pattern.

Appendix 4-Is the opportunity for sexual selection on body size influenced by OSR? Study using the measure of the first segment of the abdomen as male size proxy.

Under the different conditions of OSR, there was no evidence allowing to conclude that sexual selection acts on body size in males and females V. canescens. The intensity of sexual selection was not different from 0 in the different population cages analysed separately, except in cages 1 and 8 for males, and in cage 8 for females (Table Appendix 4).