J.-S Affagard

P Feissel

S F Bensamoun
email: sabine.bensamoun@utc.fr.

Measurement of the quadriceps muscles displacement and strain fields with ultrasound and Digital Image Correlation (DIC) techniques

Keywords: muscle, digital image correlation, ultrasound, displacement and strain fields, diffuse approximation

Objectives -In vivo characterization of muscle tissues is necessary to increase the understanding of musculoskeletal pathologies. Currently, few imaging techniques are able to characterize the in vivo mechanical properties of muscle. Thus, the purpose was to couple experimental and numerical tools to characterize the displacement and strain fields. This approach was validated on a phantom and then applied in vivo on thigh muscles.

Materials and methods -An experimental compression device was developed in order to measure the displacement field from a phantom with inclusions, and from quadriceps muscles, using compressed and uncompressed ultrasound acquisitions (9MHz). A diffuse approximation method was used to determine the optimal ZOI (Zone Of Interest) size (4, 8, 16 pixels) allowing for the realization of the DIC with a low noise and a high spatial resolution (about 100µm). Then, the strain field was represented with a diffuse approximation radius, optimized from the theoretical noise.

Results -The in vitro and in vivo displacement fields confirmed the qualitative ultrasound observations. The ZOI size was found optimal at 8 pixels. Then, the diffuse approximation radii were determined at 9 and 13 data points for the phantom and the quadriceps, respectively, providing the strain field cartographies. These cartographies validated the experimental protocol and the displacement field measurement.

Conclusion -The present study revealed the feasibility of the Digital Image Correlation technique to assess the in vivo muscle mechanical behavior, and could be extended to different soft tissues.

INTRODUCTION

The in vivo identification of the muscle mechanical behavior remained a challenge to better understand the pathologies. Medical imaging techniques such as ultrasound (US) or magnetic resonance imaging (MRI) allow for the characterization of the morphological properties in order to diagnose and follow-up treatments. To complete the structural analysis, dynamic elastography techniques using measurement of shear wave speed in function of frequencies enable the assessment of the muscle viscoelastic properties [START_REF] Debernard | In vivo characterization of the muscle viscoelasticity in passive and active conditions using multifrequency MR Elastography[END_REF][START_REF] Gennisson | Viscoelastic and anisotropic mechanical properties of in vivo muscle tissue assessed by supersonic shear imaging[END_REF][START_REF] Leclerc | Analysis of liver viscosity behavior as a function of multifrequency magnetic resonance elastography (MMRE) postprocessing[END_REF]. Moreover, it is well known that muscle tissue is a complex material with nonlinear and quasi incompressible behaviors. Thus, to further analyze these characteristics, experimental tests and numerical tools were used to measure the muscle displacement field.

Quasi-static tests such as traction [START_REF] Mofid | Mechanical skin thinning-tothickening transition observed in vivo through 2D high frequency elastography[END_REF], suction [START_REF] Delalleau | A nonlinear elastic behavior to identify the mechanical parameters of human skin in vivo[END_REF][START_REF] Kauer | Inverse finite element characterization of soft tissues[END_REF] or indentation [START_REF] Alja | Developing a bio-mechanotronic probing system for estimating soft tissue Young's modulus in vivo[END_REF][START_REF] Delalleau | Characterization of the mechanical properties of skin by inverse analysis combined with the indentation test[END_REF][START_REF] Moerman | Digital image correlation and finite element modelling as a method to determine mechanical properties of human soft tissue in vivo[END_REF] have been performed in order to measure in vivo the external displacements of the skin tissue. Subsequently, Finite Element Model Updating (FEMU) was developed to determine the skin elastic behavior. In addition, ultrasound technique was used to better characterize the geometry of the skin and the environmental tissues with suction [START_REF] Diridollou | In vivo model of the mechanical properties of the human skin under suction[END_REF][START_REF] Hendriks | A numerical-experimental method to characterize the non-linear mechanical behaviour of human skin[END_REF][START_REF] Hendriks | The relative contributions of different skin layers to the mechanical behavior of human skin in vivo using suction experiments[END_REF][START_REF] Iivarinen | Experimental and computational analysis of soft tissue mechanical response under negative pressure in forearm[END_REF] or indentation tests [START_REF] Iivarinen | Experimental and computational analysis of soft tissue stiffness in forearm using a manual indentation device[END_REF][START_REF] Koo | A mechano-acoustic indentor system for in vivo measurement of nonlinear elastic properties of soft Tissue[END_REF][START_REF] Tönük | Nonlinear elastic material property estimation of lower extremity residual limb tissues, Neural Systems and Rehabilitation Engineering[END_REF][START_REF] Tönük | Nonlinear viscoelastic material property estimation of lower extremity residual limb tissues[END_REF]. This modeling provided a more accurate identification of the soft tissues. Other studies coupled in vivo mechanical tests with MRI acquisitions to identify the hyperelastic behavior of the arm and calf muscles [START_REF] Avril | Mixed experimental and numerical approach for characterizing the biomechanical response of the human leg under elastic compression[END_REF][START_REF] Tran | In vivo characterization of the mechanical properties of human skin derived from MRI and indentation techniques[END_REF]. This identification was performed from the displacement of the contours of the muscle, bone and adipose tissues.

To further analyze their characteristics, it was proposed that experimental tests and numerical tools could be used to measure the muscle displacement field and strain fields. Thus, Zerhouni et al. [START_REF] Zerhouni | Human heart: tagging with MR imaging -A method for noninvasive assessment of myocardial motion[END_REF] used the MR tagging to measure the 2D myocardial motion during the cardiac cycle. Recently, this method was coupled with quasi-static excitations performed on the lower leg [START_REF] Fu | Motion tracking and strain map computation for quasi-static magnetic resonance elastography[END_REF] in order to obtain the displacement and strain fields. In addition, a 3D analysis was also carried out on the arm [START_REF] Moerman | Validation of SPAMM tagged MRI based measurement of 3D soft tissue deformation[END_REF][START_REF] Moerman | A novel MRI compatible soft tissue indentor and fibre Bragg grating force sensor[END_REF]. Moreover, an ultrasound technique was associated with static compression tests on phantom to image the cartography of the strain field in order to calculate the Young's modulus as the stress / strain ratio [START_REF] Céspedes | Elastography: elasticity imaging using ultrasound with application to muscle and breast in vivo[END_REF][START_REF] Ophir | Elastograph y: a quantitative method for imaging the elasticity of biological tissues[END_REF][START_REF] Ponnekanti | Axial stress distributions between coaxial compressors in elastography: an analytical model[END_REF][START_REF] Ponnekanti | Ultrasonic imaging of the stress distribution in elastic media due to an external compressor[END_REF]. Ultrasound was further improved with an algorithm based on the correlation of B-mode ultrasound signals to measure the displacement and strain fields [START_REF] Jiang | A parallelizable real-time motion tracking algorithm with applications to ultrasonic strain imaging[END_REF][START_REF] Zhu | A modified block matching method for real-time freehand strain imaging[END_REF]. This measurement was performed on phantom [START_REF] Hall | Recent results in nonlinear strain and modulus imaging[END_REF][START_REF] Pavan | Nonlinear elastic behavior of phantom materials for elastography[END_REF][START_REF] Pavan | A nonlinear elasticity phantom containing spherical inclusions[END_REF] and on breast tissue to image the hyperelastic properties and to discern tumors [START_REF] Hall | Recent results in nonlinear strain and modulus imaging[END_REF][START_REF] Barbone | Elastic modulus imaging: some exact solutions of the compressible elastography inverse problem[END_REF][START_REF] Oberai | Linear and nonlinear elasticity imaging of soft tissue in vivo: demonstration of feasibility[END_REF]. Another way to measure the displacement field was to use the digital image correlation (DIC) [START_REF] Hild | Digital image correlation: from displacement measurement to identification of elastic properties -a review[END_REF][START_REF] Hild | Correli Q4: A software for-finite-element-displacement field measurements by digital image correlation[END_REF][START_REF] Slane | The use of 2D ultrasound elastography for measuring tendon motion and strain[END_REF] that enables to measure the displacement field from images. Thus, Han et al. [START_REF] Han | Application of digital image cross-correlation and smoothing function to the diagnosis of breast cancer[END_REF] acquired ultrasound images on phantom, under several static loading, coupled with DIC method to measure the 2D and 3D displacement and strain fields. Moreover, ex-vivo tendon analysis was also performed with DIC to validate 2D ultrasound elastography results, and to demonstrate the potential of the DIC method for the characterization of in vivo biological tissue [START_REF] Slane | The use of 2D ultrasound elastography for measuring tendon motion and strain[END_REF].

Thus, the originality of this study was to use the DIC on ultrasound acquisitions of in vivo thigh muscles to obtain a local displacement measurement. Moreover, a home-made compression device was developed to perform a quasi-static loading. Thus, the purpose of this present study was to couple experimental and numerical tools to characterize the displacement and strain fields. This approach was firstly validated on a phantom and then applied in vivo on thigh muscles.

MATERIALS AND METHODS

Experimental set up

A home-made compression device (Fig. 1a-1d) was composed of three plastic plates (200 x 300 mm2).

The lower one was only used as a support for the leg and four threaded rods were fixed at each corner allowing the up and down motions of the upper plates. The two upper plates were used to compress the muscle where the level was adjusted with screws. Then, sensors (I-scan, Tekscan, 5027) (Fig. 1d) were placed between the two upper plates in order to quantify the distribution of pressure. A specific design was made on the upper plates to fit and to maintain the ultrasound probe (9MHz) during the loading and the ultrasound acquisitions.

In vitro and in vivo ultrasound acquisitions

The experimental device was applied in vitro to a phantom (CIRS, Norfolk, VA), mimicking soft tissues, composed of eight spherical inclusions (from #1 to #8) (Diameter: 10 mm and 20 mm) having different stiffnesses (Fig. 1b). The four small and large spheres were aligned and spaced (30 mm) at a depth of 15 mm and 35 mm, respectively. According to the CIRS company, the phantom background had an elasticity of 25 ± 6 kPa and the spherical inclusions had an elasticity of 8 ± 3 kPa (#1, #5), 14 ± 4 kPa (#2, #6), 45 ± 8 kPa (#3, #7), 80 ± 12 kPa (#4, #8) (Fig. 1b).

Subsequently, ultrasound images were acquired for the two central small spherical inclusions (#2, #3) (Fig. 1c).

The homemade compression device was used for in vivo experiment performed on the thigh muscles of a 33 year old male without venous pathologies (Fig. 1d). The ultrasound probe was placed in the lower third part of the thigh due to the low amount of adipose tissue allowing a better strain of the muscle tissue. Fig. 1e showed the localization of the investigated quadriceps muscles composed of the vastus medialis (VM) and vastus intermedius (VI) (Fig. 1f) acquired with the same optimized ultrasound parameters.

Measurement of the displacement field with the DIC

The DIC method (Correli_Q4 [START_REF] Hild | Correli Q4: A software for-finite-element-displacement field measurements by digital image correlation[END_REF]) is based on the grey level comparison between two B-mode ultrasound (US) acquisitions performed on uncompressed (Fig. 2a,3a) and compressed (Fig. 2b,3b) phantom and quadriceps muscle in order to measure the displacement field. The ultrasound parameters were optimized to perform efficiently the digital image correlation (DIC), as described in a previous study [START_REF] Affagard | Identification des propriétés hyperélastiques des muscles de la cuisse à l'état passif: couplage des techniques de corrélation d'images aux techniques d'imagerie médicale[END_REF][START_REF] Affagard | Development of an inverse approach for the characterization of in vivo mechanical properties of the lower limb muscles[END_REF] where the ultrasound parameters were set up to reduce the random error of the measurement.

A zone of interest (ZOI) was located inside a ROI (region of interest) defined on the ultrasound image (Fig. 1c). The DIC aimed at matching the ZOI (8 pixels) between reference and deformed images. It must be noted that the DIC performances were strongly linked to the grey level and image characteristics (i.e contrast and dynamic of the speckle). Considering f(x) and g(x), the grey level values of the reference and deformed images at position x, g(x) can be defined as:

        x b x u x f x g    (1)
u(x) is the displacement field and b(x) the noise induced by the image acquisition. The specificity of the present global DIC method is to seek the displacement field corresponding to a set of shape function. The field, u(x), can therefore be written as:

      m a m am e N x x u ,   (2)
m is the number of nodes for one element (m = 4), ναm are the degrees of freedom unknowns of the problem, Nm(x) are the shape functions and eα is a characteristic vector of the considered element.

To obtain the cartographies representing the lateral (Ux) (Fig. 2c,3c) and axial (Uy) (Fig. 2d,3d) displacement fields five ultrasound acquisitions were performed with intermediate loadings (muscle step: 15N, phantom step: 112N). Thus, the DIC convergence was facilitated between the initial (0N) and final (muscle: 60N, phantom: 450N) loadings. These loading level differences were chosen because the same external displacement levels were observed at the sample surface between the two experiments. Finally, the displacement fields were displaced to match the bone structure between the initial and final images. Thus the femur is taken as the reference.

Measurement of the strain field

The strain field was calculated using diffuse approximation from the displacement field measurement. This method based on using local least squares [START_REF] Cleveland | Smoothing by local regression: Principles and methods[END_REF] consists on performing a local diffuse approximation regression [START_REF] Avril | Estimation of the strain field from full-field displacement noisy data[END_REF][START_REF] Nayroles | La méthode des éléments diffus[END_REF]. Thus, cartographies of strain fields were represented for three different ZOI (4, 8 and 16 pixels) with three different diffuse approximation radii (4, 8 and 16 data points) (Fig. 4). These cartographies were obtained for shear strain.

The radius of the diffuse approximation is optimal when the maximum of noise is filtered without losing mechanical information (Fig. 5a-5c, 6a-6c). Thus, to determine the more appropriate diffuse approximation radius the standard deviation for each experimental strain (lateral, axial and shear strains) was represented as a function of the radius (from 4 to 20 data points).

Then, the evolution of the theoretical noise filtering, assuming a white Gaussian noise with a standard deviation (σ), can be written in function of a linear operator of the reconstruction (M) [START_REF] Avril | Estimation of the strain field from full-field displacement noisy data[END_REF]:

    T b b MM u 2 cov    (3)
This theoretical curve was shifted in order to match the same value of standard deviation with the experimental curve. Thus, a reference value (dot line) was established to represent the level of the applied filtering. Then, the standard deviation (σ) of the measured displacement field was quantified between 11 uncompressed images by considering a white Gaussian noise for each ZOI (4, 8, 16 pixels).

Finally, the intersection with the theoretical noise filtering revealed the optimal radius for which the cartographies of the lateral, axial and shear strain were represented (Fig. 5d-5f, 6d-6f).

RESULTS

Result of the in vitro and in vivo displacement fields

The ultrasound images acquired on the phantom (Fig. 2a-2b) revealed a higher deformation for the left inclusion (#3) compared to the right (#2) one. This structural result was expected due to the difference of stiffness between the two spherical inclusions.

Fig. 2c and Fig. 2d showed the cartographies of the displacement fields, represented on the uncompressed image, measured along the Y and X axis, respectively.

Along the direction of compression (Fig 2c), highest displacements (about 7mm) were observed (red color) close to surface of the phantom, and lowest displacements (blue color, about 3mm) were found in depth. In addition, a higher gradient of displacement was observed for the left inclusion (#3). This last result was in agreement with the structural deformation observed on the ultrasound acquisitions. Along the X axis (Fig. 2d), the scale indicated few displacements (from -1mm to 0.2mm) which is expected according to the applied compression test. However, a higher gradient displacement field (from -1mm to -0.4mm) was observed in the left inclusion (#3) compared to the right one (#2) (from -0.3mm to -0.1mm).

This last result confirmed the difference of mechanical properties between the two inclusions.

The result of the ultrasound images acquired for the quadriceps muscle (Fig. 3a-3b) revealed a visual sliding of both vasti muscles (VM and VI) around the aponeurosis tissue. The result of the displacement field (Fig. 3c) measured along the Y axis, showed highest displacements (about 10mm, red color) close to the fat tissue while any displacements was found around the bone tissue. Furthermore, a gradient of displacement was observed inside the muscles. Fig. 3d illustrated the cartography of the displacement field measured along the X axis. A high displacement gradient (from -5mm to 2mm), indicated by arrows, was observed around the aponeurosis located between the vastus intermedius (VI) and the vastus medialis (VM) muscles. This last result confirmed the qualitative sliding obtained after the compressive test.

Repeatability (lateral: 1.23 mm, axial: 0.96 mm) and reproducibility (lateral: 0.25 mm, axial: 0.23)

analysis showed similar quantitative and qualitative displacement fields [START_REF] Affagard | Development of an inverse approach for the characterization of in vivo mechanical properties of the lower limb muscles[END_REF].

Quantification of the experimental noise

Table 1 showed the results of the standard deviation (σ) as a function of the ZOI size (4, 8 and 16 pixels)

for the phantom and the quadriceps muscle. As expected, the value of the standard deviation decreases when the ZOI size increase. Along the ultrasound beam (axial direction), the result of the standard deviation was less important than the lateral noise. The comparison between the in vivo and in vitro experimental noises revealed a higher noise for the muscle tissue compared to the phantom.

Result of the in vitro and in vivo strain fields

Determination of the ZOI from the shear strain field

Fig. 4 illustrated the results of the shear strain fields obtained for three different ZOI and three radius of diffuse approximation.

For a ZOI of 4 pixels, the strain field was very noisy and no information can be interpreted (Fig. 4a) with a radius of 4 data points. On the contrary, the spherical inclusions were clearly represented with a radius of 16 data points (Fig. 5c). Nevertheless, the contours of the spherical insertions were not well defined due to the noise.

For a ZOI of 8 pixels, the best cartographies were obtained with a radius between 8 and 16 data points (Fig. 4e-4f). Indeed, these configurations enabled the localization of the two spherical inclusions represented with clear contours and a quasi-homogeneous strain field [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF][START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF].

Using a ZOI of 16 pixels, provided a low spatial resolution, the best compromise was represented with a radius of 4 data points (Fig. 4g). For larger radius, the strain field was strongly filtered leading to a lack of information (Fig. 4h-4i).

In agreement with all of these results, a ZOI of 8 pixels was chosen in order to determine the diffuse approximation radius allowing the interpretation of the strain fields.

Analysis of the in vitro strain field with the optimal diffuse approximation radius

The evolution of the experimental standard deviation () was represented for the lateral (Fig. 5b), axial (Fig. 5a) as a function of the radius (blue curve). As expected, the standard deviation decreased exponentially with the radius. The theoretical curve (red) showed the evolution of the standard deviation of the noise.

Fig. 5a revealed that the dotted line intersected the experimental curve (blue) for a radius value of 9 data points. As a consequence, the noise is considered filtered above a radius of 9 data points while below this value the information was degraded. Therefore, the optimal radius was chosen to 9 data points for a ZOI of 8 pixels.

The lateral standard deviation (Fig. 5b) showed that whatever the range of radius is, the experimental curve was always above the dotted line. This last result indicated that the information and the noise were in the same range. Therefore, the lateral strain was not solicited by the present mechanical test.

Thus, the diffuse approximation radius was set up to 9 data points to analyze the strain fields.

The cartographies of the strain fields were represented in the axial (Fig. 5c), lateral (Fig. 5d) and shear (Fig. 5e) directions.

Fig. 5c revealed a high strain field below the left inclusion (#3) while a weak strain field was found below the right inclusion (#2). This result was related to the environmental structure around the inclusions.

Indeed, the strain below the right inclusion (#2) was lower than the strain inside the inclusion, demonstrating that the inclusion was stiffer than the background. On the contrary, the strain below the left inclusion (#3) was higher than the strain inside the inclusion, revealing a softer inclusion than the background.

Fig. 5d was not interpretable due to the low level of information which did not allow the localization of the spherical inclusions.

Fig. 5e showed the shear strain field where the contours of the two inclusions were well identified. In addition, a homogeneous strain was observed inside both inclusions. Therefore, these results confirmed the ultrasound observations and demonstrated the feasibility of the present method to measure the displacement and strain field within the muscle tissue.

Analysis of the in vitro strain field with the optimal diffuse approximation radius

A ZOI of 8 pixels was chosen to analyze the muscle strain fields as previously presented for the phantom.

A similar process as for the in vitro study was performed for the choice of the diffuse approximation radius. A radius of 13 data points was determined as optimal for the in vivo analysis (Fig 6a) and Fig. 6b showed that the information and the noise were in the same range leading more difficult the analysis of the lateral strain field.

Fig. 6c-6e presented the axial, lateral and shear strains of the quadriceps muscle. The three cartographies revealed a sliding of both vastis around the aponeurosis area as it was observed on the ultrasound acquisitions and displacement field (Fig. 3). Moreover, Fig. 7c presented a high strain level located inside the fat tissue, demonstrating a lower stiffness for fat tissue than muscle tissue.

DISCUSSION

The originality of this study was to develop an experimental protocol to measure the muscle displacement and the strain field using a DIC. The present protocol was set up on a phantom with structural and mechanical properties know from the CIRS' company. The ultrasound structural deformation observed on the inclusions validated the developed mechanical test and also the displacement field measurement.

Furthermore, this in vitro work enabled to fix the ZOI size to 8 pixels due to the clear contours of the inclusions from the strain field. This result allowed to choose a reference ZOI size (8 pixels) for the in vivo study. Moreover, a complementary study was performed on several ZOI sizes and confirmed the use of 8 pixels for the quadriceps muscle.

The comparison of the displacement fields between the phantom (Fig. 2c) and the muscle (Fig. 3c)

showed similar maximal axial displacement levels (phantom: -7mmn, muscle: -10mm) despite different loadings were applied. Indeed, a higher loading was performed on the surface of the phantom due to the boundary conditions. In fact, the lateral direction (x axis), the boundary displacement being fixed for the phantom increasing the resulting force. This result demonstrated that the CIRS phantom could be used to mimic the muscle displacement. Obviously, the loadings were adjusted according to the respective mechanical properties. Furthermore, the cartography of axial muscle displacement field (Fig. 3c) showed a continuity of the displacement as a function of the muscle depth, and validated the developed experimental protocol (quasi static test and DIC).

On the contrary, the comparison of the lateral displacement (Fig. 2d and Fig. 3d) showed a higher range of displacement values for the muscle (from -5mm to 5mm) compared to the phantom (from -1mm to 1mm). This may be due to the difference of texture. Indeed, a speckle was present on the phantom while the presence of lateral white bands, reflecting the muscle structure, may induce measurement errors on the lateral axis. Moreover, during the mechanical test, a rotation of the muscle structures around the bone and a sliding of the tissue close to the aponeurosis may explain this difference of range.

The quantification of the experimental noise showed a higher incertitude for the lateral displacements compared to the axial ones for both samples. This result may be due to the direction of the beam which is in the compression direction inducing a higher noise in the lateral direction. The uncertainty of the lateral displacement field lead [START_REF] Oberai | Linear and nonlinear elasticity imaging of soft tissue in vivo: demonstration of feasibility[END_REF] to propose the identification of the linear and non-linear elasticity of soft tissues from a cost function only based on the axial component.

The comparison of the experimental noise between the phantom and the muscle for both directions revealed a higher error for the muscle tissue due to the difference of texture and due to muscle movements, such as involuntary contraction or blood flow, providing changes during the ultrasound acquisitions.

The method used to determine the diffuse approximation radius was similar to Feissel et al. [START_REF] Feissel | Use of diffuse approximation on DIC for early damage detection in 3D carbon/epoxy composites[END_REF] in order to obtain an optimal filtering. In the present study, the assumption of the white Gaussian noise may induce an error on the diffuse approximation radius choice. Indeed, the DIC method (Correli_Q4), revealing the displacement field as a shape function, provided another reparation of noise [START_REF] Feissel | Use of diffuse approximation on DIC for early damage detection in 3D carbon/epoxy composites[END_REF].

To our knowledge, the DIC is a technique often used in the mechanical field, this study is the first to combine the DIC method with ultrasound acquisition on in vivo muscle tissue in order to obtain the strain field. The limitation of the present study is the measurements obtained from one subject and should be realized on other participants. Indeed, the quality of images and the mechanical properties of the muscles will differ between individuals, due to the gender, age and health of muscle. Thus, the measured displacement field might be a starting point for the characterization of the in vivo elastic properties and may be extended to the measurement of different soft tissues. In addition, according the important level of sliding and strain between both vasti muscles, a way to improve the strain cartographies would be to develop an integrated-DIC method to remove the discontinuity problems.

Finally, from this measurement, the next step will be to develop an inverse method, coupling this acquired displacement field with finite element simulation, to identify the muscle mechanical behavior.

FIGURESFig. 1 .

 1 FIGURES

Fig. 2 .Fig. 3 .

 23 Fig. 2. Ultrasound images of (a) uncompressed and (b) compressed phantom. Cartographies of the displacement fields along the axial (c) and lateral (d) axis.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Evolution of the experimental (blue) and theoretical (red) standard deviation for a ZOI of 8 pixels during an axial (a) and lateral (b) strain. Cartographies of the axial (c), lateral (d) and shear (e) strain obtained for a ZOI of 8 pixels and a radius of 9 pixels.

Table 1 :

 1 Standard deviation as a function of the ZOI size for the phantom and the quadriceps muscle

	Phantom	Muscle

ACKNOWLEDGEMENT

This project is co-financed by the European Union engaged in Picardy with the European Regional Development Fund and CNRS (grant Collegium UTC CNRS INSIS).

CONFLICT OF INTEREST STATEMENT

All authors do not have conflict of interest.


Redha A Alliche
email: alliche@i3s.unice.fr

Tiago Da

Silva Barros

Ramon Aparicio-Pardo

Lucile Sassatelli

PRISMA: A Packet Routing Simulator for Multi-Agent Reinforcement Learning

Keywords: ns-3, Multi-Agent, Packet Routing, Reinforcement Learning, Network Simulation, ML tool

In this paper, we present PRISMA: Packet Routing Simulator for Multi-Agent Reinforcement Learning. To the best of our knowledge, this is the first tool specifically conceived to develop and test Reinforcement Learning (RL) algorithms for the Distributed Packet Routing (DPR) problem. In this problem, where a communication node selects the outgoing port to forward a packet using local information, distance-vector routing protocol (e.g., RIP) are traditionally applied. However, when network status changes very dynamically, is uncertain, or is partially hidden (e.g., wireless ad hoc networks or wired multi-domain networks), RL is an alternate solution to discover routing policies better fitted to these cases. Unfortunately, no RL tools have been developed to tackle the DPR problem, forcing the researchers to implement their own simplified RL simulation environments, complicating reproducibility and reducing realism. To overcome these issues, we present PRISMA, which offers to the community a standardized framework where: (i) communication process is realistically modelled (thanks to ns3); (ii) distributed nature is explicitly considered (nodes are implemented as separated threads); (iii) and, RL proposals can be easily developed (thanks to a modular code design and real-time training visualization interfaces) and fairly compared them.

I. INTRODUCTION

In the last years, Reinforcement Learning (RL) [START_REF] Sutton | Reinforcement learning: An introduction[END_REF] using Deep Neural Networks (DNNs) [START_REF] Bengio | Learning deep architectures for AI[END_REF], also called Deep Reinforcement Learning (DRL), has obtained ground-breaking results in solving highly complex tasks, such as humanlike performance results at Atari video games [START_REF] Mnih | Human-level control through deep reinforcement learning[END_REF] or beating AlphaGo [START_REF] Silver | Mastering the game of Go without human knowledge[END_REF] world champion. In communication networks, DRL has also been widely used in many networking technologies and problems. One of them is the Distributed Packet Routing [START_REF] Mukhutdinov | Multi-agent deep learning for simultaneous optimization for time and energy in distributed routing system[END_REF]- [START_REF] Manfredi | Relational deep reinforcement learning for routing in wireless networks[END_REF]. In this problem, no complete and centralized view of network topology and traffic demands is available, which poses a challenge. This is the case in multi-hop wireless networks [START_REF] Tassiulas | Stability properties of constrained queueing systems and scheduling policies for maximum throughput in multihop radio networks[END_REF] or in multi-domain optical networks [START_REF] Chen | Demonstration of distributed collaborative learning with end-to-end qot estimation in multi-domain elastic optical networks[END_REF]. More precisely, in the DPR problem, each packet can be potentially routed differently regardless of the flow they belong to. Besides, the per-packet decisions are made locally by distributed agents placed at the routing nodes. These agents exploit local

The author acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-19-CE-25-0001-01 (ARTIC project). This work was performed using HPC resources from GENCI-IDRIS (Grant 2021-AD011012577).

ISBN 978-3-903176-48-5© 2022 IFIP information, such as packet headers and neighboring nodes and link states. However, in all the above-mentioned works, the proposed DRL approaches were evaluated on ad hoc discrete-time packet-level simulation environments (typically implemented in python) tailored to the assumptions and simplifications made by each study. Namely, these works do not generally implement the RL agents as separated threads or processes (although they claim to be multi-agent studies); and, they usually assume that at most one packet per router can be transmitted at each time step, which introduces an artificial synchronization into the routers' dynamics. This has two main consequences: (i) a reduced realism of the simulated communication process; and, (ii) an obstacle to fairly comparing these DRL proposals to state-of-the-art approaches or even against each other. This lack of standardized Machine Learning (ML) tools in the networking community [START_REF] Mestres | Knowledge-defined networking[END_REF], [START_REF] Gawłowicz | ns-3 meets OpenAI Gym: The Playground for Machine Learning in Networking Research[END_REF] represents a major issue to guarantee reproducibility. We point out that ML reproducibility issues are becoming a serious concern [START_REF] Pineau | Improving reproducibility in machine learning research(a report from the NeurIPS 2019 reproducibility program)[END_REF].

As a consequence, in the last years, some tools [START_REF] Gawłowicz | ns-3 meets OpenAI Gym: The Playground for Machine Learning in Networking Research[END_REF], [START_REF] Yin | Ns3-ai: Fostering artificial intelligence algorithms for networking research[END_REF] devoted to the application of RL to networks have been developed. These tools allow python RL agents to interact with the popular ns-3 network simulator [START_REF] Nsnam | Ns-3 documentation website[END_REF]. Nevertheless, these proposals are not natively compatible with the collaborative and distributed multi-agent setting of the DPR problem, where agents collaborate to take decisions in a distributed manner.

In this paper, we propose an open-source RL simulation environment designed to overcome the aforementioned drawbacks in the application of DRL to the DPR problem. We bring the following main contributions:

1) a RL framework designed specifically for considering the distinctive characteristics of the DPR problem, serving as a playground where the community can easily validate their RL approaches and compare them. 2) more realistic modelling of the communication process based on: (i) the ns-3 [START_REF] Nsnam | Ns-3 documentation website[END_REF] network simulator; and, (ii) a multi-threaded implementation for each agent. 3) a modular code design, which allows researchers to test their own RL algorithm, without needing to work on the implementation of the environment. 4) visual tools based on tensorboard [START_REF]Tensorboard: Tensorflow's visualization toolkit[END_REF] allowing to track training and test phases. The rest of this paper is organized as follows. Section II discusses the related works. In Section III, we present the pa-per's background. Section IV presents PRISMA. The simulator usage is detailed in Section V and illustrated in Section VI. Finally, we conclude the paper in Section VII.

II. RELATED WORKS

Related works can be grouped into the next two categories. DRL for packet routing. In recent years, the first works applying DRL to the DPR problem have been published. DPR has been studied in [START_REF] Mukhutdinov | Multi-agent deep learning for simultaneous optimization for time and energy in distributed routing system[END_REF]- [START_REF] Manfredi | Relational deep reinforcement learning for routing in wireless networks[END_REF], [START_REF] Boyan | Packet routing in dynamically changing networks: A reinforcement learning approach[END_REF]. The seminal paper [START_REF] Boyan | Packet routing in dynamically changing networks: A reinforcement learning approach[END_REF] was the first to apply a RL method, the Q-learning [START_REF] Watkins | Learning from delayed rewards[END_REF], to DPR, giving rise to the Q-routing paradigm. More details about this paradigm will be given in the next section. The study in [START_REF] Mukhutdinov | Multi-agent deep learning for simultaneous optimization for time and energy in distributed routing system[END_REF] applies the Deep Q-Network (DQN) framework (i.e., a neural network fits the Q-value function) to the Q-routing algorithm, yielding to the DQN routing. In all these works, the node state only codes the destination of the packet to forward. More recently, works [START_REF] You | Toward packet routing with fully distributed multiagent deep reinforcement learning[END_REF]- [START_REF] Manfredi | Relational deep reinforcement learning for routing in wireless networks[END_REF] extended the DQN routing by adding to the destination additional inputs, particularly node buffer occupancy. DQRC (DQR with Communication) [START_REF] You | Toward packet routing with fully distributed multiagent deep reinforcement learning[END_REF] also considers action history and future destinations. Authors in [START_REF] Manfredi | Relational deep reinforcement learning for routing in wireless networks[END_REF] focus on the relational features (independent of the network graph), such as the distance to the destination. Finally, in [START_REF] Chen | Multiagent meta-reinforcement learning for adaptive multipath routing optimization[END_REF], the meta-learning framework is applied to the DQN routing to learn different tasks, where a task corresponds to a given traffic matrix. RL environments for networking. Only two tools have been developed to allow a Python agent to interact with a network simulation: first, ns3-gym [START_REF] Gawłowicz | ns-3 meets OpenAI Gym: The Playground for Machine Learning in Networking Research[END_REF]; and, more recently, ns3-ai [START_REF] Yin | Ns3-ai: Fostering artificial intelligence algorithms for networking research[END_REF]. Both of them are based on the wellknown ns-3 network simulator, but they follow different approaches to implement a RL environment. Fist, ns3-ai is not only restricted to Gym as Python environment (as ns3-gym does). Second, it makes use of a shared memory pool as mechanism to connect the ns-3 simulator with the Python framework (agent and environment), differently from the ns3-gym sockets-based approach. The ns-3, Gym and ns3-gym tools are detailed in the next section.

III. BACKGROUND

In this section, we present the elements which PRISMA is built over: the RL framework used to solve the packet routing problem, the RL toolkits, and, the ns-3 network simulator.

A. Q-routing: Reinforcement Learning for Packet Routing

The DPR problem can be formalized as Multi-Agent Partially Observable Markov Decision Process (POMDP) [START_REF] Littman | Markov games as a framework for multi-agent reinforcement learning[END_REF]: a Markov Decision Process where each agent can only locally observe its environment. The first RL method proposed to tackle the DPR was the Q-routing [START_REF] Boyan | Packet routing in dynamically changing networks: A reinforcement learning approach[END_REF], based on the classical Q-learning [START_REF] Watkins | Learning from delayed rewards[END_REF], which we will use to describe the POMDP and the RL approach used to solve it. Note that other RL methods could be also adapted to the DPR problem.

Let N be the set of routing nodes (agents), where each agent n has its own local observation space O n and its own action space A n . When a new packet arrives at time t at the node n, the node n selects as action a n the next hop node n to forward the packet. This decision is taken depending on the local observation of the router o n (e.g., the current packet destination, the node buffers' occupancy, . . . ). As a consequence, the agent n receives a reward r n : the next-hop packet delay (i.e. the packet delay to travel from n to n ), and the whole network makes a transition to a new state. When this packet (or another one) arrives to a node, this procedure is repeated. The action a n is selected to minimize an estimate of the expected end-to-end packet delay from node n to its final destination. Under the DRL scheme, this estimate, denoted as

Q n (o n , a n ; θ n ) (the Q-value),
is the output of a DNN with θ n weights. This DNN is trained to fit the target value Y Q n below:

Y Q n = r n + γ • τ • (1 -f ) (1)
where f indicates if the next hop is the packet destination, γ ∈ [0, 1] is a discount factor, r is the next-hop packet delay and τ is the remaining end-to-end delay from the next hop n to the destination computed as follows:

τ = min a n ∈A n Q n (o n , a n ; θ n ) (2)
where Q n (•; θ n ) is the output of next hop agent. The θ n weights are updated by stochastic gradient descent when minimizing the so-called Temporal Difference (TD) error:

Q(s, a; θ n ) -Y Q n 2 .

B. The ns-3 Network Simulator

The ns-3 simulator [START_REF] Nsnam | Ns-3 documentation website[END_REF] is a stable discrete-event network simulator, largely adopted by the research and educational community and accepted as a standard. It is free and opensource under GPLv2 and developed in C++ using objectoriented programming. The success of this tool is based on two main features: (i) a realistic modelling of the network systems: the protocol stack is properly abstracted by the simulator classes; (ii) a wide range of networking protocols (e.g., IP, TCP, UDP), communication technologies (e.g., Ethernet, Wi-Fi, LTE) and statistical models for channels, mobility, and traffic generation are supported by built-in classes. As a consequence, ns-3 becomes a natural choice to use as a network simulator for PRISMA.

C. Reinforcement Learning Tools

RL agents are usually implemented using Python packages (e.g., TensorFlow [START_REF]Tensorflow: An end-to-end open source machine learning platform[END_REF]) and interfaced with an environment where agents take decisions. The OpenAI™Gym [START_REF]Gym: Gym toolkit for creating reinforcement learning environments[END_REF] is one of the most popular toolkits to define RL environments in Python. Gym is used to test and compare different RL algorithms for a variety of problems, like Atari games or robotics [START_REF] Zamora | Extending the openai gym for robotics: a toolkit for reinforcement learning using ros and gazebo[END_REF]. Gym also provides a simple interface to pass the interactions between the agent and the environment, without making assumptions about the agent structure. Moreover, the Gym toolkit allows the creation of customized environments, which is very convenient since no network environment is natively available in Gym.

The first solution to this lack of networking RL environments was ns3-gym [START_REF] Gawłowicz | ns-3 meets OpenAI Gym: The Playground for Machine Learning in Networking Research[END_REF], which interfaces OpenAI ™Gym with the ns-3 network simulator. In other words, ns3-gym "transforms" a ns-3 network simulation into an RL environment by serving as gateway connecting ns-3 with Gym. Then, it provides an easy-to-use platform for developing and testing RL algorithms for networking problems. The Python process (the agent and the Gym environment) communicates with the C++ process (the ns-3 network simulation) via ZMQ sockets [START_REF]Zero mq: An open-source universal messaging library[END_REF]. ns3-gym is also open source under a GPL licence and can be easily extended. In PRISMA, we opt for this tool (instead ns3-ai) since ns3-gym is better documented.

IV. PRISMA DESCRIPTION

In this section, we present the technical details of PRISMA. First, a description of the design of the framework is given. Then, each feature is explained, ranging from the network definition to simulation visualization using tensorboard [START_REF]Tensorboard: Tensorflow's visualization toolkit[END_REF]. The PRISMA source code is publicly available at [START_REF]Prisma tool: An open marl framework for packet routing[END_REF].

A. Framework design principles

We consider the following principles: 1) Training-Forwarding separation: the agent training does not disturb the packet forwarding. 2) Modularity: code easy to reuse and modify.

3) Realistic simulation: implement the agent as close as possible to the real world. 4) Fast prototyping: easy and rapid modifications of the decision model. 5) Online simulation tracking: visualize in real-time the simulation evolution. For point 1), we choose a multi-threading approach: each node will be run on two separate threads. One is dedicated to the training process, and; the other one to the decision process. Hence, the node agents can be trained at a timestamp without disturbing the action computation at each packet arrival.

For point 2), we built node agents in an object-oriented manner, with separated functions for each task, so that a user can easily modify the behaviour of the agent without affecting the rest of the code. We document each part of the code and provide scripts for running the simulation.

Point 3) is a very important aspect of the PRISMA design, since we need reliable results as close as possible to the deployment of the agents in the real world. For that, we keep the link propagation delay and use ns-3 to simulate the network. We have also integrated the action handling to the environment and added a Poisson random traffic generator.

For point 4), PRISMA is separated into independent modules, so that a user may easily modify the code. For example, changing the neural network or the inputs of the model will not affect the other parts of the framework.

Finally, for point 5), we give to the user the ability to monitor the simulation progress in real-time. For example, the user can track some network or ML metric (e.g., average delay, buffers occupancy, error progression) to identify training issues in real-time. To do so, tensorboard has been integrated into the framework offering network monitoring. In addition, it adds the possibility to customize plots.

B. Feature description

Figure 1 depicts the PRISMA code structure, showing in green our contribution over the existing ns-3 and ns3-gym modules. We describe next these green modules. Main code. The core of the tool enables to:

• Gather the user arguments from the Arguments parser and manage the simulation parameters. • Instantiate the agents from the Agent class and create the threads for each node. • Run ns-3 simulator and tensorboard server in separated processes. • Log the simulation statistics using the Stats writer. Agent class. This class represents a node agent, containing two main methods: run_forwarder and run_trainer. The first one interacts with the ns3-gym environment as follows: (i) retrieve the observations, (ii) take action, (iii) compute the reward, (iv) store the states' transitions in the experience replay buffers, and, (v) update the information about the environment. The second method trains the model performing a gradient descent step. The run_forwarder also tracks the packets using their IDs to build the target value Y Q n at the node n, since the reward r, the next hop agent observation o n and the next hop agent Q-value Q n (•; θ n ) are associated to the packet ID.

Multi-threaded agent instances.

These are Agent class instances. Each one corresponds to a network node. They share different information among them and with the Main code using static variables. As aforementioned, we decide to allocate two threads for each instance: one calling run_forwarder method and the other one calling run_trainer method to guarantee that training and forwarding can take place separately. Stats writer. This module is called by the Main code and Agent class instances to log information about the environment or the training progress. The following variables are tracked to be visualized using tensorboard: average packet delay, loss ratio, Temporal Difference error at each agent, replay buffers' occupancy at each node, exploration ratio for each node, the number of new arriving packets, the number of delivered packets, and, the number of buffered packets. We also add the ability to define custom plots and to compare between execution's hyperparameters. Argument parser. This module is used to retrieve and parse the arguments given by the user in the script call. Poisson traffic generator. This is a ns-3 application class that generates random packets following a Poisson process at each node. The average rate is retrieved from a given traffic matrix and the packet length is fixed to a given value. The application is not installed in the node itself, but in a virtual one directly connected to the real node. Action and observation handler. As indicated in Section III-A, the action a n and the observation o n are defined as the output network interface index and the packet destination, respectively. For the local observation o n , other data such as the buffer occupancy of the outgoing interfaces Fig. 1: PRISMA code structure, separated in two parts : Python code and C++ code, highlighting our contribution in green. The figure shows the multi-threaded approach, where each node agent is composed of two threads: one forwarding the packets and another training could be considered. This module, implemented as a ns-3 handler, retrieves o n from the network simulation when a new packet arrives at the node n, and forwards the packet to the output interface a n . If the buffer interface is full, then the packet is discarded. We also provide the packet ID in the info field, so we can track each packet. MARL lib. With this module, we extend the OpenAI™ Baselines library [START_REF]Openai baselines: high-quality implementations of reinforcement learning algorithms[END_REF] to the Multi-Agent Deep Reinforcement Learning (MA-DRL) approach. The MARL lib provides the tools for defining, training, and testing an agent. It contains the following modules:

• models: a module containing the DNN models. V. USAGE In this section, we describe the parameters of PRISMA, and we present the guidelines about how to use the framework.

A. Framework parameters

PRISMA parameters are separated into the following:

• Global simulation arguments: it contains the simulation time, the base port, the seed, and whether to run training or not. • Network parameters: it contains network attributes like the path to the adjacency matrix, the maximum output buffer size, or the load factor.

• DRL agent argument: it concerns the training of the agent and contains parameters like the batch size, the learning rate, or the training time step. • Session logging arguments: it contains parameters about the session, like the name of the session and the path to store the result of the simulation. • Other parameters: some misc arguments like whether to run a tensorboard server and its port number.

B. Usage guide

First, a user must install the ns3-gym, operational system and Python dependencies. To do that, we provide installation files.

The main prisma folder is composed of four folders: source, ns-3, examples, and scripts. The source folder contains the MARL lib, the agent class and utils modules. The ns3 folder contains the ns-3 files: (i) the Poisson data generator, and; (ii) the ns-3 simulation itself (sim.cc). The latter file creates and runs the simulation scenario defined by the examples folder files, that define the network topologies and traffic matrices. Scripts for launching the simulations of Section VI are provided in the scripts folder.

The DNN model can be changed in models.py. The reward, observation, and action definitions can be modified in the methods _get_reward (in agent_class), Get_Observation (in packet-routing-gym.cc) and Get_action (in packet-routing-gym.cc), respectively. Finally, a user may test the framework by calling main.py with the corresponding arguments. Calling the latter file with "-help" argument will show all the possible parameters.

VI. ILLUSTRATIVE EXAMPLE

This section aims to show how PRISMA can be used to assist the training and testing of a DRL algorithm to solve the Distributed Packet Routing problem. Namely, we make a Deep Q-Network routing model [START_REF] Mukhutdinov | Multi-agent deep learning for simultaneous optimization for time and energy in distributed routing system[END_REF] learn the Shortest Path (SP) routing in backbone network scenarios. The tool can be applied to more challenging ad-hoc wireless networks by simply modifying the sim.cc file and properly adapting the DQN routing model. Here, we use a more simple case for illustrative purposes.

A. Simulation settings

We ran PRISMA in two different machines: (i) a Dell Precision 7920 workstation equipped with an Intel Xeon Gold 6230R Dual CPU (26 Cores, 2.1-4.0GHz Turbo, 128 GB RAM) with 2 NVIDIA RTX A5000 GPUs; and, (ii) an Intel Inspiron 14 laptop equipped with an Intel Core i7-8565U CPU (Quad Core, 1.80 GHz, 8 GB RAM) with a 1 GPU NVIDIA GeForce MX130 GPU.

The DRL agent is implemented in TensorFlow [START_REF]Tensorflow: An end-to-end open source machine learning platform[END_REF] as a three-layer neural network in models.py with a first dense layer of 32 nodes, then two layers of 128 nodes each having a Leaky ReLU activation function. This model considers only the packet destination as input and outputs the estimated endto-end packet delay based on the selected network interface. The action is retrieved by applying an argmin function to the output layer. In the examples folder, we define the two topologies considered in this experiment: Abilene (11 nodes) and Geant (23 nodes), as well as the traffic matrices (randomly generated using a uniform distribution). We fix the packet size to 512 KB, the link propagation delay to 2 ms, and the maximum output buffer length to 30 packets. We chose the reward r n (the next hop packet delay) to be 1, which means that the agent policy will aim to minimize the total number of hops from source to destination These parameters along with the model hyperparameters are set as arguments when calling main.py. We expect that DQN-routing will have a similar performance to Shortest Path routing, since both algorithms minimize the number of hops taken by a packet to reach its destination.

B. Model training

In this section, we show how the training process is performed successfully. The model is trained during 1 minute at each 7 ms (network simulation time) with a learning rate of 0.001, batch size of 512, load factor of 50% and γ of 1. The duration of the simulation time for training was 30 seconds. Moreover, we use an -greedy approach ( decays from 1 to 0.01) to move from exploration to exploitation.

Thanks to tensorboard visualization tools (see Fig- cost represents the average hop count per packet. We can observe an early cost increase due to the initial exploration: the network is flooded with packets needing many hops to arrive at their destinations, since the forwarding decisions are mainly random. Afterward, we move progressively to the exploitation phase, where the decision comes from the Deep Neural Network (DNN) model. As we see, the packets reduce the hop count to reach their destination, since model decisions improve. Figs 3c and 3dshow the Temporal Difference (TD) error for each node over the simulation time. We can observe that this error is decreasing, which means that all the agents are converging (learning) to a feasible routing policy (packets are not forwarded indefinitely).

C. Model testing

In this part, we evaluate the performances of the trained DQN-routing agents at different traffic intensities (low, medium, high, and very high traffic, in which the load factor was 50%, 100%, 150% and 200%, respectively). In Figure 4, we compare the results with SP routing in terms of the average end-to-end delay and the packet loss ratio. For low and medium traffic rates, DQN-routing presents the same performance as SP since both methods are using the shortest path decision policy, which is sufficient to handle all the packets at this rate. For high traffic rate, DQN-routing outperforms SP routing in terms of the end-to-end delay for both topologies; and, in terms of loss ratio, for Abilene. On the contrary, SP routing is slightly better in loss ratio in Geant. Anyway, these differences are not significant since we intend to learn a routing close to the SP, which is the case.

Low Medium High Very High

Traffic Rate Intensity

VII. CONCLUSION

In this paper, we have presented the PRISMA tool, which is, to the best of our knowledge, the first Deep Reinforcement Learning (DRL) framework specifically devoted to the Distributed Packet Routing (DPR) problem. In this problem, routing agents take packet forwarding decisions based only on local information. When the network state is uncertain or partially hidden (e.g., a global view is not available) or changes very dynamically, DRL appears as a solution to discover an efficient routing policy. In this vein, PRISMA aims to provide a playground for researchers interested in the application of this machine learning paradigm to this challenging problem, allowing fast prototyping and bench-marking. We illustrated its main functionalities by applying the tool to learn in a distributed manner a Shortest Path routing policy in two backbone networks.

ure 2 )

 2 , we can watch the training progress. Figures 3 show in detail the most relevant metrics to watch the learning progress for Abilene and Geant. Figure3a and 3bshow the average cost over the simulation time. The cost is computed by dividing the sum of the rewards over time by the total number of packets. Since the reward is 1 at each hop, the

Fig. 2 :

 2 Fig. 2: Screenshot of tensorboard interface. The x axis represents the time in ns, and the y axis represents the variables tracked during the simulation. For the training metrics (bottom subfigures), each curve corresponds to an agent.

Fig. 3 :

 3 Fig. 3: The average cost per packet and TD errors at each node during training for Abilene (rightmost subfigures) and Geant (leftmost subfigures) topology.

Fig. 4 :

 4 Fig. 4: The average end-to-end delay per arrived packet and packet loss ratio for different traffic rate intensities for Abilene (right row) and Geant (left row) topology

  agent: a class defining an DQN agent in the context of DPR. It contains methods like step, train and sync neighbor target neural network • utils: a set of utility functions like save_model and load_model.

• Replay buffer: a class handling the experience replay buffer. It contains methods like: add, sample or save.

•