Lando Bosma
email: lbosma@umcutrecht.nl

Cornel Zachiu

Mario Ries

Baudouin Denis De Senneville

Bas Raaymakers

Quantitative investigation of dose accumulation errors from intra-fraction motion in MRgRT for prostate cancer

Keywords: Motion-compensated dose accumulation, intra-fraction plan adaption, MRguided radiotherapy, dose reconstruction

Accurate spatial dose delivery in radiotherapy is frequently complicated due to changes in the patient's internal anatomy during and in-between therapy segments. The recent introduction of hybrid MRI radiotherapy systems allows unequaled soft-tissue visualization during radiation delivery and can be used for dose reconstruction to quantify the impact of motion. To this end, knowledge of anatomical deformations obtained from continuous monitoring during treatment has to be combined with information on the spatio-temporal dose delivery to perform motion-compensated dose accumulation (MCDA).

Here, the influence of the choice of deformable image registration algorithm, dose warping strategy, and MR image resolution and SNR on the resulting MCDA is investigated. For a quantitative investigation, four 4D MRI-datasets representing typical patient observed motion patterns are generated using finite element modeling and serve as a gold standard. Energy delivery is simulated intra-fractionally in the deformed image space and, subsequently, MCDA-processed. Finally, the results are substantiated by comparing MCDA strategies on clinically acquired patient data.

It is shown that MCDA is needed for correct quantitative dose reconstruction. For prostate treatments, using the energy per mass transfer dose warping strategy has the largest influence on decreasing dose estimation errors.

Introduction

Radiotherapy for cancer treatment aims to deliver a high dose to the tumor while keeping the dose to the surrounding healthy tissue as low as possible. Biological motion (e.g. respiration, bladder and rectal filling, peristalsis) during treatment might result in geometric uncertainties and, thus, dose delivery errors to the tumor. Margins around the tumor are typically used to optimize dose coverage [START_REF] De Luca | Prescribing, recording, and reporting intensity-modulated photon-beam therapy (imrt)(icru report 83)[END_REF]. However, such margins generally lead to over-exposing healthy tissue or organs-at-risk (OAR) situated within the margins and are often dose-limiting.

Image-guided radiotherapy offers the possibility to reduce these uncertainties, allowing a reduction of these margins and alleviating the dose limitations, potentially improving treatment outcomes [START_REF] Jaffray | Image-guided radiotherapy: from current concept to future perspectives[END_REF][START_REF] Mundt | Image-guided radiation therapy: a clinical perspective[END_REF]. For example, on-board imaging such as cone-beam computed tomography (CBCT) can be used to compare the daily treatment anatomy to the anatomy at the time of planning, allowing the correction of translations and rotations of the anatomy on an inter-fraction basis [START_REF] Guckenberger | Image-guided radiotherapy based on kilovoltage cone-beam computed tomography-a review of technology and clinical outcome[END_REF][START_REF] Thilmann | Correction of patient positioning errors based on in-line cone beam cts: clinical implementation and first experiences[END_REF]). CBCT-guidance can thus be used to decrease geometric uncertainties to a particular extent; however, the poor soft tissue contrast renders the estimation of more complex deformations of individual anatomical areas challenging [START_REF] Mcbain | X-ray volumetric imaging in image-guided radiotherapy: the new standard in on-treatment imaging[END_REF]. Additionally, CBCT-guidance also leads to the delivery of imaging-related radiation, which is undesirable. Recently, magnetic resonance (MR) imaging has also been integrated with radiotherapy delivery systems [START_REF] Lagendijk | The magnetic resonance imaging-linac system[END_REF][START_REF] Mutic | The viewray system: magnetic resonance-guided and controlled radiotherapy[END_REF], allowing MR guided radiotherapy (MRgRT). MRgRT offers high soft-tissue contrast suitable for tumor tracking at high frame-rates without additional radiation exposure to the patient. This enables direct tumor visualization and OAR localization from the treatment table, not only prior to the daily treatment but also during radiation delivery. An example of dose reconstruction for prostate treatments at the 1.5T MRI radiotherapy system is presented by [START_REF] Kontaxis | Delivered dose quantification in prostate radiotherapy using online 3d cine imaging and treatment log files on a combined 1.5 t magnetic resonance imaging and linear accelerator system[END_REF], where the information on the anatomical deformations is coupled with knowledge of the delivered dose under the assumption of correct image registration of the different time points. Thus, having knowledge of the delivered dose can aid more precise reporting for treatment response assessment. Additionally, the time-resolved dose reconstruction can also be used as input to adapt the treatment, i.e. for intrafraction plan adaption.

In particular, for radiotherapy for prostate cancer, several therapy innovations have been suggested that rely on an accurate and precise motion-compensated dose accumulation (MCDA), or motion-compensated dose reconstruction, as an input. Hypofractionation has recently gained traction as it combines positive therapy results with cost benefits and decreased patient burden [START_REF] Den Hartogh | Planning feasibility of extremely hypofractionated prostate radiotherapy on a 1.5 t magnetic resonance imaging guided linear accelerator[END_REF][START_REF] Widmark | Ultra-hypofractionated versus conventionally fractionated radiotherapy for prostate cancer: 5-year outcomes of the hypo-rt-pc randomised, non-inferiority, phase 3 trial[END_REF]. However, when increasing the dose-per-fraction while decreasing the number of fractions, there is less to none averaging out of the dose error over fractions, making precise dose delivery and repeated adaptation even more of an imperative. Furthermore, since fractions in a hypofractionated plan are generally of longer duration, there is a higher chance for (more extensive) biological motion to occur (de Muinck Keizer, Pathmanathan, Andreychenko, Kerkmeijer, van Zyp, Tree, van den Berg & de Boer 2019, de Muinck Keizer, Kerkmeijer, Maspero, Andreychenko, van Zyp, Van den Berg, Raaymakers, [START_REF] De Muinck Keizer | Soft-tissue prostate intrafraction motion tracking in 3d cine-mr for mr-guided radiotherapy[END_REF]. Note that the possibility of large prostate motion is somewhat dependent on the treatment method used as for modern non-MR-hybrid linacs, the duration of a standard fraction can be as short as 2 minutes, decreasing the possibility of large intra-fraction prostate motion.

Moreover, due to the relatively slow and episodic component of the changes in the anatomy due to bladder filling and peristalsis, intra-fraction dose reconstruction and plan adaptation are conceptually feasible for prostate MRgRT, as the temporal MRresolution is about one image per 10 to 20 seconds. For example, if the reconstructed dose deviates from the planned dose distribution, the treatment might be interrupted and continued with an adapted treatment plan, which accounts for the dose delivered so far. This ensures convergence of the delivered dose to the planned dose, irrespective of motion [START_REF] Kontaxis | Towards adaptive imrt sequencing for the mr-linac[END_REF], Kontaxis et al. 2015a[START_REF] Kontaxis | Towards fast online intrafraction replanning for free-breathing stereotactic body radiation therapy with the mr-linac[END_REF], thus increasing the quality of the treatment.

Since MCDA can potentially serve as an input for such motion correction strategies, its reliability, precision, and accuracy become paramount, as any errors may end-up accumulating over time. Moreover, such errors could even amplify due to feed-back into the (re-)planning algorithm, up to a point where the (erroneous) corrections become detrimental for therapy. The quality of MCDA (i.e. accuracy and precision) depends on multiple factors such as (a) the deformable image registration (DIR) algorithm, (b) the image resolution, (c) the signal-to-noise-ratio (SNR), and (d) the employed dose warping strategy.

Here we investigated the impact of these factors on the accumulated dose of MCDA for intra-fraction MRgRT for prostate cancer. To this end, a quantitative comparison of the different MCDA strategies is performed on finite-element simulated motion patterns representing 'typical' motion patterns observed in patients. A quantitative comparison is then possible as the ground truth deformations are known, and the different MCDA strategies can be compared to a gold standard. Subsequently, all MCDA strategies are qualitatively evaluated on a single patient clinical dataset. The aim of this study is to give an overview of the impact of the most relevant methodological choices on the quality of the resulting MCDA and to provide guidance for making well-considered choices for clinical implementations.

Materials and methods

In the following section, the experimental setup used for analyzing the impact of the choice of registration algorithm, image resolution and SNR, and dose warping method on the dose accumulation process will be described. Four finite element models (FEM) of different typical prostate motion patterns are simulated. From these, a set of cine-MR images depicting the motion trace is generated for a quantitative assessment of the DIR algorithms. A single radiotherapy fraction is simulated, using a typical clinical treatment plan and incorporating the simulated motion patterns.

The data

As a basis for the FEM-generated images used here for the MRgRT simulation, we have used a clinical 3D MR scan acquired on a prostate cancer patient treated on the 1.5T MR-Linac Unity system (Elekta Unity AB, Stockholm, Sweden) installed at the UMC Utrecht, The Netherlands. The current study was conducted in agreement with the required standards and regulatory requirements. Ethical approval was provided by the Ethics Board of the University Medical Center Utrecht. The original spatial dimensions of the image are 448 × 448 × 63, with a reconstructed voxel spacing of 0.93 × 0.93 × 2 mm 3 . The image was acquired using 3D bTFE with fat suppression (relaxation time T R = 4.7 ms, echo time T E = 2.3 ms, flip angle = 50 • , field strength B 0 = 1.5 T). The delineations of the bladder, prostate, and rectum performed on the image by an experienced radiation oncologist were then imported into the FEM software FEBio [START_REF] Maas | FEBio: finite elements for biomechanics[END_REF]. The prostate was modeled as an isotropic elastic material with a density ρ = 1000 g/L, Young's modulus E = 21 kPa [START_REF] Barr | Shear wave ultrasound elastography of the prostate: initial results[END_REF]) and Poisson's ratio ν = 0.4 [START_REF] Bharatha | Evaluation of three-dimensional finite element-based deformable registration of pre-and intraoperative prostate imaging[END_REF]. The rectum and bladder were modeled as biphasic materials, allowing them to expand. Volumetric changes and movements of the bladder and rectum were subsequently induced and used as actuators within the FEM simulation, which in turn led to movements and deformations of the surrounding anatomy.

We modeled four different anatomical deformations, mimicking typical patient observed motion patterns, which we have outlined below (from largest to smallest average prostate motion amplitude).

• Rectal filling motion. The rectum expands and moves in the superior direction, modeling the passage of a gas-bubble. Under this effect, the prostate moves in the anterior-superior direction, following an individual patient motion pattern observed and reported within a previous study ( • Residual motion only. The most common anatomical motion in the abdomen during treatment is the absence of large scale deformations or displacements, with only residual motion due to the digestive tract's peristalsis. In this simulation, the bladder expands only slightly, pushing the prostate marginally in the inferior direction. This dataset is added mostly as a control in order to investigate the effect of MCDA on dose accumulation errors when minimal motion is present.

The resulting average displacements of the prostate for all motions are shown in figure 1 and the resulting anatomical deformations for the sixth time point of the simulation for the rectal filling motion are shown in figure 2. The simulated motion was extracted in ten equally distanced time-steps, resulting in ten 3D deformation vector fields (DVF), which relate each voxel's position in 3D space to their original position. This set of vector fields serves as the gold standard for all DIR algorithms investigated in the scope of this study. The 11 time-steps correspond to motions observed during 6 to 10 minutes.

The FEM-generated deformations were then used to warp the aforementioned base 3D MR scan, resulting in a series of eleven MR images (the original scan was also included), which depict the simulated anatomical changes under the influence of the simulated motion. Due to hardware limitations, the MR images and gold standard vector fields were re-sampled onto an image matrix of size 256 × 256 × 128 (corresponding to a resolution of 1.63 × 1.63 × 0.98 mm 3 ).

In addition to the synthetic MR sequence, a clinical cine-MR series acquired on the MR-Linac consisting of sixteen images with a temporal resolution of 16.9 seconds acquired during beam-on time (patient in the supine position, same acquisition parameters as above) was also submitted to the analysis.

Simulating a radiation therapy fraction

For the systematic dose error analysis, we simulated the delivery of a single fraction for a prostate tumor. We established a treatment plan using MatRad [START_REF] Wieser | Development of the open-source dose calculation and optimization toolkit matrad[END_REF] with clinical constraints that are frequently used at the Radiotherapy department of the UMC Utrecht for prostate cancer therapy (8 Gy to the prostate, maximally 4 Gy to the rectum and bladder, using step-and-shoot IMRT as is clinically performed on the MR-Linac with 7 photon beams at 0 • , 50 • , 100 • , 150 • , 210 • , 260 • and 310 • angles). The simulated treatment time is 10 minutes and 14 seconds. The electron density and Houndsfield units were extracted from the clinical planning CT-image of the patient, pre-registered to the first MR-image in the cine series. Subsequently, both electron density and Houndsfield units were warped using the FEM-generated DVFs, which in turn allowed the calculation of the absorbed radiation dose for each of the deformed anatomical states. To simulate the aggregate dose absorption under the influence of motion, MatRad was modified to allow the simulation of sequential dose absorption for all different gantry and MLC positions with an updated anatomical situation and to store the corresponding dose fragments separately. In summary, this strategy provided for each of the simulated time points a complete set of (I) anatomical MRIs, (II) electron density and Houndsfield units, and (III) the resulting dose absorption maps for each delivery configuration of the Linac in the known frame-of-reference of the deformed anatomy. Finally, both the gold standard DVFs and those estimated by the DIR algorithms were used to accumulate the dose onto the original anatomy, as detailed in section 2.3.3. This regularization represents an assumption on the pattern of estimated deformations, which essentially reduces the dimensionality of the solution space. The regularization parameter α determines the relative weight of the regularization. One of the most frequently used regularization terms for DIR algorithms constrains the estimated deformations to be spatially smooth, and is defined as

R s = ∇u( r) 2 2 , (1)
where u = (u 1 , u 2 , u 3 ) is the 3D DVF and • 2 is the Euclidean norm. Aside smoothness, previous studies have also employed incompressibility as a constraint on the estimated deformations [START_REF] Yang | Singular features in sea surface temperature data[END_REF][START_REF] Haber | Numerical methods for volume preserving image registration[END_REF][START_REF] Haber | Image registration with guaranteed displacement regularity[END_REF][START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF]). This can be applied when, due to their high water content, biological soft tissues (e.g. liver, kidney, prostate etc.) are not expected to undergo volumetric changes at the time scale of a radiotherapy fraction. The corresponding mathematical formulation of this incompressibility regularization is

R i = J( r + u( r)) -1 2 2 , (2)
where we have used the Jacobian determinant of the deformation:

J( r + u( r)) = 1 + ∂u 1 ( r) ∂x ∂u 1 ( r) ∂y ∂u 1 ( r) ∂z ∂u 2 ( r) ∂x 1 + ∂u 2 ( r) ∂y ∂u 2 ( r) ∂z ∂u 3 ( r) ∂x ∂u 3 ( r) ∂y 1 + ∂u 3 ( r) ∂z . (3)
The five DIR algorithms used in this project are listed below. Their specifics are detailed in the Appendix A.

• The original optical flow (OF22), introduced in the seminal paper by [START_REF] Horn | Determining optical flow, Techniques and Applications of Image Understanding[END_REF]. Assumes that moving voxels conserve their brightness, combined with a smooth regularisation. The implementation was adapted from (Zachiu, Denis de Senneville, Moonen & Ries 2015a) ‡.

• Improved optical flow (OF21), introduced in (Zachiu, Papadakis, Ries, Moonen & Denis de Senneville 2015). Substitutes the L 2 -norm of the original optical flow by an L 1 -norm, making the algorithms less sensitive to local grey-level intensity variations.

• EVOlution (EVO), introduced by Denis de [START_REF] Denis De Senneville | EVolution: an edge-based variational method for non-rigid multi-modal image registration[END_REF]. This algorithm aims to align similar contrast patterns by matching image gradients. Also uses a smooth regularisation.

• Incompressible EVOlution (EVI), introduced in [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF]. Using the same data-fidelity term as EVO, combined with an incompressibility-regularisation.

• Spatially adaptive EVOlution (AEVO), introduced by Zachiu, Denis de Senneville, Willigenburg, de Boer, Raaymakers, Ries et al. ( 2020). Uses a spatially-adaptive ‡ See http://bsenneville.free.fr/RealTITracker/ combination of both regularisations, allowing incompressibility in certain regions (e.g. bones and the prostate) and smoothness in others, combined with the same data-fidelity term as EVO.

In this study, an exhaustive search in the parameter space of the regularization parameter(s) was employed, minimizing the endpoint error (further discussed in section 2.4) between the estimated and gold standard deformations. The results are shown in section 3.1.

2.3.2.

Resolution & signal-to-noise-ratio of the MR-images MRI acquisitions typically imply a compromise between spatial resolution, temporal resolution and SNR. Since both accuracy and precision of image registration algorithms depend on both the spatial resolution and SNR of the underlying images, we evaluate the impact of both factors on the MCDA results.

To assess the influence of the spatial resolution, the synthetic cine-MR data were downsampled with an isotropic factor of 2, 3, and 4. A sagittal slice for the four resolutions is shown in figure 3. After performing the registration for each downsampling (12,9,6,4,3) or PSNR-values of (30,26,22,19,17). A sagittal slice for the resulting images is shown in figure 4.

Figure 4: A sagittal slice of the reference MR-image for the different levels of SNR (excluding the original SNR already shown in figure 3).

Dose warping strategies

For a motion-compensated dose accumulation, the delivered dose for all moving images is warped to the reference image using their respective DVF. In this study, two different strategies are used for this warping: direct dose mapping (DDM) and energy per mass transfer (EMT), as detailed in [START_REF] Li | Direct dose mapping versus energy/mass transfer mapping for 4d dose accumulation: fundamental differences and dosimetric consequences[END_REF]. Due to the deformable character of the motion, tissues from multiple voxels in the moving anatomy might end up, for example, in a single reference voxel. Using DDM in this situation, the reference voxel dose is found by averaging the doses from the moving voxels. EMT on the other hand, is motivated by the physical definition of dose as the energy absorbed by a mass. This method warps both the energy and mass from the moving voxels and then divides the energy by the mass to get the reference dose. The EMT method thus requires a mass-map for the anatomy, which we created from the CT-image, differentiating between air (a density of ρ = 1.225 g/L was assumed) and soft tissue (density ρ = 1000 g/L).

Error measures

To assess the quality of an estimated DVF, it is compared with the gold standard DVF in terms of absolute, component-wise differences, or the endpoint error [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF]:

∆DV F (x, t) = DVF gs (x, t) -DVF est (x, t) 2 , (4)
i.e. the Euclidean distance between the 3D deformation vectors of the gold standard and the estimated vector field. This results in a scalar for each voxel, providing a distribution of errors. In a clinical setting, however, the gold standard DVF is not available. Therefore, an alternative surrogate quality assurance criterion is required for evaluations of the clinical data included in this study. Evaluated on our simulated datasets, we found the structural similarity index (SSI) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]) and image intensity differences (IID) of the MR-data inside a contour or set of contours to have a Pearson correlation coefficient with the mean endpoint error of 0.90 and 0.95, respectively. This value is much higher than found for the Dice similarity coefficient (Dice 1945, Sorensen 1948) (0.79) or Hausdorff distance [START_REF] Hausdorff | Measurement of signal intensities in the presence of noise in mr images[END_REF]) (0.57). Therefore, the former two criteria will be used to evaluate the registration performances on the clinical dataset. Additionally, we use the Jacobian determinant in incompressible regions to investigate the physiological plausibility of a DVF, as suggested in [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF], Zachiu, Denis de Senneville, Raaymakers & Ries 2020). For the intra-fraction investigation, only minor volume variations of the prostate are expected (and therefore modeled), such that large deviations of the Jacobian determinant from one are an indication of implausible deformations.

To assess the performance of MCDA, the accumulated dose using an estimated DVF and the dose warping method of choice (i.e. DDM or EMT) is compared to the accumulated dose resulting from the gold standard DVF in combination with the EMT dose warping method. This choice is motivated by the physical definition of dose (i.e. amount of absorbed radiation energy per unit mass). The relative difference in grays is what we will refer to as the dose (accumulation) error. As for the endpoint error, this gives a value for each voxel, providing a distribution of errors.

Results

Deformable image registration algorithm configuration

The approach used to investigate and optimize the regularization parameters for the five DIR algorithms is described in section 2.3.1. The optimal values for the regularization parameters used in this study are reported in table 1. Of the four single-parameter algorithms, EVO and EVI are most robust with respect to a choice of α, while OF22 is the most sensitive. Furthermore, OF22 has its optimum with respect to the endpoint error at the lowest analyzed regularization level. For EVO and EVI, the error is rather stable for 0.3 ≤ α ≤ 1, but EVI does not converge in all voxels for α ≤ 0.2. For the interested reader, the complete results of the endpoint error averaged over the bladder, AEVO has two free parameters: α for the smooth regularization and β for the strength of the incompressible regularization. Their influence on the endpoint error can be seen in the appendix in figure A2. To reduce computation time, only the last time point for the rectal filling motion is evaluated. The regularization parameter space displays a global minimum for the endpoint error at 0.1 ≤ α ≤ 0.2, 0.2 ≤ β ≤ 0.8.

Influence of the deformable image registration algorithm & dose warping method

In figure 5, a boxplot of the distribution of the absolute dose accumulation errors on the prostate in [Gy] (grays) for the four different motions are shown. The central red line in the boxplots indicates the median while the bottom and top edges of the box mark the 25 th and 75 th percentiles, respectively. The indicators connected by the dotted lines (the 'whiskers') give the most extreme data points not considered outliers, while the outliers are shown as black dots (for a normally distributed variable they consist of 0.7% of the data). The results for performing no registration and the five DIR algorithms are presented, using both DDM and EMT to accumulate the dose. From the figure, it is clear that the choice of dose mapping method generally has a much larger effect on the error than choosing the optimal DIR algorithm. When performing no MCDA, the median dose error on the prostate is 0.42 Gy for the rectal filling motion. Using DDM, this is reduced to 0.24-0.28 Gy for the different algorithms. When using EMT, however, the median dose error decreases by a factor 6-8 to 0.06-0.07 Gy. Furthermore, using EMT outperforms DDM in every case.

A noteworthy observation is that for the residual motion pattern (bottom right graph in figure 5), OF22 and OF21 coupled with a DDM dose warping strategy lead to larger dose errors compared to the case where no motion compensation is performed. On the other hand, when using EMT or one of the EVOlution-based algorithms for MCDA, the resulting dose error is slightly decreasing.

It should be noted that all y-axes are scaled to the maximum dose error of performing no registration. The (out-of-range) maxima are OF22 (7.1 Gy) for the rectal filling motion, OF22 (3.3 Gy) for the bladder filling motion, OF22 (6.0 Gy) and OF21 (3.0 Gy) for the average prostate motion, and OF22 (4.3 Gy) and OF21 (2.0 Gy) for only residual motion, all associated with a DDM dose warping strategy.

Figure 5: Absolute dose accumulation errors for the four motion patterns on the prostate in grays. From left to right in each graph are shown the results for performing no registration, the original optical flow, the improved optical flow, EVOlution, incompressible EVOlution, and spatially adaptive EVOlution, both using direct dose mapping (DDM) and energy per mass transfer (EMT). Noticeable is that EMT (greatly) outperforms DDM in every case and that MCDA using EMT can decrease dose errors by a factor of 2.3 to 8.

In figure 6, the dose accumulation errors on the rectal wall are shown for the rectal filling motion. Here, the relative (i.e. not using the absolute value) errors are shown. This offers little additional information for the prostate region as all errors are fairly symmetric around zero for this anatomy. For the rectal wall, however, the noncompensated dose error is not centered around zero. For this particular deformation pattern, where the rectum would move into the high dose region, performing no MCDA underestimates the dose received by the rectal wall, by a median of almost 2 Gy, ranging up to 6 Gy, on a plan with a dose constraint of 4 Gy to the rectum. Performing any MCDA centers the dose error around zero.

The median absolute dose accumulation error on the rectal wall for no MCDA is 1.9 Gy. This error is reduced to 0.2-0.4 Gy when using DDM. When using EMT, the error is reduced by a factor of 12-16 to 0.12-0.16 Gy. For the other three motion patterns, the results are qualitatively similar but of a lower magnitude. It is also worth noting that MCDA always decreases the maximal and median dose error, except when using an optical flow algorithm combined with DDM. Moreover, similar to the results obtained on the prostate (see figure 5), EMT outperforms DDM in every case.

For the interested reader, we have also included figure A3 in the appendix, showcasing an example of the spatial distribution of the dose accumulation errors for one of the motion patterns (i.e. rectal filling motion). The resulting observations are well in line with the ones made above. For the reader interested in comparing the quality of the registration only, the distribution of endpoint errors within the prostate, bladder and rectal wall for all time point of the four different simulated motion patterns are shown in figure A5 in the appendix. The lowest maximum and median endpoint errors are achieved by either OF21 or EVO, while the highest outliers arise from OF22.

Figure 6: Relative dose accumulation errors for the rectal filling motion on the rectal wall in grays. Shown from left to right are the results for performing no registration and for the five DIR algorithms using direct dose mapping (DDM) and energy per mass transfer (EMT). Without MCDA, there is an underestimation of the dose delivered to the rectal wall of 2 Gy. This asymmetry vanishes for all algorithms, and the median absolute dose error decreases by a factor of 12-16 when using EMT.

Figure 7 shows the Jacobian determinants within the prostate of the simulated and estimated DVFs for the bladder filling motion. The results for the other three motion patterns (not shown) are qualitatively similar. The deformations from all the simulated time points were included in the analysis. OF21, OF22, and EVO have a wider range of Jacobian determinants (0.93 -1.09, 0.91 -1.10, and 0.93 -1.07, excluding outliers) compared to EVI and AEVO (0.94 -1.04 and 0.96 -1.04), which are closer to the gold standard values (1.00 -1.00). For OF22, the outliers even range from -0.25 to 3.0, indicating tissue folding and extreme expansion.

Figure 7: Jacobian determinants on the prostate for the bladder filling motion pattern for the gold standard, the original optical flow, the improved optical flow, EVOlution, incompressible EVOlution, and spatially adaptive EVOlution (from left to right). Note that the gold standard shows only slight deviations from 1. EVI and AEVO produce a narrow range around 1, while OF21, OF22, and EVO show a wider spread.

For completion, the endpoint error distribution on the prostate, bladder and rectal wall is shown in figure A5 in the appendix, showing that performing image registration considerably reduces the endpoint error with errors after registration mostly under half a voxel.

The influence of image resolution and signal-to-noise-ratio

Figure 8 shows the influence of resolution on the dose accumulation errors within the prostate and rectal wall for the rectal filling and bladder filling motion patterns. The dose accumulation errors for the other two motion patterns (not shown) are qualitatively in-between those shown, with slight relative differences between the optical flow and the EVOlution-based algorithms. Given the observations made in section 3.2, we only report the results obtained by using EMT, for better graph readability. Note that the six leftmost boxes correspond to registrations performed at the original image resolution and are the same as in figure 5.

Except for OF22 on the lowest resolution on the prostate for the bladder filling motion, the maximum and median absolute dose accumulation error is decreased by performing MCDA with any DIR algorithm on any resolution, for both anatomies and shown motion patterns. For the bladder filling motion, the median absolute dose accumulation error on the prostate decreases with a factor 5-6 on the highest resolution, 4-5 on the second-highest, 3-4 on the second-lowest, and 1.3-2.6 on the lowest resolution. On the rectal wall, this becomes a factor 10-17 on the highest, 2-9 on the second-highest, 3-9 on the second-lowest, and 1.4-3.2 on the lowest resolution. EVO and AEVO are the most robust algorithms for lower resolutions, but the differences between algorithms are relatively small. Figure 8: Absolute dose accumulation error on the prostate and rectal wall for different algorithms at no, two-fold, three-fold and four-fold down-sampled images for the rectal filling and bladder filling motions, respectively. The resolutions in millimeters are provided in the graphs. The DIR algorithms shown from left to right are no registration, the original optical flow, the improved optical flow, EVOlution, incompressible EVOlution, and spatially adaptive EVOlution, all combined with energy per mass transfer. The effect of lowering the resolution is multiple factors higher than the differences between algorithms. Even on the lowest resolution performing MCDA decreases the dose error notably, while a resolution below 2 mm 3 is optimal.

Figure 9 shows the influence of SNR on the absolute dose accumulation errors within the prostate and rectal wall, again for the rectal filling and bladder filling motions. The results for the other two motion patterns (not shown) are qualitatively similar to those of the bladder filling motion. It can be observed that all DIR algorithms decrease the median as well as the maximum absolute dose accumulation error on both anatomies for all levels of the SNR. For the bladder filling motion, the median absolute dose accumulation errors on the prostate decrease by a factor of 2.0 to 4.3, even for the lowest SNR. On the rectal wall, this becomes a factor of 3.4 to 6.8.

Evaluation on the clinical data

For the clinical dataset, the prostate experiences a sudden displacement between the 8 th and 10 th time point of about 4-5 millimeters in the anterior direction and about 3 millimeters in the superior direction. The prostate then only slowly drifts towards its original position for about 0.5-1.5 and 1 millimeters, respectively, during the final 6-7 time-steps. For the interested reader, an assessment of the motion as estimated by the five used DIR algorithms is shown in figure A4 in the appendix.

As discussed in section 2.4, for the clinical cine-MR series, we use the SSI and IID to evaluate the registration performances since no gold standard DVF is available. Note that we evaluate 1-SSI (as the SSI gives a value between 0 and 1) such that it increases for larger registration errors, as does the endpoint error. Figure 10 showcases the SSI (in (a)) and the IID (in (b)) quantified jointly inside the prostate, bladder, and rectum wall. Note that the same regularization parameters as for the simulated data were used. On a qualitative level, the results are similar in nature to those obtained for the simulated data illustrated in figure A5: the differences between algorithms are relatively small, and the registration error decreases notably for all algorithms compared to performing no registration. In terms of both SSI and IID, it can be observed that EVI performs worst of all algorithms. Excluding EVI, the median registration error decreases by a factor of 1.4 to 3.0, depending on the DIR algorithm and quality assurance criterion used. Figure 10 (c) displays the Jacobian determinants on the prostate for all time points of the clinical dataset. As for the simulated data, OF22, OF21, and EVO show a wider range of Jacobian determinants for voxels inside the prostate (0.64 -1.41, 0.72 -1.31, and 0.86 -1.13, excluding outliers), compared to EVI and AEVO (0.95 -1.04 and 0.93 -1.07, respectively). For OF22, the outliers even range from -0.5 to 6.2, indicating tissue folding and extreme expansion.

Due to its competitive registration errors (figure 10 (a) and (b)) and better anatomical plausibility of the estimated deformations in terms of Jacobian determinant values (figure 10 (c)), the accumulated dose provided by the AEVO algorithm was chosen as a point of comparison to analyze dose accumulation errors on the clinical data (the results on the simulated data also supported this choice). Note that for dose warping, the EMT strategy was employed as before.

The resulting dose accumulation errors are shown in figure 11, showing great similarities to the results for the simulated data. Without registration, the dose accumulation errors on the prostate range towards 2.5 Gy, with outliers towards 5.6 Gy. As for the simulated data, the dose warping method has a larger impact than the choice of the DIR algorithm. Moreover, as for the simulated rectal filling motion, there is a median underestimation of the dose delivered to the rectal wall when not performing any MCDA of almost 2 Gy, ranging up to 4.5 Gy while the rectum had a dose constraint of 4 Gy in the therapy plan. This asymmetry is again addressed by all algorithms and warping methods.

Figure 11: Dose accumulation errors for the clinical data (AEVO with EMT was used as a point of comparison). The absolute error on the prostate (left) and the relative error on the rectal wall (right) are shown for the five different DIR algorithms, as well as no registration, for energy per mass transfer (EMT) and direct dose mapping (DDM). The difference between the dose warping methods is much larger than between DIR algorithms, with EMT outperforming DDM substantially. All methods resolve the underestimation of the dose received by the rectal wall when not performing MCDA.

Discussion

All presented dose accumulation results show that performing MCDA using variational algorithms for motion mitigation can considerably reduce errors in the estimated delivered dose for prostate cancer therapy. Figure 5 shows that using an EVOlutionbased algorithm, or EMT, always results in a decrease in the dose error. When using EMT, the median dose error on the prostate is reduced by a factor of 1.5 to 7, depending on the motion pattern and DIR algorithm. For the rectal wall, the median dose error decreases with a factor of 2.5 to 16 for the non-residual motion patterns. Additionally, MCDA resolves any asymmetry or bias in the uncorrected dose error. Of similar importance is the observation that this type of MCDA in the absence of largescale biological motion does not degrade dosimetry compared to the uncorrected dose estimate.

Dose warping strategy

The choice of dose warping strategy has the most substantial impact on the dose accumulation error of all factors considered in this project, with EMT always outperforming DDM. This reconfirms the findings of [START_REF] Siebers | An energy transfer method for 4d monte carlo dose calculation[END_REF] on a deformable phantom (1.1% difference between EMT and DDM along the beam and up to 25% in the penumbra), [START_REF] Zhong | Monte carlo dose mapping on deforming anatomy[END_REF] on a lung plan (7.3% difference) and [START_REF] Li | Direct dose mapping versus energy/mass transfer mapping for 4d dose accumulation: fundamental differences and dosimetric consequences[END_REF] on a lung plan (11.3% difference in the planning target volume and 4.4% in the internal target volume minimum doses). We observe median dose differences of approximately 2.5% on the prostate and 10% on the rectal wall.

These earlier works evidenced that the typical situations where DDM becomes increasingly problematic are deformations of boundary interfaces displaying high tissue density variations and tissue areas that display internal density fluctuations. This trend can clearly be seen in figure 5, where the difference between using DDM and EMT is found to be larger for algorithms for which the Jacobian determinants are unequal to one (figure 7), i.e. where compression or expansion is estimated. This is due to the fact that EMT intrinsically compensates for density modifications introduced by the realignment process, while DDM is very susceptible to this effect. As a consequence, algorithms that have been optimized to minimize density variations (such as EVI and AEVO) display a considerably better performance when paired with DDM, as shown in figure 5 and6. Nevertheless, even in this configuration, EMT resulted in a lower dose error compared to DDM. It is noteworthy that the usage of EMT compared to DDM introduces only a small level of additional complexity and thus computing time.

Different deformable image registration algorithms

First and foremost, the selection of algorithms in the scope of this paper has been motivated by the established track record of variational algorithms concerning robust registration performance even for image data of low SNR and resolution [START_REF] Glitzner | On-line 3d motion estimation using low resolution mri[END_REF], Zachiu, Denis de Senneville, Moonen & Ries 2015b[START_REF] Zachiu | An improved optical flow tracking technique for real-time mr-guided beam therapies in moving organs[END_REF][START_REF] Denis De Senneville | EVolution: an edge-based variational method for non-rigid multi-modal image registration[END_REF]. Since variational algorithms also employ 'surrounding' information of respective voxels, accuracy and precision can, for coherent deformations, typically exceed the image resolution. We observed that the 75 th percentile of the endpoint error for all motion patterns and all algorithms is around half a voxel.

When using DDM, clear differences in dose error between DIR algorithms are observed. EVI and AEVO with biophysical plausibility (in the form of an incompressibility regularization) built-in outperform the other algorithms, in particular OF22 and OF21. The biophysical plausibility is observed in the narrow range of Jacobian determinants (figure 7), as described in [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF], Zachiu, Denis de Senneville, Raaymakers & Ries 2020). Exemplary for the implausible deformations is for example the tendency of the OF22 algorithm to estimate tissue folding in problematic regions. This phenomenon is most likely caused by violations of the brightness constancy assumption. Using all gray-level information, the optical flow algorithm is more likely to be biased by image artifacts to converge into local minima. Although the effect of these unrealistic deformations on the accumulated dose is reduced by using the EMT dose warping strategy, such an approach may not be physically applicable for warping other quantitative information such as Houndsfield units, diffusion/perfusion values, mass maps, etc.

Furthermore, an asset of the EVOlution-based algorithms is that they perform well for a wide range of regularization parameters (figures A1 and A2 in the appendix). This can be of considerable clinical value in particular for online scenarios, where the tuning of regularization parameters with the patient on the treatment table is either impossible or for a smooth work-flow highly undesirable.

Based on these grounds, we consider AEVO the optimal choice for the DIR algorithm, as the clinical data shows its superiority over EVI in terms of registration errors for sudden deformations.

Resolution and SNR

Generally, a resolution of about 1 mm 3 seems good, 2 mm 3 seems adequate, and from 3 mm 3 onward, we start to under-sample. This corresponds to approximately half the size of some relevant structures like the rectal wall (8 by 24 mm) and the dose gradient (11 to 2 Gy) between the prostate and rectal wall (7-11 mm).

The EVOlution-based algorithms prove to be the most robust with respect to lower resolutions. Interestingly, for these algorithms, even at the lowest resolution, the median dose error is lower for all simulated motions and anatomies compared to the one obtained by performing no motion mitigation. For the rectal filling and bladder filling motion, the error on these low-resolution images is still reduced by a factor of two. It is thus overall beneficial to perform MCDA even when only compensating for bulk motion; however, the optimal resolution value resulting from the present study is above 2 mm 3 . All variational methods are fairly robust and still perform well for a low SNR. The difference in median absolute dose error for the rectal filling motion between the lowest and highest SNR is only a factor of 1.5 to 3.5. EVI is the optimal DIR algorithm when working with images of very low SNR due to incompressibility being a strong-natured regularization. For all DIR algorithms, the step in dose errors from an SNR of 12 to an SNR of 9 is about as big as from an SNR of 9 to an SNR of 3. Using an SNR above 10 thus gives the optimal results, but even for an SNR of 3, performing MCDA (using EMT) is beneficial.

Clinical data

The dose errors found for the clinical data (figure 11) confirm the findings of the simulated data. First, performing MCDA again proves to be able to decrease errors in the estimated dose considerably. Secondly, the choice of dose warping strategy again impacts the dose errors the most. Finally, as for the simulated rectal filling motion, when not performing MCDA, there is a significant underestimation of the dose delivered to the rectal wall that is resolved by performing MCDA.

The dose accumulation errors for the clinical data reported in figures 11 are found using AEVO with EMT as a the point of comparison. AEVO was used as it produces particularly low registration differences (figure 10 (a) and (b)) and a narrow Jacobian determinant range within the prostate (figure 10 (c)) as well as competitive dose errors for the simulated data (figures 5 and 6).

Combined with the simulated data results, there is a strong case for performing MCDA in scenarios where the delivered dose is of interest. This will provide another step in the direction of hypofractionation and intra-fractionally adapted treatments that count on reliable MCDA for further development and implementation.

The results presented here for the prostate are expected to carry over to other tumor sites with motion patterns having a slow motion component (like peristaltic events, respiratory baseline, muscle relaxation or digestive activity). For inter-fraction dose accumulation, deformations are generally more extreme than for intra-fraction, and it is to be expected that using EMT over DDM has an even larger impact.

Conclusion

The presented study analyzed the influence of the choice of DIR algorithm, dose warping strategy, MR image resolution, and MR image SNR on the resulting motioncompensated accumulated dose. For a quantitative investigation, deformations based on clinical observations were simulated using finite element modeling, which provided a gold standard for both the deformations estimated by the registration algorithms and the accumulated dose. The evaluation was also extended to include clinical patient data.

The analysis conducted on both the simulated and the clinical data demonstrates that, in general, performing MCDA leads to a notable improvement of the geometric uncertainty with respect to the radiation dose delivery. We have found that, for an optimal selection of the investigated influence factors, for large-scale biological motion such as the passage of a rectal filling, the dose error when performing MCDA decreases by a factor of 7 and 15 within the anatomical structures of interest. While not to the same extent, MCDA continued to provide improvements in geometric uncertainty even for the lower end of the image resolution and SNR investigated within the scope of this study. Moreover, MCDA did not introduce any relevant errors for the case in which only minimal motion is present.

Out of all the investigated factors, we conclude that the most important one impacting the accumulated dose is the dose warping strategy. The differences between the two dose warping strategies were considerably larger than the differences between the registration algorithms. Furthermore, using energy per mass transfer over direct dose mapping decreases dose errors in all cases, at no extra scan time, and only implies limited extra computing power.

In terms of the employed registration algorithm, it was found that, due to its accuracy, anatomical plausibility, and robustness to deviations of the regularization parameters, AEVO is a marginally better choice than the rest of the investigated algorithms. Finally, an isotropic image resolution of 2 mm and an image SNR above 10 were found to be optimal for an accurate MCDA.

Overall, the present work provides an extensive analysis of the dependency of MCDA on multiple contributing factors, providing a guideline for their choice in accordance with the specific technical and functional specifications of ongoing and future MRgRT workflows.

anatomy. Knowledge of the underlying anatomy can thus be used to enforce incompressibility within e.g. biological soft tissues, while allowing other regions to expand or compress, using the smoothness regularization. In effect, this algorithm employs two regularization parameters: one for the incompressible region and one for the smooth region.

Figure 1 :Figure 2 :

 12 Figure 1: Average displacements of the prostate in millimeters for the four different motion patterns exported from the finite element models. The motion-patterns are modelled after motions observed in patients during 6 to 10 minutes.

2. 3 .

 3 Considered factors influencing motion-compensated dose accumulation 2.3.1. Deformable image registration algorithms In the scope of this study, we have selected five different variational DIR algorithms previously proposed in the context of MRgRT (Zachiu, Papadakis, Ries, Moonen & Denis de Senneville 2015, Zachiu, Denis de Senneville, Willigenburg, de Boer, Raaymakers, Ries et al. 2020). A variational registration algorithm typically implies the minimization of a two-part cost function: E = D + αR. The data fidelity part D is a function quantifying the dissimilarity between the images. Since the minimization of the data fidelity term alone is usually an under-determined problem, a regularization term R is added to ensure well-posedness.

Figure 3 :

 3 Figure 3: A sagittal slice of the reference MR-image for the four evaluated resolutions. The top-left image displays the original resolution.

Figure 9 :

 9 Figure 9: Absolute dose accumulation errors on the prostate and rectal wall for different values of the SNR and different algorithms. The rectal filling and bladder filling motion patterns are shown. The SNR is provided within the graph. No registration, original optical flow, improved optical flow, EVOlution, incompressible EVOlution, and spatially adaptive EVOlution are shown from left to right, all using energy per mass transfer. The difference in dose errors between an SNR of 12 and 9 is about as big as the difference between an SNR of 9 and 3.

Figure 10 :

 10 Figure 10: (a) the time-averaged structural similarity index (shown as 1-SSI) and (b) image intensity differences (IID) within the bladder, prostate, and rectal wall for the clinical data. The results for performing no registration, the original optical flow, the improved optical flow, EVOlution, incompressible EVOlution, and spatially adaptive EVOlution are shown, from left to right. For both error measures, the differences between algorithms are relatively small, and the error decreases considerably compared to no registration. (c) the Jacobian determinants for the original optical flow, the improved optical flow, EVOlution, incompressible EVOlution, and spatially adaptive EVOlution on the prostate for the clinical data. OF22, OF21, and EVO show a wider range of values than EVI and AEVO.

Figure A1 :

 A1 Figure A1: Endpoint error averaged over the voxels of the bladder, prostate, and rectal wall for different values of the regularization parameter α for the original optical flow (OF22), the improved optical flow (OF21), EVOlution (EVO), and incompressible EVOlution (EVI), respectively. The different colors correspond to the four simulated motion patterns. The vertical green line depicts the optimal choice for the regularization parameter used in this paper. OF22 in particular is very sensitive to the correct value of the regularization parameter. TODO: change legend to rectal filling motion and adapt scale for optical flows α → α 2 !.

Figure A2 :

 A2 Figure A2: Average endpoint error on the bladder, prostate, and rectal wall in voxels for different values of the two regularization parameters for spatially adaptive EVOlution. Results for the last time point of the rectal filling motion are shown. The arrow points at the optimal configuration used in this paper.

Figure A3 :

 A3 Figure A3: Sagittal slice of delivered dose distribution and the dose accumulation errors in grays for the rectal filling motion pattern. Shown are the results for no MCDA, original optical flow, improved optical flow, EVOlution, incompressible EVOlution, and spatially adaptive EVOlution combined with direct dose mapping (DDM) or energy per mass transfer (EMT) to warp the dose. The prostate (center) and rectum (right) contours are shown in black. A large influence of the choice of dose warping method can be observed.

Figure A4 :

 A4 Figure A4: Average displacement of the prostate as estimated by the five DIR algorithms for the clinical dataset. There is a sudden shift in the position of the prostate in the anterior (4-5 mm) and superior (3 mm) directions between the 8 th and 10 th time point. Thereafter it slowly drifts back in the posterior (0.5-1.5 mm) and inferior (1 mm) directions.

Figure A5 :

 A5 Figure A5: Image registration performance in terms of endpoint errors in millimeters on the prostate, bladder, and rectal wall for the different DIR algorithms and motion patterns. From left to right in each graph are shown the results for performing no deformable image registration, the original optical flow, the improved optical flow, EVOlution, incompressible EVOlution, and spatially adaptive EVOlution. All algorithms lower the median endpoint error considerably, from a factor of 5-13 for the rectal filling motion to a factor of 1.4-2.1 for residual motion. Note that the y-axis is scaled differently for the four different motion patterns.

Table 1 :

 1 Optimal values for the regularization parameters. , and rectal wall as a function of α for the four single-parameter algorithms for the four different simulated motion patterns is shown in figureA1in the appendix.

	DIR algorithm	α	β
	OF22	0.05
	OF21	0.35
	EVO	0.30
	EVI	0.35
	AEVO	0.10	0.30

prostate

Acknowledgments

The authors wish to thank Daan de Muinck Keizer and Hans de Boer for helpful discussions. This work has been supported in part by the Topconsortia for Knowledge and Innovation -LifeSciences & Health, project no. LSHM19031 (POPCORN).

Appendix A. DIR algorithms

The specifics of the five DIR algorithms used in this project are outlined below.

• The original optical flow (OF22). It relies on the assumption that voxels conserve their gray-level intensity as they move. This implies minimizing a cost function which uses the square Euclidean distance between the images as a data fidelity term, coupled with a smoothness regularization:

where (I x , I y , I z ) are the image spatial partial derivatives, I t is the temporal derivative.

• Improved optical flow (OF21). Replaces the quadratic data fidelity term from the original optical flow by a linear L 1 norm, reducing the impact of local gray-level inconsistencies between the images to be registered. In effect, this implies the minimization of the following cost function:

• EVOlution (EVO). This algorithm aims to align similar contrast patterns between the registered images. This is achieved by the minimization of the function

where Γ is a neighborhood around r and I 1,2 are the reference and moving images, respectively.

• Incompressible EVOlution (EVI). This algorithm uses the same data fidelity term as EVO, but replaces the smoothness regularization term with a global penalty on incompressibility.

• Spatially adaptive EVOlution (AEVO). Employs the same data fidelity term as EVO and EVI. However, it uses the smoothness and the incompressibility regularization selectively, depending on the physical properties of the underlying


Julien Vidalie
email: julien.vidalie@irt-systemx.fr

Michel Batteux

Jean-Yves Choley

Faïda Mhenni

Mohamed-Sami Kendel

Julien Vidalie Quartz

Typology of the differences Between Model-Based System Engineering (MBSE) and Safety Assessment (MBSA) models: Analysis of a Reference System

Keywords: MBSE, MBSA, Model consistency, differences typology, AltaRica 3.0, SysML

  L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Developing complex systems is a multidisciplinary process that requires many different models to be used to represent and simulate one same system for different purposes. As those different models are made by different actors, it is obvious that there is a high risk that they present some inconsistencies. In this context there is a need for verification of consistency of the MBSA model with the MBSE model it is derived from. System Engineering and Safety Assessment are two disciplines that are deeply correlated. The system engineer and safety analyst need to work together to prove that the system is safe, especially in industrial fields that require safety certification.

This paper aims at providing a typology of differences that can occur between a Model Based System Engineering (MBSE) model and a Model Based Safety Assessment (MBSA) model. It is a basis for further work centered on MBSE/MBSA synchronization. This typology does sort differences between models based on their causes and allows engineers to get a better understanding over what should or shouldn't be corrected in the models. This is the key for a future formal definition of what is an inconsistency.

To establish this typology we modeled a reference system, a landing gear study case [START_REF] Boniol | The Landing Gear System Case Study[END_REF], from both the MBSE and MBSA point of view, and analyzed the differences that occurred between those models.

The remainder of the paper is organized as follows. Section II describes the landing gear study case that we used for this study. Section III presents the modeling of the system in both MBSE and MBSA points of view. Section IV presents some previous work about model synchronization, comparison of the models and foundations of our typology of differences. Finally the conclusion is given in the last section, along with some perspectives for future work.

II. THE LANDING GEAR STUDY CASE

A. Study case

The landing gear study case was described in [START_REF] Boniol | The Landing Gear System Case Study[END_REF] and served as a benchmark for techniques and tools for the assertion of system behavior. This system is a standard aircraft landing gear composed of three gears (front, rear-left and rearright). It describes the system pilot interface, its mechanical and hydraulic parts, and its digital control part.

This system is also relevant to a MBSA study, since aeronautic systems are required to be compliant with CS25 regulations, an aeronautic recommended practice describing safety analyses that are authorized to be completed on aircraft equipment for certification is the ARP4761a [START_REF] International | Guidelines and methods for conducting the safety assesment process on civil airborne systems and equipment ARP4761[END_REF]. This study case will also be used in further work for the definition of a mathematical framework around the S2ML language [START_REF] Batteux | From Models of Structures to Structures of Models[END_REF], which will be referred to in the IVth section of this paper, and for MBSE/MBSA synchronization.

For this work, two models of this system were created, the first one aiming at modeling the system architecture, made with the Cameo System Modeler tool with the SysML language [START_REF] Friedenthal | A Practical Guide to SysML[END_REF], and the second one is a safety analysis view that was created using the OpenAltaRica tool with the AltaRica 3.0 [START_REF] Batteux | AltaRica 3.0 in 10 Modeling Patterns[END_REF] modeling language. To comply with the aeronautic certification requirement, these models have to be separated and to be made by separated people. This allows safety analysis to independently verify the compliance of the architecture described in the MBSE model. In this work we reproduce a realistic workflow, with models written by two different people and did not eliminate differences before the final review. By this protocol we want our workflow to present realistic differences, and we want not to avoid inconsistencies by not having independency between both models. The creation of the MBSA model was based on the MBSE model, and we aim at detecting differences that occurred in this creation.

B. Introduction to the Landing Gear System

As depicted in Fig. 1, the system is composed of three main parts :

1) The Pilot Interface

This part allows the system to communicate with the pilot. It is composed of a Handle which is used by the pilot to order the system to be up or down. This handle communicates its position to the digital part.

Three lights indicate to the pilot the status of all three gears using the following code:

 Green light: "Gear locked down"

 Orange light: "Gear maneuvering"

 Red light: "Landing gear system failure"

 No light: "Gear locked up"

2) The Mechanical and Hydraulical Parts

The structure of the Hydraulic part is described in Fig. 2. The system is composed of three landing sets (front, rear-left and rear-right). Each set has a box (containing all the components, to be fitted in the aircraft landing gear well), door which is opened and closed by a cylinder and a landing gear that is extended and retracted by another cylinder.

The hydraulic power is provided to the cylinders from the aircraft hydraulic circuit by a set of electro-valves:  One general electro-valve supplies all the system from the aircraft hydraulic circuit  One electro-valve provides pressure to the portion of the system related to door opening  One electro-valve provides pressure to the portion of the system related to door closing  One electro-valve provides pressure to the portion of the system related to gear extending  One electro-valve provides pressure to the portion of the system related to gear retracting Each valve is controlled by electrical order from the digital part. 3

) The Digital Part

The digital part is composed of two redondant computing modules that run the same control software. It is used to receive and analyze data from the system sensors, and to command the system's components.

III. SYSTEM MODELING

A. MBSE modeling 1)

Methodology In our study, the MBSE modeling was realized following the SysML methodology described in [START_REF] Mhenni | A SysML-based methodology for mechatronic systems architectural design[END_REF]. This methodology first focuses on a black box analysis of the system describing requirements, system context, lifecycle and operational scenarios. Then some white box views of the system represent its functional and physical structure in addition to its behavior.

For synchronization of the MBSA model we only focus on the white box views of the system. In fact, our interest is about the structural and behavioral features of our models.

2)

Modeling The system architecture is modeled around its 3 main subsystems which can be observed in the Block Definition Diagram (BDD) shown in Fig. 3. This BDD shows the system's breakdown structure. Arrows in the diagram represents composition links, meaning that one block (or system) is composed of the components that arrows are pointing to. This view of the system shows us its structure but does not define its architecture for we lack connections and flows between the components, which are specified in the Internal Block Diagram (IBD) in Fig. 5. In this IBD, we can observe the relations between components. The names of the connections between them give us information over their type of flows. Finally they can also contain typed variables. It is important to note that this model does not intend to simulate our system but rather to give us a communication tool over it, and means for traceability. Therefore naming on this model is of high importance and carries a lot of information compared to usual simulation models where the content of variables, flows, and other quantitative values are the most important information carried by the model.

The two views we previously described are the ones that will be interesting in the context of MBSE/MBSA synchronization for the scope of this work. In further work we will also investigate the state machines diagrams and sequence diagrams that describe the functional behavior of the system.

B. MBSA Modeling 1) Methodology

The MBSA modeling of a system can be part of its safety analysis and is a method that is accepted by authorities as a safety analysis method with its integration in the protocol described in ARP 4761a [START_REF] International | Guidelines and methods for conducting the safety assesment process on civil airborne systems and equipment ARP4761[END_REF] regulations.

Although most safety analysis nowadays is carried through analysis such as Fault Tree Analysis (FTA), MBSA modeling is a good tool for safety analysis thanks to its high expressivity. It is easier to understand and to communicate with this model, making it easier for the safety analyst to shows problems in the architecture to the system engineer.

The MBSA model that we made for this work was created using the OpenAltaRica platform, based on the AltaRica 3.0 modeling language. It represents the system through its structure and dysfunctional behavior. Unlike the MBSE tool SysML, which is a graphical notation, AltaRica is a formal language, meaning a textual syntax and a semantic, even though some AltaRica tools such as Simfia provide graphical interface to design parts of the AltaRica model.

2) Modeling

The MBSA modeling was achieved using the article [START_REF] Boniol | The Landing Gear System Case Study[END_REF] presenting the system as a reference document and based on the MBSE modeling it aims at verifying. We here considered the MBSE model as a specification document of the system and expect to verify its compliance to safety requirements through MBSA modeling. Except for a few differences which will be talked over in section IV, the model has a very similar structure to the one presented in the SysML IBD and BDD in Fig. 3 and Fig. 5 respectively. In AltaRica 3.0, we represent the system by a main "block" which is a container that will be considered and simulated by our tool. The components of our system are described in classes that are instantiated in the main block. Fig. 7 describes the main block representing the system and Fig. 8 and Fig. 9 describes the class ElectroValve and Cylinder which are components of the landing gear system. Instances of those components are linked in the system to allow for retraction and extension of the landing gear with assertions such as the one presented in Fig. 6, this allows for the input of the cylinder to be at all time equal to the value of the output of the electrovalve. All components in our system are extending the NonReparaibleComponent class presented in Fig. 4. This class describes the state machine for failure of a component. Those events are characterized by delays. This allows the execution to compute the time after which the transition shall be fired. Delays are described using probability distributions such as, in our study case, an exponential distribution. They also allow for computation of probability of event happening for the generation of fault trees in the case of static systems.

From this general class we derive all components of the system, specializing this class by adding new variables that are ports of our components and assertions that represents connections between these variables. As an example Fig. 6 represents such a connection. It means that the output value of the udHandle component of the pilot interface is given to the input of the CPIOM1 component of the digital part. This means that during the execution, the value of the output of the udHandle component, which is the handle used by the pilot to actuate the landing gear, will be given to the input of the first CPIOM unit of the digital part. Assertions can also be used to give values to variables based on component state or other information.

The interest of having this formal representation of the system rather than using a notation such as SysML is that it allows for formal computation over the system safety. Thanks to this model we are able to compute minimal cut sets of the system with their probabilities. We can also execute stochastic simulation of the system with failures, and identify propagation paths of the failures. This wouldn't be possible if there was any ambiguity in the representation of the system. Whereas for human communication with the MBSE model this isn't an issue. MBSE and MBSA are two important parts of the design of the system, but as we explained it before, they serve two different purposes. Therefore, differences occur between those models, should it be for lack of communication between the teams or for more fundamental reasons linked to the nature of those models. These differences could lead to models presenting two distinct and different systems instead of representing the same model. Synchronization between MBSE and MBSA models is thus necessary to ensure consistency. Work is in progress to create consistency methods that will be described in IV.A. We think that in order to improve upon this work, there is a need for a formal definition of inconsistencies. In order to prepare this definition we established the typology of differences that is described in IV.B.

A. Existing methods for MBSE/MBSA synchronisation

As synchronization between heterogeneous models is an important concern for researchers in the modeling field, some tools were already developed to help model synchronization, even in the case of MBSE/MBSA synchronization.

In [START_REF] Rauzy | Foundations for model-based systems engineering and model-based safety assessment[END_REF], the authors suggest a synchronization methodology based on three phases: Abstraction, Comparison and Concretization, illustrated in Fig. 10. The different models are first translated to intermediary models written in a same formalism that will allow comparison, this is the abstraction phase. Those intermediary models are then compared to detect the differences that exist between them, this is comparison. Finally the source models are annotated with the differences that are detected, and corrective actions are proposed to the designers, this is concretization.

The SmartSync synchronization framework [START_REF] Batteux | Synchronization of System Architecture and Safety Models : a Proof of Concept[END_REF] does provide a methodology and tools for synchronization between MBSE and MBSA derived from this methodology. This is achieved by the abstraction of both models to a common formalism, the S2ML language [START_REF] Batteux | From Models of Structures to Structures of Models[END_REF], then computational comparison of the two abstracted models and finally a concretization step, where inconsistencies are taken into account and the models are adjusted to match each other.

The work proposed in [START_REF] Berriche | Towards Model Synchronization in Model Driven Engineering of Mechatronic Systems[END_REF] has a similar perspective as the ones in [START_REF] Rauzy | Foundations for model-based systems engineering and model-based safety assessment[END_REF] proposal. It focuses on the topological aspect of models, and performs the comparison relying on graph theory after abstracting the different models into graphs. This methodology carries a flat view of the model. This methodology is to be applied by an "interface expert" that communicates with specific model designers. The MOISE project provided the system engineer and safety analyst with a method to create "synchronisation points" that are linked on one side to elements of the MBSE model and to elements of the MBSA model on the other side. This allows the users to review consistency by reviewing each consistency point and attributing it a status and a rationale to justify that elements of MBSE and MBSA models represent the same component.

All those solutions highlight the need for an intermediate representation of the models that enables the comparison. However we feel that there is still a need for a typology and a formal definition of the inconsistencies between MBSE and MBSA models.

B. Comparison of models

In order to identify the types of differences, we reviewed both MBSE and MBSA models together and searched for all differences that occurred between them. We will not cite all of them here because some are the same on different elements of the model but we will list all different ones. The first differences that occur when reviewing the models are the names of the elements of the model. The names can vary between both models, for example the handle from the pilot interface is called "Handle" in the MBSE model and "udHandle" in the MBSA model. Naming differences are linked to the system engineer and the safety analyst calling the same elements in different ways, which can happen because of their different technical backgrounds, or their naming practices (upper/lower cases, shortcuts…). Other naming inconsistencies may come from the naming rules (that can be either restrictions or only recommendations) in both tools. For example the pilot interface subsystem is called "Pilot interface" in the MBSE model and "pilotInterface" in the AltaRica 3.0 model, this difference is caused by the differences of uses of the two models: AltaRica 3.0, as a formal modeling language, is used to compute reliability or safety indicators, thus objects shall be represented in one word; whereas SysML is used to declare and communicate and as a consequence, there are no restrictions on naming objects.

It also happens that some elements of both models can sometimes not be named, whereas their counterpart in the other model is. Assertions in AltaRica 3.0 are unnamed, they serve a pure purpose of simulation, defining rules for the calculation of variables values. On the contrary connections in the SysML model are sometimes named using the type of flow or actions they convey, when the ports they connect are often unnamed by the engineer and automatically named p1, p2, p3… by the modeler. As the previous one, this difference is also caused by the differences of uses of the two models and by the modeling habits of both engineers.

Most variables in the AltaRica model are typed with discrete values, this is explained by the fact that we are mainly interested in knowing whether they are nominal or dysfunctional and not by their accurate values. The system engineer and the safety analyst do not have the same point of view on the system. Therefore they do not represent system variables in the same way. Those variables type are sometimes not affected to the ports in the SysML model but to the connections, depending on the engineer modeling habits. Because the system engineer wants to prove that the system answers requirements that can be linked to those values, so they are usually typed as their physical unit.

We also observe different structural differences between our models.

The first one is that some subsystems are specified at different levels of abstraction. For example, the arborescence for the electro-valve that brings hydraulic pressure to the door extension cylinders is "Landing gear system\mecahydraulic part\Hydraulic part\Open door Electro-Valve" in the MBSE model and "LandingSys\mechahydraulicPart\hydraulicSys\ DoorHydraulicSys\extensionDoorElectroValve". Although the naming differs, subsystem levels between those paths match apart for the "DoorHydraulicSys" level in the MBSA model that does not match any subsystem in the MBSE model. This originates from the Safety analyst regrouping parts in a different way, which may be due to his wish to only consider failure of a group of parts rather than every unique part, and it could also happen that the MBSA model only specified "DoorHydraulicSys" without modeling the electro-valve inside itSome components have also been specified at a different place of the Product Breakdown Structure. This is the case of the cylinders in our comparison. The system engineer considered them to be part of the hydraulic set of parts, whereas the safety analyst did put them in the landing sets along with the gears/doors they are connected to. Such a difference could be caused either by a different point of view over the system as it is here, or by a modeling error.

Finally, in the Internal Block Diagram we observe some connections to the outside of the system that are not considered in the AltaRica model, for example the "Electric Power" input, this is due to the need for the system engineer to represent all interactions within the system and with its environment. However, even though this connection has a real impact on the system, it was not considered relevant for safety analysis of the landing gear system and thus not modeled with AltaRica. Such a difference could be either considered a modeling error or not depending on whether that connection has an impact or not on safety analysis.

We note that although it is not the case in our models, some connections could have been placed between ports that do not necessarily exist in the MBSE or MBSA representations of the system. This could either occur by modeling error, or because those values are not relevant to one or the other modeling intent.

C. Typology of differences

From this comparison we deduce three main types of differences in our model that are in fact related to the cause of the differences between the models. Some of those differences are caused by modeling errors and lead to the models describing different systems, those are inconsistencies, but we also note differences that are due to modeling practices and tools.

The first type of differences we encountered is related to differences that are caused by the different modeling tools and practices. Examples of this are different naming rules or connections between ports or assertions that relate variables carrying names in different ways. This type of differences could be handled by modeling practices or rules in certain cases, for example the implementation of naming rules in SysML similar to the ones that exist in other modeling languages. It could also be taken into account in the comparison by not taking into account names that have no counterpart, or by cleverly comparing them, for example in our case, comparing connections names in SysML to variables names in AltaRica. The second type of differences that we denote are differences linked to modeling intent. This is the case of the difference in abstraction that could occur between both modeling as we showed it with the electro-valve, or in the case of the different value types that are observed. These differences are necessary for both the system engineer and safety analyst to work correctly. Their existence is the reason for having two separate models rather than modeling all information in one global model.

Our third type of differences are the ones caused by modeling errors. The aim of model synchronization is to eliminate those differences. They can be different naming, wrong values, incorrect links between components, etc. These inconsistencies will probably be more difficult to recognize from the second type of differences (modeling intent) than the first (modeling practices and tools).

We also raise another interesting way to classify differences. We encountered differences between the models that were either related to a particular element of the system, such as a naming difference or a variable type difference, or differences that were related to the structure of the model, such as abstraction differences, different placement of an element in the Product Breakdown Structure, or wrong connections between component ports/variables. These two kinds of differences are also interesting because we intuit that they should translate very differently in a mathematical framework around the models.

Type

Observed Since system engineering and safety analysis are serving different purposes, they need to be using two different models and two different formalisms. This results in a risk of inconsistency between those models.

In this paper, we proposed a typology of the differences that may exist between MBSE and MBSA models. This is made possible by MBSA and MBSE modeling of a reference system, and analyzing the differences noticed between the two models.

This typology shows that differences can be sorted by their causes. These causes are the use of different modeling tools and practices, different modeling intents, and finally modeling errors. We think that differences due to modeling intents and standards are important to the models since they are the reason for two models being used instead of one, and they should be preserved. Whereas inconsistencies due to modeling errors should be eliminated from the model, and differences due to different modeling tools and practices should be reduced as much as possible by unifying naming practices for example.

On the basis of this typology, we want to write a formal definition of what is an inconsistency between MBSE and MBSA models. This will be helpful in the formalization of a consistency assessment method and model reconciliation. We think that reconciliation requires to develop strategies to eliminate unacceptable inconsistencies, i.e. ones caused by modeling errors, and reduce other differences to an acceptable threshold.

Moreover, in further work, we believe that the formalization of a mathematical framework supporting those models will help us with consistency assessment, and therefore we intend in formalizing the S2ML language used in the SmartSync methodology with the category theory formalism, to allow better comparison between models, and easier translation from SysML and AltaRica to S2ML.

Fig. 1 :

 1 Fig. 1 : Global Architecture of the system

Fig. 2 .

 2 Fig.2. : Architecture of the hydraulic part[START_REF] Boniol | The Landing Gear System Case Study[END_REF]

Fig. 4 :

 4 Fig. 4 : Class NonRepairableComponent in AltaRica 3.0

  domain nrpState {OK, KO} class NonRepairableComponent nrpState s (init = OK); parameter Real lambda = 1.0e-5; event failure (delay = exponential (lambda)); transition failure: s == OK -> s := KO; end

Fig. 3 :

 3 Fig. 3 : Product Breakdown Structure of the Landing Gear System (BDD)

  digitalPart.CPIOM1.input:= pilotInterface.udHandle.output; Fig. 6 : Example of assertion in AltaRica

Fig. 5 :

 5 Fig. 5 : System Physical Architecture (IBD)

Fig. 10 :

 10 Fig.10: Model Synchronization approach[START_REF] Rauzy | Foundations for model-based systems engineering and model-based safety assessment[END_REF]

Fig. 11 :

 11 Fig. 11 : Typology of differences associated to examples from VI.B

ACKNOWLEDGMENTS

The authors want to thank the S2C project at IRT SystemX and IRT Saint-Exupéry and all its industrial and academic partners for the funding of this research and the trust they give us.