Nathan Decaux
email: nathan.decaux@imt-atlantique.fr

Pierre-Henri Conze

Juliette Ropars

Xinyan He

Frances T Sheehan

Christelle Pons

Ben Salem

Sylvain Brochard

François Rousseau

Semi-automatic muscle segmentation in MR images using deep registration-based label propagation

Keywords: semi-automatic segmentation, musculoskeletal system, label propagation, deep registration

Fully automated approaches based on convolutional neural networks have shown promising performances on muscle segmentation from magnetic resonance (MR) images, but still rely on an extensive amount of training data to achieve valuable results. Muscle segmentation for pediatric and rare diseases cohorts is therefore still often done manually. Producing dense delineations over 3D volumes remains a time-consuming and tedious task, with significant redundancy between successive slices. In this work, we propose a segmentation method relying on registration-based label propagation, which provides 3D muscle delineations from a limited number of annotated 2D slices. Based on an unsupervised deep registration scheme, our approach ensures the preservation of anatomical structures by penalizing deformation compositions that do not produce consistent segmentation from one annotated slice to another. Evaluation is performed on MR data from lower leg and shoulder joints. Results demonstrate that the proposed semi-automatic multi-label segmentation model outperforms state-of-the-art techniques.

Introduction

Pediatric musculoskeletal segmentation is a key pre-processing stage for clinical decision-making, as the generated 3D muscle and bone models enable the extraction of key biomarkers related to disease progression, which, in turn, optimizes identifying therapeutic interventions. In this context, manual segmentation remains the gold standard to achieve valuable delineations [START_REF] Muzic | Reproducibility of manual segmentation in muscle imaging[END_REF]. Indeed, the variability and scarcity of musculoskeletal data, especially in pediatric cohorts, is an obstacle to the generalization of automatic segmentation models. Thus, 3D atlas-based label propagation methods, whether based on registration [START_REF] Miller | Mathematical textbook of deformable neuroanatomies[END_REF][START_REF] Sambaturu | Scribblenet: Efficient interactive annotation of urban city scenes for semantic segmentation[END_REF] or local patch similarity [START_REF] Rousseau | A supervised patch-based approach for human brain labeling[END_REF], are not always well suited. To reduce the processing-time while maintaining a quality close to the interexpert variability, semi-automatic techniques requiring few annotated images appeared to be the best trade-off to date [START_REF] Ogier | Overview of MR image segmentation strategies in neuromuscular disorders[END_REF].

Fully supervised methods based on deep learning (DL) have been investigated with sparse annotated 2D slices [START_REF] Conze | Deep convolutional encoder-decoders for deltoid segmentation using healthy versus pathological learning transferability[END_REF][START_REF] Conze | Healthy versus pathological learning transferability in shoulder muscle MRI segmentation using deep convolutional encoder-decoders[END_REF]. Despite superior performance to atlas-based methods [START_REF] Chen | Deep learning vs. atlas-based models for fast autosegmentation of the masticatory muscles on head and neck CT images[END_REF], these DL-based approaches require a significant number of training examples [START_REF] Cheng | Automatic quadriceps and patellae segmentation of MRI with cascaded U2-Net and SASSNet deep learning model[END_REF]. Learning from few examples remains an open issue. Interactive approaches have also been explored in the field of DL for medical image segmentation purposes. Point or scribble-based methods allow for the least possible interaction, but suffer from a lack of reproducibility and thus require complex optimization procedures [START_REF] Feng | Interactive few-shot learning: Limited supervision, better medical image segmentation[END_REF][START_REF] Zhang | Interactive medical image segmentation via a point-based interaction[END_REF][START_REF] Sakinis | Interactive segmentation of medical images through fully convolutional neural networks[END_REF]. Chanti et al. explored the use of long short-term memory (LSTM) units to propagate an annotated sub-volume through 3D ultrasound images, but still requires hundreds of annotated slices for training [START_REF] Al Chanti | IFSS-Net: Interactive few-shot siamese network for faster muscle segmentation and propagation in volumetric ultrasound[END_REF].

As an alternative, intra-volume semi-automatic segmentation enables to deal with the limited size of annotated data [START_REF] Ogier | Overview of MR image segmentation strategies in neuromuscular disorders[END_REF]. In recent medical imaging literature, semi-automatic segmentation by label propagation, which involves spreading annotated slices throughout a volume or sequence, has not been widely studied. This technique is more commonly seen in the field of video segmentation, where the first annotated frame is propagated throughout the sequence [START_REF] Jabri | Space-time correspondence as a contrastive random walk[END_REF][START_REF] Sun | Munet: Motion uncertaintyaware semi-supervised video object segmentation[END_REF]. In the context of medical images, morphological-based interpolation of distant annotated slices provides a fast strategy to propagate labels over a 3D volume, but may also require an important number of interactions and often fails on multi-structure objects [START_REF] Azimbagirad | Towards the definition of a patient-specific rehabilitation program for TKA: A new MRI-based approach for the easy volumetric analysis of thigh muscles[END_REF][START_REF] Albu | A morphology-based approach for interslice interpolation of anatomical slices from volumetric images[END_REF]. In addition, this method does not capture local structure variations, as it relies solely on segmentation contours. To address this issue, Ogier et al. investigated a 2D registration framework for propagating annotated slices towards 3D segmentation maps [START_REF] Ogier | Individual muscle segmentation in MR images: A 3D propagation through 2D nonlinear registration approaches[END_REF]. In the same spirit as [START_REF] Feng | A dual propagation contours technique for semi-automated assessment of systolic and diastolic cardiac function by cmr[END_REF][START_REF] Khalvati | Inter-slice bidirectional registration-based segmentation of the prostate gland in MR and CT image sequences[END_REF] for static or dynamic images, the approach proposed in [START_REF] Ogier | Individual muscle segmentation in MR images: A 3D propagation through 2D nonlinear registration approaches[END_REF] allowed to obtain contour delineations for a sequence of slices delimited by two annotated slices. The anatomical coherence is reinforced by exploiting both distant deformations between annotations and successive deformations between images, weighting each contribution according to their distance to the closest annotated slices.

Our approach differs from these previous works by using an unsupervised deep learning framework for deformation generation. Numerous studies have been conducted in the topic of deep learning registration, with promising outcomes, especially for its generalizability and its lower processing time compared to direct optimization methods [START_REF] Fu | Deep learning in medical image registration: a review[END_REF][START_REF] Tustison | Learning image-based spatial transformations via convolutional neural networks: A review[END_REF][START_REF] Wei | Deep morphological simplification network (ms-net) for guided registration of brain magnetic resonance images[END_REF]. Among these methods, Voxelmorph reaches similar performances of popular 3D free-form registration algorithms between subjects such as ANTs SyN [START_REF] Balakrishnan | Voxelmorph: a learning framework for deformable medical image registration[END_REF][START_REF] Avants | Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain[END_REF]. In this study, we propose to use a similar approach to perform 2D registration. One key distinction between this study and the work of Voxelmorph is the focus on 2D registration within a single volume, rather than on 3D registration between multiple volumes.

In this work, our objective is to develop an intra-subject 3D segmentation method for pediatric muscles from magnetic resonance (MR) images based on very few manually annotated 2D slices. To this end, we investigate the use of a learning-based registration framework to propagate labels through the full volume. More precisely, the 3D segmentation problem is modeled as a 2D label propagation problem based on unsupervised DL-based registration. The registration approach relies on intensity similarity between successive slices and on muscle shapes from annotated slices. A regularization term is introduced via the definition of a dedicated loss from combined deformation fields. Furthermore, propagated masks from different manual segmentations are merged through a novel weighting method based on image similarity measures. The proposed approach is evaluated on two clinical datasets and compares favorably to state-of-the-art methods. The PyTorch implementation of this work is publicly available1 .

Methods

Let us consider a volumetric MR image as a stack of 2D slices. The use-case considered in this work consists of a clinical expert who minimally delineates muscles of interest in few slices only. The objective is to provide a full 3D segmentation map of the muscles by propagating the 2D annotations. We model this problem into an intra-patient registration-based propagation framework. We denote a full volume X as a set of stacked 2D images X = (x 1 , x 2 , ..., x n ) and the corresponding segmentation map Y as a set of stacked sparsely annotated 2D images Y = (y 1 , y 2 , ..., y n ), i.e. the majority of the slices in the set are not annotated. In such a registration-based propagation framework, the key assumption is that muscle shapes in two consecutive 2D images are highly similar. As recently highlighted in [START_REF] Ogier | Overview of MR image segmentation strategies in neuromuscular disorders[END_REF], methods based on intra-subject non-linear registration currently provide state-of-the-art results. The purpose of the proposed approach is to propagate a small amount of annotated 2D slices through the volume X to estimate a full 3D muscle segmentation map. To this end, we propose a DL-based approach that simultaneously registers consecutive 2D MR slices and propagates annotated muscle contours.

Unsupervised 2D intensity-based registration

The proposed label propagation framework relies on non-linear registration of consecutive 2D slices. Let ϕ k→k+1 be a non-linear mapping that maps coordinates of x k to coordinates of x k+1 . Such a mapping ϕ k→k+1 is usually estimated by solving the following optimization problem: where L sim k→k+1 (X , ϕ) is an image dissimilarity measure between the consecutive slices x k and x k+1 , L smooth a regularization term and λ 1 a trade-off weight between both terms. In this work, the non-linear mapping ϕ is modeled as ϕ = Id + u where Id is the identity transform and u is the displacement field. We consider a diffeomorphic mapping ϕ through the integral of a stationary velocity vector field v such that ϕ preserves the topology and is invertible [START_REF] Ashburner | A fast diffeomorphic image registration algorithm[END_REF]. Moreover, the inverse mapping ϕ -1 k→k+1 is obtained using the negative velocity field -v, computed with a similar process as in forward integration.

φk→k+1 = arg min ϕ L sim k→k+1 (X , ϕ) + λ 1 L smooth (ϕ), (1)
L sim k→k+1 (X , ϕ) penalizes the dissimilarity between the slice x k deformed using ϕ k→k+1 and the slice x k+1 , and reciprocally using ϕ -1 k→k+1 . To be robust to any intensity variations inside the full volume X , the dissimilarity loss used in this work is the normalized local cross-correlation (NCC):

L sim k→k+1 = -NCC(ϕ k→k+1 (x k ), x k+1 ) -NCC(ϕ -1 k→k+1 (x k+1 ), x k ). (2)
To ensure smooth deformation fields, the regularization term L smooth penalizes for each pixel p the spatial derivatives of the velocity field v [START_REF] Balakrishnan | Voxelmorph: a learning framework for deformable medical image registration[END_REF]: L smooth (ϕ k→k+1 ) = p ∥∇(v k→k+1 (p))∥ 2 . Similarly to VoxelMorph [START_REF] Balakrishnan | Voxelmorph: a learning framework for deformable medical image registration[END_REF], the deformation field is computed in an unsupervised way, where the velocity field is modeled using a neural network.

Propagation constraints

Our objective is to segment a full MR volume by propagating 2D annotations using estimated deformation fields between consecutive slices. We propose to consider a loss related to segmentation to make the registrationbased label propagation more accurate and anatomically plausible. Let consider two annotated slices x i and x j and their corresponding annotations y i and y j . The slice x i can be mapped to x j using ϕ i→j , the combination of nonlinear mappings ϕ i→j = ϕ j-1→j • ... • ϕ i→i+1 . To improve the anatomical correspondences, we consider a overlap-based segmentation loss L seg i→j between the deformed annotated slice ϕ i→j (y i ) and y j : L seg i→j = -DSC(ϕ i→j (y i ), y j ) where DSC is the Dice similarity coefficient quantifying the surface overlap for a given muscle. The combination of non-linear mappings can lead to an accumulation of registration errors between annotated sections. Therefore, to regularize the estimated displacement fields through the propagation, we propose to investigate the use of the following intensity-based criterion:

L comp i→j = -NCC(ϕ i→j (x i ), x j ).

Label propagation framework

Considering two slices (x i , x j ) and their annotations (y i , y j ), the corresponding loss function L i→j is defined as follows:

L i→j = 1 j -i j-1 k=i L sim k→k+1 + λ 1 L smooth (ϕ k→k+1 ) + λ 2 L seg i→j + λ 3 L comp i→j , (3)
with weighting factors λ 2 and λ 3 . For more robustness, label propagation can be performed in both directions by considering a bidirectional loss L = L i→j + L j→i . The generated dense annotations are then merged by giving more weight to the pixels that are most likely to be correctly registered:

ŷk = β k (X , ϕ).(ϕ i→k (y i )) + (1 -β k (X , ϕ)).(ϕ j→k (y j )), (4)
with β k the weighting function. As an approximation of the DSC evolution along the propagation axis, [START_REF] Ogier | Individual muscle segmentation in MR images: A 3D propagation through 2D nonlinear registration approaches[END_REF] suggested using

β k = 1-arctan(C(k-(j-i)/2))
π as a weighting function, where C is a smoothing factor. The approach therefore assumes a symmetry in the propagation quality in both directions, which may not be the case. This function is denoted as distance-based weighting (DW) in Sect. 3.7. Instead of distance, we hypothesize a correlation between image similarity and segmentation quality by considering now β k (X , ϕ) as a function of the propagated slices ϕ i→k (x i ) and ϕ j→k (x j ), such that :

β k (X , ϕ) = sim(ϕ i→k (x i ), x k ) sim(ϕ i→k (x i ), x k ) + sim(ϕ j→k (x j ), x k ) . (5)
Three sim(ϕ i→k (x i ), x k ) functions as composition of the pixel-level function NCC(., .) have been investigated: (1) the average cross-correlation map p:

1 N N p NCC(ϕ i→k (x i ), x k )(p)
, denoted as correlation-based weighting (CW);

(2) the average of the masked map

N p ϕ i→k (y i )(p)NCC(ϕ i→k (x i ),x k )(p) N p ϕ i→k (y i )(p)
, denoted as masked correlation-based weighting (MCW); (3) the pixel-level map NCC(ϕ i→k (x i )) where each pixel is weighted independently, denoted as pixel-wise correlationbased weighting (PCW) in Sect. 3.7.

Experiments and results

MR muscle datasets

We evaluate the proposed approach on two muscle MR datasets, acquired from two distinct anatomies: shoulder and thigh.

The first one (referred to as MR Shoulders) is a shoulder dataset collected from a previous IRB approved study [START_REF] Pons | Shoulder muscle atrophy and its relation to strength loss in obstetrical brachial plexus palsy[END_REF] originally established to investigate the muscle volume-strength relationship in 12 children with unilateral obstetrical brachial plexus palsy (averaged age of 12.1±3.3 years). Informed consents from a legal guardian and assents from the participants were obtained for all subjects. Data from children with atrophic muscles exhibit a large variability in muscle structure and texture between subjects, which represents a challenge for learning-based segmentation algorithms. For each subject, 3D axial-plane T1-weighted gradient-echo MR images were acquired for the affected shoulder. For each MR volume, a sparse set of 2D axial slices were selected to delineate four rotator cuff muscles: deltoid, infraspinatus, supraspinatus and subscapularis (Fig. 4, top). Manual annotations were performed by an expert in physical medicine and rehabilitation. Size for axial slices is constant for each subject (416×312 pixels). Resolution varies from 0.55×0.55 to 0.63×0.63mm. The number of axial slices fluctuates from 192 to 224 while slice thickness is fixed to 1.2mm. The average axial slice number of interest is 106±28 for deltoid, 63±16 for infraspinatus, 24±8 for subscapularis and 71±16 for supraspinatus. The second dataset employed to evaluate our method is the publicly available MyoSegmenTUM dataset [START_REF] Schlaeger | Thigh muscle segmentation of chemical shift encoding-based water-fat magnetic resonance images: the reference database myosegmentum[END_REF], consiststing of Dixon water MR scans of thighs, with dense segmentation of four muscles: quadriceps femoris, sartorius, gracilis, hamstring, differentially marked for left and right thigh (Tab.4, bottom). The dataset is composed of a population of 15 healthy adults (80% of males) with age, weight, and height of resp. 29±8 years old, 86±13kg, and 180±9cm. Image size varies from 400×400 (63 axial slices) to 560×560 (65 axial slices), while resolution is 0.93×0.93 with slice thickness of 4mm. The average number of annotated axial slices is 48±10, all muscles included.

Implementation details and evaluation

A U-Net architecture is used for the velocity field generator, designed as a contraction/expansion path with skip-connections. It takes as input two successive images of the volume and estimates as output two maps representing the velocity fields in both directions of the plane. The displacement field is obtained by repeatedly integrating the velocity field over small-time steps, following the method described in the DARTEL [START_REF] Ashburner | A fast diffeomorphic image registration algorithm[END_REF] paper. The velocity field is first scaled by a factor of 1 2 N where N is the number of squaring steps. Then, the scaled field is integrated by self-application N times, which is considered an approximation of the computation of the exponential of the velocity field. The number of squaring steps is set to 7, as suggested in DARTEL method [START_REF] Ashburner | A fast diffeomorphic image registration algorithm[END_REF]. The implementation used is a publicly available U-Net network from the MONAI library2 . The number of features is set to [START_REF] Azimbagirad | Towards the definition of a patient-specific rehabilitation program for TKA: A new MRI-based approach for the easy volumetric analysis of thigh muscles[END_REF]32,32,32) and the number of residual units to 2, other parameters being kept by default. Final convolution layer, velocity field integration, deformation field mapping function and original losses are based on the official Voxelmorph Pytorch implementation3 . As in [START_REF] Balakrishnan | Voxelmorph: a learning framework for deformable medical image registration[END_REF] when using NCC, the weighting parameter λ 1 (Eq. 1) is set to 1. λ 2 and λ 3 is empirically set to 1. NCC sliding-window size is set to 9 × 9 in L sim and L comp while set to 31 × 31 in Eq. 5. The total number of parameters is similar to the one involved in Voxelmorph (114k versus 110k).

For all experiments, the model is trained using the Adam optimizer with a learning rate fixed at 10 -3 and a weight decay of 10 -8 . The best trained model over 200 epochs is selected by keeping the one that provides the highest Dice score between propagated labels and ground-truth annotated training slices. The prediction evaluated is the segmentation map resulting from the fusion of the propagation in both directions, according to Eq. 4, using PCW function.

To provide a comprehensive assessment, the metrics employed for evaluation purposes are the Dice coefficient (DSC), the average surface distance (ASD) as well as the Hausdorff distance (HD), computed in 2D according to the axial plane.

Existing methods

The proposed approach is compared with two segmentation methods: a morphology-based interpolation method [START_REF] Albu | A morphology-based approach for interslice interpolation of anatomical slices from volumetric images[END_REF] as well as a supervised U-Net deep segmentation model [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF]. The morphology-based approach is an iterative method that can explicitly handle inter-slice topology changes by decomposing a many-to-many correspondence. This method is very suitable in the scenario considered in this work, namely a semi-automatic segmentation strategy which propagates annotated slices within the same volume. This can thus be considered as a reference method. A publicly available code is provided by Kitware in ITK 4 . This approach is denoted as Morph in Tab. 1 and Tab. 2.

The second method is a supervised DL technique, relying on the U-Net architecture, similar to previous studies [START_REF] Conze | Healthy versus pathological learning transferability in shoulder muscle MRI segmentation using deep convolutional encoder-decoders[END_REF][START_REF] Ahmad | Semantic segmentation of human thigh quadriceps muscle in magnetic resonance images[END_REF][START_REF] Ding | Deep learning-based thigh muscle segmentation for reproducible fat fraction quantification using fat-water decomposition MRI[END_REF]. The network is trained using a cross-entropy loss and data augmentation with random affine transformations (rotations between -20 and 20 degrees, scaling from 0.8 to 1.2 and shearing from -20 to 20). In the context of semi-automatic segmentation, the U-Net model is trained during 200 epochs and tested in a supervised learning fashion for each subject independently. The implementation is based on a publicly available Pytorch implementation 5 . It follows the original architecture and number of filters, with the addition of a batch normalization operation after each convolution. The number of parameters is approximately 14 million. This approach is denoted as U-Net in Tab. 1 and Tab. 2.

The final method presented in the paper is a semi-automatic deep learning approach for video label propagation tasks. It uses space-time graphs to represent video sequences, where patches from each frame are the nodes and transitions between patches in adjacent frames are the edges. The model learns a representation that defines the transition probability of a random walk on the graph. The representation is optimized to place high probability along paths of similar nodes, and is trained using cycle-consistency without supervision. At test time, the first frame annotation is propagated to the entire sequence using the most likely transitions for each patch. To remain fair to our method, we apply the same label propagation framework during inference, described in Sec. 2.3, propagating for each subsequence (i.e. between two annotated slices) in both directions, and merging the predicted segmentations as a function of distance from the nearest annotated slice (see DW function in Sec. 2.3). A model is trained for each dataset separately, using every sequence without annotation. Models are initialized using pretrained weights provided in the public implementation6 of the method, and trained during 200 epochs. During training, sequences are divided into clips of 8 frames, dropout is set to 0.1, softmax temperature to 0.05, and learning-rate and batch size are respectively set to 10 -4 and 32. For data augmentation, we use the random flip function provided in the code, which provides the best performance in our case. Other hyperparameters are set to default. The number of parameters is approximately 2.7 million. This approach is denoted as VideoWalk in Tab. 1 and Tab. 2.

Intra-Subject Segmentation with Minimal Supervision

In this section, we consider a difficult intra-subject segmentation scenario: only three slices are annotated for each MR volume to segment. The middle slice allows a fine delineation of the muscle while the two other slices provide its spatial extent (in z-direction). In this case, the typical distance between two consecutive annotations is 32 ± 17 for shoulders and 24 ± 5 for thighs. These three slices are the only annotation seen for each methods, as we perform the experiment for each subject separately.

Tab. 1 reports the DSC, HD and ASD scores for the shoulder dataset averaged over each slice and muscle. With DSC scores of 76.6, 81.0, 80.7, and 72.7% for deltoid, infraspinatus, subscapularis, and supraspinatus muscles, respectively, our method outperforms interpolation-based techniques [START_REF] Albu | A morphology-based approach for interslice interpolation of anatomical slices from volumetric images[END_REF], U-Net [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF] and VideoWalk [START_REF] Jabri | Space-time correspondence as a contrastive random walk[END_REF] for all muscles. Significant improvements in ASD stability (3.0±4.6 versus 5.3±16.8 mm with [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF] for supraspinatus) and HD (18.8±16.1 versus 32.8±26.4 mm with [START_REF] Albu | A morphology-based approach for interslice interpolation of anatomical slices from volumetric images[END_REF] for infraspinatus) are achieved. Overall, VideoWalk shows the lowest performance on this shoulder dataset, presumably because the approach is not suitable for the propagation of fine structures on low contrast images such as these T1w images. . The DSC scores for supraspinatus are overall lower than for the other muscles due to its thin and elongated shape. For a more complete assessment, we also include in Tab. 1 the results from [START_REF] Conze | Healthy versus pathological learning transferability in shoulder muscle MRI segmentation using deep convolutional encoder-decoders[END_REF], based on a U-Net model trained in a leave-one-out (LOO) manner, i.e., using every subject except the one tested. Compared to [START_REF] Conze | Healthy versus pathological learning transferability in shoulder muscle MRI segmentation using deep convolutional encoder-decoders[END_REF], our method shows better results for all muscles despite a much smaller amount of annotated data. A strong improvement in DSC is moreover visible for infraspinatus (from 71.4 to 81.0% in DSC) although only 3 slices are used for label propagation.

Similar conclusions can be noted for MyoSegmenTUM in Tab. 2, where the proposed approach consistently outperformed other methods in DSC, ASD, and HD metrics. Specifically, we find that small muscles such as gracilis and sartorius challenge both U-Net [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF] and interpolation-based technique [START_REF] Albu | A morphology-based approach for interslice interpolation of anatomical slices from volumetric images[END_REF] in this semi-automatic intra-subject segmentation case (10.2 for [29] and 30.6% for [START_REF] Albu | A morphology-based approach for interslice interpolation of anatomical slices from volumetric images[END_REF] in average DSC for gracilis and sartorius muscles and left/right legs versus 81.4% for our approach). Thanks to the good contrast of the water sequences in this dataset, the VideoWalk [START_REF] Jabri | Space-time correspondence as a contrastive random walk[END_REF] approach achieves better performance on these small muscles in DSC, but remains significantly less accurate compared to our method (47.0% for gracilis and sartorius muscles and left/right legs versus 81.4% for our approach).

The proposed method shows similar ASD scores with [START_REF] Albu | A morphology-based approach for interslice interpolation of anatomical slices from volumetric images[END_REF] on larger muscles, but shows significant improvements in HD (12.9±12.0 versus 30.7±27.2 mm for [START_REF] Albu | A morphology-based approach for interslice interpolation of anatomical slices from volumetric images[END_REF] on right quadriceps femoris). The VideoWalk method [START_REF] Jabri | Space-time correspondence as a contrastive random walk[END_REF] shows similar performance for the global quadriceps femoris in terms of DSC, with a score of 82.9% for their method compared to 86.7% for our method. However, their method did not perform as well in terms of ASD and HD, with scores of 5.3 in ASD and 26.7 in HD compared to our scores of 3.1 in ASD and 13.0 in HD.

Moreover, Fig. 5 shows the evolution in the z-direction of the averaged DSC of both datasets. It appears that the proposed label propagation is less sensitive to distance from the nearest annotated slice, compared to interpolationbased approaches [START_REF] Albu | A morphology-based approach for interslice interpolation of anatomical slices from volumetric images[END_REF], U-Net [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF] and VideoWalk [START_REF] Jabri | Space-time correspondence as a contrastive random walk[END_REF]. Results also demonstrate a higher overall robustness with our method, as evoked by the standard deviation (colored area around the mean curves).

Ablation study

We now evaluate the contribution of L seg and L comp losses and the propagation schemes: forward (F) model, backward (B) model (i.e. propagated in the opposite direction as F) and the fusion of these two predictions as described in Eq. 4). To this end, 7 annotated slices per subject are used for training. Results are reported in Tab. 3 for the two datasets. First, it can be seen that, as expected, adding a segmentation-based loss improves the label propagation quality, especially while propagating in a single direction (82.3±14 versus 71.1±21.2 in DSC for MR Shoulders in backward direction and 90.5±7.1 versus 83.8±15.4 for MyoSegmenTUM in forward direction). Second, adding a composition loss to regularize the estimation of the deformation fields does not lead to a significant improvement in the segmentation results. Finally, the bidirectional fusion strategy allows in a large majority of cases an improvement of the segmentation results, and more particularly when L seg is not used (DSC scores of 82.2±13.9 versus 73.7±19.2 and 71.1±21.2 in single directions for MR Shoulders, 90.1±7.4 versus 83.8±15.4 and 86.5±10.3 for MyoSegmenTUM).

Pre-training

The training stage of the proposed method can be performed in a subjectspecific label propagation context for each volume independently, similarly to [START_REF] Ogier | Individual muscle segmentation in MR images: A 3D propagation through 2D nonlinear registration approaches[END_REF]. However, it is possible to train the neural network from the available unannotated data set as well. We focus in this section on the influence of such a pre-training step. More precisely, we evaluate performance and convergence time with or without pre-training. For a given dataset and a given propagation direction, the pre-training stage corresponds to a joint learning of every subject successive deformations, without any annotation. This corresponds then to minimizing Eq. 1 for each slice of each subject, in both directions. We consider the minimal supervision learning scenario described as in Sect. 3.4. Results for both datasets are reported in Tab. 4. The experiment shows the interest of pre-training in a difficult context of learning from only 3 slices, with improvements in mean DSC score from 2.7 to 7.2% for MR Shoulders, and from 4.0 to 22.6% for MyoSegmenTUM, as well as for other metrics.

DSC convergence over epochs on the MR shoulder dataset is studied in

Study of weighting strategies

In this experiment, we compare performances of each weighting strategies described in above Sect. 2.3. We study their influence on both datasets, and with a different number of annotated slices per subject : 3, 5 and 7.

Results for DSC, ASD and HD scores are reported in Tab. 5. Pixel-wise fusion (PCW) shows the best average improvements compared to the baseline distance fusion, with DSC gains from 0.1 to 1.2% and up to -1.3 in HD. In MR Shoulders, the local similarity strategy (MCW) reaches the lowest ASD even though it gives the lowest DSC as well.

Visual comparison between weighting strategies based on distance and pixel-wise similarity is provided in Fig. 9. The PCW approach shows more robustness to structure splitting of the quadriceps femoris (blue label) and hamstring. slices are used, while the number of annotated slices is set to 7 in Sec. 3.5. In Sec. 3.7, we have studied the influence of the number of annotated slices per volume: 3, 5, and 7 slices. In the previous sections, the experiments were performed with a fixed number of annotated slices, rather than an interval between annotations, to demonstrate the application of our approach for a fixed interaction time, invariant to the subject morphology. The advantage of this approach is that the user knows how long it takes to segment a dataset using the proposed propagation method. However, the same anatomical region of interest may occupy a different number of slices depending on the subject's height, weight, or age. Longer muscles, for a given slice thickness, therefore imply larger propagation distances for a given number of annotated slices. In this section, we examine the performance of the proposed approach as a function of the distance between annotated slices, with the hypothesis that a greater propagation induces an accumulation of registration errors and thus a lower segmentation accuracy. Figure 10 shows the results of the experiment described in Sec. 3.7 with the pixel-wise correlation-based weighting PCW strategy for the three configurations of the number of annotated slices (3, 5 and 7) according to the propagation distance. Each point corresponds to the average DICE for a given subject and muscle as a function of the distance between two annotated slices used for propagation. As expected, the smaller the distance between the annotated slices, the better the segmentation is propagated. For example, it can be noted that for the deltoid, the DICE remains above 75% of DSC despite large distances between annotated slices (up to 69 slices). Figure 10 shows the variability of performance depending on the muscle to be segmented. The choice of the optimal distance between the annotated slices depends on the anatomy of the subject, the size of the structures to be segmented but also the contrast in the images. These results require further investigation. Indeed, these results should be interpreted with caution, as the DSC coefficient is sensitive to the size of the structure to be segmented. The determination of the optimal distance for each muscle and each subject remains an avenue to be explored for future work, for example with a view to using the approach in a routine clinical case.

Conclusion

In this paper, we describe a deep registration-based label propagation method for intra-subject muscle segmentation from MR images. The propagation process is guided by image intensity, muscle shape and registration consistency. This operation is performed in both directions, resulting in two segmentation predictions that are merged with image registration quality estimation as weighting guidance. To improve performance and speed up convergence, an unsupervised pretraining stage is suggested. To obtain smooth and accurate fusions of segmentations from bidirectional propagation, several weighting functions have been studied. Experiments on clinical data show that 3D full scan segmentation can be achieved from very limited manual annotations using the proposed approach. One limitation of using diffeomorphic registration for anatomical structure segmentation is that it does not adequately handle topological discontinuities. In such cases, the annotated slices must accurately describe these changes in topology. A potential solution to this problem is to incorporate an active learning mechanism that guides the annotator in selecting the most informative slices for manual segmentation. On the other hand, our approach relies on context information in the form of annotated extreme slices. This requirement can be challenging for annotators because extreme slices are often difficult to segment manually and provide little information about the label structure. One possible way to address this limitation is to investigate the use of weak annotations, such as clicks or scribbles, for these slices. In our study, we focused on single-axis propagation, which is more applicable to MR images where only one plane is likely to be isotropic. However, a multi-view propagation approach could be considered for isotropic volume segmentation, such as in CT or 3D ultrasound imaging. In the context of musculoskeletal medical images analysis, for which there is very little annotated data, the described approach can serve as a basis for the development of interactive user-friendly segmentation tools. It is also very likely adaptable to other modalities that benefit from high spatial resolution in all planes (CT, 3D US) or organs.

Acknowledgment

This work was partially funded by ANR (AI4Child project, grant ANR-19-CHIA-0015) and Innoveo from CHRU de Brest.

Compliance with ethical standards

Data acquisition was performed in line with the principles of the Declaration of Helsinki. Informed consents from the participants were obtained for all subjects. Authors declare that they do not have any conflicts of interest.

Figure 1 :

 1 Figure 1: Deep registration-based label propagation from an annotated slice {x i , y i } to the next one {x j , y j } towards dense muscle segmentation in MR images. Please refer to the text for complete description of notations.

Figure 2 :

 2 Figure 2: Unsupervised 2D intensity-based registration model between two successive slices x k and x k+1 . Please refer to the text for complete description of notations.

Figure 3 :

 3 Figure 3: Label propagation framework scheme. Labels are propagated in both direction, resulting in two prediction for intermediate slices, that are fused using a weighting function β. Please refer to the text for complete description of notations.

Figure 4 :

 4 Figure 4: Sample images from the two evaluated shoulder MR dataset in each plane, with corresponding labels. (Top) 3D MR scan of a pediatric shoulder of a child with obstetrical brachial plexus palsy from the MR Shoulders dataset. Four muscles of the rotator cuff are represented in color: deltoid, infraspinatus, subscapularis and supraspinatus muscles. (Bottom) 3D MR scan of a healthy adult thigh extracted from the MyoSegmenTUM dataset. Four muscle groups are represented in color: gracilis, hamstring, quadriceps femoris and sartorius. Each muscle is differentiated by leg (left-right).

Figure 5 :

 5 Figure 5: Averaged Dice (DSC) for each annotated slice over the MR Shoulders (a) and MyoSegmenTUM (b) datasets, in a minimal supervision setting. DSC is displayed with respect to the normalized axial slice number obtained by linearly scaling slice number from [z min , z max ] to [0, 1] where z min (z max ) is the minimal (maximal) slice index displaying a muscle. Vertical dotted black lines represents the location of annotated slices used for training. Colored areas deal with standard deviation.

Fig. 6 .

 6 Fig. 6. It shows that the pre-training stage provides a good initialization (58.1 versus 41.1 at epoch 0). Moreover, the pre-trained models reach on

  average per epoch performance of non pre-trained models from 14 epochs (70.5 at epoch 14 versus 70.0 at epoch 175).

Figure 6 :

 6 Figure 6: Performance evolution over training epochs on MR shoulders dataset, with or without pretraining. Shown metric is the 2D DSC averaged over all subjects and muscles. Colored areas deal with standard deviation.

Figure 7 :

 7 Figure 7: Visual comparison (first row : axial, second : sagittal, last : coronal) on MyoSeg-menTUM dataset of existing methods and our approach in the minimal supervision setting (Sect. 3.3).

3. 8 .Figure 8 :

 88 Figure 8: Visual comparison (first row : coronal plane, last rows : axial) on MR Shoulder dataset of pretraining stage influence in the minimal supervision setting.

Figure 9 :

 9 Figure 9: Visual comparison of axial slice predictions for different weighting strategies on MyoSegmenTUM dataset, in the three annotated slices setting. Images surrounded by solid lines yellow lines are annotated slices used for training and propagation.

Figure 10 :

 10 Figure 10: Average Dice (DSC) per subject as a function of the propagation distance, over the MR Shoulders (a) and MyoSegmenTUM (b) datasets. Dotted lines show linear regressions per label, with 95% confidence interval in colored area.

Table 1 :

 1 Quantitative assessment over the MR Shoulder dataset in the minimal supervision setting. Best scores are highlighted in bold.

			no learning	intra-subject learning	LOO learning
	metric	muscle	Morph [17] U-Net [29] VideoWalk [14] Proposed	[7]
		deltoid	45.4±29.2 56.8±28.7	49.7±32.6	76.3±15.5	68.9±29.9
	DSC ↑	infraspinatus 53.4±27.6 67.0±27.7 subscapularis 55.6±30.5 66.3±25.5	46.1±31.6 56.7±25.9	81.0±11.5 80.7±12.1	71.4±24.7 78.1±18.1
		supraspinatus 38.6±25.1 50.2±29.0	39.1±29.1	72.7±18.3	64.9±28.0
		deltoid	3.4±5.7	3.2±9.0	4.8±4.9	2.9±3.2	-
	ASD ↓	infraspinatus subscapularis	3.6±7.5 2.1±4.2	1.8±4.3 2.3±7.8	7.3±8.8 7.8±5.9	1.9±1.3 2.4±1.2	--
		supraspinatus	4.0±5.5	5.3±16.8	7.4±6.6	3.0±4.6	-
		deltoid	39.6±32.3 42.7±44.2	37.8±35.2	27.6±27.4	-
	HD ↓	infraspinatus 32.8±26.4 25.7±29.6 subscapularis 18.6±12.9 20.6±21.4	36.1±28.8 31.2±20.6	18.8±16.1 12.8±9.3	--
		supraspinatus 33.9±19.5 35.9±33.3	33.7±21.5	22.5±17.9	-

Table 2 :

 2 Quantitative assessment over the MyoSegmenTUM dataset in the minimal supervision setting. Best scores are highlighted in bold.

				no learning	intra-subject learning
	metric	muscle	leg Morph [17] U-Net [29] VideoWalk [14] Proposed
		gracilis	R L	27.2±29.6 27.8±30.9	8.7±21.5 2.5±7.5	40.8±30.4 38.0±30.3	83.6±12.2 80.6±16.6
	DSC ↑	hamstring quadriceps femoris	R L R L	78.6±14.3 78.0±15.8 76.8±18.2 75.6±20.0	36.6±31.3 34.7±29.1 38.5±31.8 37.0±35.2	76.9±14.5 77.6±13.2 83.1±9.7 82.7±10.8	88.0±10.2 86.7±11.0 87.4±8.8 86.0±10.1
		sartorius	R L	36.1±31.0 31.5±27.4	12.1±20.6 17.7±27.0	56.9±25.5 52.5±27.4	82.9±16.8 78.5±21.4
		gracilis	R L	7.8±3.7 7.9±3.9	152.4±70.5 174.7±60.3	2.0±1.9 2.1±2.0	0.9±0.4 1.1±0.9
	ASD ↓	hamstring quadriceps femoris	R L R L	3.1±1.3 3.1±1.7 3.3±2.7 3.2±1.6	90.1±71.1 94.9±64.5 113.3±73.3 112.5±80.6	4.8±4.4 5.1±4.2 5.2±4.6 5.4±4.7	2.1±0.9 2.2±1.0 2.9±2.5 3.3±3.0
		sartorius	R L	9.3±6.2 9.5±5.8	137.2±76.0 116.4±83.7	2.8±2.9 2.7±2.1	1.5±2.7 1.7±1.7
		gracilis	R L	15.2±5.3 15.6±5.9	263.6±104.7 295.3±69.7	7.8±7.4 8.1±7.8	3.1±1.9 3.8±3.8
	HD ↓	hamstring quadriceps femoris	R L R L	16.5±12.0 210.8±131.6 16.9±12.8 229.2±114.2 30.7±27.2 251.6±114.7 28.4±26.8 249.8±141.7	20.3±15.1 21.3±14.9 24.9±22.4 28.6±28.3	10.3±8.4 11.3±10.1 12.9±12.0 13.1±12.8
		sartorius	R L	24.7±14.9 259.7±122.3 25.4±13.1 218.9±141.1	9.8±8.5 10.5±8.7	4.8±5.7 5.1±4.3

Table 3 :

 3 Ablation study. Influence of L seg and L comp losses for different propagation schemes. Best results for a given loss combination is in bold. The asterisk (*) indicates the best results over all loss combinations.

					MR Shoulders			MyoSegmenTUM
	direction L seg L comp	DSC ↑	HD ↓	ASD ↓	DSC ↑	HD ↓	ASD ↓
	F			73.7±19.2	22.4±23.0	5.1±8.2	83.8±15.4 8.0±10.7	1.6±1.0
	B	✗	✗	71.1±21.2	23.3±25.1	3.4±4.2	86.5±10.3 7.6±10.6	2.0±3.1
	Fused			82.2±13.9	14.0±15.8	2.3±4.3	90.1±7.4	5.4±8.1	1.3±1.8
	F			83.8±11.1	15.5±19.0	2.7±4.8	90.5±7.1	5.9±9.6 1.2±0.7*
	B	✓	✗	82.3±14.0	16.5±22.1	1.9±2.9	90.6±6.7	5.8±9.3	1.3±1.5
	Fused			84.4±11.6	13.1±16.9	2.0±3.7 91.1±6.4* 5.0±7.9	1.2±1.2
	F			76.4±18.4	21.8±23.0	4.8±7.9	84.5±15.6 8.1±10.7	1.7±1.1
	B	✗	✓	74.2±20.0	22.6±24.6	3.3±4.6	87.4±9.7	7.5±10.6	1.9±2.8
	Fused			83.1±13.3	14.0±15.8	2.3±4.2	89.9±8.3	5.6±8.5	1.4±1.8
	F			84.6±10.8	15.3±18.3	2.7±4.7	90.3±7.2	6.0±9.6 1.2±0.7*
	B	✓	✓	83.2±13.8	16.5±22.4	1.9±4.0	90.6±6.4	5.8±9.1	1.3±1.6
	Fused			85.3±11.2* 12.3±15.4* 1.8±3.4* 91.0±6.4 4.9±7.5* 1.2±1.2

Table 4 :

 4 Evaluation of pre-training step on MR Shoulders (a) and MyoSegmenTUM (b) datasets. Best results are highlighted in bold.

		DSC ↑		ASD ↓		HD ↓
	pre-training	No	Yes	No	Yes	No		Yes
	muscle
	deltoid	69.5±16.4 76.3±15.5 4.4±4.4 2.9±3.2 34.2±29.7 27.6±27.4
	infraspinatus 78.3±12.4 81.0±11.5 2.6±2.2 1.9±1.3 21.9±18.9 18.8±16.1
	subscapularis 77.8±13.3 80.7±12.1 3.2±1.8 2.4±1.2 15.2±9.7	12.8±9.3
	supraspinatus 65.5±21.9 72.7±18.3 4.7±6.1 3.0±4.6 25.4±19.5 22.5±17.9
				(a)
			DSC ↑		ASD ↓		HD ↓
		pre-training	No	Yes	No	Yes	No	Yes
	muscle	leg
	gracilis	L R	64.8±24.3 80.6±16.6 2.4±2.5 1.1±0.9 7.8±7.3 61.0±30.7 83.6±12.2 2.1±2.2 0.9±0.4 7.0±5.6	3.8±3.8 3.1±1.9
	hamstring	L R	82.5±12.6 86.7±11.0 3.0±1.4 2.2±1.0 13.0±7.8 11.3±10.1 84.0±11.6 88.0±10.2 2.9±1.5 2.1±0.9 12.7±7.2 10.3±8.4
	quadriceps femoris	L R	81.6±13.6 86.0±10.1 4.9±4.3 3.3±3.0 18.4±14.0 13.1±12.8 83.4±10.8 87.4±8.8 4.8±4.3 2.9±2.5 19.8±14.2 12.9±12.0
	sartorius	L R	61.7±28.2 78.5±21.4 4.4±6.3 1.7±1.7 12.9±13.1 61.9±29.8 82.9±16.8 4.2±5.7 1.5±2.7 12.4±12.9	5.1±4.3 4.8±5.7
				(b)

Table 5 :

 5 Performance table summarizing average metrics per dataset, number of annotated slices, and by weighting strategy. Please refer to the Sect. 2.3 for details about weighting strategies

	annotated slices strategy	DSC ↑	HD ↓	ASD ↓	DSC ↑	HD ↓	ASD ↓
		DW	75.9±16.7	22.8±24.7 2.6±3.1 83.8±14.6 9.3±12.5 1.9±1.3
	3	MCW CW	75.3±19.9 76.6±16.3	22.3±24.0 2.3±2.8 83.1±17.8 8.8±11.3 2.1±2.9 22.4±24.2 2.5±3.0 84.5±14.3 8.7±11.6 1.8±1.2
		PCW	77.1±15.5 21.5±20.4 2.6±3.2 84.3±14.3 8.2±9.3 2.0±2.1
		DW	82.5±13.1	16.1±20.8 1.9±2.8	90.7±7.6	5.1±7.6 1.1±1.0
	5	MCW CW	82.5±15.0 83.0±13.1	15.2±19.6 1.8±2.8 15.7±20.7 1.8±2.8	90.5±9.8 90.9±7.4	4.8±6.9 1.1±1.2 4.9±7.2 1.2±1.3
		PCW	83.3±11.9 15.2±17.5 2.2±3.7 90.9±7.3	4.9±7.4 1.2±1.4
		DW	84.8±11.6	13.4±18.2 1.7±3.0	92.0±6.8	4.3±6.8 1.0±1.1
	7	MCW CW	84.6±13.9 85.3±11.2 13.1±17.9 1.7±2.9 12.8±17.5 1.6±2.8	92.0±6.9 92.0±6.9	4.3±7.1 1.0±1.2 4.3±6.8 1.0±1.1
		PCW	85.3±11.2 12.3±15.4 1.8±3.4 92.1±6.7	4.4±7.4 1.0±1.3

https://github.com/nathandecaux/labelprop

https://docs.monai.io/en/stable/networks.html#unet

https://github.com/voxelmorph/voxelmorph

https://github.com/KitwareMedical/ITKMorphologicalContourInterpolation

https://github.com/milesial/Pytorch-UNet

https://github.com/ajabri/videowalk


Mathieu Andraud
email: mathieu.andraud@anses.fr

Pachka Hammami

Brandon Hastings Hayes

Jason Ardila Galvis

Timothée Vergne

Gustavo Machado

Nicolas Rose

Modelling African swine fever virus spread in pigs using time-respective network data: Scientific support for decision makers

Keywords: Transmission model, Infectious Diseases, Epidemiology, Surveillance, Control

la diffusion de documents scientifiques 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

Infectious diseases represent a major concern for livestock production due to economic consequences and related societal impacts as well as veterinary and public health issues. In the last decade, African swine fever (ASF) became the main concern for the European and Asian swine industry and more recently the North American one because of its reintroduction to the Caribbean [START_REF] Amat | Peste porcine africaine (PPA) en République dominicaine : déclaration le 28/07/2021 de deux premiers foyers -Actualisation[END_REF][START_REF] Gonzales | African swine fever in the Dominican Republic[END_REF][START_REF] Wang | Advance of African swine fever virus in recent years[END_REF]. Following its introduction into Georgia in 2007, ASF has gradually spread in Europe and is now considered endemic in some countries [START_REF] Chenais | Epidemiological considerations on African swine fever in Europe 2014-2018[END_REF][START_REF] Pejsak | Four years of African swine fever in Poland. New insights into epidemiology and prognosis of future disease spread[END_REF]. To date, only two countries that had been infected with the ASF virus have regained their free status: Belgium and the Czech Republic [START_REF] Danzetta | African Swine Fever: Lessons to Learn From Past Eradication Experiences. A Systematic Review[END_REF][START_REF] Emond | Doing the "dirty work": how hunters were enlisted in sanitary rituals and wild boars destruction to fight Belgium's ASF (African Swine Fever) outbreak[END_REF]. Recently, an introduction of the virus into the wild boar population at the German-Polish border was detected (Sauter-Louis et al.). At the moment of the study, Novembre 2021, this epidemics was still ongoing with geographical expansion in Germany despite the implementation of control measures. Although limited to wildlife for about 10 months after introduction, the first case in domestic pigs was registered in July 2021 (Federal Ministry of Food and Agriculture, 2021). Understanding the spreading patterns of infectious agents at local and global scales is necessary to design relevant control measures [START_REF] Beaunée | Modelling of paratuberculosis spread between dairy cattle farms at a regional scale[END_REF]. The identification of pathogen transmission routes and the assessment of their relative role in the transmission process are essential [START_REF] Andronico | Highly pathogenic avian influenza H5N8 in south-west France 2016-2017: A modeling study of control strategies[END_REF][START_REF] Ezanno | How mechanistic modelling supports decision making for the control of enzootic infectious diseases[END_REF]. Although contacts between individuals and contributions from the environment are critical for pathogen transmission at farm level, additional contact structures -such as animal movements and indirect contacts through transportation or geographical/social proximity between actors -were highlighted as pivotal in the study of transmission of animal diseases at a large geographical scale [START_REF] Büttner | Epidemic Spreading in an Animal Trade Network -Comparison of Distance-Based and Network-Based Control Measures[END_REF][START_REF] Galvis | Modelling the transmission and vaccination strategy for porcine reproductive and respiratory syndrome virus[END_REF].

Social network analysis has the potential to shed light on contact structures, which, joined to epidemiological data, provided insights into the role of specific actors in the transmission process for livestock infections [START_REF] Hardstaff | Livestock trade networks for guiding animal health surveillance[END_REF][START_REF] Salines | Pig movements in France: Designing network models fitting the transmission route of pathogens[END_REF][START_REF] Schulz | Network analysis of pig movements: Loyalty patterns and contact chains of different holding types in Denmark[END_REF]. Metapopulation models of epidemics were also developed to study the role of animal movements on infection spread based on either real movements or summaries derived from movement datasets [START_REF] Beaunée | Modelling of paratuberculosis spread between dairy cattle farms at a regional scale[END_REF][START_REF] Nickbakhsh | Implications of within-farm transmission for network dynamics: consequences for the spread of avian influenza[END_REF]. Focusing on transmission models in pigs, most of those were dedicated to regulated diseases such as foot-and-mouth disease, classical swine fever and more recently ASF [START_REF] Andraud | Modelling infectious viral diseases in swine populations: A state of the art [Review][END_REF]. The development of generic modelling frameworks (NAADSM, ISP, DTU-DADS, Be-Fast) was highly promoted to assess the impact of control strategies using unified frameworks (Bates et al., 2003b;[START_REF] Boklund | Comparing the epidemiological and economic effects of control strategies against classical swine fever in Denmark[END_REF][START_REF] Harvey | The North American Animal Disease Spread Model: a simulation model to assist decision making in evaluating animal disease incursions[END_REF][START_REF] Martínez-López | Evaluation of the risk of classical swine fever (CSF) spread from backyard pigs to other domestic pigs by using the spatial stochastic disease spread model Be-FAST: The example of Bulgaria [Article][END_REF].

These models account for direct between-farm (animal movement) and indirect transmission processes between epidemiological units such as farms, with anthropogenic, vectorial or environmental origins such as material sharing, feeding, vehicles used in the transportation of animals or feed, or airborne dissemination (Bates et al., 2003a;[START_REF] Boklund | Simulating the spread of classical swine fever virus between a hypothetical wild-boar population and domestic pig herds in Denmark[END_REF][START_REF] Halasa | Simulating the epidemiological and economic effects of an African swine fever epidemic in industrialized swine populations[END_REF]. However, the integration of movement data requires a pre-treatment of network data to summarize the contact rates between epidemiological units according to their intrinsic characteristics (e.g. herd types and sizes) and relational characteristics (e.g. distances and frequency of contacts according to herd types) [START_REF] Büttner | Epidemic Spreading in an Animal Trade Network -Comparison of Distance-Based and Network-Based Control Measures[END_REF][START_REF] Halasa | Simulating the epidemiological and economic effects of an African swine fever epidemic in industrialized swine populations [Article][END_REF][START_REF] Hu | Bayesian inference of epidemiological parameters from transmission experiments[END_REF].

The conjuncture of ASF in Europe demonstrates the permanent risk of introduction of the virus into disease-free areas and the need for predictive tools to assess both the consequences such an introduction would have and the responsiveness of the surveillance protocols. Based on social network structure, a previous study demonstrated the interest of a mathematical modelling approach for targeted surveillance towards at-risk herds among the swine commercial network in France, where the network was summarized based on movement frequency and distances between production sites (Andraud et al., 2019). The present study comes to complete this modelling framework through the integration of realized movement data, represented on a realistic time scale. An original metapopulation model was developed, fed with raw movement data, and parameterized according to ASF epidemiological specifications, to assess the extent to which ASF would spread among domestic pig holdings in France, as well as the effectiveness of surveillance and control strategies.

Material and Methods

Data

Data were extracted from the national swine movement database BDPORC [START_REF] Rautureau | Structural vulnerability of the French swine industry trade network to the spread of infectious diseases[END_REF][START_REF] Salines | Pig movements in France: Designing network models fitting the transmission route of pathogens[END_REF]. Briefly, notification of all shipments of live animals is compulsory in France since 2010, allowing for tracing of contacts between swine production units. The origin, the destination, the date, the number of animals that were exchanged and the truck identification number are collected at every shipment. Farm specific data are also available, allowing for geolocation and informing on the herd types and capacities (total number of animals). The data used hereafter cover the period from 2017 to 2019, over which two types of networks were designed, namely the animal introduction model (AIM) and the transit model (TM) (Figure 1) [22]. While the AIM approach considers effective contacts through the movements of live animals between production sites, the transit model accounts for complex rounds consisting of successive loadings and unloadings in different sites. The order of visited farms during the rounds was informed by the database, allowing for reproducing the exact contact tracing. TM links all herds that were visited by the same truck during a round. In the AIM, movements to slaughterhouses and rendering plants were neglected as they were considered as epidemiological endpoints [START_REF] Rautureau | Structural vulnerability of the French swine industry trade network to the spread of infectious diseases[END_REF]. However, if several herds were visited to collect slaughter-aged animals during the same round and using the same truck, they were fully considered in the TM due to the possibility of viral transmission from one herd to the other.

Epidemiological model

A multiscale model was designed to account for the within-herd transmission process, which governs in turn the transmission to commercial partners assuming three different transmission routes: local transmission affecting neighbouring herds (farm-to-farm proximity), transmission through the AIM, and indirect transmission related to TM.

Within-herd transmission dynamics

A stochastic Susceptible-Exposed-Infected-Removed (SEIR) model was used to represent the dynamics of infection in every infected herd. The initial number of susceptible animals was determined from the total pig capacity provided in the database, under the assumption of homogenous mixing. Running the simulations over a relatively short period of time (c.f. Simulations subsection described later), no birthdeath process was considered except for infection-related mortality. Therefore, for agents with a high lethality rate such as ASF virus, the within-herd population decreases as the infection progresses, thereby affecting the transmission process, due to the frequency dependent formulation of the force of infection. At each time step 𝑡 (day), the number of exposed animals in herd 𝑗, E j (𝑡) , is updated using a binomial process acting on the number of susceptibles 𝑆 𝑗 (𝑡) with probability 𝜋 𝑗

(𝑡) = 1 - 𝑒𝑥𝑝 (-𝛽 𝐼 𝑗 (𝑡) 𝑁 𝑗 (𝑡) )
, where 𝛽 denotes the transmission rate and 𝐼 𝑗 (𝑡) and 𝑁 𝑗 (𝑡) , the numbers of infectious and total number of animals in herd 𝑗 at time 𝑡. Transitions from exposed to infectious, and from infectious to removed stages are modelled in the same way, with an exponentially distributed infectious period of mean 𝛿 and 𝜎, respectively (transition rates: 1/𝛿, 1/𝜎). When lethal infections are considered, mortality was assumed to occur at the end of the infectious period, with a proportion 𝜇 of deaths among the number of newly removed animals. A herd was defined as infected when at least one pig was infected, initiating the within-herd transmission process.

Between-herd transmission dynamics

Three transmission routes were considered for the spread of the infectious agent between production sites.

Local transmission

The transmission in the neighborhood of infected herds may have different origins, e.g. airborne transmission, indirect contact through shared material, passage of vehicles, human factors, as well as fomites or vector transmission (e.g. hematophagous flies), which can be non-negligible routes of infection for ASF [START_REF] Bonnet | Putative Role of arthropod vectors in African swine fever virus transmission in relation to their bio-ecological properties [Review][END_REF][START_REF] Mauroy | Semi-quantitative risk assessment by expert elicitation of potential introduction routes of African swine fever from wild reservoir to domestic pig industry and subsequent spread during the Belgian outbreak (2018-2019)[END_REF][START_REF] Vergne | Mechanical transmission of African swine fever virus by Stomoxys calcitrans: Insights from a mechanistic model [Article[END_REF]. A distance-dependent transmission rate 𝛽 𝐿𝑜𝑐 (𝑑 𝑖,𝑗 ) was considered using a Gaussian kernel centered on every infected farms within a radius of 𝐷 = 2𝜎, 𝜎 being the standard error of the Gaussian kernel (expressed in meters).

The force of infection was represented as proportional to the number of infected individuals in the infected herds (𝐼 𝑖 ). Therefore, the animals housed in a susceptible farm 𝑗, located at a distance 𝑑 𝑖,𝑗 from an infected herd 𝑖, were exposed to a force of infection 𝜆 𝐿𝑜𝑐 (𝑗) (𝑡) = 𝛽 𝐿𝑜𝑐 (𝑑 𝑖,𝑗 ) * 𝐼 𝑖 (𝑡) , leading to a probability of infection at each time step equal to 𝜋 𝐿𝑜𝑐 (𝑗) (𝑡) = 1 -𝑒𝑥𝑝 (-𝜆 𝐿𝑜𝑐 (𝑗) (𝑡)). As described for within-herd transmission, the introduction process is governed by a binomial process at each time step, leading or not to the infection of some animals in the surrounding herds. Once a herd is infected via this route, meaning that at least one animal is infected and initiates the within-herd transmission process, local transmission is then assumed negligible for this herd.

Network driven transmission

In this study, raw data were directly integrated to parameterize the movements respecting their origin and destinations, as well as their actual timing, daily contact networks.

Animal introduction model (AIM)

A reaction-diffusion approach was considered to represent the spread of the agent between farms.

For each recorded shipment sourced from an infected farm, the transferred individuals are randomly drawn from the total population of the farm (herd size) according to the number of animals actually moved recorded in the database, and regardless of their infection status. They are then integrated into the target population (replacing former susceptible animals to keep the population size constant) and replaced by susceptible individuals in the source population. A susceptible farm was considered infected when at least one infectious individual was introduced, initiating the within-herd transmission process.

Transit model

In addition to the reaction-diffusion approach considered for animal movements, when a transport vehicle passes through an infected farm, a probability 𝜋 𝑇𝑀 of transmission of the agent to farms downstream in the round (dashed lines in Figure 1), in the absence of unloading of animals, was considered. This probability was assumed constant, independent from herd sizes or the numbers of farms visited during the round, using a Monte-Carlo algorithm to characterize the occurrence of infection events. When a farm was infected through this route, only one individual was shifted from susceptible to an exposed epidemiological state, its transition to the infectious stage initiating the within-farm transmission process.

Monitoring and control measures

Detection

The detection module, due to passive surveillance, was based on two indices reflecting the evolution of the infection at the herd level: a prevalence threshold, defined by the proportion of infected individuals in the farm, and a mortality threshold, defined by the cumulative number of dead (or removed) individuals in the farm (Table 1). In the absence of demographics, only pathogen-related mortality was considered. When a positive herd was detected, control measures were triggered and surveillance was enhanced by lowering the detection thresholds to notification of every dead animal and a prevalence of infected individuals of 1%.

Surveillance

In case of regulated diseases, the model evaluates the number of herds to be monitored based on the French emergency plan (DGAL/SDSPA/2019-41; DGAL/SDSPA/2019-41). When a case is detected, surveillance and protection zones (SZ and PZ) are set up for 40 days within a 10-and 3-km radius, respectively. All herds in these two zones are monitored during that period. Mortality notification is compulsory and systematically induces clinical inspection leading to further laboratory analyses (virology and serology). Downward and upward tracings are also performed to identify the contacts due to animal movements that occurred up to 30 days prior to detection.

Control measures

Implemented control measures were based on movement restrictions. Detected herds were depopulated and considered as not infectious anymore by any transmission route. A complete containment was implemented for herds identified in surveillance and protection zones as well as those traced as potential contacts of the infected herd, similarly to what would be implemented in the event of an outbreak.

Simulations General Considerations

The parameter set used for simulation purposes is provided in table 2. For a matter of illustration, these parameters were selected in ranges matching to ASF epidemiological knowledge and were mainly derived from the literature [START_REF] Guinat | Experimental pigto-pig transmission dynamics for African swine fever virus, Georgia 2007/1 strain [Article[END_REF][START_REF] Halasa | Simulating the epidemiological and economic effects of an African swine fever epidemic in industrialized swine populations [Article][END_REF]. Based on the pyramidal structure of the swine industry in France [START_REF] Rautureau | Structural vulnerability of the French swine industry trade network to the spread of infectious diseases[END_REF][START_REF] Salines | Pig movements in France: Designing network models fitting the transmission route of pathogens[END_REF], the herds were distributed in three categories corresponding to nucleus and multiplier herds (further called "Breeders"), farrowing-only and post-weaning production farms ("Farrowers"), and farms having a finishing sector ("Finishers"). For each of these farm categories, 1,000 simulations were run, each with a randomly selected seeding herd within the given category. Each simulation could last up to 40 days, but was stopped as soon as there were no more infected herds in the system. To minimize the risk of stochastic extinction in the first days of simulations, the number of infectious animals in the index herd was set to five. Thus, simulations were initialized at a random date (in the 3 years of data) through the introduction of 5 infected individuals in a random index herds, selected among those having at least one animal movement in the 40 days of simulation.

Scenarios

The baseline scenario consisted of no-control spread. Although unrealistic for a regulated disease such as ASF, this scenario demonstrated the extent to which an infectious agent could spread over the network when considering realized shipment data.

Detection and movement restriction

The initial detection thresholds were set to 5% of infected individuals at the herd level, along with the death of 5 individuals (Table 2). Based on that scenario, the prevalence threshold was then decreased to 1% in scenarios 3, followed by a decrease of the mortality threshold down to 1 animal in scenario 4.

After the first detection, the threshold were set to scenario 4 values to represent the onset of enhanced surveillance protocol.

Outcomes

For each scenario and seeded introduction at each herd category, the number of simulations where secondary infections occurred was recorded. For those simulations, the number of secondary infected herds -along with the transmission route and the distance to the source of infection -were registered. Where applicable, the delay between infection and detection, otherwise stated as the 'at risk' period, for each infected and detected herd was also analyzed.

The model was developed in R software, using the igraph package to represent the time-respective networks, rgeos, rgdal and sf to deal with spatial features and shiny for user-friendly implementation interface. The model is available on GitHub, along with a data sample for testing (https://github.com/mandraud/EpiNet/tree/main).

Results

Network descriptions

Animal introduction model

The network accounts for 12,397 active herds located in France and involved in 268,672 shipments over the three years of data (2017-2019). Considering shipments with loadings and unloadings of animals between production sites (figure 3A and 3B), the average number of animals per shipment was 150 pigs (2.5 -97.5 percentiles: 1 -672) and were dependent on the category of the source herd.

Indeed, shipments from nucleus farms mostly involved small groups of animals, with a median size of 34 animals (range [1 -280]), while those from farrowing and post-weaning herds concerned full batches with a median group size about 208 piglets ([8 -830], Figure 3B). However, the highest out degree was reached by nucleus herds with a median number of 195 outgoing contacts ([2 -908], Figure 3A) during the three years. Only 32% of herds with finishing barns were involved in outgoing shipment with unloadings of post-weaned animals in other finishing herds, and they had very limited outgoing contact rate (median outdegree 12 [1 -124]). Accounting for truck transit in rounds dramatically densifies the contact structure with up to 916,271 links in the network involving 12,795 production sites in 238,159 rounds. This increase is obviously due to extensive rounds visiting high numbers of farms. The role of herds with finishing sectors was found essential, covering 75% of sites involved in the TM network due to sequential loadings from different farms to slaughterhouses (Figure 3C and3D).

Based on network characteristics, three groups of herd types with similar characteristics were designed for further analysis: Farrowers corresponding to herds with farrowing-only and/or post-weaning sectors (represented in yellow on boxplots in figure 3), Finishers, to herd having a finishing sectors (orange on boxplots in figure 3), and Breeders, to multiplier and nucleus herds (purple on boxplots in figure 3). The color scheme remains consistent throughout the paper.

Simulations Scenario 1: "No-control" Baseline situation

The theoretical results in the absence of control measures are described in Figure 4. On average, about 70% of the simulations resulted in secondary cases. Transmission was quasi-systematic when the introduction was set in breeders and farrowers (79 and 82% of simulations; Figure 4A) compared to finishers (49%). Introductions in Breeder herds led to major outbreaks with a median of 12 secondary cases, but potentially rising up to 52 cases (Figure 4B). In contrast, introductions in other production herds (farrowers or finishers) seeded the infectious agent with a lower intensity, affecting a median number of 3 herds (95% interval: [3 -12]); the distance ranges were also two-times lower reaching 200 km (Figure 4C). Direct contact represented the major route of transmission of the virus, between 68% and 85% of transmission events. This was followed by transmission through indirect contact in the neighbourhood which, despite a fairly small radius of interaction (500 m), represented 14 to 29% of events depending on the category of farm of introduction. Finally, transit movements (TM) were implicated as being responsible for 2% of the cases recorded during the simulations (Figure 4D).

Transmission through the movement of animals or vehicle transit occurred mainly on a relatively small geographical scale, with median distances ranging between 10 and 30 km. Some movements nevertheless generated more distant infections, which could reach 400 km when the virus was introduced in breeders.

Control measures

Scenario 2. Detection criteria: 5% Prevalence threshold and 5 dead animals

The detection process, based on a prevalence of 5% of infected animals and a disease-related mortality of 5 animals per farm, induced a median time to detection of 14 to 21 days post-infection depending on the category of infected unit (Figure 5A-E). All detection events were effective within 30 days following infection, supporting the duration of the contact-tracing period. Farrowers were detected earlier than breeders and finishers, due to lower herd sizes leading to earlier fulfillment of detection criteria. Finishers were the most affected farm category with up to 4743 detection events over the whole set of simulations (3000 iterations, Figure 5E). During the 'at-risk' period, defined as the delay between infection and detection; the occurrence of secondary infections was reduced by 17% when breeders were selected as the initial seeders (675 simulations) with transmission chains limited to 5 herds ([2 -20], Figures 5A and5B). Detection and subsequent control measures reduced the number of secondary infections by 28% and 32% when seeding infection in finishers and farrowers, respectively, but did not affect the length of transmission chains when transmission occurred, with a median of 2 to 3 infected farms per run. The overall impact of the control measures induced a reduction of the total number of infections due to animal transfer, making the relative contribution of local transmission route more important, representing 26% to 27% of transmission events when the virus was introduced in breeders and farrowers, and becoming the major transmission route when introduction occurred in finishers (50.6%, Figure 5D). The median distance between new cases and their source of infection was evaluated to 66.5 km for direct transmission when seeding in breeders (43 and 32 km for farrowers and finishers, respectively, Figure 5C).

Scenario 3. Detection criteria: 1% Prevalence threshold and 5 dead animals

Decreasing the prevalence of infected pigs' threshold for detection to 1% led to a faster response to the transmission process (Figure 6A-E). The median time to detection was evaluated from 11 to 14 days for all farm types. The earliest detections were still observed when introduced in farrowers, starting on day 6 post-infection, but differences between herd categories were reduced (Figure 6E).

Although the detection time was significantly decreased when primary infection was seeded in farrowers and breeders, the impact on the number of simulations with effective transmission when starting from finishers was not significant. In contrast, breeders were detected much earlier when lowering the prevalence detection threshold (21 days in scenario 2, 13 days in scenario 3), confirming the importance of the out-degrees in regards to the transmission process. The proportion of simulations with effective transmission to secondary cases fell to 21% when the agent was introduced in finishers, and was limited to 43% in the alternative introduction schemes (Figure 6A). A median of one secondary outbreak was observed wherever the introduction occurred, but breeders could spread the infection up to 13 herds due to an extended "at-risk" period and important outgoing contact chains (Figure 6B). A movement ban after the first detection drastically reduced the number of infected herds, which were mostly due to local transmission thereafter (Figure 6D).

Scenario 4. Detection criteria: 1% Prevalence threshold and 1 dead animal

Decreasing the mortality threshold to 1 animal improved the detection in simulations when the infection was seeded in farrowers (median time to detection estimated to 7 days, Figure 5A-E). It remained unchanged in breeders and finishers. This was directly reflected by the proportion of simulations leading to secondary infections which fell to 33% when infection was introduced in farrowers (Figure 7A). None of the evaluated control measures were found to inflect the distance pattern for transmission processes (Figures 5-7C), despite the higher proportion of local transmission events.

Discussion

African swine fever represents a major threat for the swine industry and related economic sectors.

Preparedness to face such a threat is therefore a need in countries that are still free of the disease.

The development of predictive tools, based on risk assessment or mechanistic models, offers scientific support to decision makers. Several mechanistic models were developed to tackle potential introduction and dissemination of ASF given the increased risk of global spread -as highlighted with its expansion throughout Europe and Asia and now into the Caribbean [START_REF] Gonzales | African swine fever in the Dominican Republic[END_REF]. A recent review identified four main objectives for these studies [START_REF] Hayes | Mechanistic modelling of African swine fever: A systematic review[END_REF]. The present study comes at the crossroad of three of those through the evaluation of the contribution of the transmission routes, the assessment of the epidemiological consequences of a hypothetical outbreak and the assessment of effectiveness of control measures. Based on farms' characteristics and movement of animals, the results highlighted the relationship between the transmission potential of the herds and their position in the swine pyramidal production system. A quantitative evaluation of the consequences of pathogen introduction in swine trade network was conducted. Although live pig movements were identified as the main transmission route, the results evidenced the need of compliance with external biosecurity, to limit the transmission in the neighborhood of infected farms as well as disinfection of vehicles.

Here we developed a metapopulation model to represent the on-farm spread of the pathogen combined with a transmission process between production sites based on movement data. Such an approach was already described in several studies and revealed the huge potential of multiscale modelling [START_REF] Andraud | Modelling infectious viral diseases in swine populations: A state of the art [Review][END_REF][START_REF] Bastard | Dynamics of livestock-associated methicillin resistant Staphylococcus aureus in pig movement networks: Insight from mathematical modeling and French data[END_REF][START_REF] Beaunée | Modelling of paratuberculosis spread between dairy cattle farms at a regional scale[END_REF]. The work of Beaunée et al. on Mycobacterium avium subsp. paratuberculosis (Map) spread in cattle clearly illustrated the advantage of integrating information at different scales to have an exhaustive view of the epidemiological context [START_REF] Beaunée | Modelling of paratuberculosis spread between dairy cattle farms at a regional scale[END_REF]. Including data at the farm level along with records of animal movements allowed analyzing the conditions for introduction and persistence of Map on farms and at a regional scale. Dedicated to ASF in the domestic pig compartment, Halasa et al. proposed a multiscale approach, combining a within-herd transmission model with a reconstruction of the contact network based on the frequency of contacts registered between herd types and distances between production sites [START_REF] Halasa | Simulating the epidemiological and economic effects of an African swine fever epidemic in industrialized swine populations [Article][END_REF]. More recently, Machado et al. designed a model representing ASF spread in Brazil, taking 3-by-3km cells of a geographical raster as epidemiological units and representing movements between these units (Gustavo Machado et al., 2021).

Three transmission routes were considered for between-sites transmission: the animal movement pathway remains dominant but as soon as movement restriction occurs through the control measures, the relative contribution of local transmission increases. This is under the assumption that there is no modification of local biosecurity measures in case of an outbreak, which may be somewhat unrealistic.

Hence, an evolution of local biosecurity, through awareness of all local actors, appears as a key element for disease control. Although playing a minor role in the transmission process, truck transit should not be neglected and complex rounds serving several production sites have to be considered in the tracking surveillance system.

The delay from infection to detection, or the 'at-risk' period, was analyzed depending on the sensitivity of detection in infected farms. Using different prevalence and mortality thresholds, the time to detection varied between 8 and 21 days, in line with results obtained from Halasa's model [START_REF] Andraud | Threat to the French Swine Industry of African Swine Fever: Surveillance, Spread and Control Perspectives[END_REF][START_REF] Halasa | Simulating the epidemiological and economic effects of an African swine fever epidemic in industrialized swine populations [Article][END_REF]. During this period infected herds may transmit the virus to epidemiologically related herds: the shorter this at-risk period, the lower the number of infected herds. The detection delay was also highly dependent on the herd types (breeding, farrowing or finishing herds). Due to their position in the pyramidal structure of the swine production system, the probability of first detecting finishing herds is rather low and the modification of detectiontriggering parameters did not affect their detection delay (15 to 17 days). This suggests that such a threshold should be adapted according to both the herd characteristics (e.g. the type of herds defining their contact structure) and their neighbouring context, based on the local density of farms as well as their connectivity to other sites. Although the awareness of the local actors (farmers, veterinarians) is essential to limit disease spread, high levels of vigilance -in regards to social and societal aspectswould be difficult to maintain for an extremely long period of time, requiring adjustments according to the level of risk.

Highly connected herds were identified as potential key spreaders in this study, especially nucleus and multiplier herds, which have the longest outgoing contact chains. Conversely, finishing herds -with a high number of ingoing connections -appeared more vulnerable but were not at a high risk of spread because of their low number of outgoing connections. It is worth noting all herds were assumed to have similar internal and external biosecurity levels. This shortcut may create an overestimation of the transmission process, especially originating from nucleus and multiplier herds which are likely to observe strict biosecurity protocols limiting the risk of introduction of pathogens [START_REF] Silva | Assessment of biosecurity practices and development of a scoring system in swine farms using item response theory[END_REF].

The incorporation of data-based classification of herds' biosecurity levels would improve the model, providing for better guidance for surveillance targets. In the absence of such biosecurity qualifications, the choice was made to set a baseline level for all herds representing the most severe epidemiological situation.

In contrast with models based on summary statistics of movement data (Andraud et al., 2019;[START_REF] Halasa | Simulating the epidemiological and economic effects of an African swine fever epidemic in industrialized swine populations [Article][END_REF], raw movement data were directly used to build a dynamic network on a realistic time-scale. This approach allows for a rapid update of input data without any pre-treatment, which could be important in terms of reactivity, should an introduction occur. However, the results remain based on historical data, which may raise the question of prediction quality in the event of an actual introduction. One solution to this issue may rely on the generation of simulated networks, with similar characteristics. In that view, Exponential Random Graph Models (ERGMs) offer a promising perspective [START_REF] Kukielka | Modeling the live-pig trade network in Georgia: Implications for disease prevention and control [Article][END_REF][START_REF] Phoo-Ngurn | Modeling the spread of porcine reproductive and respiratory syndrome virus (PRRSV) in a swine population: transmission dynamics, immunity information, and optimal control strategies[END_REF][START_REF] Relun | Spatial and functional organization of pig trade in different european production systems: Implications for disease prevention and control [Article[END_REF]. This statistical tool allows for characterizing the main factors governing the network structure based on its structural features and the node characteristics, from which simulated networks can be generated. This way, the model would get free from the past, improving outcome predictiveness.

Our model focuses on transmission dynamics of an infectious agent on swine production networks after the introduction of ASF in the domestic pig compartment. The introduction in France is nevertheless very likely to occur in wild boars, as was the case in several countries [START_REF] Chenais | Epidemiological considerations on African swine fever in Europe 2014-2018[END_REF][START_REF] Podgórski | Spatial epidemiology of African swine fever: Host, landscape and anthropogenic drivers of disease occurrence in wild boar [Article][END_REF][START_REF] Sauter-Louis | Joining the club: First detection of African swine fever in wild boar in Germany[END_REF]. To analyze the risk of introduction into the domestic pig reservoir, the interface between wild-boar and domestic pig compartments needs to be represented [START_REF] Hayes | Mechanistic modelling of African swine fever: A systematic review[END_REF]. This could induce an external force of infection, especially on free-range herds which were already highlighted as critical point for surveillance [START_REF] Andraud | Threat to the French Swine Industry of African Swine Fever: Surveillance, Spread and Control Perspectives[END_REF].

In this study, we focused exclusively on control measures based on movement restriction and immediate depopulation of infected herds after detection. In the meantime, movement ban was implemented in surveillance and protection zone. In that view, the modelling framework could be used to evaluate alternative surveillance and control measures (e.g. sizes and duration of zones) to optimize their epidemiological and economic impact (G. Machado et al., 2021). The identification of 'at-risk' herds from farms and network characteristics could help designing preventive strategies, reducing so the infectious pressure in the neighborhood of infected farms, especially on high density areas [START_REF] Jurado | Relevant measures to prevent the spread of African Swine Fever in the European Union Domestic Pig Sector [Review[END_REF].

In conclusion, a transmission model of infectious agents on the swine commercial network was developed to analyze the effectiveness of surveillance and control measures. Fed with raw movement data, the modelling framework is flexible with the transmission process represented on a timerespective schedule. Applied to ASF, this model identified animal movements as the pivotal transmission route. However, local transmission and mechanical transmission through truck transit should not be neglected, and requires a global awareness of all actors from the production system.

Figure 1 .

 1 Figure 1. Illustration of the Animal introduction Model (AIM; left) and the Transit Model (TM; right)

Figure 2 .

 2 Figure 2. Multi-scale epidemiological units and transmission processes. Within-farm a stochastic

Figure 3 .

 3 Figure 3. Network characteristics for animal movements (AIM, panels A and B) and transit (TM, panels C and D) from 2017 to 2019 for the different herd types in France (FA: Farrowing, PW: Post-Weaning, FPW: Farrowing-Post-Weaning, PWF: Post-Weaning-Finishing, FI: Finishing, FF: Farrow-to-Finish, MU: Multipliers, NU: Nucleus). A. Number of outgoing contacts (loadings); B. Number of animals exchanged per movement; C. Number of transit movement from each farm type; D. Number of transit movement to each farm type.

Figure 4 .

 4 Figure 4. Simulation results in the absence of control measures according to the category of index farm by the end of 40 days of disease propagation: A. number of simulations leading to secondary transmissions (among 1000 simulations); B. number of infected farms; C. distances between source farms and secondary cases; D. relative contribution of different transmission routes (direct for pig movements, local and transit related transmissions)

Figure 5 .

 5 Figure 5. Analysis of detection probability according to time since infection depending on the farm category (scenario with prevalence and mortality detection thresholds of 5% and 5 animals per farm, respectively). A. Number of simulations leading to secondary transmissions (among 1000 simulations); B. Number of infected farms; C. Distances between source farms and secondary cases; D. Relative contribution of different transmission routes (direct for pig movements, local and transit related transmissions); E. Survival analysis of the delay of detection.

Figure 6 .

 6 Figure 6. Simulation results under control scenario 3 (for prevalence and mortality detection thresholds of 1% and 5 animals) according to the type of index farm. A. number of simulations leading to secondary transmissions (among 1000 simulations); B. number of infected farms; C. distances between source farms and secondary cases; D. relative contribution of different transmission routes (direct for pig movements, local and transit related transmissions); E. survival analysis of the delay of detection (3000 simulations).

Figure 7 .

 7 Figure 7. Simulation results under control scenario 3 (for prevalence and mortality detection thresholds of 1% and 1 animal) according to the type of index farm. A. number of simulations leading to secondary transmissions (among 1000 simulations); B. number of infected farms; C. distances between source farms and secondary cases; D. relative contribution of different transmission routes (direct for pig movements, local and transit related transmissions); E. survival analysis of the delay of detection (3000 simulations).

Table 1 .

 1 Parameters of the epidemiological model.

	Parameter	Definition (unit)	Value (Reference)
	𝛽	Between-individuals Transmission rate (day -1 )	0.6 (Guinat et al., 2016)
	𝛽 𝑙𝑜𝑐	Maximum Local Transmission rate (day -1 )	0.0001 (Halasa, Botner,
			et al., 2016)
	𝐷	Radius for local transmission kernel (m)	500
	𝛽 𝑇𝑀	Transit transmission probability	0.05 (Expert opinion)
	𝛿 -1	Latent period (day)	4 (Halasa, Botner, et al.,
			2016)
	𝛾	Infectious period (day)	7 (Halasa, Botner, et al.,
			2016)

Table 2 .

 2 Parameters for detection of the first herd initiating control measures and enhanced surveillance. After the first detected case, the parameters are set to those of scenario 4.

			Detection parameters
			Prevalence of infected	Mortality Threshold (number
	Scenario	detection	animals Threshold	of dead pigs)
	1	FALSE	-	0
	2	TRUE	0.05	5
	3	TRUE	0.01	5
	4	TRUE	0.01	1

Acknowledgements

We would like to thank the members of BDPORC for making holdings and animal movement data available.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à

Data availability statement

Code and example datasets are made available on GitHub: https://github.com/mandraud/EpiNet/tree/main.

Conflict of interest statement

The authors declare no conflict of interest.


Shaohua Yu
email: shaohua.yu@njust.edu.cn

Ann Melissa Campbell
email: ann-campbell@uiowa.edu

Jan Fabian Ehmke
email: jan.ehmke@univie.ac.at

Jakob Puchinger
email: jpuchinger@em-normandie.fr

Last-mile Robot Deliveries with Access Restrictions

Keywords:

Introduction

In recent years, the improvement of quality of life in cities through the reduction of local pollutant emissions and traffic has been an increasingly important topic for city planners and researchers. This leads to classical practices such as using vans for deliveries to be restricted and other modes preferred in some areas such as campuses, pedestrian zones or inner-city districts. Recent technological developments and experiments point to the potential of automated lastmile deliveries using small battery-powered robots. We therefore investigate a two-echelon hybrid van-based robot delivery system for last-mile logistics operations.

In existing research, using robots for second-echelon delivery has been considered. For example, a two-echelon delivery model (2EVRP) is proposed in [START_REF] Bakach | A two-tier urban delivery network with robot-based deliveries[END_REF]. At the first echelon, vans transport packages from a large depot to small local robot hubs. At the second echelon, robots deliver items to customers. A van-based robot delivery model (VRD) concept with time windows is studied in [START_REF] Yu | Two-echelon urban deliveries using autonomous vehicles[END_REF]. Here, robots are transported in the vans to perform second-level route deliveries. According to a sensitivity analysis performed in [START_REF] Yu | Van-based robot hybrid pickup and delivery routing problem[END_REF], each of these delivery models have their advantages. Depending on fixed costs of vans, robots or satellites, as well as other factors, one of the models may perform better than the other.

We propose a hybrid delivery model to select the optimal distribution strategy considering conventional vans, inner city hubs with robots, and delivery vans carrying robots, potentially improving distribution efficiency. We model the problem as mixed-integer program and propose an adaptive large neighborhood search (ALNS) to solve large instances efficiently. We conduct computational experiments to evaluate how our proposed model performs and how different characteristics (size, time) of the access restrictions impact its performance. We also evaluate how different characteristics impact solution structure.

Problem Description

We propose to combine the models in the literature to build a hybrid van-based robot delivery model (HVRD). This model lets the distribution system optimize the distribution strategy (using hubs, parking nodes or a combination of the two) to improve distribution efficiency.

A central depot, vans, vans carrying robots, robot hubs, and parking nodes are all given. The vans or vans carrying robots move along the first-level routes, serve customers, or visit parking nodes to drop off/pick up their robots, or visit hubs to unload goods. Robots start from parking nodes or hubs to handle deliveries along second-level routes. Two kinds of customers are considered : van customers can be visited by either the van or the robot, whereas robot customers can only be visited by a robot. They correspond to the customers in access-restricted van no-go areas.

Customers must be visited exactly once and have a service time and a time window. We assume parking nodes are used for the van to drop off/pick up and replenish its onboard robots, and hubs are used to store robots and goods. Each hub can only be visited once and contains a maximum number of robots ; it can also act as a parking node. Parking nodes can be visited multiple times. We also assume the depot is equipped with a sufficient number of vans.

Robots have a maximum battery capacity, which determines the distance they can travel, but there is no set maximum travel distance restriction for vans. There is a maximum capacity for vans and robots. We assume that when the robot is on board the van, the van and its robot's total load cannot exceed the van's maximum capacity. The robots can visit multiple customers during their trips. The onboard robots cannot leave the depot to serve customers directly and cannot go back to the depot independently.

For recharging, we assume every time a van picks up its robots, the battery is swapped. We assume an onboard robot can be replenished with goods from its corresponding van at parking nodes or hubs, and a hub robot can be replenished with goods at hubs. Note that hubs can be replenished with goods from vans, and vans' freight must be loaded at the depot.

The objective is to minimize the total cost, consisting of travel costs, equipment usage costs, and labor costs.

Results and Conclusions

We solve the HVRD model using CPLEX and an ANLS approach. We compare our HVRD model with a two-echelon vehicle routing model (2EVRP) and a van-based robot delivery model (VRD) to evaluate how our proposed hybrid model performs. Our research shows that the number of solvable cases of the HVRD model is larger than that of VRD and 2EVRP models with our instances. Even for the instances where all the customers are feasible for three models, the outputs of the HVRD model outperform those of the 2EVRP and VRD model. Therefore, allowing the system to select the appropriate distribution mode automatically can help the enterprise's distribution accessibility and efficiency.

For sensitivity experiments, we explore the impact of van no-go areas and van no-go area time windows. Our experimental results show that as the van no-go area shrink, the cost of the distribution system is significantly reduced. Besides, as the length of the access time window increases, their impact gradually decreases.

With regards to the increasing importance of an on-demand economy and same-day deliveries dynamic versions of the HVRD problems might also be of interest for future research.