Mateo Villa

Julien Bert

Antoine Valeri

Ulrike Schick

Dimitris Visvikis

Fast Monte Carlo-based Inverse Planning for Prostate Brachytherapy by Using Deep Learning

Keywords: Brachytherapy, Dosimetry, Deep Learning, Inverse Planning

Inverse planning is an essential tool for optimizing the delivered radiation dose on low-dose-rate (LDR) prostate brachytherapy. Clinical inverse planning systems use the TG-43 dose computation formalism in order to perform a fast optimization. However, this method is an approximation that often leads to a dose overestimation, resulting on sub-optimal plans. Alternatively, Monte Carlo simulation (MCS) can be used to obtain an accurate dose distribution, but considerably increasing the estimation time. We propose a fast inverse planning method for LDR prostate brachytherapy that uses a deep convolutional neural network (DCNN) trained on a GPU-based MCS generated database to estimate the dose distribution on the prostate and organs at risk. Segmentations of the organs and seeds' positions are given as the DCNN input. The mean percent error on the test set was -1.19±0.94% within the prostate. The DCNN was used to estimate the dosimetric parameters in each organ for every configuration of the optimization loop. The dosimetric parameters of the final DCNN-based brachytherapy plans were in good agreement compared to the same plans recalculated with a full MCS. The proposed inverse planning based on DCNN was capable to reach an equivalent level of accuracy with Monte Carlo with a runtime in less than 1 min using conventional GPU card.

I. INTRODUCTION

B RACHYTHERAPY is an internal radiotherapy technique widely used on the treatment of the prostate cancer. It involves the application of radioactive sources through the prostate on the tumor proximity, with the aim of destroying the tumor while preserving healthy distant organs and tissues from the radiation exposure. Low-dose-rate (LDR) brachytherapy is one of the most popular for the case of low-stage prostate cancer. In this method, multiple small seeds of radioactive material such as iodine 125 ( 125 I) or paladium 103 ( 103 Pd) are implanted and remain permanently into the affected organ region. The seeds implantation is done by using hollow needles which are inserted through the perineum, with the guide of trans-rectal ultrasound (TRUS) imaging. The intervention is thus less invasive than in other kind of procedures such as external radiotherapy or high-dose-rate (HDR) brachytherapy.

Nevertheless, the spatial distribution of the inserted seeds during LDR brachytherapy is highly determinant on the resulting radiation dose. Therefore, an improper configuration plan may result on harmful side effects caused by high and All the authors are with LaTIM, INSERM, Faculté de Médicine, IBRBS 22, Av. Camille Desmoulins, 29200 Brest, France. (email: ingmatvillaa@gmail.com)

This work was funded by the French National Research Agency through the FOCUS project (ANR-16-CE19-0011), the OptimiX project (ANR-18-CE45-0011), and within the Investissements d'Avenir program (LabEx CAMI 2) under reference ANR-11-LABX-0004. heterogeneous zones of radiation. For avoiding such issues, inverse planning is an essential tool for guiding LDR prostate brachytherapy routine. In this technique, an optimization algorithm searches the most suitable seeds locations by minimizing a cost function, which is composed by pre-defined dose objectives. The final evaluation of the procedure is made by using dose-volume histograms (DVHs), which give information about the dose homogeneity and intensity on the prostate and organs at risk (OARs), specifically the urethra and rectum. The DVHs must fill different criteria specified by the American Association of Physicist in Medicine (AAPM) on their Task Group (TG) TG-137 report [START_REF] Nath | Executive Summary of Report of Task Group 137[END_REF]. Fast simulated Annealing (FSA) [START_REF] Lessard | Inverse planning anatomy-based dose optimization for HDR-brachytherapy of the prostate using fast simulated annealing algorithm and dedicated objective function[END_REF] is one of the most popular optimization algorithms used in practice. Other kind of approaches such as genetic algorithms [START_REF] Yu | A genetic algorithm for the optimization of prostate implants[END_REF], [START_REF] Mcgeachy | An opensource genetic algorithm for determining optimal seed distributions for low-dose-rate prostate brachytherapy[END_REF] and compressed sensing inspired solvers [START_REF] Guthier | A new optimization method using a compressed sensing inspired solver for real-time {LDR}-brachytherapy treatment planning[END_REF] have also been proposed. All of these methods allow to achieve plans in less than one minute. However, in the current clinical applications, the dose objectives are not optimized on the entire organs' volume, but only on points of the surface. As a result, clinical planning systems tend to obtain plans that does not fulfil the DVHs criteria entirely, and a manually correction by an expert radiologist should often be made, consequently increasing the time and risk of the intervention. DVH-based optimization has been successfully implemented on inverse planning systems for external radiotherapy, but they are not yet used in real practice on prostate LDR brachytherapy due to its elevate time computation.

Moreover, clinical inverse planning systems compute the dose distribution based on the TG-43 formalism, an analytical model established by the AAPM [START_REF] Rivard | Update of AAPM Task Group No. 43 Report[END_REF], [START_REF] Rivard | Supplement to the 2004 update of the AAPM Task Group No. 43 Report[END_REF]. This method, regardless of allowing a fast computation, it does not consider the organs densities and the interaction between the sources, causing an overestimation of the dose [START_REF] Rivard | The evolution of brachytherapy treatment planning[END_REF]- [START_REF] Carrier | Postimplant Dosimetry Using a Monte Carlo Dose Calculation Engine: A New Clinical Standard[END_REF], that could reach up to 12 % for the D 90 parameter (the dose absorbed by at least the 90 % of the prostate volume) [START_REF] Carrier | Impact of interseed attenuation and tissue composition for permanent prostate implants[END_REF]. Monte Carlo (MC) dosimetry can be used to obtain a more precise estimation. Nevertheless, computation time is not suitable for a clinical application. Different Monte Carlo simulation (MCS) frameworks based on Graphics Processing Units (GPU) have been proposed [START_REF] Sempau | DPM, a fast, accurate Monte Carlo code optimized for photon and electron radiotherapy treatment planning dose calculations[END_REF]- [START_REF] Lemaréchal | GGEMS-Brachy: GPU GEant4-based Monte Carlo simulation for brachytherapy applications[END_REF], and they proved a reduction on computation time from multiple hours to seconds for a single dose map. Even though, the inverse planning algorithms usually perform a thousand of iterations before obtaining a satisfactory solution, leading to a very long-time calculation to get the final optimized plan. A DVH-based optimization algorithm using MCS dosimetry for LDR prostate brachytherapy has previously proposed [START_REF] Mountris | DVH-Based Inverse Planning Using Monte Carlo Dosimetry for LDR Prostate Brachytherapy Radiation Oncology[END_REF]. In this work, MCS is used to compute dose kernels for each single seed position, followed by the use of FSA. At each optimization step, the current dose kernel is obtained by adding the dose distribution of the selected seeds. This method achieved treatment plans within less than one minute, with most of the DVH critical values matching the criteria. Still, this method does not include the effect of seeds' interaction, and the need of storing a large data set of dose kernels may rise issues on memory management and input/output time access. Additionally, the number of simulated particles used to estimate each dose kernel had to be reduced in order to achieve a fast planning, which was reflected in an increase of uncertainty.

On the other hand, deep convolutional neural networks have been broadly used in the last years in many medical imaging applications due to their ability to accurately learn a mapping between an input image and corresponding output. Recent works in the field of radiotherapy and medical imaging dosimetry have been released, with the aim to predict the resulting dose in to the patient's body. The general approach consist to train a deep convolutional neural network (DCNN) with a medical image (CT, PET, MRI) or a segmentation map of the organs as the input and the desired dose map as an output. In a recent work [START_REF] Mao | RapidBrachyDL: Rapid Radiation Dose Calculations in Brachytherapy via Deep Learning[END_REF], a DCNN was trained to learn the dosimetry prediction in HDR brachytherapy for a given treatment plan, based on the patient's image and organs' segmentation. However, due to time computation to build a large training data set from MCS, they pre-train a network using water-based dose. The dose prediction time was 1.7 seconds. In HDR there is only one source at the time within the patient and the source only interacts with human tissues. Such method cannot be used in LDR, since the interactions between the seeds have to be considered. The water-based pre-training method may introduce bias because the dose in water and seed are very different. In addition, the optimization in LDR needs a large number of iterations to explore every possible seed configuration. A much faster dose prediction is required in LDR.

We propose in this paper a fast inverse planning method for LDR prostate brachytherapy that uses a deep convolutional neural network (DCNN) trained on a GPU-based MCS generated database to estimate the dose distribution on the prostate and organs at risk (OARs). Segmentation maps of the concerned organs and seeds' positions are given as the input of the DCNN, and the dose distribution obtained by Monte Carlo simulation is the output to be learned. The DCNN was used to estimate the dosimetric parameters in each organ for every configurations test by the optimization loop.

II. METHODS AND MATERIALS

We propose the use of a DCNN which learns to compute the dose into the prostate and OARs for a given patient anatomy and set of seeds' positions. Figure 1 show a slice example of the input and output. A fast GPU-based MCS is used to generate the dose maps database. The trained model is then used to obtain the dose at each step of the inverse planning. The following sections explains the approach in detail.

A. Patient's Database

Prostate and OARs contours from 273 ultrasound (US) image patients who were undergone with LDR brachytherapy at University of Brest Hospital were extracted. The whole data set contains prostate volumes with sizes varying between 17 cm 3 and 67 cm 3 .

The original ultrasound data had a size of 550 x 550 x 88 voxels with a pixel spacing of 0.15 mm x 0.15 mm and a slice thickness of 1 mm .

For each US image patient, 50 layout of seed positions were generated randomly, with a number of seeds varying randomly between 48 and 70. This lower and upper values were obtained according the number of seeds usually implanted to the patient, which depends on the prostate volume size. The space of possible positions follows the standard grid of 5 mm spacing used during the needles' insertion. The resulting number of configurations was 13650 for the entire data set.

B. MC Dose Computation

Since MCS usually involves long time computation, a generation of a large dose map data set would be difficult with conventional MCS software. Therefore, dose distributions on the prostate and OARs were generated by using GGEMS, a validated GPU-based MC simulation platform for dosimetry applications including LDR brachytherapy [START_REF] Bert | Geant4-based Monte Carlo simulations on GPU for medical applications[END_REF], [START_REF] Lemaréchal | GGEMS-Brachy: GPU GEant4-based Monte Carlo simulation for brachytherapy applications[END_REF]. The simulation was performed by using 125 I sources at the desired positions on a voxelized phantom generated from the patient contours. The complete geometry of the seeds (STM1251, Bard Medical Division, Covington, GA, USA) was taken into account, similarly to [START_REF] Lemaréchal | GGEMS-Brachy: GPU GEant4-based Monte Carlo simulation for brachytherapy applications[END_REF]. For each organ label (prostate, urethra, rectum) the corresponding material composition was assigned, based on the recommendations of the International Commission on Radiation Units and Measurements Report No. 46 [START_REF] Bethesda | Icru report n. 46: Photon, electron, proton and neutron interaction data for body tissues[END_REF]. The surrounding area of the aforementioned organs is considered as a soft-tissue and its material is also included in order to fully take into account tissue heterogeneities. GGEMS software was based on Geant4 physics effects. The GPUbased MC simulation was setup to use only photon physics effects without considering electrons. Electron ranges for the 125 I sources are very small compared to the size of a voxel [START_REF] Chibani | MCPI©: A sub-minute Monte Carlo dose calculation engine for prostate implants[END_REF]. Therefore, their doses were deposited locally within each voxel they were emitted from. The dose calculation was performed using the track-length estimator [START_REF] Williamson | Monte Carlo evaluation of kerma at a point for photon transport problems[END_REF] which is a variance reduction technique allowing to reduce time computation without introducing any approximations. A total of 5 × 10 7 particles were simulated, allowing to achieve approximately 0.5 % mean statistical uncertainty per dose map on the prostate volume, 0.6 % on the urethra and 1 % on the rectum. Simulations took on average 20 seconds approximately per configuration on a NVIDIA GTX 1080 GPU.

C. Deep Convolutional Neural Network (DCNN)

We propose a DCNN model that uses as input the patient's 3D segmentation map, composed by the prostate, urethra, and the rectum. Additionally, the different seeds positions were directly encoded on the image by defining a 1-pixel label for each seed and included on the organs' segmentation map, as shown in the figure 1. In order that the network considers organ structures and seed positions, all information was numerically encoded by using a label value inside the input image. For example, all voxels containing prostate tissue provided by the segmentation were set to the label value 1. The same for the urethra and the rectum with label value 2 and 3 respectively. Voxels containing soft tissue were filled with the value zero. The voxel that corresponds to a seed position was filled with a label value 4. An example of final labelled image is illustrated on figure 1 (see the input image). Contour of each organ was achieved by the urologist using the transrectal US image while the patient is under general anesthesia in operating room. Then an inverse planning is run to determine an optimal plan, this is the standard clinical workflow on prostate brachytherapy. The seeds label will be changed during the optimization process, when a new seed position is explored, the current seed position is cleared by filling the corresponding voxel this prostate label, and the voxel that will receive the new seed position is filled with the label 4.. The output of the network will be the prediction of 3D dose distribution according organ contours and a given seed layouts.

1) DCNN Architecture: We designed a DCNN based on the 3D-Unet [START_REF] Abdulkadir | 3D U-net: Learning dense volumetric segmentation from sparse annotation[END_REF], which was originally intended for volumetric semantic segmentation. The architecture, which is shown in Figure 2 consist of two main stages: the contracting path (encoder), which extracts different features of the images by applying 3D convolutional filters and reduces their dimensionality using 3D max-pooling layers; and the expanding path (decoder), which gradually up-samples the images by using 3D deconvolutions (transposed convolutions) while concatenating them with the different features of the contracting path, until arriving to an output of the same dimension of the input. The main difference with the original architecture lies in the number of filters per convolution layer, which in our case were reduced by a factor of 8, in order to avoid over-fitting and being able to have a fast prediction. Additionally, the use of a Rectified-Linear-Unit (ReLU) activation layer instead of the softmax layer (which is more suitable for classification) at output layer enabled the network to be trained for this regression task.

The original images were re-sampled to a resolution of 1 mm x 1 mm x 1 mm (the slice thickness did not change) in order to allow the network to perform a fast inference, leading to image with input dimensions of 88 x 88 x 88 voxels. These images were centered at the prostate isocenter and a zeropadding was performed in order to match the input size of the neural network.

2) Training and Validation: The data set was divided patient-wise in a stratified fashion according to 6 ranges of prostate sizes. Training, validation, and test sets were created by randomly taking 60 %, 20 % and 20 % of the patients of each interval, as shown in Figure 3. This strategy allows to have a representative sample of the data set on each of the splits, thus reducing the risk of over-fitting. The data split results on 163, 55 and 55 patients for each set. With the 50 generated plans per patient (and corresponding dose maps), the final split consists of 8150 data for train, 2750 for validation and 2750 for test.

The training phase was made using adaptive momentum algorithm (ADAM) [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] for minimizing the mean squared error (MSE) between the predicted dose and the MCS dose. During the network training, both training and validation loss were monitored. Consequently, the train process was stopped at the epoch when the validation loss ceased showing improvement with respect to training loss (i.e., the 17 th epoch), as shown in figure 4. The architecture was developed and trained on Keras, using Tensorflow back-end, using a Nvidia GeForce RTX 2080Ti GPU.

D. DVH-based FSA for Inverse Planning

Inspired by the work of [START_REF] Mountris | DVH-Based Inverse Planning Using Monte Carlo Dosimetry for LDR Prostate Brachytherapy Radiation Oncology[END_REF], we performed an inverse planning based on dose-volume-histograms (DVHs) by using the Fast Simulated Annealing (FSA) algorithm [START_REF] Lessard | Inverse planning anatomy-based dose optimization for HDR-brachytherapy of the prostate using fast simulated annealing algorithm and dedicated objective function[END_REF]. The cost function to minimize, given by the equation 1, includes the different DVH dosimetric criteria given by the AAPM TG-137 recommendations [START_REF] Nath | Executive Summary of Report of Task Group 137[END_REF]:

C = w v100 * Θ(V 100 L -V 100 )+ i w vi Θ(V i -V i H ) + j w dj Θ(D j -D j H )+ k w dk Θ(D k -D k H ) + w N needles * N needles (1)
where V i represents the percent of the prostate volume that receives at i% of the prescribed dose, D j represents the percent of the prescribed dose absorbed by j % of the urethra's volume, and D k the amount of dose received in kcm 3 of the rectum's volume. The subscripts H and L represent the critical values (L for lowest, H for highest) of the aforementioned metrics. N is the number of needles used to insert the seeds. In the same way w vi , w dj , w dk and w N are weight factors given to each metric. Here, in order to give the same importance to all parameters, each weight was equally distributed i.e. with a value of 1/N w , where N w is the number of weights uses in the cost function. Θ(X) represents the heavyside step function that overrides the cost contribution when its value is negative. The indices values are i = [150, 200] and j = [START_REF] Carrier | Postimplant Dosimetry Using a Monte Carlo Dose Calculation Engine: A New Clinical Standard[END_REF]30]

and k = [2cc, 0.1cc].
Before optimization, a fixed number of seeds is selected for the current patient, using the equation introduced by [START_REF] Pouliot | Optimization of permanent 125I prostate implants using fast simulated annealing[END_REF]:

N seeds = 4 + 4.674 × V 0.562 A ( 2
)
where V is the volume of the prostate, and A represents the strength of the source. The FSA outline starts with a random configuration of seeds. The resulting dose map is estimated using our trained DCNN and then the corresponding DVHs from the prostate, urethra and rectum are extracted. The cost function is calculated with the derived DVH dosimetric values. For the following steps, each configuration is obtained by randomly changing a single seed position. If the cost of the new configuration is lower than the previous one, then the new one is maintained. Differently, the new configuration is accepted with the following probability:

P (∆C) = exp(-∆C/T (k)) ( 3
)
Where k is the current optimization iteration, ∆C is the difference between the current and previous cost function, and T is the annealing temperature at the k th iteration. The temperature is changed at each step by applying a pre-defined cooling-rate factor between 0 and 1. The optimization takes end after a pre-defined number of iterations. This schedule makes that at the first iterations, non-conform configurations could be chosen, allowing to explore a wide range of solutions and avoiding the algorithm to get stuck on a local minima. In our study, we empirically set T = 10 5 and C= 1% and max 2000 iterations. We set a fixed source strength of 0.645 U for all the optimizations, which corresponds to the average of the sources used in the real plans. The DVHs were calculated with the same prescribed dose of the real treatments, which is equal to 160 Gy.

E. DCNN Dose Prediction Evaluation

The trained model was used to estimate dose maps on the 2750 images of the test set, which for recall, correspond to 50 random configurations of 55 patients. The quality of the predicted dose maps was assessed in terms of the mean percent error (MPE) between the prediction and the ground-truth for each organ, which is defined as follows:

M P E organ = 100% N organ i D pred i -D M C i D M C i ( 4
)
Where D pred i and D M C i represents the absorbed dose on the i voxel of the organ for the prediction and the MC dose map, and N organ is the total number of voxels of the organ. We report the mean and the standard deviation of the MPE across the entire test set for each organ, as well as the graphical distribution in terms of the median, first and third quartiles, and first and ninth deciles.

F. Inverse Planning Quality Evaluation

The quality of the inverse planning has been evaluated on the 55 patients of the test set. For each patient, using the segmented organs, an inverse planning was performed using the FSA algorithm with our trained DCNN model as the dose calculation model. For each final optimized plan, the seed configuration was exported and used to recalculate a dose map using Monte Carlo simulation with the same parameters presented in II-B. DVH parameters were estimated on both dose map, the one directly obtained by the DCNN-based inverse planing after optimization and on the one recalculated by Monte Carlo simulation. The aim is to assess if after a DCNNbased optimization the equivalent MCS DVH parameters in terms of the TG-137 criteria are fulfilled. The DVH parameters are always calculated with respect to a prescribed dose of 160 Gy, which is the standard dose applied on the real treatments of the database.

III. RESULTS

A. DCNN Evaluation

Training of the DCNN was run on approximately 12 hours and dose prediction was achieved in approximately 38 ms per dose map on a Geforce RTX 2080 Ti GPU, which is highly faster than the 20 seconds of the MCS.

The distributions of the MPE of the DCNN dose prediction for the prostate and OARs along the test set are shown on figure 6. The mean ± STD MPE for the entire test set are -1.19% ±0.94 %, 0.52 % ± 1.42 %, and -1.47 % ± 1.90 % on the prostate, urethra and rectum respectively. According to this results, the DCNN shows a general tendency to underestimate the dose in the prostate and the rectum, while in the urethra the dose is generally over-estimated. Nevertheless, the mean error magnitude does not exceed the limit of 2 % for the prostate and urethra, and 4% for the rectum.

In general it was found that the dose predictions were visually in good agreement with the ground-truth. The figure 5 shows a representative middle-slice example of dose distribution of one of the test set configurations, in comparison to the MCS dose map. We can observe again the general tendency to over-estimation of dose in the urethra, which appears to be predicted with similar values of the prostate neighboring region. In the case of the prostate, we confirm the general trend of under-estimation of the dose, especially in zones of low and middle dose. For the rectum, most of the doses are under-estimated. Despite the fact that some pixels of the image may have errors reaching the 8 %, the absolute difference is still relatively low, and the accuracy the overall image and different organs is not significantly affected.

B. Inverse Planning Quality

The FSA algorithm achieved acceptable configurations in 2000 iterations on average. Adding to the DCNN prediction time (38 ms per dose map) to the time taken by the DVHs extraction and cost function calculation, the total optimization time was on average 1 minute. Which was very fast considering that dose maps accuracy were close to a fully Monte Carlo simulation.

The table I summarizes the different DVH parameters obtained with our DCNN inverse planning approach and their equivalent by running MCS on the resulting plans, as well as the TG-137 criteria for each parameter, with respect to reference prescribed dose of 160 Gy. We can observe an underestimation of the dose for the V 100 and the V 150 which does not exceed the 1% and a very reduced over-estimation of the V 200 of 0.2%. For the urethra, the D 10 and D 30 does not exceed 2 % of difference with the resultant MCS. For the case of the rectum, we can observe an under estimation of the D 0.1cc of 3%, whereas the D 2cc is predicted above 0.4 % of the MCS value, which is still non-significant. In spite of those differences, it can be seen that all the DVH parameters recalculated by MCS on the final configuration complied the TG-137 criteria.

The figure 7 shows a representative example of the DVHs of an optimized plan by using the DCNN dose calculation, compared with the DVHs obtained after doing MCS on the final optimized configuration. The DCNN DVHs fit almost the entire MCS DVHs despite of having the over and underestimation behaviors highlighted previously. Finally, figure 8 shows a comparison of the iso-doses of (contours of the dose map containing the same value of dose) of both DCNN and the resultant MCS, for another typical case. The levels of the iso-doses were chosen based on the usual values analyzed by the radiologists. The DCNN iso-doses match accurately with the MCS, which confirms again that the DCNN has the capability to match the MCS dose distribution in terms of dose homogeneity.

IV. DISCUSSION

The generalization capability of a DCNN model depends primarily on the network model parameters, training hyperparameters, the quality of the data set, and the way it is split for training, validation and test. Therefore, we set the number of convolutional filters per layer as small as possible without losing accuracy with the aim to prevent over-fitting and perform fast predictions. On the other hand, we needed to split the patients in a way that the three ensembles had a good representation of the entire data set. Despite the fact that the geometry of the prostate is relatively consistent, we found that the volume size was substantially variable from one patient to another, and we thus decided to use it as the criteria for splitting. This stratified approach contributed significantly to avoid to over-fitting, and moreover made easy the hyperparameter tuning.

According to the dose prediction results, both errors on the prostate and the urethra were below 2%, and rectum error was below 5%, which is a good indicator of the performance of our method. The lower accuracy on the rectum's dose prediction may be caused by different reasons. First of all, the rectum is generally more distant to the sources than the prostate (which contains all of them) and the urethra (which is located within the prostate region), which is reflected in less radiation dose. Since the DCNN optimizer tends to minimize the mean square error, then differences in low doses would be less relevant on the weights updating than the high doses errors. Another phenomenon frequently observed is that the DCNN has a smoothing effect in regions where the dose level changes in a considerable proportion, as observed in the figure 5 for the dose within the urethra, which appears to be clearly distinguishable in the MCS dose map but in the DCNN prediction the dose it is homogenized to the dose level of the neighboring region on the prostate. This is probably caused by the convolution operations that works on a filtering principle, and thus tend to smooth the changes of intensity. Despite the results on the error distribution and the qualitative observations, the dose prediction of the DCNN and the ground-truth are highly similar. Moreover, the final results on the inverse planning are not significantly affected by this prediction error, as could be stated in the DVH results. Concerning the inverse planning quality results on the DVH parameters, our method tends to find optimal configurations with an under-estimated dose with respect to MCS on the prostate and rectum, whereas the DVH parameters on the urethra are slightly above the MCS values. However, these errors are below 1%, 2% and 3% for the prostate, urethra and rectum respectively. Yet, these errors are acceptable and still much lower than the TG-43 formalism which could reach the 12 %. Moreover, these errors do not affect the final outcome of the dose and all the DVH criteria are fulfilled on the recalculated MCS, which was the main purpose of our study.

Prediction time was highly reduced with respect to the MCS. For the same low level of uncertainty, we reached a dose map computation around 710 times faster than the GPU-based One of the advantages of our DCNN approach compared to others of the state-of-the-art is that it is not dependent of the imaging system, due to the fact that it directly uses the annotated organ contours. However, a good segmentation must be guaranteed in order to avoid inaccuracies on the delivered dose, which is not the focus of our study. Finally, the most important innovation of our approach is that all the seeds' positions are encoded directly into the patient's image, which compared to previous approaches, prevents from estimating each seed effect separately and thus improving computation time and accuracy.

As a future work, we will study the applicability of our approach on focal brachytherapy. Additionally, we will study the feasibility of a DL technique for decrease the number of iterations needed by the FSA optimization.

V. CONCLUSION

We proposed a new DCNN-based method for performing DVH-based inverse planning dosimetry for LDR prostate brachytherapy. The DCNN model is capable to learn accurately from a MCS generated dose map database. Our approach is a potential solution for being adapted on the commercial inverse planning systems as an alternative calculation of the existing TG-43 formalism.

Fig. 1 .

 1 Fig. 1. Example slice of input (left) and output (right) of the DCNN-based method. The light blue label represents the prostate, the light green the urethra, the yellow the rectum, and the red points the seeds of the current configuration. The output used for training is the dose map obtained from MCS.

Fig. 2 .

 2 Fig. 2. Architecture of the proposed DCNN model.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Patients' prostate volume distribution along the training, validation and test sets

Fig. 5 .

 5 Fig. 5. Dose distribution comparison on an intermediate slice for a test set patient. From left to right: The input of the network (segmented organs and seeds), the ground-truth dose map obtained through MCS, the dose map predicted by our DCNN model, the absolute and the relative difference between the predicted and the groun-truth dose.

Fig. 6 .

 6 Fig. 6. Distribution of the mean relative error between DCNN predicted dose and MCS over the entire test set. The middle line represents the median, the central rectangle spans from the first quartile to the third quartile, and the vertical line extends from the first decile to the ninth decile

Fig. 7 .

 7 Fig. 7. DVH obtained after inverse planning on a test set patient by using our DCNN dose map calculation method.

Fig. 8 .

 8 Fig. 8. Iso-dose distribution comparison on a mid-slice between the DCNN and the computed MCS after optimization. The values are represented as % of the prescribed dose, which is equal to 160 Gy.

TABLE I EVALUATION

 I OF THE DVH PARAMETERS RESULTING FROM INVERSE PLANNING USING THE DCNN DOSE CALCULATION METHOD WITH A REFERENCE PRESCRIBED DOSE OF 160 GY.

	Organ	DVH Metric Criteria	DCNN	MCS (Ground-Truth)	Dif (%)
		V 100 (%)	≥ 95	95.1 ± 0.3	95.8 ± 1.0 -0.8 ± 1.0
	Prostate	V 150 (%)	≤ 50	47.4 ± 2.2	48.1 ± 4.4 -0.6 ± 3.1
		V 200 (%)	≤ 20	19.9 ± 1.5	19.7 ± 2.4	0.2 ± 1.3
	Urethra	D 10 (%) D 30 (%)	< 150 < 130	135.2 ± 7.8 123.2 ± 5.3	133.5 ± 8.0 122.3 ± 5.5	1.7 ± 3.5 0.9 ± 2.7
	Rectum	D 2cc (%) D 0.	< 100	71.2± 7.8	70.8 ± 7.6	0.4 ± 2.6

1cc (%) < 150 112.1 ± 15.8 114.9 ± 21.1 -2.8 ± 9.2

ACKNOWLEDGMENT

This work was funded by the French National Research Agency through the FOCUS project (ANR-16-CE19-0011), the OptimiX project (ANR-18-CE45-0011), and within the Investissements d'Avenir program (LabEx CAMI 2) under reference ANR-11-LABX-0004.


Guillaume Billet
email: guillaume.billet@mnhn.fr

Jeremie Bardin

ALLOMETRY OF MOLAR PROPORTIONS IN PLACENTALS Segmental series and size: clade-wide investigation of molar proportions reveals a major evolutionary allometry in the dentition of placental mammals

Keywords: Size, Dentition, Evolution, Segmentation, Phylogeny, Model

Iterative segments such as teeth or limbs are a widespread characteristic of living organisms. While their proportions may be governed by similar developmental rules in vertebrates, there is no emerging pattern as regards their relation to size. Placental mammals span eight orders of magnitude in body size and show a wide spectrum of dietary habits associated with size and reflected in their dentitions, especially molars. Although variation in size constitutes an important determinant for variation in biological traits, few major allometric trends have been documented on placental molars so far. Molar proportions have been intensively explored in placentals in relation to developmental models, but often at a small phylogenetic scale. Here, we analyzed the diversity of upper molar proportions in relation to absolute size in a large sample of placental species (n = 286) encompassing most of the group's dental diversity. Our phylogenetically informed analyses revealed a twofold pattern of evolutionary integration among upper molars: while molars covary in size with each other, their proportions covary with the absolute size of the entire molar field. With increasing absolute size, posterior molars increase in size relative to anterior ones, meaning that large-sized species have relatively large rear molars while the opposite is true for smallsized species. The directionality of proportional increase in the molar row exhibits a previously unsuspected allometric patterning among placentals, showing how large-scale variations in size may have influenced variation in dental morphology. This finding provides new evidence that processes regulating the size of individual molars are integrated with overall patterns of growth and calls for further testing of allometric variation in the dentition and in other segmental series of the vertebrate body.

Placental mammals repeatedly evolved a wide range of body sizes in their history and span over eight orders of magnitude in size today [START_REF] Alroy | Cope's rule and the dynamics of body mass evolution in North American fossil mammals[END_REF][START_REF] Baker | Adaptive evolution toward larger size in mammals[END_REF][START_REF] Price | The macroevolutionary relationship between diet and body mass across mammals[END_REF][START_REF] Bokma | Testing for Depéret's rule (body size increase) in mammals using combined extinct and extant data[END_REF]. Allometry, the size-related changes of traits, constitutes a major and pervasive pattern in morphological evolution [START_REF] Gould | Allometry and size in ontogeny and phylogeny[END_REF][START_REF] Hallgrímsson | Integration and the developmental genetics of allometry[END_REF], and yet it remains only sporadically documented in placentals. Scaling relationships were probably important along their diversification in size and impacted several aspects of their skull [START_REF] Nummela | Scaling of the mammalian middle ear[END_REF][START_REF] Cardini | Larger mammals have longer faces because of size-related constraints on skull form[END_REF]. Placentals also evolved a wide spectrum of dietary habits in association to body size [START_REF] Price | The macroevolutionary relationship between diet and body mass across mammals[END_REF] and dental morphology [START_REF] Ungar | Mammal teeth: origin, evolution, and diversity[END_REF]. However, although mammalian teeth have been studied intensively for their phylogenetic and functional signals, allometric relationships within their dentition have only been sparsely explored (e.g., [START_REF] Kurtén | Observations on allometry in mammalian dentitions: its interpretation and evolutionary significance[END_REF][START_REF] Gould | On the scaling of tooth size in mammals[END_REF][START_REF] Singleton | Allometric and metameric shape variation in Pan mandibular molars: a digital morphometric analysis[END_REF].

Within the differentiated dentition of placentals [START_REF] Ungar | Mammal teeth: origin, evolution, and diversity[END_REF], molars represent an outstanding model to study allometry in a segmental series. Placentals generally present three molars per jaw quadrant, which typically grade in size anteroposteriorly (e.g., [START_REF] Butler | Studies of the Mammalian Dentition. Differentiation of the Post-canine Dentition[END_REF][START_REF] Polly | Development with a bite[END_REF]). Like other teeth, molars are serially repeated organs. Together with limbs, vertebrae and phalanges, they constitute segmental series of the body whose proportions may be governed by similar developmental rules [START_REF] Kavanagh | Predicting evolutionary patterns of mammalian teeth from development[END_REF][START_REF] Young | Shared rules of development predict patterns of evolution in vertebrate segmentation[END_REF][START_REF] Green | Positional information and reaction-diffusion: two big ideas in developmental biology combine[END_REF]. The size of a given segment within these series tends to covary with that of adjacent segments, as proposed by the Inhibitory Cascade (IC) model for lower molars in mice [START_REF] Kavanagh | Predicting evolutionary patterns of mammalian teeth from development[END_REF]). This experimental model is based on an activator-inhibitor logic controlling sequential tooth development, with each molar inhibiting the formation of the following one. Assuming that this cascade system is the sole determinant of molar size differences along the jaw, the second molar will make up one third of the total molar area [START_REF] Kavanagh | Predicting evolutionary patterns of mammalian teeth from development[END_REF][START_REF] Polly | Development with a bite[END_REF]. The diversity of lower molar proportions in placentals was much studied following this finding. These subsequent studies reported varying degrees of adhesion to the IC model across different taxa, and found associations between molar proportions, diet and phylogeny (e.g., [START_REF] Hlusko | Modularity in the mammalian dentition: mice and monkeys share a common dental genetic architecture[END_REF][START_REF] Wilson | Testing a developmental model in the fossil record: molar proportions in South American ungulates[END_REF]Bernal et al. 2013, Halliday and[START_REF] Halliday | Testing the inhibitory cascade model in Mesozoic and Cenozoic mammaliaforms[END_REF][START_REF] Asahara | Unique pattern of dietary adaptation in the dentition of Carnivora: its advantage and developmental origin[END_REF][START_REF] Carter | The evolution of anthropoid molar proportions[END_REF][START_REF] Evans | A simple rule governs the evolution and development of hominin tooth size[END_REF][START_REF] Roseman | The Inhibitory Cascade Model is Not a Good Predictor of Molar Size Covariation[END_REF]. However, the existence of potential allometric patterns within the large diversity of molar proportions in placentals has rarely been addressed. Two studies suggested that lower molar proportions may vary jointly with body size in a few isolated placental clades, but this was never investigated at a large phylogenetic scale within placentals [START_REF] Wilson | Testing a developmental model in the fossil record: molar proportions in South American ungulates[END_REF][START_REF] Evans | A simple rule governs the evolution and development of hominin tooth size[END_REF]).

Here, we explore upper molar proportions in a large and phylogenetically diverse sample of placental mammals and test for an allometric pattern accompanying this diversity.

By focusing on the upper molars, we also assess how well the IC model, which was originally proposed on lower molars, fits variation observed in the upper jaw. Our results highlight a previously overlooked evolutionary link between the absolute size of molars and their proportions that adds to our knowledge of the sequential competition faced by adjacent teeth [START_REF] Kavanagh | Predicting evolutionary patterns of mammalian teeth from development[END_REF]. The discovery of this allometric gradient in the tooth row strongly emphasizes the value of exploring scaling patterns among other segmental series within the body of vertebrates and other living organisms.

MATERIAL & METHODS

Specimen Selection

We measured the size of the three upper molars (M1, M2, M3) in complete molar rows of 311 extant and fossil specimens belonging to 199 extant and 100 extinct placental species, in 277 genera. 286 species were represented by one specimen, and 13 species by 2 specimens. The intraspecific variation of molar proportions, generally low in placental species [START_REF] Vitek | Mammal molar size ratios and the inhibitory cascade at the intraspecific scale[END_REF], was assumed to be negligible in comparison to the magnitude of the interspecific variation related to allometry for Placentalia (see Supplementary Material available on Dryad, http://dx.doi.org/10.5061/dryad.dfn2z34xq). The sample covered most extant placental orders, some being excluded for being edentulous or enamel-free taxa (Cingulata, Pilosa, Tubulidentata), having high-crowned teeth (Lagomorpha) (see Gomes [START_REF] Rodrigues | Ontogenetic variations and structural adjustments in mammals evolving prolonged to continuous dental growth[END_REF] for reasoning), or supernumerary molars (Sirenia). The sample was mostly composed of original specimens or casts present in the collections of the Muséum national d'Histoire naturelle, Paris (MNHN), while images for the remaining specimens were collected from the literature (see Supplementary Material available on Dryad).

Measurements

Molar size was quantified by measuring area in occlusal view using the software Fiji [START_REF] Schindelin | Fiji: an open-source platform for biological-image analysis[END_REF]) on scaled images of upper molar rows. For MNHN specimens (original or cast), images of molar occlusal surfaces were taken with a Canon® EOS 60D camera fitted with a Canon® EF Macro-Lens 100 mm, f/2.8. The camera was attached to a Giottos® camera tripod (VT-901). The occlusal surface and the lens of the camera were oriented parallel to each other, in the horizontal plane (controlled by a leveling device). A scale was placed parallel to, and at the same distance from the lens as, the basal plane of the crown (i.e., cementoenamel junction), as in previous studies [START_REF] Gómez-Robles | A geometric morphometric analysis of hominin upper first molar shape[END_REF].

Measurements for non-MNHN specimens were made on scaled illustrations available in the literature. On the scaled images, molar outlines were delineated at the cementoenamel junction using the tool Freehand Selection in Fiji. The area of each molar was measured three times following this method and then averaged. Three variables were obtained: M1area, M2area, and M3area. From these, three variables were computed for further analyses: M2area/M1area; M3area/M1area; log (size) = log (M1area + M2area + M3area). The latter variable thus measured the absolute size of the molar field (total molar area). The three variables are noted M2/M1, M3/M1 and log (size) hereafter.

In some taxa, molars experienced some degree of labial overlap that partially hid small parts of the parastylar region or mesial edge of the tooth. In such case, the molar outline in this region was drawn by following a straight path between the non-hidden parts.

Tree building

A fully dichotomous composite cladogram was constructed for phylogenetically-informed statistical analyses (see below). It is largely overlapping with that of [START_REF] Billet | Serial homology and correlated characters in morphological phylogenetics: modeling the evolution of dental crests in placentals[END_REF] and its construction is based on the same publications and rationale. We used the phylogenetic patterns obtained by [START_REF] Meredith | Impacts of the Cretaceous Terrestrial Revolution and KPg extinction on mammal diversification[END_REF] and [START_REF] Wible | The eutherian mammal Maelestes gobiensis from the Late Cretaceous of Mongolia and the phylogeny of Cretaceous Eutheria[END_REF] for interordinal relationships of extant taxa and for the placement of successive fossil outgroups to crown Placentalia, respectively. We used a combination of sources for dating the divergences on the reconstructed composite cladogram and followed the same rationale as [START_REF] Billet | Serial homology and correlated characters in morphological phylogenetics: modeling the evolution of dental crests in placentals[END_REF] for placing fossils (see Supplementary Material available on Dryad).

Statistical analyses

Statistical analyses were performed on the complete sample (Placentalia), as well as on selected clades within Placentalia, with n ≥ 15 sampled species. The relationships between molar proportions, i.e. M2/M1 and M3/M1, were first analyzed using both ordinary Reduced Major Axis (RMA) and Phylogenetic Reduced Major Axis (PRMA) regressions respectively with the R packages rgr [START_REF] Garrett | The 'rgr' package for the R Open Source statistical computing and graphics environment -a tool to support geochemical data interpretation[END_REF] and phytools [START_REF] Revell | phytools: An R package for phylogenetic comparative biology (and other things)[END_REF]. RMA and PRMA were used as none of the variables are known to have a direct influence on each other. Confidence of RMA is estimated with 95% confidence intervals, p-values of f-tests and Pearson correlation coefficient. PRMA confidence is based on the same metrics except for confidence intervals that are not available.

The influence of the size of the molar field (log (size)) on M2/M1 and M3/M1 was quantified using both Generalized Least Squares (GLS) and Phylogenetic GLS (PGLS) regressions with the R package caper (Orme et al. 2013). Using maximum likelihood, we simultaneously optimized Pagel's lambda to integrate a variable phylogenetic signal. GLS and PGLS robustness was estimated with standard error of estimates, p-values of corresponding t-tests and R²; their relative adequacy to the data was compared using the Akaike Informative Criterion (AIC; [START_REF] Akaike | A new look at the statistical model identification[END_REF]). It should be noted that the R² of a PGLS is hardly comparable with that of an Ordinary Least Squares regression as their covariance matrices are respectively the covariance matrix derived from the model of trait change and the phylogeny and the identity matrix. In order to visualize the correlated evolution of log (size) and molar proportions on the cladogram, we used the function contMap of phytools [START_REF] Revell | phytools: An R package for phylogenetic comparative biology (and other things)[END_REF]) that optimizes ancestral node values that are the most likely under a Brownian motion evolution.

Finally, to investigate the variation of the regression outcomes due to phylogenetic uncertainty in placentals, we performed PGLS on altered topologies. Using phangorn [START_REF] Schliep | phangorn: Phylogenetic analysis in R[END_REF], we ran 100 times 50 topological changes by Nearest Neighbour Interchange rearrangement.

RESULTS

The distribution of specimens in the morphospace is concentrated mainly in one direction, suggesting that molar sizes largely covary (Fig. 1a), but this direction is shifted towards larger M2 values relative to the proposed IC model [START_REF] Kavanagh | Predicting evolutionary patterns of mammalian teeth from development[END_REF][START_REF] Polly | Development with a bite[END_REF]) and shows a weaker slope than the latter (see Supplementary Material available on Dryad). In the morphospace of molar proportions, taxa with the largest molar field generally have relatively large posterior molars (M1<M2<M3 and M1<M2>M3) while smaller taxa show a relatively larger M1 (M1>M2>M3) (Fig. 1a).

The phylogenetic co-distributions of the molar proportions and of the absolute size of the molar field (log (size)) confirm this pattern: when the absolute size of the molar field increases in a clade, both ratios M3/M1 and M2/M1 most often increase as well (Fig. 2 and see Supplementary Material available on Dryad). These co-distributions are characterized by multiple events of increases and decreases of M2/M1, M3/M1 and log (size) during placentals' evolutionary history. Generalized Least Square analyses, whether phylogenetically informed (PGLS) or not (GLS), show a strong statistical support for correlated evolution of the absolute size of the molar field and molar proportions (Table 1). Bivariate plots of the proportions and size variables (Fig. 1b-c) show that specimens belonging to various clades cluster in the extremal parts of the morphospace, and thus discard the possibility of a Felsenstein's worst case scenario (singular evolutionary events; [START_REF] Uyeda | Rethinking phylogenetic comparative methods[END_REF]) at the Placentalia level. The Akaike scores indicate that the favored regression model for both ratios is the PGLS one (Table 1). The slope of the regression of M3/M1 on log (size) is more than twice the slope of the regression of M2/M1 on log (size), in both GLS and PGLS analyses (Fig. 1b-c & Table 1). This shows that M3/M1 increases more than twice as fast as M2/M1 with an increasing log (size). The intercept values, one of which has poor statistical support (PGLS, M3/M1 vs log (size)), are much higher for the regressions of M2/M1 on log (size) than for those of M3/M1 on log (size). Besides, the slope values are lower in the PGLS than in the GLS analyses of the same variables while it is the opposite for the intercept values, indicating that the consideration of phylogenetic relationships pulls the regression line higher and makes it less steep (Fig. 1b-c & Table 1). At a lower taxonomic level, the Afrotheria, Primates and Glires also show a well-supported correlated evolution for the ratios M2/M1 & M3/M1 with the absolute size of the molar field, while Chiroptera show an opposite trend for M3/M1 (see further comments on variations among clades in Supplementary Material available on Dryad).

The analyses on altered topologies also show a robust statistical support, for both regressions of M2/M1 on log (size) and M3/M1 on log (size). Both the slope and intercept values of the randomized analyses remain very close to the ones found by the PGLS analyses on the baseline cladogram (Table 1 and see Supplementary Material available on Dryad). The distribution of the p-values of t-tests are very low except for intercepts of a few runs for the regression of M3/M1 on log (size). The results of the PGLS analyses on the entire Placentalia sample are thus not prone to change when the phylogenetic relationships are slightly changed.

DISCUSSION

Evolutionary Allometry of Molar Proportions: Robustness and Scope

The large diversity of upper molar proportions within placentals is accompanied by an allometric pattern along their favored direction of covariation: placental mammals with large dentitions present larger posterior molars relative to the anterior one. The area of M2 and even more that of M3 increases faster than that of M1 with increasing size at the macroevolutionary level (Figs. 1b-c and see Supplementary Material available on Dryad). This is congruent with the well-known high variability in size of third molars, which are the last to develop (Gould and Garwood 1966;[START_REF] Bernal | Developmental processes, evolvability, and dental diversification of New World monkeys[END_REF]). An allometric pattern of molar proportions was very strongly supported at the level of Placentalia, and well-supported in clades like Primates and Afrotheria. The taxonomic scale for such analyses may be important as the evolution of several independent and large-magnitude changes in size might be more likely in longexisting and diverse lineages. Clades with a lower taxonomic diversity and/or lower magnitude of size differences in our sample often show no (or less) well-supported allometric relationship (see Supplementary Material available on Dryad). Our results seem robust to phylogenetic variations, which is important given that placental phylogeny especially suffers from uncertainties concerning the position of fossil lineages [START_REF] Springer | Technical comment on "The placental mammal ancestor and the post-K-Pg radiation of placentals[END_REF][START_REF] Beck | Improvements in the fossil record may largely resolve current conflicts between morphological and molecular estimates of mammal phylogeny[END_REF].

Our study shows that variation in upper molar proportions across placentals deviates from the line predicted by the IC model, and maybe even slightly more than for lower molars (e.g., [START_REF] Kavanagh | Predicting evolutionary patterns of mammalian teeth from development[END_REF][START_REF] Bernal | Developmental processes, evolvability, and dental diversification of New World monkeys[END_REF]. Antagonist molars may partly differ in their development (e.g., [START_REF] Hayden | Developmental variability channels mouse molar evolution[END_REF] and probably experience different spatial constraints on the palate compared to the narrow lower jaw [START_REF] Renvoisé | Mechanical constraint from growing jaw facilitates mammalian dental diversity[END_REF]. Whether these factors could cause differences in proportional relationships among the upper and lower molar rows remains to be evaluated. Despite these differences, antagonist molars are often strongly integrated in both shape and size [START_REF] Renaud | Developmental constraints revealed by co-variation within and among molar rows in two murine rodents[END_REF][START_REF] Gómez-Robles | Morphological integration in the hominin dentition: evolutionary, developmental, and functional factors: morphological integration in the hominin dentition[END_REF] and molar proportions in the upper jaw often resemble those of the lower jaw [START_REF] Rodrigues | Ontogenetic variations and structural adjustments in mammals evolving prolonged to continuous dental growth[END_REF][START_REF] De Muizon | New remains of kollpaniine "condylarths" (Panameriungulata) from the early Palaeocene of Bolivia shed light on hypocone origins and molar proportions among ungulate-like placentals[END_REF]. The genetic determinants of molar size variation also appear partially similar for the upper and lower jaws (Navarro & Maga, 2019). For all these reasons, we postulate that the evolutionary allometry detected here on upper molars may be expressed similarly, but not identically, on lower molars in placentals, which echoes previous findings in isolated placental clades [START_REF] Wilson | Testing a developmental model in the fossil record: molar proportions in South American ungulates[END_REF][START_REF] Bernal | Developmental processes, evolvability, and dental diversification of New World monkeys[END_REF][START_REF] Evans | A simple rule governs the evolution and development of hominin tooth size[END_REF].

While choosing a local measure of size (molar field), we assumed that molar proportions were more closely linked to local processes than to those determining variation in overall body size or size of other organs [START_REF] Hallgrímsson | Integration and the developmental genetics of allometry[END_REF]). Yet, the absolute size of molars is also roughly indicative of body sizes in mammals (e.g., Copes and Schwarz 2010).

When combined with our results, this suggests that large-bodied placentals will tend to show much larger M3 and M2 than M1 and the opposite for small-bodied ones. Most extant placental taxa that lost the upper and/or lower third molar count as small-to medium-sized organisms (e.g., some macroscelids, chrysochlorids, bats, callitrichines monkeys, otochonid lagomorphs, carnivorans, murid and bathyergid rodents; taxa with dental displacement and ever-growing teeth not considered; [START_REF] Ungar | Mammal teeth: origin, evolution, and diversity[END_REF][START_REF] Rodrigues | The great disparity of dental structures and dynamics in rodents: new insights into their ecological diversity[END_REF]. This seems in good agreement with our results, because if small mammals tend to have relatively reduced posterior molars, they could be more prone to a loss of their upper and lower third molars.

Similarly, cases of small placentals with much enlarged posterior molars relative to M1 should be rare. Apart from taxa with reduced or homodont dentitions, major exceptions to the allometric pattern are found within rhinoceroses, which have relatively small posterior molars while having a large absolute size, and in the large-sized carnivorans presenting relatively large M1, and small-to-absent M2 and M3 (e.g., felids, ursids). Stem and early diverging carnivorans with reduced-to-absent M3 and M2 were however rather small-to medium-sized placentals [START_REF] Bokma | Testing for Depéret's rule (body size increase) in mammals using combined extinct and extant data[END_REF], as predicted by the allometric pattern. Later increased specialization towards a flesh-eating diet and the anterior position of the carnassials might have caused the deviation of late diverging carnivorans from this pattern (see below).

Although it is robustly supported at the placental level, the strength and direction of the allometric pattern vary within placental groups and exceptions exist. This suggests that size may represent one among several determinants of molar proportions, each one having a fluctuating influence along evolutionary branches.

Potential Correlates and Mechanisms behind the Allometric Pattern

Our work suggests a biased macroevolution of molar proportions towards relatively larger posterior molars in large-sized taxa. Potential mechanisms for this phenomenon could implicate a weaker intermolar inhibition with increasing absolute size of the molar field and/or a higher dental mesenchymal activation [START_REF] Kavanagh | Predicting evolutionary patterns of mammalian teeth from development[END_REF][START_REF] Wilson | Testing a developmental model in the fossil record: molar proportions in South American ungulates[END_REF][START_REF] Evans | A simple rule governs the evolution and development of hominin tooth size[END_REF]. Overall, this would mean a higher ratio of activator to inhibitor in taxa with a larger molar field but it is unclear if the mouse-derived IC model applies to organisms with very different life histories (e.g., larger and longer-lived organisms; Speakman 2005; Carter and Worthington 2016). Another question that arises from this is whether the allometric pattern of molar proportions detected here at the macroevolutionary scale is also expressed at a microevolutionary one, and how. Further research on the size-related variation of molar proportions will be needed within species and along ontogeny to test for expression of this allometric patterning across biological levels (Klingenberg and Zimmerman 1992;[START_REF] Pavlicev | Genetic variation in pleiotropy: differential epistasis as a source of variation in the allometric relationship between long bone lengths and body weight[END_REF][START_REF] Hlusko | Modularity in the mammalian dentition: mice and monkeys share a common dental genetic architecture[END_REF][START_REF] Navarro | Genetic mapping of molar size relations identifies inhibitory locus for third molars in mice[END_REF]. Among other aspects, this allometric pattern may be worth investigating in relation to the genetic correlation between body size and molar width previously detected in baboons [START_REF] Hlusko | Evolution of genetically correlated traits: tooth size and body size in baboons[END_REF]. At the macroevolutionary level, our study suggests a more complex pattern of integration among molars than previously thought, with molars covarying in size with each other, and their proportions covarying with the absolute size of the entire molar field. Like any pattern of covariation, this twofold pattern might involve complex combinations of developmental processes, genetic and environmental effects [START_REF] Hallgrímsson | Deciphering the palimpsest: studying the relationship between morphological integration and phenotypic covariation[END_REF][START_REF] Roseman | The Inhibitory Cascade Model is Not a Good Predictor of Molar Size Covariation[END_REF].

Potential correlates of molar proportions and thus of the detected evolutionary allometry include the configuration of neighboring premolars and the dietary habits of the taxa. In rodents, molar proportions were linked to the presence or absence of the fourth premolar [START_REF] Labonne | When less means more: evolutionary and developmental hypotheses in rodent molars[END_REF]) which might constitute a spatial constraint on the development of the first molar. The fourth permanent premolar is absent in only one clade composed of four rodents in our sample and thus was not an important parameter. The proportions of posterior premolars and their relation to the allometry shown by molars remain to be investigated, although, for instance, species with large posterior premolars do not seem to be associated with a given range of molar proportions in our sample (e.g., rhinos vs pyrotheres).

The timing of premolar initiation and replacement could also be an important parameter in relation to the space available for molar formation in the jaw (see below).

Molar proportions were reported to reflect dietary habits in various groups of placental mammals [START_REF] Kavanagh | Predicting evolutionary patterns of mammalian teeth from development[END_REF][START_REF] Halliday | Testing the inhibitory cascade model in Mesozoic and Cenozoic mammaliaforms[END_REF][START_REF] Carter | The evolution of anthropoid molar proportions[END_REF] but see also Monson et al., 2019a). Herbivorous and faunivorous forms were then found at opposite extremes of the morphospace with the former showing relatively large posterior molars (e.g., upper right corner of Fig. 1a). Diet might also be correlated to body mass with large terrestrial mammals being mostly herbivorous and small ones mostly omnivorous and carnivorous [START_REF] Price | The macroevolutionary relationship between diet and body mass across mammals[END_REF]. It is therefore not surprising that both diet and size may be reflected in the diversity of molar proportions among placentals. Interestingly, the extinct hyainalourine creodonts included giant carnivorous forms with relatively large posterior molars [START_REF] Borths | Simbakubwa kutokaafrika, gen. et sp. nov. (Hyainailourinae, Hyaenodonta,'Creodonta,'Mammalia), a gigantic carnivore from the earliest Miocene of Kenya[END_REF], which agrees better with our results and contrasts with the case of carnivorans. Besides allometry, the chewing biomechanics and the various positions of the carnassial teeth in carnivorous forms [START_REF] Asahara | Unique pattern of dietary adaptation in the dentition of Carnivora: its advantage and developmental origin[END_REF]) may play a role in molar proportions.

One important aspect of this study is that it may imply differential growth dynamics of the anterior and posterior molars depending on the size of the molar field. It is currently unclear whether and how a different timing of initiation and/or different growth rates of posterior molars could explain diverging molar proportions [START_REF] Kavanagh | Predicting evolutionary patterns of mammalian teeth from development[END_REF][START_REF] Carter | The evolution of anthropoid molar proportions[END_REF], and even more so in relation to size. With respect to their skull, mammals generally present well-marked differential growth between their face and neurocranium, known as craniofacial allometry. This allometry is present both at the evolutionary and ontogenetic levels in mammals with larger individuals having relatively larger faces than smaller ones [START_REF] Cardini | Larger mammals have longer faces because of size-related constraints on skull form[END_REF][START_REF] Cardini | Craniofacial allometry is a rule in evolutionary radiations of placentals[END_REF]. Large mammals could therefore be characterized by both relatively large faces and large posterior molars, which raises the question of shared mechanisms between these two allometric patterns. Upper molars grow sequentially, with the M3 forming last, in the posterior region of the palate formed by the maxillary bone. As a main facial component, the maxillary bone strongly increases in size during mammalian ontogeny and is much affected by the evolutionary craniofacial allometry [START_REF] Cardini | Larger mammals have longer faces because of size-related constraints on skull form[END_REF][START_REF] Cardini | Craniofacial allometry is a rule in evolutionary radiations of placentals[END_REF]. The forming molars and surrounding jaw are tightly integrated during development, with mechanical aspects of the growing jaw actively constraining the shape of molars [START_REF] Renvoisé | Mechanical constraint from growing jaw facilitates mammalian dental diversity[END_REF]) and with dynamic bone remodeling around forming molars [START_REF] Butler | The ontogeny of molar pattern[END_REF][START_REF] Lungová | Tooth-bone morphogenesis during postnatal stages of mouse first molar development[END_REF]. Knowing the spatiotemporal association between jaw and molar formations and the late completion of the posterior part of the palate in its growth sequence, one could speculate that jaws growing an overall larger size might be developmentally and mechanistically associated with a greater growth of lateforming posterior molars in comparison to anterior ones. Speculation about an association between molar proportions and jaw size is not new (e.g., [START_REF] Ford | Callitrichids as phyletic dwarfs, and the place of the Callitrichidae in Platyrrhini[END_REF]Monson et al., 2019b), and links between jaw dimensions and third molar development are for instance welldocumented in humans [START_REF] Kömerik | Skeletal facial morphology and third molar agenesis[END_REF][START_REF] Marchiori | Initial third molar development is delayed in jaws with short distal space: An early impaction sign?[END_REF]. However, more research is needed in that direction as the validity of an association between jaw growth and molar proportions and the potential mechanisms behind it remains extremely obscure.

Overall, the allometric relationship found in our study is congruent with the idea that processes regulating the size of organs, such as individual molars, and overall patterns of growth are integrated (e.g., [START_REF] Shingleton | Evolution and the regulation of growth and body size[END_REF]. This finding adds complexity to our Afrotheria; Art., Artiodactyla; Car., Carnivora; Chi., Chiroptera; Eua., Euarchontoglires; Eul., Eulipotyphla; Lau., Laurasiatheria; P., Perissodactyla; Pri., Primates; Rod., Rodentia. Silhouettes from phylopic.org TABLE 1. Results of GLS and PGLS regressions with molar proportions as the dependent variables and the size of the molar field (log (size)) as the independent variable in Placentalia.

  understanding of the covariation of molar proportions and prompts further questions about the mechanisms generating this structured diversity. The current lack of mechanistic explanations should stimulate further research along with further testing of allometric patterns among other segmental systems of living organisms.

FIGURE 2 .

 2 FIGURE 2. Phylogenetic co-distribution of molar proportions (M3/M1) and size of the molar

ACKNOWLEDGMENTS

SUPPLEMENTARY MATERIAL

Supplementary material, including data files and/or online-only appendices, can be found in the Dryad data repository at http://datadryad.org, http://dx.doi.org/10.5061/dryad.dfn2z34xq. 1). Abbreviation: SANUs = South American native ungulates.

FIGURES & TABLES