Nour El

Houda Bouzouita
email: nour-el-houda.bouzouita@ens-lyon.fr

Anthony Busson
email: anthony.busson@ens-lyon.fr

Herve Rivano
email: herve.rivano@inria.fr

Exploiting Frame Aggregation to Enhance Access Point Selection

Keywords: Networks → Wireless local area networks; Network performance modeling IEEE 802.11 networks, aggregated MAC protocol data unit (A-MPDU), Performance evaluation, Network load, Markov chain

IEEE 802.11 (known as Wi-Fi) has emerged as a vital wireless network access technology for mobile devices. By providing the potential for high connectivity speeds, this technology has led to a dramatic rise in the number of access points (APs). In such environments, mobile devices have the choice to join several Wi-Fi networks. Despite its importance to user Quality of Experience (QoE), the AP selection is still trivial since it focuses at best on the received signal strength if not only the user's history. Crucial metrics that capture the overall dynamics of the AP load condition, such as the network load, are not taken into account. In this paper, we propose to use the Busy Time Fraction (BTF) as a metric to choose the best AP to attach to. The BTF level of a given channel is inferred based on the frame aggregation scheme proposed since the 802.11n standard. In this regard, we propose an analytical model based on a Discrete-Time Markov chain that discerns the theoretical DownLink aggregation levels for probe traffic concurrent to cross traffic. We validate the accuracy of our proposed approach against ns-3 simulations under several scenarios.

CCS CONCEPTS

•

technology for accessing the Internet. Nearly, all companies, universities, and homes use Wi-Fi to satisfy their connectivity needs, which has led to the explosion of wireless data usage and the colossal rise of Access Points (APs), smartphones, and various mobile devices. According to [START_REF]Annual Internet Report[END_REF], Internet users represented 3.9 billion in 2018, and they will exceed 5.3 billion by 2023, while the number of Wi-Fi hotspots will grow four-fold from 169 million hotspots in 2018 to nearly 628 million public Wi-Fi hotspots by 2023. With this huge proliferation, selecting the AP that offers the required performance is still a striking concern and a critical ongoing challenge.

In dense environments, a device seeking connectivity must choose among all the APs that are within its radio range. This association is initiated by the device's operating system based on the best Received Signal Strength Indicator (RSSI) or the user's history or preferences. RSSI or equivalent metrics may help to estimate the physical transmission rate that will be used. But, these simple criteria do not take into account the load of the channels, APs, or any information relative to the network performance imposed by the other stations and the current traffic.

In contrast to these metrics, the network load seems to be an excellent fit since it conditions the throughput and gives a broader view of the network state independently of the transmission conditions of the device. This latter can be expressed in many ways. In this work, we consider the Busy Time Fraction (BTF), defined as the fraction of time the wireless medium is sensed busy due to successful or unsuccessful transmissions. It captures concurrent transmissions, constituting the AP load, as well as inter-network interference.

This paper presents a method allowing a device to estimate the BTF of the channels used by the AP in its vicinity. We consider an active probing approach to estimate this quantity between two devices connected to the same AP. These two devices are not the only recipients of these measures, but they are supposed to be made available in a public database or through a crowd-sensing application in which participant mobile devices measure and share the load of their nearby networks. By building collective knowledge, other devices using the same application will be able to choose their AP as a function of the load and produce their own measures.

The classical channel load estimation tools based on time dispersion cannot be used in Wi-Fi [START_REF] Farshad | On the impact of 802.11n frame aggregation on end-to-end available bandwidth estimation[END_REF][START_REF] Song | Leveraging Frame Aggregation for Estimating WiFi Available Bandwidth[END_REF]. Indeed, in the most modern Wi-Fi versions (802.11n, ac, ax, and the most recent roll-up [START_REF]Standard for Information Technology-Telecommunications and Information Exchange between Systems -Local and Metropolitan Area Networks-Specific Requirements -Part 11: Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) Specifications[END_REF]), a frame aggregation scheme has been introduced. This scheme allows an AP or a station (STA) to send several frames through a unique access to the medium and thus making inter-packet arrival null. Hopefully, the load of the channel conditions the aggregation level. In this paper, we propose a Markov chain that describes the aggregation level, defined as the mean number of sub-frames within the same aggregated frame, as a function of the load. A device can then compare the actual aggregation level to the analytical model to estimate the network load. Our model is compared to ns-3 simulations [START_REF]The Network Simulator[END_REF] allowing us to consider realistic scenarios. For the considered scenarios, we show that the mean aggregation level is an accurate metric to estimate the load of an AP.

The paper is organized as follows. Section 2 describes in detail the background of this work. Section 3 presents our proposed analytical model to estimate the DownLink frame aggregation levels of the probe traffic. Section 4 is dedicated to the performance evaluation and numerical results. Papers that relate to the studied problem are discussed in Section 5. We conclude this paper in Section 6.

BACKGROUND

Before plunging into the details of our proposed system, we present several key background concepts. First, we give an overview of the IEEE 802.11 frame aggregation scheme in Section 2.1. We then give a brief description of how the AP queue handles packets under the frame aggregation.

Frame aggregation overview

Ever since their first introduction in 1997, IEEE 802.11 network has considerably enhanced its throughput by coming with new technological advances to improve the physical (PHY) and the Medium Access Control (MAC) layers. At the MAC layer, one of the significant steps forward was the introduction of the frame aggregation scheme by the 802.11n standard. Since then, this latter has become the default mechanism of sending data on modern Wi-Fi (802.11n, ac, ax). The general concept of frame aggregation is to allow several packets, that are destined to the same destination, to be clustered and sent out as one aggregated frame, and thus reducing the access time (Distributed Inter Frame Space (DIFS), Back-off, etc.) and the PHY overhead as all the aggregated packets share one PHY header. Frame aggregation can be achieved at two levels: at the Aggregate Mac Protocol Data Unit (A-MPDU) and the Aggregate MAC Service Data Unit. Since A-MSDU is rarely implemented, we focus on the A-MPDU aggregation.

Figure 1 depicts the structure of an aggregated frame. The A-MPDU aggregation consists in sending several sub-frames (A-MPDU) with a common Wi-Fi physical header (PHY Header). Each subframe is composed of an MPDU delimiter, followed by the MPDU, consisting of its own MAC header, MAC payload (a packet or any layer 3 Protocol Data Unit), and Frame Check Sequence (FCS). Padding bytes may be appended at the end.

AP queuing system under frame aggregation

Without frame aggregation and for a given enhanced distributed channel access (EDCA) class (an extension of the Distributed Coordination Function scheme that includes Quality of Service traffic prioritization for data packets), the AP follows a First In First Out (FIFO) per-packet scheduling method. With frame aggregation, the aggregated transmissions act as a batch scheduler which alters the timing characteristics of received packets. To explain this approach, we show an example in Figure 2. It consists of an AP queue with two interleaved sequences of packets sent to two different clients, A and B, on a Wi-Fi link. 𝐴𝑖 defines the 𝑖 𝑡ℎ packet targeted to client A while 𝐵𝑖 is the 𝑖 𝑡ℎ packet destined to client B. In the case depicted in Figure 2, the packet at the head of the AP queue is addressed to A. All A's packets (within the limit of the maximum number of aggregated sub-frames) will be aggregated and sent out as an A-MPDU. Then, an aggregated frame with B's packets is sent.

A1 A2 A3 A4 A5 B1 B2 B3 B4 A6 B1 B2 B3 B4

SYSTEM DESCRIPTION

We consider a general Wireless Local Area Network (WLAN) based on IEEE 802.11 Distributed Coordination Function (DCF). Figure 3 depicts the scenario considered in our model. An AP and three stations operate on the same Wi-Fi channel. The AP is connected to a distribution system (a wired network) which is not modeled here. Each station or AP has its own physical transmission rate. It is a combination of the following parameters: the modulation type and the coding rate which is defined by the Modulation and Coding Scheme index (MCS), the number of spatial streams, the channel width, and the guard interval length. In our model, we distinguish three types of traffics:

• The UpLink (UL) probe traffic sent from the probe traffic client to the AP. It is a Constant Bit Rate (CBR) traffic, sent at regular interval 𝑑 𝑝 . • The DownLink (DL) probe traffic forwarded to the probe server.

• The cross traffic coming from the distribution system (not represented in the figure) and sent to the cross traffic server.

It is modeled by a constant bit rate source sending packets at regular interval 𝑑 𝑐 and managed by a unique queue. The AP and the probe traffic client buffers are assumed to have a finite size. More precisely, we assume that when an aggregated frame is sent for a given destination, the buffer becomes empty for this destination. This assumption simplifies drastically the model and its computation. However, the simulation results show that it does not introduce a significant error. Table 1 summarizes the principal notation used in the paper. In this section, we derive a mathematical model that estimates the frame aggregation levels of the DL probe traffic in this scenario. This frame aggregation level is a function of the cross traffic. The model aims to be compared with the aggregation levels as measured by the probe traffic server to deduce the cross traffic load (load of the AP or channel).

The frame aggregation level of the DL probe traffic depends on the competition for medium access. In particular, when the cross traffic increases, the DL probe traffic has to wait longer. More packets may hence accumulate in the queue of the AP and be aggregated when the next DL probe transmission occurs.

The proposed analytical model is a Discrete-Time Markov Chain (DTMC). We consider the Markov chain defined as (𝑋 𝑛 , 𝑌 𝑛 , 𝑍 𝑛 , 𝑆 𝑛 ) 𝑛 ≥0 . 𝑋 𝑛 describes the number of probe sub-frames contained in the AP queue at the moment of the 𝑛 𝑡ℎ frame transmission departure.

Transition probabilities

Having defined the set of possible states for each process, we shall now derive the transition probabilities. The transition probabilities are denoted 𝑃 (𝑖,𝑗,𝑘,𝑢) (𝑙,𝑚,𝑞,𝑣) , and represent the probability to go from state (𝑋 𝑛 , 𝑌 𝑛 , 𝑍 𝑛 , 𝑆 𝑛 ) = (𝑖, 𝑗, 𝑘, 𝑢) to state (𝑋 𝑛+1 , 𝑌 𝑛+1 , 𝑍 𝑛+1 , 𝑆 𝑛+1 ) = (𝑙, 𝑚, 𝑞, 𝑣).

The transition probabilities depend on the time between two successive transmissions. As both UL probe and UL cross traffics are deterministic, this time sets the number of packets that arrived in the AP buffer and the probing client buffer between two transmissions and thus the number of frames that will be sent in the next aggregated frame. We therefore analyze the events that may occur between two successive transmissions.

The next stage of our modeling is to decide when a transition from one state to another is allowed and to compute its probability. Note that impossible transitions have zero probability and that non-null transition probabilities are computed assuming that all concerned nodes are equally likely to access the channel. For ease of illustration, we categorize the state transitions into the following three classes.

Class I:Transition from state APP. For the sake of clarity, we recourse to Figure 4a to describe the possible transition probabilities and their computations. It illustrates an example of possible events between two successive transmissions when the current state is APP (an aggregated frame containing probe packets is transmitted by the AP). In this example, at time 𝑡 = 0, we start from state (𝑖, 𝑗, 𝑘, 𝐴𝑃𝑃) (i.e (3,2,2,APP) in Figure 4a). The AP sends the 𝑖 corresponding probe frames currently in its buffer (𝑖 = 3 in this example), aggregating them in a single A-MPDU using its current MCS. Its transmission duration is denoted 𝑓 (𝑖) which is given in Appendix A. Note that the AP has a single buffer that contains at once the probe and the cross traffics but in order to clarify the explanation of the model, we distinguish between them in Figure 4a. ⌋ of cross traffic packets ((4,4) in the Figure 4a). At the end of this transmission, the buffer of the AP does not have probe frames to send (𝑋 𝑛+1 = 0 almost surely), and another APP transmission is impossible (𝑆 𝑛+1 ≠ 𝐴𝑃𝑃 almost surely). So from this state, only two transitions are allowed: to APC or SP with 𝑋 𝑛+1 = 𝑙 = 0. It can occur only if 𝑌 𝑛+1 = 𝑚 > 0 and 𝑍 𝑛+1 = 𝑞 > 0 respectively.

We derive the non-null transition probabilities as follows. If the AP gains access to the channel, it will send the cross traffic currently in its buffer. For 𝑚 > 0, the next transmission will be APC with probability:

𝑃 (𝑖,𝑗,𝑘,𝐴𝑃𝑃 ) (0,𝑚,𝑞,𝐴𝑃𝐶) = 𝑃 𝑆 𝑛+1 = 𝐴𝑃𝐶 |𝑆 𝑛 = 𝐴𝑃𝑃, 𝑍 𝑛+1 = 𝑞, 𝑌 𝑛+1 = 𝑚 > 0 • 1 𝑘+ ⌊ 𝑓 (𝑖 ) 𝑑𝑝 ⌋=𝑞 • 1 𝑗+ ⌊ 𝑓 (𝑖 ) 𝑑𝑐 ⌋=𝑚
where 1 𝑐𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛 is the indicator function that equals to 1 if 𝑐𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛 is true and 0 otherwise. 𝑃 𝑆 𝑛+1 = 𝐴𝑃𝐶 |𝑆 𝑛 = 𝐴𝑃𝑃, 𝑍 𝑛+1 = 𝑞, 𝑌 𝑛+1 = 𝑚 > 0 denotes the probability that the event APC will occur given that the event APP has already occurred. In the interest of brevity, we postpone the computation of such probabilities to Appendix B. Now, if the probe station gains access, the next event will be SP and the transition probability from APP to SP for 𝑞 > 0 is given by:

𝑃 (𝑖,𝑗,𝑘,𝐴𝑃𝑃 ) (𝑙,𝑚,𝑞,𝑆𝑃 ) = 𝑃 𝑆 𝑛+1 = 𝑆𝑃 |𝑆 𝑛 = 𝐴𝑃𝑃, 𝑌 𝑛+1 = 𝑚, 𝑍 𝑛+1 = 𝑞 > 0 • 1 𝑘+ ⌊ 𝑓 (𝑖 ) 𝑑𝑝 ⌋=𝑞 • 1 𝑗+ ⌊ 𝑓 (𝑖 ) 𝑑𝑐 ⌋=𝑚
Class II:Transition from state APC. Once again, when presenting the transition probabilities of this class, we resort to Fig- ure 4b. It depicts a timeline of feasible events between two successive transmission departures when we start from the state APC.

The AP sends the 𝑗 cross traffic frames currently in its buffer as an A-MPDU using the transmission rate associated to the cross traffic server. Its transmission duration is denoted ℎ( 𝑗). During the period ℎ( 𝑗), the probe traffic client buffer receives ⌊ Conversely to the previous class of transitions where only two possible transitions are allowed from the state APP, there are here three possible transitions under some conditions. 𝑆 𝑛+1 can be APP, APC, or SP.

First, we suppose that there will be another APC, the transition from APC to APC is deemed possible if and only if 𝑋 𝑛 = 𝑖 = 0 and 𝑌 𝑛 = 𝑙 = 0: it is impossible to have two successive APC transmissions if 𝑖 > 0 or 𝑙 > 0 due to the AP queuing system detailed in Section 2.2. The transition probability from APC to APC with 𝑌 𝑛+1 = 𝑚 > 0 is:

𝑃 (0,𝑗,𝑘,𝐴𝑃𝐶) (0,𝑚,𝑞,𝐴𝑃𝐶) = 𝑃 𝑆 𝑛+1 = 𝐴𝑃𝐶 |𝑆 𝑛 = 𝐴𝑃𝐶, 𝑋 𝑛 = 𝑋 𝑛+1 = 0, 𝑌 𝑛+1 = 𝑚 > 0, 𝑍 𝑛+1 = 𝑞 • 1 ⌊ ℎ ( 𝑗 ) 𝑑𝑐 ⌋=𝑚 • 1 𝑘+ ⌊ ℎ ( 𝑗 ) 𝑑𝑝 ⌋=𝑞
Now if we suppose that the AP gains the medium access to transmit the DL probe flow, the transition probability from APC to APP is possible only if 𝑋 𝑛 = 𝑖 = 𝑋 𝑛+1 = 𝑙 since the AP cannot receive probe frames in its buffer during the transmission of the cross traffic and 𝑋 𝑛+1 = 𝑙 > 0. The corresponding probability is given by:

𝑃 (𝑖,𝑗,𝑘,𝐴𝑃𝐶) (𝑙,𝑚,𝑞,𝐴𝑃𝑃 ) = 𝑃 𝑆 𝑛+1 = 𝐴𝑃𝑃 |𝑆 𝑛 = 𝐴𝑃𝐶, 𝑋 𝑛+1 = 𝑋 𝑛 = 𝑙 > 0, 𝑍 𝑛+1 = 𝑞 • 1 ⌊ ℎ ( 𝑗 ) 𝑑𝑐 ⌋=𝑚 • 1 𝑘+ ⌊ ℎ ( 𝑗 ) 𝑑𝑝 ⌋=𝑞
If we assume that the probe traffic client gains the medium access, the transition from APC to SP is also allowed only if 𝑋 𝑛+1 = 𝑙 = 𝑋 𝑛 = 𝑖 and 𝑍 𝑛+1 = 𝑞 > 0. We have:

𝑃 (𝑖,𝑗,𝑘,𝐴𝑃𝐶) (𝑙,𝑚,𝑞,𝑆𝑃 ) = 𝑃 𝑆 𝑛+1 = 𝑆𝑃 |𝑆 𝑛 = 𝐴𝑃𝐶, 𝑍 𝑛+1 = 𝑞 > 0, 𝑋 𝑛+1 = 𝑋 𝑛 = 𝑙, 𝑌 𝑛+1 = 𝑚 • 1 ⌊ ℎ ( 𝑗 ) 𝑑𝑐 ⌋=𝑚 • 1 𝑘+ ⌊ ℎ ( 𝑗 ) 𝑑𝑝 ⌋=𝑞
Class III:Transition from state SP. In order to derive the transition probabilities of the last class, we apply the same principle. Due to space limitations, we do not describe the details of the derived transition probabilities.

From the SP state, there are three possible transitions APP, APC, or SP depending on the competition for the channel resource. First, if the next transition is APP with 𝑋 𝑛+1 = 𝑙 > 0, the transition probability is given by:

𝑃 (𝑖,𝑗,𝑘,𝑆𝑃 ) (𝑙,𝑚,𝑞,𝐴𝑃𝑃 ) = 𝑃 𝑆 𝑛+1 = 𝐴𝑃𝑃 |𝑆 𝑛 = 𝑆𝑃, 𝑋 𝑛+1 = 𝑙 > 0, 𝑍 𝑛+1 = 𝑞, 𝑌 𝑛+1 = 𝑚 • 1 𝑙=𝑖+𝑘 • 1 𝑗+ ⌊ 𝑔 (𝑘 ) 𝑑𝑐 ⌋=𝑚 • 1 ⌊ 𝑔 (𝑘 ) 𝑑𝑝 ⌋=𝑞
Second, if the next transition is APC with 𝑌 𝑛+1 = 𝑚 > 0, the transition probability is defined as: 𝑃 (𝑖,𝑗,𝑘,𝑆𝑃 ) (𝑙,𝑚,𝑞,𝐴𝑃𝐶) = 𝑃 𝑆 𝑛+1 = 𝐴𝑃𝐶 |𝑆 𝑛 = 𝑆𝑃, 𝑌 𝑛+1 = 𝑚 > 0,

𝑍 𝑛+1 = 𝑞, 𝑋 𝑛+1 = 𝑙 • 1 𝑙=𝑖+𝑘 • 1 𝑗+ ⌊ 𝑔 (𝑘 ) 𝑑𝑐 ⌋=𝑚 • 1 ⌊ 𝑔 (𝑘 ) 𝑑𝑝 ⌋=𝑞
Finally, if the next transition is SP with 𝑍 𝑛+1 = 𝑞 > 0, the transition probability is thus formulated as:

𝑃 (𝑖,𝑗,𝑘,𝑆𝑃 ) (𝑙,𝑚,𝑞,𝑆𝑃 ) = 𝑃 𝑆 𝑛+1 = 𝑆𝑃 |𝑆 𝑛 = 𝑆𝑃, 𝑍 𝑛+1 = 𝑞 > 0, 𝑋 𝑛+1 = 𝑙, 𝑌 𝑛+1 = 𝑚 • 1 𝑙=𝑖+𝑘 • 1 𝑗+ ⌊ 𝑔 (𝑘 ) 𝑑𝑐 ⌋=𝑚 • 1 ⌊ 𝑔 (𝑘 ) 𝑑𝑝 ⌋=𝑞

Stationary probability

So far, we have evaluated all the transition probabilities. The last phase consists in deriving the stationary probability and computing the frame aggregation levels. Let us remind that a Markov chain is irreducible if and only if every state can be reached by any other state through one or several transitions. As our Markov chain is irreducible and has a finite number of states, hence it exists a unique stationary distribution. We solve this Markov chain through a numerical method. We compute the stationary probabilities and we denote by 𝜋 the vector containing the corresponding values. The mean DL probe traffic aggregation levels, denoted 𝑀𝑒𝑎𝑛𝐴𝑔𝑔 are then given by: 𝑀𝑒𝑎𝑛𝐴𝑔𝑔 = 𝐴𝑀𝑃𝐷𝑈 𝐴𝑃 𝑛=1 𝑛 • 𝜋 𝑛 .

PERFORMANCE EVALUATION

In this section, we assess our model's accuracy by comparing its outcomes with those delivered by the discrete-event simulator ns-3 (version 3.30).

Simulation settings

The simulation setup consisted of the one described in Figure 3 where a probe traffic client associated with an AP transmits User Datagram Protocol (UDP) flow to the probe server. Another server station is associated with the AP to receive a UDP cross traffic. The cross traffic generation rate is varied over a wide range of values covering saturated and unsaturated network conditions. The simulated cross traffic is set in order to have six different BTFs: 0, 0.125, 0.25, 0.375, 0.5, and 0.625. All the nodes use the 802.11n standard and operate on the 2.4 GHz frequency band. We provide three different scenarios to evaluate the performance of our model. The different scenarios share the same parameters but differ in the transmission rate used by each node.

Evaluation

We start by examining the performance details of our model when the same MCS index (MCS 15 with a physical transmission rate of 144.4 Mbps) is used for the probe traffic client and the AP. The maximum number of aggregated sub-frames 𝐴𝑀𝑃𝐷𝑈 𝐴𝑃 = 𝐴𝑀𝑃𝐷𝑈 𝑃 =36. Figure 5 shows the mean aggregation level for the probe traffic as a function of the probe packet gap for the model and simulations. The probe packet gap gradually varies from 50𝜇s to 1000𝜇s.

According to these results, it appears that the model performs well since it follows closely the pattern of the ns-3 simulations for the three levels of the network loads. It is also the case for the three other BTFs (0, 0.125, and 0.625) that are not shown by lack of space. We can observe that the aggregation levels depend on loads of the network. Indeed, when the BTF increases, the probe traffic has to wait longer, and more packets are received between two successive probe transmissions.

Figure 6 shows the model and ns-3 simulations outcomes when we use different MCS for the probe traffic client and the AP. Figure 6a plots the results when we use the MCS15 (144.4 Mbps) for the AP to send the DL probe traffic and the cross traffic and the MCS11 (57.8 Mbps) for the probing client. Figure 6b provides the corresponding results when we use MCS15 for the AP and MCS13 (115.6 Mbps) for the probe traffic client. The results show that considering different MCS for the AP and the probe traffic client does not impact the accuracy of our approach.

Figure 7 illustrates the mean probe aggregation levels as a function of the probe packet gap for all the levels of BTF used in our work for the case where all the nodes (AP and stations) use the same transmission rate (144.4Mbps). These results nicely highlight the fact that varying the level of the network load significantly affects the probe aggregation levels. Consequently, the results are sufficiently separated to be used to infer the load level.

Inspired by the correlation between the aggregation level and the network load, we propose the following method that aims to infer the network load (BTF). For doing so, we send a batch of probe streams with an increasing inter-packet arrival from a client to a server. Upon receiving the probe packets, the server computes the mean aggregation levels. By relying on the model's results, the server then computes the error between the measured and theoretical (derived by the model) aggregation levels. The estimated BTF is the one that minimizes this error. Table 2 and Table 3 present the estimated BTFs for the ns-3 simulations when we use the same physical transmission rate (MCS15 with 144.4 Mbps) for all the nodes or different transmission rates (MCS15 with 144.4 Mbps for the AP and MCS13 with 115.6 Mbps for the probe traffic client) respectively. The 'ground truth column' gives the real value of the BTF. 'S' indicates the BTF used in simulations and 'M' presents the value returned by our approach. This approach finds the good result when the two letters ('S'/'M') are in the same box.

Results presented in Table 2 show that the method finds the good results for all the levels of BTF. We note that our model takes into account the time to send block ACK requests while computing the frame transmission duration. These requests mean that a response is requested upon transmission of a frame whose sequence number is distant at least by a given threshold multiplied by the transmit window size from the starting sequence number of the transmit window. Table 3 shows that the method finds the good results for the five BTFs 0, 0.25, 0.375, 0.5, and 0.625 and slightly overestimates the BTF for the case 0.125. This error is only of 0.125 (approximately 10%) and is due to the fact that we neglect some protocol aspects, such as the beacons sent by the AP. We obtain similar results as Table 3 for the other case where MCS11 is used for the probe traffic client. However, by lack of space, we do not present its outcomes in this paper.

Overall, the results show that the proposed analytical model accurately captures the mean aggregation level across different network settings and predicts with reasonable accuracy (sometimes with a shift of one BTF level) the BTF level.

RELATED WORK

Performance evaluation of Wi-Fi networks has been intensely studied in the literature to offer relevant information and methods to help its configuration and improvement. In this paper, we focus on three categories of work: Available bandwidth estimation, crowdsensing of wireless networks, and Analytical models.

Available bandwidth estimation. Conversely to achievable throughput estimation tools, there exists another pivotal metric for estimating the overall network performance which is the available bandwidth. Active Available bandwidth measurements that involve the usage of probing packets can be generally divided into two categories: Packet Rate Model (PRM) and Probe Gap Model (PGM). The works based on PRM, such as TOPP [START_REF] Melander | A new end-to-end probing and analysis method for estimating bandwidth bottlenecks[END_REF], pathload [START_REF] Jain | Pathload: A Measurement Tool for End-to-End Available Bandwidth[END_REF], pathChirp [START_REF] Ribeiro | PathChirp: Efficient Available Bandwidth Estimation for Network Paths[END_REF] and DietTOPP [START_REF] Johnsson | An Analysis of Active End-toend Bandwidth Measurements in Wireless Networks[END_REF], focus on comparing the probe packet and received packet rates. They inject increasing probe traffic in the network path and try to find the point where link congestion starts occurring to discern the available bandwidth. These methods differ according to the probing rate settlement.

Conversely to PRM, PGM techniques (e.g., Spruce [START_REF] Strauss | A Measurement Study of Available Bandwidth Estimation Tools[END_REF], Initial Gap Increase/Packet Transmission Rate (IGI/PTR) [START_REF] Hu | Evaluation and characterization of available bandwidth probing techniques[END_REF], Delphi [START_REF] Ribeiro | Multifractal cross-traffic estimation[END_REF]...) leverage the dispersion in time between two received packets to infer the cross traffic. All these approaches were primitively focused on wired networks and unfortunately failed with the complexity of CSMA/CA of Wi-Fi networks. To copy with this latter, more complex schemes were developed, such as Exact [START_REF] Shah | Available bandwidth estimation in IEEE 802.11-based wireless networks[END_REF], IdleGap [START_REF] Lee | Bandwidth Estimation in Wireless Lans for Multimedia Streaming Services[END_REF] and WBest [START_REF] Li | WBest: A bandwidth estimation tool for[END_REF]. These techniques intensely depend on the fact that the received packet rate or packet gap should be an effective indicator to reveal the cross traffic. However, the dissimilar packet scheduling system under frame aggregation, proposed by modern Wi-Fi versions, prevents the use of these techniques as inter-packet arrival becomes null when they were embedded in the same aggregated frame. For monitoring the network conditions of these new versions of Wi-Fi, WBest+ [START_REF] Farshad | On the impact of 802.11n frame aggregation on end-to-end available bandwidth estimation[END_REF] and AIWC [START_REF] Song | Leveraging Frame Aggregation for Estimating WiFi Available Bandwidth[END_REF] use the characteristics of frame aggregation to discern the available bandwidth by relying on a threshold-based method: if an inter-arrival time is less than a given threshold, the two packets are considered as aggregated.

In [START_REF] Bouzouita | Analytical study of frame aggregation level to infer IEEE 802.11 network load[END_REF], the mean aggregation level is compared to the one of an analytical model to estimate the network load. The network scenario used in these papers consists of a server with a wired Internet connection that sends [START_REF] Farshad | On the impact of 802.11n frame aggregation on end-to-end available bandwidth estimation[END_REF][START_REF] Song | Leveraging Frame Aggregation for Estimating WiFi Available Bandwidth[END_REF] or receives [START_REF] Bouzouita | Analytical study of frame aggregation level to infer IEEE 802.11 network load[END_REF] probing traffic along the network path to/from a client with a last-hop wireless connection that assumed to be the bottleneck. The scenario considered here is different as both the probe traffic generator and the receiver are associated with the same Wi-Fi AP. The behavior of the mean aggregation is then much more complex as it may be different on the uplink and downlink links due to their different parameters (e.g., MCS). Moreover, the probe traffic accesses twice the channel that decreases the level of probe traffic for which the network becomes saturated and the size of the aggregated frames. While evaluating the probe aggregation levels of the UL probe traffic, the model presented in [START_REF] Bouzouita | Analytical study of frame aggregation level to infer IEEE 802.11 network load[END_REF] assumes that the connection between the AP and the server is ideal. It is therefore modeled as if the server were implemented on the AP. However, the model presented here appraises the aggregation levels of the DL probe traffic by considering that the server is embedded on a second wireless device which is more realistic and practical.

Crowd-sensing. Alternatively, some techniques focus on crowdsensing applications where mobile devices become active sensors to actively estimate the end-to-end network performance. The authors of [START_REF] Rosen | MCNet: Crowdsourcing wireless performance measurements through the eyes of mobile devices[END_REF] proposed MCNet, a tool based on active smartphone measurements to estimate the user-perceived performance in enterprise wireless networks. Another work based on channel scans (spectrum sensing) was also proposed in [START_REF] Shi | A walk on the client side: Monitoring enterprise Wifi networks using smartphone channel scans[END_REF]. It is based on passive measurements while our work investigate the effectiveness of active ones. The authors in [START_REF] Farshad | Urban WiFi characterization via mobile crowdsensing[END_REF] proposed a measurement crowdsensing study in the city of Edinburgh to characterize urban Wi-Fi. It revealed several problems in Wi-Fi deployments in public spaces. However, none of these works involves the network load while crowd-sensing the network performance.

Analytical models. Wi-Fi Frame aggregation has attracted a lot of attention due to the significant efficiency enhancements brought by this feature. Consequently, several frame aggregation analytical models have been studied in recent years. [START_REF] Hajlaoui | An accurate two dimensional Markov chain model for IEEE 802.11n DCF[END_REF][START_REF] Mohammad | Modeling and Analyzing MAC Frame Aggregation Techniques in 802.11n Using Bi-dimensional Markovian Model[END_REF] predict throughput while incorporating frame aggregation under saturated conditions. A Discrete-Time Markov Chain (DTMC) that considers the frame aggregation has been proposed to show its performance under various aspects under unsaturated conditions [START_REF] Kim | Effect of Frame Aggregation on the Throughput Performance of IEEE 802[END_REF]. Unfortunately, these models are not applicable or extensible to our work since they consider dissimilar scenarios and traffic patterns.

CONCLUSION

In this paper, we have proposed an analytical model that predicts the DownLink probe aggregation levels concurrent to cross traffic. Based on that, we proposed an approach to estimate the expected network load via the Busy Time Fraction. This information could be then used by a crowd-sourcing application, in which participant devices measure and share the load of their surrounding networks. They would build a collective knowledge so that when a device arrives in some area and wants to connect to a Wi-Fi network, it could choose the less loaded AP. Our model was validated through ns-3 simulations across a variety of scenarios.

Overall, the model estimates with a satisfactory degree of precision the aggregation levels. We showed that this frame aggregation can effectively capture the AP Busy Time Fraction with at most 10% of errors. Our future efforts include addressing some of the open challenges in network load estimation by an unmodified mobile device.

Figure 1 :

 1 Figure 1: Illustration of A-MPDU frame aggregation.

Figure 3 :

 3 Figure 3: Scenario considered in the model. The probe traffic client and server estimate the current load. The cross traffic server aims to represent the current load of the system.

  The possible states are {0, .., 𝐴𝑀𝑃𝐷𝑈 𝐴𝑃 }. 𝑌 𝑛 defines the number of cross traffic sub-frames in the AP queue. The possible states are {0, .., 𝐴𝑀𝑃𝐷𝑈 𝐴𝑃 }. 𝑍 𝑛 represents the number of packets at the probe traffic client. The possible states are {0, .., 𝐴𝑀𝑃𝐷𝑈 𝑃 }. 𝑆 𝑛 describes the 𝑛 𝑡ℎ transmission. It takes three possible values {𝐴𝑃𝐶, 𝐴𝑃𝑃, 𝑆𝑃 }. 𝐴𝑃𝑃 (Access Point Probe) is a DL transmission of probe traffic from the AP. 𝐴𝑃𝐶 (Access Point Cross) denotes a cross traffic transmission from the AP to the cross server, while 𝑆𝑃 (Station Probe) corresponds to a UL transmission from the probe traffic client.

Figure 4 :

 4 Figure 4: Possible transitions

⌋

  cross traffic packets. At the end of the APC state, the probe traffic client buffer and the AP will thus contain respectively, 𝑘 + ⌊ ℎ ( 𝑗) 𝑑 𝑝 ⌋ of probe packets and ⌊ ℎ ( 𝑗) 𝑑 𝑐 ⌋ of cross traffic packets.

5 Figure 5 :Figure 6 :

 556 Figure 5: Mean aggregation level versus Busy Time Fraction

Figure 7 :

 7 Figure 7: Mean Aggregation level for all BTFs for the simulations.

Table 1 :

 1 Principal notation

	Parameter (unit)	Definition
	𝑇 𝐷𝐼 𝐹 𝑆 (𝜇s)	Distributed Inter Frame Space duration
	𝑇 𝑆𝐼 𝐹 𝑆 (𝜇s)	Short Inter Frame Space duration
	𝑇 𝑃𝐻𝑌 (𝜇s)	The preamble and PHY header duration
	𝐹𝐶𝑆 (𝐵𝑦𝑡𝑒𝑠)	Frame Check Sequence
	𝑇 𝐵𝑙𝑜𝑐𝑘𝐴𝐶𝐾 (𝜇s)	Required time to send the block acknowledgment
	𝑇 𝑏𝑎𝑐𝑘𝑜 𝑓 𝑓 (𝜇s)	Average backoff time
	𝑇 𝑠𝑙𝑜𝑡 (𝜇s)	Slot time
	𝐶𝑊 Required time to send k UL probe aggregated packets
	ℎ (𝑘) (𝜇s)	Required time to send k cross traffic sub-frames
	𝑑 𝑐 (𝜇s)	Inter-arrival time of cross traffic sub-frames
	𝑑 𝑝 (𝜇s)	Inter-arrival time of probe traffic sub-frames
	𝐴𝑀𝑃𝐷𝑈 𝐴𝑃	Maximum A-MPDU size for the AP
	𝐴𝑀𝑃𝐷𝑈 𝑃	Maximum A-MPDU size for the probe traffic client

𝑚𝑖𝑛 Minimum size of the contention window 𝑃 (𝑖,𝑗,𝑘,𝑢) (𝑙,𝑚,𝑞,𝑣) Transition probability from state (𝑖, 𝑗, 𝑘, 𝑢) to state (𝑙, 𝑚, 𝑞, 𝑣) 𝑓 (𝑘) (𝜇s) Required time to send k DL probe aggregated packets 𝑔 (𝑘) (𝜇s)

Table 2 :

 2 BTF estimations for ns-3 simulations where the MCS15 (144.4 Mbps) is used for all the nodes

	Ground Truth BTF		Estimated BTF
		0	0.125 0.25 0.375	0.5	0.625
	0	S/M
	0.125		S/M
	0.25		S/M
	0.375		S/M
	0.5			S/M
	0.625				S/M

Table 3 :

 3 BTF estimations for ns-3 simulations where different MCS (MCS15 and MCS13) are used for the AP and the probe traffic client

	Ground Truth BTF			Estimated BTF
		0	0.125 0.25 0.375	0.5	0.625
	0	S/M
	0.125		S	M
	0.25			S/M
	0.375			S/M
	0.5				S/M
	0.625				S/M

A FRAME TRANSMISSION DURATION

In this paper, the frame transmission duration takes into account the data packet transmission itself and all protocol overheads (backoff, SIFS, DIFS and ACK or block ACK transmission). The transmission duration of 𝑙 sub-frames is expressed as :

𝑃ℎ𝑦𝑠𝑖𝑐𝑎𝑙 𝑡𝑟𝑎𝑛𝑠𝑚𝑖𝑠𝑠𝑖𝑜𝑛 𝑟𝑎𝑡𝑒 Where freq denotes the frequency of sending of block ACK requests. 𝑇 𝑏𝑎𝑐𝑘𝑜 𝑓 𝑓 denotes the mean duration of the backoff period for frames without retransmissions. It is given by:

𝑇 𝐵𝑙𝑜𝑐𝑘𝐴𝐶𝐾 and 𝑇 𝐴𝐶𝐾 denote respectively the time required to send the block acknowledgment frame and the acknowledgment frame. They also count their physical header. Note that the only difference between the functions 𝑓 (.), ℎ(.) and 𝑔(.) is the physical transmission rate and the packet size that depend on the local properties of the nodes.

B ADDITIONAL PROBABILITIES


K Pantzas

G Beaudoin

M Bailly

A Martin

A Grisard

D Dolfi

O Mauguin

L Largeau

I Sagnes

G Patriarche

Relaxation mechanism of GaP grown on 001 Si substrates: influence of defects on the growth of AlGaP layers on GaP/Si templates

HAL is

Introduction

The (Al, Ga)P system possesses a significant number of properties highly desirable for nanophotonics integrated on Si. A non-exhaustive list of these properties includes: lattice mismatches of only 0.36 percent between GaP and Si [START_REF] Soga | Dislocation generation mechanisms for GaP on Si grown by metalorganic chemical vapor deposition[END_REF], the smallest among all III-V compounds and Si, and 0.3 % between AlP and GaP, allowing for the growth of thick, high-quality layers and Bragg mirrors [START_REF] Beyer | Influence of crystal polarity on crystal defects in gap grown on exact Si (001)[END_REF][START_REF]Ioffe NSM archive[END_REF]; a refractive index difference between AlP and GaP as large as 0.4 at 700 nm and 0.57 at 500 nm. Only fifteen AlP/GaP pairs are, therefore, required to yield a peak reflectivity of 98 % for wavelengths between 530 nm to 690 nm [START_REF] Hestroffer | Transport properties of doped AlP for the development of conductive AlP/GaP distributed bragg reflectors and their integration into light-emitting diodes[END_REF]; band-gap energies of 2.5 eV and 2.24 eV for AlP and GaP, respectively, making both materials transparent to the green-to-red region of the visible spectrum and allowing active structures to be optically pumped by conventional green lasers; high second-order nonlinearity, Pockels, and piezoelectric effects [START_REF]Ioffe NSM archive[END_REF][START_REF] Aitchinson | The nonlinear optical properties of AlGaAs at the half band gap[END_REF]; high thermal conductivity, 110 W K -1 m -1 , twenty times higher than that of active region material InGaP [START_REF] Colman | Temporal solitons and pulse compression in photonic crystal waveguides[END_REF], and close to that of Si.

Despite both the low lattice and thermal mismatch between Si and the (Al, Ga)P system, their integration is not straightforward. Indeed, Si crystallizes in the diamond, non-polar structure, while (Al,Ga)P materials crystallize in the polar zinc-blende structure. This difference in polarity leads to the apparition of anti-phase domains (APDs). The presence of such APDs is well-documented in the literature [START_REF] Beyer | GaP heteroepitaxy on Si(001): Correlation of Si-surface structure, gap growth conditions, and Si-III/V interface structure[END_REF][START_REF] Guillemé | Antiphase domain tailoring for combination of modal and 4-quasi-phase matching in gallium phosphide microdisks[END_REF][START_REF] Lucci | Universal description of iii-v/Si epitaxial processes[END_REF], and significant effort has been directed towards their suppression [START_REF] Beyer | Influence of crystal polarity on crystal defects in gap grown on exact Si (001)[END_REF][START_REF] Kunert | Si (001) surface preparation for the antiphase domain free heteroepitaxial growth of GaP on Si substrate[END_REF]. As a result, GaP on quasi-exact (0 0 1) Si templates are now commercially available. According to the manufacturer's website, these templates are free of dislocations, stacking-faults and possess an APD-free surface with excellent morphology [START_REF]Official NAsP website[END_REF]. 1 In the present contribution, the use of these templates for the growth of GaP and AlGaP using metalorganic chemical vapor deposition (MOCVD) is investigated. Regrowth on these templates typically commences with a hightemperature annealing step in a phosphine ambient to remove the oxide that forms on top of the template. In the conditions developed by the present group of authors for the growth of high-quality, defect-free GaP and Al 0.7 Ga 0.3 P layers [START_REF] Xie | Nonlinear, response of a gallium phosphide nanopatterned photonic waveguide in the CW regime[END_REF], this annealing step as well as the subsequent growth of AlGaP are carried out at 850 • C. At this temperature, the lattice mismatch between GaP and Si increases to 0.53 %. While the mismatch is still low, the critical layer thickness at this temperature according to the model of Matthews and Blakeslee [START_REF] Matthews | Defects in epitaxial multilayers: I. misfit dislocations[END_REF] is a little over 100 nm and is bound to be exceeded during growth. The GaP layer is, therefore, expected to plastically relax strain through the generation of misfit dislocations.

The objectives of the present paper are, therefore, to assess the mechanical stability of the GaP/Si template during high-temperature annealing and subsequent growth, and to evaluate the impact of any defects induced by plastic relaxation of misfit strain on the growth of AlGaP layers. The paper is structured as follows: a succinct description of the GaP/Si templates and the experiments that were carried out is given first. The stability of GaP/Si templates during annealing and the efficiency of the annealing in removing the native oxide are then assessed for an array of temperatures. Finally, the impact that defects induced by the annealing have on the subsequent growth of GaP is discussed.

Experiment

A 200 mm GaP on Si template was purchased from NAsP III/V GmbH [START_REF]Official NAsP website[END_REF].

This template consists of 50 nm of GaP on 001 Si. The Si substrate is 1 mm thick and is almost exact, with a small intentional miscut of 0.18 • in [1 1 0]. According to Reference [START_REF] Beyer | Influence of crystal polarity on crystal defects in gap grown on exact Si (001)[END_REF], this miscut, in conjunction with specific growth conditions, allows one to favor one crystal polarity over the other, resulting in a GaP/Si template with no observable APDs emerging at the GaP surface. Furthermore, given the small mismatch between GaP and Si -only 0.36 % at room temperature -the GaP layer is pseudomorphically accommodated on the Si template and, therefore, is under tensile strain and contains very few dislocations. AFM images of the surface of GaP as received are shown in Figure 1. The surface exhibits step flow, with a root mean square (RMS) roughness of 0.2 nm and no visible defects in areas as large as 10 µm by 10 µm.

Cross-sections were prepared for observation in a transmission electron microscope (TEM). The cross-sections for transmission electron microscopy were prepared from some of the samples using Focused Ion Beam (FIB) ion milling and thinning. Prior to FIB ion milling, the sample surface was coated with 50 nm of carbon to protect the surface from the platinum mask deposited used for the ion milling process. Ion milling and thinning were carried out in a FEI SCIOS dual-beam FIB-SEM. Initial etching was performed at 30 keV, and final polishing was performed at 5 keV. All TEM cross-sections were observed in an aberration-corrected FEI TITAN 200 TEM-STEM operating at 200 keV. The grid has an approximate period of 1.5 µm. The presence of this grid of dislocations, indicates the onset of plastic relaxation in the GaP/Si template, despite the small mismatch and relatively thin GaP layer. No threading dislocations are visible in the lamella, however. Given that the cross-sections are 25 µm long by 80 nm wide, the threading-dislocation density is less than 5 × 10 7 cm -2 .

[Figure 1 about here.] The 200 mm template was diced into 50 mm wafers using laser cutting. Annealing and growth was carried out in a Veeco D180 Turbodisk reactor. The annealing experiments were carried out at three different annealing temperatures of: 645 • C, 745 • C, and 845 • C, for 30 min and under a phosphine flow of 300 sccm. The lowest temperature is close to the one disclosed in Reference [START_REF] Beyer | Influence of crystal polarity on crystal defects in gap grown on exact Si (001)[END_REF], discussing the fabrication of the templates. The highest temperature is similar to one used for the obtained high quality GaP and AlGaP in Reference [START_REF] Xie | Nonlinear, response of a gallium phosphide nanopatterned photonic waveguide in the CW regime[END_REF]. In addition to the annealing experiments, two samples were grown on the GaP/Si templates using MOCVD: a 200 nm GaP layer to test GaP growth, and, separately, the following structure to test AlGaP growth: a nominally 325 nm thick GaP buffer followed by two periods consisting of a 47 nm thick AlGaP layer, containing 73 at. % aluminum, and a 64 nm thick GaP layer. Additionally, the 80 nm of the GaP buffer contain five 1 nm thick AlP markers, spaced every 16 nm. Both samples were grown at a reactor temperature of 845 • C and a reactor pressure of 20 Torr, using trimethylaluminum, trimethylgallium, and phopshine as precursors to elementary aluminum, gallium, and phosphorous, respectively.

[Figure 2 about here.]

Results

The four AFM images in Figure 3 show the evolution of the surface morphology of the GaP/Si template as a function of temperature. At the lowest annealing temperature, 645 • C, the GaP layer preserves its surface morphology, exhibiting step-flow and no visible defects. Nevertheless, the surface still exhibits spheres of oxide, indicating that the annealing temperature was not sufficient to completely remove the 3 nm oxide present at the surface of the template. Annealing for an hour at the same temperature did not yield any different results (not shown here). At higher temperatures, the oxide is completely removed and the step-flow morphology is still present, but the surface also exhibits stripes that are between 250 nm and 500 nm in width. At the highest annealing temperature, 845 • C, these stripes are more markedly present, creating a 5 nm high step. Additionally, several dislocations emerge at the surface of GaP, proof that the GaP layer plastically relaxes during the thermal annealing. The estimated threading-dislocation density is 2 × 10 8 cm -2 , at least an order of magnitude higher than that present in the template initially.

[Figure 3

about here.]

The defects observed on the annealed samples persist in the 200 nm GaP layer grown on the GaP template, as shown the AFM images of the sample surface in Figure 4. The threading dislocation density, estimated from the density of pits in the AFM images, is 2 × 10 8 cm -2 , no different than that observed for the sample at 845 • C. The GaP surface is also rougher owing to the onset of cross-hatch, with an RMS roughness of 0.9 nm, although step-flow is still preserved.

[Figure 4 about here.] The defects observed in the GaP sample significantly perturb the growth of AlGaP, as revealed by the dark-field TEM images of the AlGaP sample, shown in Figure 5. These images reveal that the stripes observed already during annealing are bordered by pairs of micro-twins. When these cross the AlGaP layer, they induce severe faceting, mainly [1 1 2] and [1 1 4] planes. This faceting is mainly visible at the AlGaP/GaP interface. Subsequent growth of GaP tends to smooth this roughness, as seen at the second GaP/AlGaP interface in Figure 5, but the surface becomes again rougher during the growth of the second AlGaP layer.

The RMS roughness at the surface is 5 nm, more than an order of magnitude higher than the initial roughness of the template.

[Figure 5 about here.] Finally, the GaP/Si interface now exhibits a dense misfit-dislocation grid, with an approximate period of 30 nm. Upon close inspection, this dislocation grid is composed up to 60 % of edge dislocations, and 40 % partial dislocations.

Half of the partial dislocations at the interface combine with a Lomer-Cottrell dislocation in Si to form a stair-rod dislocation [START_REF] Hull | Introduction to Dislocations[END_REF]. This type of defect is more common in highly mismatched systems, such as Ge/Si that has a lattice mismatch of 4 %, an order of magnitude higher than GaP/Si [START_REF] Arroyo | Effect of thermal annealing on the interface quality of Ge/Si heterostructures[END_REF], and to the extent of our knowledge has not been reported for lowly mismatched systems.

The stair-rod itself appears as an inverted triangle at the GaP/Si interface. One such defect is shown in Figure 6 (a). Two apexes of the the triangle are partial misfit dislocations at the GaP/Si interface. Tracing the Burgers' circuits of the partial dislocations -Figures 6 (b) and (c) -reveals that their Burgers' vector is b = 1 6 a 112 , in keeping with theory. Stacking faults extend from them and slip along [1 1 1] planes in Si until they meet at about 30 nm from the GaP/Si interface. There they combine to form the Lomer-Cottrell dislocation, that has a Burgers' vector of b = 1 2 a 110 and is sessile.

[Figure 6 about here.]

Discussion

Both AFM and TEM results show that the 50 nm GaP layer in the template does not contain any emerging APDs, and presents a smooth surface with a stepflow morphology and a roughness of only 0.2 nm. Nevertheless, the TEM images reveal that even at this thickness, which is lower than the critical thickness for plastic relaxation of GaP, reported at 90 nm [START_REF] Soga | Initial stage of GaP/Si heteroepitaxial growth by metalorganic chemical vapor deposition[END_REF], a loose misfit-dislocation grid is already present at the GaP/Si interface. In the literature, this has been associated with thermal-stress induced relaxation of GaP during the cooling down phase of the GaP/Si template following GaP growth at temperatures higher than 400 • C. The reported threading dislocation densities in optimized conditions were of the order of 1 × 10 7 cm -2 [START_REF] Tachikawa | Dislocation generation of GaAs on Si in the cooling stage[END_REF][START_REF] Olson | MOCVD Growth and Characterization of GaP on Si[END_REF], in agreement with the estimate obtained from TEM on the reference sample.

In two samples grown by MOCVD, the total GaP thickness is more than 250 nm, exceeding the critical thickness. In the literature full plastic relaxation of misfit strain that is due to lattice mismatch between GaP and Si is not achieved until a thickness of 200 nm [START_REF] Soga | Initial stage of GaP/Si heteroepitaxial growth by metalorganic chemical vapor deposition[END_REF][START_REF] Olson | MOCVD Growth and Characterization of GaP on Si[END_REF][START_REF] Emmer | Fabrication of Single Crystal Gallium Phosphide Thin Films on Glass[END_REF]. At the reported thickness the samples are, therefore, fully plastically relaxed, and, accordingly, exhibit dislocation densities in the 1 × 10 8 cm -2 range. The plastic relaxation of the misfit strain is further supported by the appearance of Lomer-Cottrell locks in Si. Indeed, the appearance of misfit dislocations in the Si substrate has been previously associated with the relaxation of strain that is due to the lattice mismatch, whereas thermal stress tends to generate defects at the GaP/Si interface or in GaP [START_REF] Tachikawa | Dislocation generation of GaAs on Si in the cooling stage[END_REF]. Although the RMS roughness of GaP sample increases, it still remains within reasonable limits, i.e. less than 1 nm. This roughness could be improved by increasing the V/III ratio and/or the reactor pressure.

The most cumbersome defects, however, are the micro-twins that already appear during high-temperature annealing of the GaP/Si template. These defects induce heavy faceting in the AlGaP layer, rendering these layers inappropriate for the growth of AlGaP/GaP Bragg mirrors. The micro-twins reported here are denser when the annealing temperature is increased. This result is contrary to what has been observed in the case of GaP grown on 4 • off-cut Si using migration-enhanced molecular-beam epitaxy [START_REF] Yamane | Growth of low defect density GaP layers on Si substrates within the critical thickness by optimized shutter sequence and post-growth annealing[END_REF]. Indeed, in that reference, micro-twins are eliminated by implementing an annealing step after the initial nucleation of a 20 nm GaP layer. It is unclear whether the miscut is solely responsible for this difference, or whether the micro-twins could be annihilated by further tweaking of the growth and annealing conditions on the current templates. Alternatively, growth could be performed in a temperature range closer to what has been reported by NAsP III/V Gmbh., i.e. between 450 • C and 600 • C

[2], using precursors that are better suited to low-temperature growth, such as triethylgallium and tri-tert-butylphosphine. Nevertheless, incorporation of high amounts of aluminum at low temperatures may prove challenging.

Conclusion

Despite excellent surface properties and the absence of emerging APDs, the GaP/Si templates used in the present contribution are revealed not to be compatible with the growth conditions developed for high-quality AlGaP and GaP growth on GaP substrates. Indeed, the templates are unstable at high temperatures, plastically relaxing to form micro-twins that detrimentally affects the growth of AlGaP in particular. Furthermore, the oxide layer present at the surface of the templates is not readily removed through annealing in a phosphine ambient at lower temperatures. A combination of wet etching to remove the oxide and growth with alternative metalorganic precursors that are better suited to low-temperature growth may be the only solution towards growth of high quality GaP and AlGaP on these templates. . This undulation is due to crosshatch. Crosshatch in GaP has been shown to occur anisotropically, more dense along the 0 1 0 direction than the 1 0 0 direction [START_REF] Lucci | Universal description of iii-v/Si epitaxial processes[END_REF]. The surface is dotted by small pits where threading dislocations emerge, pointed out by white arrows in image (b). The threading dislocation density is estimated at 1 × 10 8 cm -2 , at least four times higher than the estimated density of the initial template.

List of Figures

In addition to the emergence pits and the crosshatch, at least one 500 nm wide band delimited by two lines is pointed out by arrows in image (a), indicating that the bands that appear during annealing persist throughout subsequent growth of GaP. . . . . . 17 . This undulation is due to crosshatch. Crosshatch in GaP has been shown to occur anisotropically, more dense along the 0 1 0 direction than the 1 0 0 direction [START_REF] Lucci | Universal description of iii-v/Si epitaxial processes[END_REF]. The surface is dotted by small pits where threading dislocations emerge, pointed out by white arrows in image (b). The threading dislocation density is estimated at 1 × 10 8 cm -2 , at least four times higher than the estimated density of the initial template. In addition to the emergence pits and the crosshatch, at least one 500 nm wide band delimited by two lines is pointed out by arrows in image (a), indicating that the bands that appear during annealing persist throughout subsequent growth of GaP.

A

  high-angle annular dark-field image of the GaP/Si template is shown in Figure 2 (a). The image shows the presence of a 3 nm thick oxide layer on the surface of the 50 nm thick GaP layer. At this scale no defects are present in the GaP layer. Furthermore, 002 Dark-field TEM images of the samples confirm that antiphase domains are indeed buried and do not emerge at the surface of the template -see Figure 2 (b). Nevertheless, misfit dislocations can be seen at the GaP/Si interface, as shown in the 220 Dark-field TEM image in Figure 2 (c).

1 (a) 100 µm 2

 2 and (b) 4 µm 2 AFM images showing the surface of GaP on a GaP/Si template procured from NAsP III/V GmbH.Step-flow morphology is observed, with the step-flow occurring along the 001 direction. The RMS roughness is 0.2 nm. . . . . . . 14 2 (a)HAADF-STEM image of the as-received GaP/Si template. The image shows that the 50 nm thick GaP layer is covered by a 3 nm thick oxide layer. (b) 002 Dark-field TEM image showing the buried antiphase domains. In all the samples tested none of the domains emerge on the surface. (c) 220 Dark-field TEM image showing dislocations (black arrows) spaced by 1.5 µm along the GaP/Si interface, indicating that despite the low mismatch and the relatively thin GaP layer, plastic relaxation has set on in the GaP/Si template . . . . . . . . . . . . . . . . . . . . . . . . . 15 3 25 µm 2 AFM images of the surface of the GaP/Si template (a) before, and after annealing at (b) 645 • C, (c) 745 • C, and (d) 845 • C for 30 min. At 645 • C the atomic steps are preserved, but the surface is riddled with small pockets of oxide that was not removed. At higher annealing temperatures, these pockets disappear and the step-flow morphology is preserved. At 745 • C and 845 • C, however, bands that are on average 375 nm wide appear to rise above from the background. These are particularly visible in the sample annealed at 845 • C. At this temperature the elevation of these bands is 5 nm. These bands are indicative of plastic relaxation taking place in the GaP layer of the template during high temperature annealing. . . . . . . . . . . . . . . . . . . . . . 16 4 (a) 100 µm 2 and (b) 4 µm 2 AFM images of the GaP surface after the growth of a 200 nm of GaP on the GaP/Si template using MOCVD. Although step flow is still observed, the surface undulates along the [1 1 0]

  (a) 220 dark-field TEM image of the AlGaP/GaP structure grown on the GaP template. A dense grid of misfit dislocation is present at the GaP/Si interface, indicating that the structure is fully plastically relaxed. In addition to the dislocation grid, micro-twins appear in the structure. These often appear in pairs spaced by 365 nm, as revealed in (b) showing a 220 dark-field TEM image close up of one such pair. These micro-twins induce faceting in AlGaP, significantly roughening the GaP/AlGaP interfaces. This faceting is highly detrimental to the optical properties of AlGaP/GaP pairs. As shown in Figure3, these micro-twins already appear during the annealing of the GaP/Si template. . . . 18 6 (a) HAADF-STEM image of a the GaP/Si interface in the Al-GaP sample showing a stair-rod dislocation. The stair-rod forms at the junction of the stacking faults (SF 1 and SF 2 ) that extend from two partial dislocations (PD 1 and PD 2 ) at the GaP/Si interface and meet in the Si substrate, forming a sessile Lomer-Cottrell dislocation. (b,c) Insets showing magnified portions of the HAADF-STEM image in (a), centered around PD 1 and PD 2 , and their Burgers' circuits. The polarity of GaP is not inverted within the domain of the Lomer-Cottrell lock. . . . . . . . . . . . 19

Figure 1 :

 1 Figure 1: (a) 100 µm 2 and (b) 4 µm 2 AFM images showing the surface of GaP on a GaP/Si template procured from NAsP III/V GmbH. Step-flow morphology is observed, with the stepflow occurring along the 001 direction. The RMS roughness is 0.2 nm.

Figure 2 :Figure 3 :

 23 Figure 2: (a)HAADF-STEM image of the as-received GaP/Si template. The image shows that the 50 nm thick GaP layer is covered by a 3 nm thick oxide layer. (b) 002 Dark-field TEM image showing the buried antiphase domains. In all the samples tested none of the domains emerge on the surface. (c) 220 Dark-field TEM image showing dislocations (black arrows) spaced by 1.5 µm along the GaP/Si interface, indicating that despite the low mismatch and the relatively thin GaP layer, plastic relaxation has set on in the GaP/Si template

Figure 4 :

 4 Figure 4: (a) 100 µm 2 and (b) 4 µm 2 AFM images of the GaP surface after the growth of a 200 nm of GaP on the GaP/Si template using MOCVD. Although step flow is still observed, the surface undulates along the [1 1 0]. This undulation is due to crosshatch. Crosshatch in GaP has been shown to occur anisotropically, more dense along the 0 1 0 direction than the 1 0 0 direction[START_REF] Lucci | Universal description of iii-v/Si epitaxial processes[END_REF]. The surface is dotted by small pits where threading dislocations emerge, pointed out by white arrows in image (b). The threading dislocation density is estimated at 1 × 10 8 cm -2 , at least four times higher than the estimated density of the initial template. In addition to the emergence pits and the crosshatch, at least one 500 nm wide band delimited by two lines is pointed out by arrows in image (a), indicating that the bands that appear during annealing persist throughout subsequent growth of GaP.

Figure 5 :

 5 Figure 5: (a) 220 dark-field TEM image of the AlGaP/GaP structure grown on the GaP template.A dense grid of misfit dislocation is present at the GaP/Si interface, indicating that the structure is fully plastically relaxed. In addition to the dislocation grid, microtwins appear in the structure. These often appear in pairs spaced by 365 nm, as revealed in (b) showing a 220 dark-field TEM image close up of one such pair. These micro-twins induce faceting in AlGaP, significantly roughening the GaP/AlGaP interfaces. This faceting is highly detrimental to the optical properties of AlGaP/GaP pairs. As shown in Figure3, these micro-twins already appear during the annealing of the GaP/Si template.

Figure 6 :

 6 Figure 6: (a) HAADF-STEM image of a the GaP/Si interface in the AlGaP sample showing a stair-rod dislocation. The stair-rod forms at the junction of the stacking faults (SF 1 and SF 2 ) that extend from two partial dislocations (PD 1 and PD 2 ) at the GaP/Si interface and meet in the Si substrate, forming a sessile Lomer-Cottrell dislocation. (b,c) Insets showing magnified portions of the HAADF-STEM image in (a), centered around PD 1 and PD 2 , and their Burgers' circuits. The polarity of GaP is not inverted within the domain of the Lomer-Cottrell lock.

A detailed investigation of these claims is presented in Section

Acknowledgements

The authors would like acknowledge J. Nagle, A. De Rossi, and S. Combrié for their support and fruitful discussions. The authors thankfully acknowledge funding from the CNRS Renatech network and the ANR Labex TEMPOS (ANR-10-EQPX-0050).


Sandrine Charles
email: sandrine.charles@univ-lyon1.fr

Aude Ratier

Aude Ratier # Christelle Lopes

Generic Solving of One-compartment Toxicokinetic Models

Keywords: Ordinary differential equations, Environmental risk assessment, Living organisms, Active substances, Toxicokinetic models, Bioaccumulation metrics ADME, absorption, distribution, metabolism and excretion, EFSA, European food safety authority, ERA, environmental risk assessment, ODE, ordinary differential equation, PBPK, physiologically based pharmacokinetic, PBTK, physiologically based toxicokinetic, TK, toxicokinetic

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

Introduction

In this paper, we consider a generic one-compartment toxicokinetic (TK) model describing uptake and elimination processes within living organisms when exposed to chemical substances. From a general point-of-view, the term TK refers to the movement and fate (also called disposition) of chemical substances that are potentially toxic. In particular, this term is used when describing the time course of absorption, distribution and elimination (accounting for both biotransformation and excretion) of substances within organisms (ADME). Toxicokinetics is also related to pharmacokinetics and is often considered as the same scientific field, the main difference being linked to the type of compound (toxicants or pharmaceuticals). In comparison with pharmaceutical exposure, exposure to toxicants is often uncontrolled and variable over time predominantly at larger concentrations. [START_REF] Gupta | General Principles of Toxicology[END_REF] Mathematical TK models are particularly useful to predict the adverse effects of xenobiotics or to prevent undesired residues within animal tissues from entering the human food chain. [START_REF] Côté | A novel toxicokinetic modeling of cypermethrin and permethrin and their metabolites in humans for dose reconstruction from biomarker data[END_REF] Indeed, TK models allow the prediction of tissue concentrations over time, involving several kinetic parameters, such as uptake and elimination rates, from which useful metrics for assessing risk may be derived. [START_REF] Paini | Assessment of the predictive capacity of a physiologically based kinetic model using a read-across approach[END_REF] Additionally, TK models can be used to facilitate a better understanding of the underlying physiological mechanisms driving chemical ADME processes. 4, [START_REF] Lautz | Generic physiologically based kinetic modelling for farm animals: Part I. Data collection of physiological parameters in swine, cattle and sheep[END_REF] A wide variety of TK models exists, chosen from both the available data and the use intended for the model (e.g., describe and/ or predict); simple or more complex models may be preferred accordingly. [START_REF] Riviere | Veterinary Pharmacology and Therapeutics[END_REF] Basically, we distinguish classic TK models from physiologically-based (PB) TK models, referred to as PBTK, or PBPK (standing for PB Pharmaco-Kinetic) models in the field of pharmacology. As fully detailed in Gehring et al., [START_REF] Gehring | Veterinary-specific applications of pharmacokinetic modelling and theory: Opportunities and challenges[END_REF] PBPK models are a crucial tool in Environmental Risk Assessment (ERA) for regulatory bodies, such as the US Environmental Protection Agency (US EPA) and European Food Safety Authority (EFSA). 8, [START_REF] Mclanahan | Physiologically based pharmacokinetic model use in risk assessment-Why being published is not enough[END_REF] In fact, PBPK models are an ethical and scientifically sound method to predict exposure to toxic xenobiotics in humans through animalto-human extrapolation or based on human biomonitoring data. (PB)TK models are also of particular interest for ERA, when the calculation of bioaccumulation metrics is required by regulators to evaluate dossiers for substance marketing authorizations. [START_REF]Test No. 305: Bioaccumulation in Fish: Aqueous and Dietary Exposure[END_REF][START_REF] Nyman | Importance of toxicokinetics for interspecies variation in sensitivity to chemicals[END_REF][START_REF] D'argenio | Analytical solution of linear multi-compartment models with non-zero initial condition and its implementation with R[END_REF] The generic TK model we unravel in this paper is a one-compartment model that considers organisms as a whole, in which chemical compounds may enter and from which these compounds can be eliminated. More complex TK models (namely PBTK models) refine the description of contamination pathways within organisms, distinguishing organs and tissues from physiological hypotheses on potential targets of exposure compounds; PBTK models have been mainly developed for humans 13,14 and big mammals. 15, [START_REF] Mavroudis | Development and validation of a physiology-based model for the prediction of pharmacokinetics/toxicokinetics in rabbits[END_REF] Recently, a PBTK model was proposed by Gestin et al. [START_REF] Gestin | One and multi-compartments toxico-kinetic modeling to understand metals' organotropism and fate in Gammarus fossarum[END_REF] for invertebrates, pushing the immense potential of these models a step further.

Our TK model, even if not fully refined as a PBTK model, has the great advantage of accounting for several exposure sources (e.g., water, sediment, and/or food), several elimination processes (e.g., direct elimination, dilution by growth and/or biotransformation), and several potential metabolites of the parent chemical compound to which organisms are exposed. [START_REF] Ratier | A Bayesian framework for estimating parameters of a generic toxicokinetic model for the bioaccumulation of organic chemicals by benthic invertebrates: Proof of concept with PCB153 and two freshwater species[END_REF] The next section details the dynamic system of ordinary differential equations (ODEs) and the mathematical exact solution, followed by a summary of the full set of ODE in section 3. Then, section 4 illustrates the use of model based on simulations from parameter values obtained by the study of real datasets.

One-compartment TK set of ODE

Our generic one-compartment TK model is composed of two sets of ODE: accumulation and depuration sets. The accumulation phase includes both absorption and elimination processes, as those related to the biotransformation of the parent compound into metabolites, during which organisms are exposed to a given compound. The depuration phase concerns only elimination processes, including biotransformation, during which organisms are transferred into a clean medium. The transition from one set of ODE to the other takes place at time t c corresponding to the duration of the accumulation phase (see Table 1).

In practice, the exposure concentration to which organisms are exposed may vary over time in real environments. This means

= U R E M + NA = ( ) e D k E M -+   NA, ∀ℓ = 1…L ( )
( ) 0 = 1 E M t c Q C R e + - - NA
# stands for numbers, while NA means 'Not Applicable'

there is no analytical solution of the TK model, so that only a numerical solution can be obtained with an appropriate algorithm. Therefore, our work assumes that the exposure concentration remains constant over time regardless of the exposure source. Such an experimental condition can be ensured for most chemical compounds when performing laboratory experiments. In addition, this assumption provides the exact solution of the full one-compartment TK model by considering as many routes of exposure and as many elimination processes as desired, as well as an infinite number of phase I metabolites, i.e., directly derived from the parent compound to which organisms are exposed. Kuo et al. [START_REF] Kuo | Deriving in vivo biotransformation rate constants and metabolite parent concentration factor/stable metabolite factor from bioaccumulation and bioconcentration experiments: An illustration with worm accumulation data[END_REF] previously suggested a partially resolved one-compartment TK model but only for the accumulation phase and one exposure route. In a similar manner, D'Argenio et al. [START_REF] D'argenio | Analytical solution of linear multi-compartment models with non-zero initial condition and its implementation with R[END_REF] published an analytical solution for a linear multi-compartments TK model with a non-zero initial condition.

Accumulation phase (0 ≤ t ≤ t c )

The set of ODEs describing the accumulation phase (0 ≤ t ≤ t c ) is as follows:

dC t dt U E M C t dC t dt k t k C t p p m m p e m ( ) ( ) ( ) ( ) ( ) ( ) = -+ × = × -× ∀ (a) , C = (b) 1 L              (1)
All parameters and variables, with their meaning and units when applicable, are provided together in Table 1.

Equation (1.a) for the parent compound is a linear first-order ODE with constant coefficients and a second member. Equation (1.a) admits

( ) = = part U C t R E M +
as a particular solution. Equation (1.a) without its second member is:

( ) ( ) ( ) ( ) = 0 ( ) = , p E M t p p dC t E M C t C t K e K dt  -+ × + -+ × ⇔ × ∈ (2)
Given the initial condition C (t = 0) = C 0 (C 0 ≥ 0), we can obtain the full analytical solution of Equation (1.a), providing that the internal concentration of the parent compound over time during the accumulation phase (0 ≤ t ≤ t c ) is as follows:

( )

( ) 0 ( ) = E M t p C t C R e R -+ × - × + (3)
See Table 1 for the definition of parameter R. Equation (1.b) is also a linear first-order ODE with constant coefficients and a second member, with the following analytical solution when the second member is removed:

( ) ( ) = 0 ( ) = , , = 1
k t m e e m m dC t k C t C t K e K L dt - × + -× ⇔ × ∈ ∀         (4)
The method of variation of a constant consists of writing the general solution of Equation (1.b) as:

( ) = ( ) k t e m C t K t e - × ×   ( 5
)
which is used to find function K(t) by deriving ( )

m C t 
then insert the result into Equation (5). The derivative of Equation (1.b) then becomes:

( ) ( ) = ( ) k t k t m e e e dC t dK t e K t k e dt dt - × - × × - × ×     (6)
Inserting the result from the previous equation into Equation (1.b)

leads to: 0 ( ) = ( ) k t D t e m dK t k C R e R e dt × ×   × -× + ×        (7) which integrates into: 0 ( ) = , k t D t e m e C R R K t k e e C C D k × ×   -   × × + × + ∈           (8)
See Table 1 for the definition of parameter D ℓ . The general solution of Equation (1.b) finally becomes:

( ) 0 ( ) = , k t e E M t m m e C R R C t k e C e C D k - × -+ ×   -   × × + + × ∈           (9)
When considering the initial condition

( = 0) = 0 m C t 
when the accumulation phase starts, organisms are only exposed to the parent compound, so that there are no metabolites. With this condition, we can finally acquire the full analytical solution of Equation (1.b), providing the internal concentration of metabolite ℓ over time during the accumulation phase:

( ) 0 ( ) = 1 k t e E M t m m k t e e C R C t k e e D R e k - × -+ × - ×  -   × × -           + × -             ( 10
)

Depuration phase (t > t c )

The set of ODEs describing the depuration phase (t > t c ) is as follows:

dC t dt E M C t dC t dt k t k C t p p m m p e m ( ) ( ) ( ) ( ) ( ) ( ) ( ), = -+ × = × -× ∀ = a C 1 1 b L ( )              (11)
The definitions and units (when applicable) of the parameters and variables are described together in Table 1.

Equation (11.a) is a linear first-order ODE without a second member, so that it has a general solution of the following form:

( ) ( ) = , E M t p C t K e K  -+ × + × ∈ ( 12
)
For the depuration phase and the parent compound, the initial condition comes from the calculation of the internal parent compound concentration at the end of the accumulation phase (i.e., at time t = t c ) using Equation ( 13):

( ) 0 ( ) = ( ) E M t c p c C t C R e R -+ × -× + ( 13
)
From the general analytical solution given by Equation ( 12), we get

( ) ( ) = E M t c p c C t Ke -+
, from which constant K in Equation ( 12) is determined by:

( ) 0 = E M t c K C R R e + × -+ × (14)
Then, the final analytical solution of Equation (11.a), providing the internal concentration of the parent compound over time during the depuration phase (t > t c ), becomes:

( ) ( ) ( ) ( ) 0 ( ) = 1 E M t E M t c p C t C R e e + × -+ × -× - × (15)
For simplicity reasons, Equation ( 15) can be written as C p (t) = Qe -(E+M)×t , where Q is defined in Table 1.

Equation (11.b) is a linear first-order ODE with constant coefficients and a second member. The analytical solution of Equation (11.b) without its second member is:

( ) = , k t e m C t K e K - × + × ∈    ( 16
)
As previously stated, the method of the variation of a constant provides the general solution of Equation (11.b) as

( ) = ( ) k t e m C t K t e -  
, from which function K(t) must be determined.

The derivative of C m ℓ (t) from Equation ( 16) is:

( ) = ( ) k t k t e e e dK t e k K t e dt - × - × -× ×    ( 17
)
Inserting Equation (17) into Equation (11.b) leads to:

( ) = D t m dK t k Q e dt × × ×   ( 18
)
which integrates into:

( ) = , D t m Q K t k e C C D × × × + ∈     ( 19
)
finally leading to the general analytical solution of Equation (11.b):

( ) ( ) = , k t e E M t m m Q C t k e C e C D - × -+ × × × + × ∈      (20)
Constant C can be determined from the initial condition, i.e., from the internal concentration of metabolite ℓ at t = t c both at the end of the accumulation phase and at the beginning of the depuration phase. From Equation (20), we get

( ) ( ) = k t E M t e c c m c m Q C t k e Ce D - -+ +    
, and from Equation (10), we obtain

( ) 0 ( ) = 1 k t k t E M t e c e c c m m e C R R C t k e e e D k - - -+   -      - + -                  
. Finally, the expression for constant C can be determined using:

0 = 1 k t k t e c e c m e R C R R C k e e k D D × ×   -     × × -- -×               (21)
Replacing constant C in Equation (20) gives the final analytical solution of Equation (11.b), which yields the internal concentration of metabolite ℓ for the depuration phase (t > t c ) as follows:

( ) ( )
( ) 0 ( ) = ( ) k t t k t e c e E M t m m e k t k t t e e c Q R C t k e e e D k C R R e e D D - × - - × -+ × - × - × -   × × + × -     - - × -×           (22)
Replacing constant Q by its own expression (Table 1) leads to:

( ) 0 ( ) ( ) ( ) ( ) ( ) = ( k t e E M t
+ × - × - - × - × - -+ × - -   × × -       + × -       + × -              (23)

Generic set of solutions

Considering the following intermediate notations = U R E M +

and D ℓ = k e ℓ -(E + M), we can obtain the full set of analytical solutions corresponding to the whole one-compartment TK set of ODE describing the time-course for the accumulation and depuration phases of the parent compound and its potential metabolites.

• The analytical solution for the internal concentration of the parent compound during the accumulation phase, previously referred as Equation ( 3), becomes:

( ) 0 ( ) =
, 0

E M t p c U U C t C e t t E M E M -+ ×   - × + ∀   + +    
• The analytical solution for the internal concentration of metabolite ℓ during the accumulation phase, previously referred as Equation (10), is:

( ) 0 ( ) = 1 ( ) , 0 k t m m e m e e k t e E M t c k k U C t e k E M k E M U C e e t t E M - × - × -+ ×   × × - +   + -+       × - × - ∀     +             
• The analytical solution for the internal concentration of the parent compound during the depuration phase, previously referred as Equation ( 15), is written as:

( )

( ) ( ) 0 ( ) = 1 , > E M t E M t c p c U C t C e e t t E M + × -+ ×   - × - × ∀   +  
• The analytical solution for the internal concentration of metabolite ℓ during the depuration phase, previously referred as Equation ( 23), is determined using:

( ) 0 ( ) ( ) ( ) ( ) ( ) = ( ) , > ( ) k t m e E M t
- × -+ × - × - - × - × - -+ × -     × - × -     -+ +       + × × -   +     + × × - ∀   -+ +             
Then, the final generic analytical solution of the one-compartment TK model can be achieved to estimate all parameters from observed data using an inference process by replacing constants U, E and M with

=1 =1 =1 = , = and = I J L u i e m i j i j U k c E k M k × ∑ ∑ ∑   , respectively.

Model simulations based on real case-studies

Currently, the above exact mathematical solution of the generic onecompartment TK model is used in support of the development of the R-package rbioacc 20 as well as its web interface MOSAIC bioacc , which is freely available at https://mosaic.univ-lyon1.fr/bioacc. [START_REF] Ratier | New perspectives on the calculation of bioaccumulation metrics for active substances in living organisms[END_REF] Both tools allow fitting the TK model to bioaccumulation data in a user-friendly way, providing all necessary outputs to check for goodness-of-fit and obtain bioaccumulation metrics, regardless of the considered speciescompound combination. The database associated to the web interface is also available online at http://lbbe-shiny.univ-lyon1.fr/mosaicbioacc/data/database/TK_database.html, with a wide collection of accumulation-depuration data in support of toxicokinetic modelling. [START_REF] Ratier | Accumulation-depuration data collection in support of toxicokinetic modelling[END_REF] From this database, we selected three real case-studies from which we got fitting results (especially the joint posterior probability distribution of all kinetic parameters) in order to perform the subsequent model simulations, thus illustrating the use of the generic TK model. All model simulations of the generic set of mathematical solutions (Equations ( 3), ( 10), ( 15) and ( 23)) were performed using R software with the rbioacc package and function predict(). [START_REF] Ratier | rbioacc: an R-package to analyse toxicokinetic data[END_REF] Note that the same simulations can be performed directly online with the MOSAIC bioacc platform and its prediction module (https://mosaic.univ-lyon1.fr/bioacc). We used 500 time points in [0; t f ], where t f stands for the final time of each simulation. These simulations illustrate the use of the generic TK model in three case studies where different compounds are bioaccumulated by different species. In order to proceed, the time duration of the accumulation phase (parameter t c ) is required, as well as the exposure concentrations in the media and the model parameter values (Table 2); all associated numerical values were taken from fitting results performed with real data sets. Figures 1 to 5 show the simulations of internal concentrations over time for the three speciescompound combinations, as subsequently described. For each case study, TK model parameters were varied one-at-a-time with 20, 50, and 80% of increase from their original value.

Elementary one-compartment TK model

Equations ( 3) and (15) were simulated for a simple case study in which fish are exposed to a highly hydrophobic chemical contaminated water and only excretion is considered (Fig. 1 and Table 2). [START_REF] Lautz | Generic physiologically based kinetic modelling for farm animals: Part I. Data collection of physiological parameters in swine, cattle and sheep[END_REF] The corresponding inputs are the exposure concentration c w (referring to c i and I = 1), uptake rate from water u w k (referring to u i k and I = 1), and excretion rate e e k (referring to e j k and J = 1). When parameter u w k increases, internal concentrations become higher (e.g., blue curve in Fig. 1a) than with the original value (orange curve in Fig. 1a). Biologically speaking, the higher the uptake rate is for a given substance, the more it is bioaccumulated by organisms. Conversely, an increase in parameter e e k leads to a lower internal concentration (e.g., blue curve in Fig. 1b), which is consistent with the known underlying biological mechanism: the faster a contaminant is eliminated, the quicker its concentration decreases in organisms.

One-compartment TK model with several exposure routes, no metabolites

Equations ( 3) and ( 15) are simulated for case, in which freshwater shrimp are exposed to an organic chlorine compound through contaminated sediment and food, and only natural excretion is considered as the elimination process (Fig. 2 and Table 2). [START_REF] Crookes | Estimation of fish bioconcentration factor (BCF) from depuration data[END_REF] The corresponding inputs are the exposure concentration via sediment c s and food c f (referring to c i and I = 2), the two uptake rates from sediment and food, 2a) than with the original value (orange curve in Fig. 2a). Biologically, the higher the uptake rate from sediment is for a given substance, the stronger it is bioaccumulated by organisms. Regarding the exposure from food, as parameter u f k is low in this case, its variation has minimal influence on the internal concentration of the contaminant (Fig. 2b). This means that the exposure via sediment is the major route of contamination for these organisms. Besides, as previously shown in section 4.1, an increase in parameter e e k leads to a faster decreasing concentration (e.g., blue curve in Fig. 2c).

One-compartment TK model with one exposure route, several metabolites

Equations ( 3), ( 10), ( 15) and ( 23) are simulated for case, in which freshwater shrimp are exposed to an organic biocide by contaminated water. Three metabolites derived from the parent compound are considered [START_REF] Ashauer | Significance of xenobiotic metabolism for bioaccumulation kinetics of organic chemicals in Gammarus pulex[END_REF] together with the natural excretion (Figs. 3 to 5 and Table 2). The corresponding inputs are the exposure concentration via water c w (referring to c i and I = 1), the uptake rate u w k (referring to u i k and I = 1), excretion rate e e k (referring to e j k and J = 1), three metabolization rates k increases, internal concentrations of metabolite 1 increase (e.g., blue curve in Fig. 3b) greater than with the original value (orange curve in Fig. 3b). Conversely, internal concentrations are lower for the parent compound (e.g., blue curve in Fig. 3a). In other words, the more the biotransformation rate for a given metabolite increases, the higher its concentration becomes within organisms due to the highly biotransformation of the parent compound. This leads to a lower internal concentration of the parent compound than with the original value of 1 m k (e.g., blue curve, Fig. 3a). An increase in 1 m k also induces a decrease in the internal concentrations of the other metabolites (Figs. 3c andd). Besides, an increase in parameter 1 e k only affects the internal concentration of metabolite 1 (Fig. 4b). In addition, as previously viewed (sections 4.1 and 4.2), an increase in parameter u w k will induce high internal concentrations for both the parent compound and its metabolites (Fig. 5). Indeed, the more the internal concentration of the parent compound increases, the more the biotransformation process will intensify, leading to high internal concentrations for each metabolite.

Future directions

Based on the results of this work, we suggest future research aims to acquire a better understanding of the variation of TK parameters according to the chemicals and their exposure concentrations in the environment. In particular, the influence of concentration dependency may be investigated through the TK parameter values. For example, if the uptake rate remains on the same order of magnitude rather than increasing two-fold when the exposure concentration doubles, this would most likely indicate a concentration dependency. A second interesting research direction would be the elaboration of an improved modelling framework that includes phase II metabolites within the TK model; some experimental data are already available to support such research but no model exists yet.

Conclusions

As supported by complete mathematical and exact calculations, our work introduces a practical and generic one-compartment TK model to deal with different species-compound combinations of interest. To-date, such a model is required by regulatory bodies to perform risk assessment of chemical substances. The generic feature of this model is, thus, in full compliance with the numerous situations that regulators are expected to be faced. The associated ready-to-use tools we developed to perform different actions, including calibration of experimental data, validation of simulations compared to observed data, and predictions in untested situations, should be of pertinent use to such regulators.

  and

  to u i k and I = 2), and excretion rate e e k (referring to e j k and J = 1). When parameter u s k increases,

Fig. 1 .

 1 Fig. 1. Example of model simulations for a simple TK model (exposure via water and elimination by excretion) and the influence of variations of parameters (a) u w k (uptake rate) and (b) e e k (excretion rate). TK, toxicokinetic.

Fig. 2 .

 2 Fig. 2. Example of model simulations for a TK model with multiple exposure routes (via sediment and food) and influence of variations of parameters (a) u s k (uptake rate from sediment), (b) u f k (uptake rate from food) and (c) e e k (excretion rate). TK, toxicokinetic.

Fig. 3 .

 3 Fig. 3. Example of model simulations for a TK model with biotransformation (three metabolites) and influence of the variations of parameter m 1 (biotransformation rate for metabolite 1) on (a) the parent compound, (b) metabolite 1, (c) metabolite 2, and (d) metabolite 3. TK, toxicokinetic.

Fig. 4 . 1 (

 41 Fig. 4. Example of model simulations for a TK model with biotransformation (three metabolites) and influence of the variations of parameter e k 1 (elimination rate of metabolite 1) on (a) the parent compound, (b) metabolite 1, (c) metabolite 2, and (d) metabolite 3. TK, toxicokinetic.

Fig. 5 .

 5 Fig. 5. Example of model simulations for a TK model with biotransformation (three metabolites) and influence of the variations of parameter u w k on (a)

Table 1 . Symbols, meaning, and units of all parameters and variables involved in the full set of ordinary differential equations defining the generic one- compartment toxicokinetic model Symbol Meaning and unit

 1

	t								time (in time units)
	t c								duration of the accumulation phase (in time units)
	C p (t)		internal concentration of the parent compound at t (in μg•g -1 )
	( ) C t m 		internal concentration of metabolite ℓ at t (in μg•g -1 )
	U	=	∑	=1 I i k c u i i	sum of all uptake terms
	c i								concentration via the exposure route i (in μg•mL -1 )
	E	=	∑	J j	=1 k	e	j	sum of all elimination terms for the parent compound
	M	=	=1 L k ∑   m	sum of all elimination terms for metabolite ℓ
	i								index of exposure sources, i = 1…I
	j								index of elimination processes, j = 1…J
	ℓ								index of metabolites, ℓ = 1…L
	I								total number of exposure sources (#)
	J								total number of elimination processes (#)
	L								total number of metabolites (#)
	k	u i				uptake rate of exposure source i (per time units)
	k	e j				excretion rates of elimination process j (per time units)
	k	e 				excretion rates of metabolite ℓ (per time units)
	k	m 				metabolization rate of metabolite ℓ (per time units)

Table 2 . Sets of inputs for model simulations according to Figures 1 to 5

 2

	Input	Units	Simple (section 4.1)	Multiple exposures (section 4.2)	Biotransformation (section 4.3)
	c w	μg•mL -1	0.0044	-	15.53
	c s			μg•g -1	-	56.6	-
	c f						μg•g -1	-	1.46	-
	t c						days	49	7	1
	k	u w	days -1	10.46	-	16,740
	k	u s	days -1	-	0.071	-
	k	u f	days -1	-	0.013	-
	k		1 m	days -1	-	-	73.27
	k		2 m	days -1	-	-	0.5166
	k		3 m	days -1	-	-	0.1957
	k	e e	days -1	0.04	0.178	4.164
	k	1 e	days -1	-	-	561
	k		2 e	days -1	-	-	0.123
	k	3 e	days -1	-	-	0.7808
	(-) Hyphens stand for no required inputs.

Acknowledgments

The authors would like to express their sincere thanks to Miléna Kaag, Yacout Lahlou and Nino Molin who designed the MOSAI-C bioacc prediction tool as part of their 4 th year study project at the National Institute of Applied Sciences (INSA) in Lyon (France). This work was made under the umbrella of the French GDR "Aquatic Ecotoxicology" framework, aiming at fostering stimulating scientific discussions and collaborations for more integrative approaches.

Funding

This work was financially supported by the Graduate School H2O'Lyon (ANR-17-EURE-0018) and "Université de Lyon" (UdL), as part of the program "Investissements d'Avenir" run by "Agence Nationale de la Recherche" (ANR).

Data sharing statement

No additional data are available.

Conflict of interest

The authors have no conflicts of interest related to this publication.

Author contributions

Contributed to study concept and design (SC, AR and CL), generic solving (SC, CL), acquisition of the data (AR), assay model simulations (AR), drafting of the manuscript (SC, AR and CL), supervision (SC, CL).