Richard Ríos

Renaud De Crevoisier

Juan D Ospina

Frédéric Commandeur

Caroline Lafond

Antoine Simon

Pascal Haigron

Jairo Espinosa

Oscar Acosta

Richard Rios

Frederic Commandeur

Population model of bladder motion and deformation based on dominant eigenmodes and mixed-effects models in prostate cancer radiotherapy

Keywords: principal component analysis, motion/deformation models, mixed-effects models, prostate cancer radiotherapy

A C C E P T E D M A N U S C R I P T

Highlights

• We propose a population model to estimate the probability of the bladder presence on a given region during treatment using only the planning CT scan as input information.

• We train a motion/deformation model, based on longitudinal data, to predict bladder motion and deformation between fractions.

• We propose a longitudinal analysis using mixed-effect models to separate intra-and inter-patient variability in order to control confounding.

• We reduce, in a factor of 10, the number of variables required to represent bladder surface using spherical harmonics (SPHARM).

A C C E P T E D M A N U S C R I P T

Prostate cancer is one the most common cancer types worldwide. One of the standard treatment methods is external radiotherapy, which involves delivering ionizing radiation to a clinical target, namely the prostate and seminal vesicles. Due to the uncertain location of organs during treatment, consisting around 40 radiation fractions delivering a total dose ranging from 70 to 80Gy, safety margins are defined around the tumor target upon treatment planning. This also leads to portions of healthy organs neighboring the prostate or organs at risk -the bladder and rectum -to be included in the target volume, potentially resulting in adverse events affecting patients'urinary (hematuria and cystitis, among others) or rectal (rectal bleeding, fecal incontinence, etc.) functions. Several studies have shown that increasing dose delivery to the prostate leads to improved local cancer control, up to approximately 80Gy [START_REF] Zietman | Randomized trial comparing conventional-dose with highdose conformal radiation therapy in early-stage adenocarcinoma of the prostate: long-term results from proton radiation oncology group/american college of radiology 95-09[END_REF]. However, such dose increases are limited by their associated risks of treatment-related toxicity involving the organs at risk [START_REF] Fiorino | Dose volume effects for normal tissues in external radiotherapy: Pelvis[END_REF][START_REF] Fonteyne | Intensity-modulated radiotherapy as primary therapy for prostate cancer: Report on acute toxicity after dose escalation with simultaneous integrated boost to intraprostatic lesion[END_REF].

The bladder is particular for presenting the largest inter-fraction shape variations during treatment, caused by continuous volume changes [START_REF] Roeske | Evaluation of changes in the size and A C C E P T E D M A N U S C R I P T location of the prostate, seminal vesicles, bladder, and rectum during a course of external beam radiation therapy[END_REF][START_REF] Turner | Bladder movement during radiation therapy for bladder cancer: implications for treatment planning[END_REF]. These shape variations introduce geometric uncertainties that render assessment of the actual dose delivered to the bladder during treatment difficult, thereby leading to dose uncertainties that limit the possibility of modeling dose-volume response for late genitourinary (GU) toxicity [START_REF] Fiorino | Dose volume effects for normal tissues in external radiotherapy: Pelvis[END_REF]Thor et al., 2013b;[START_REF] Nassef | Quantification of dose uncertainties in cumulated dose estimation compared to planned dose in prostate IMRT[END_REF]Palorini et al., 2016a;[START_REF] Landoni | Predicting toxicity in radiotherapy for prostate cancer[END_REF]. The Quantitative Analysis of Normal Tissue Effects in the Clinic (QUAN-TEC) project has stated that a similar dose-response to that of late gastrointestinal (GI) toxicity is far from being established. The dosimetric variables obtained from the planning CT prove to be very poor surrogates for the real-delivered dose

A C C E P T E D M

A N U S C R I P T [START_REF] Viswanathan | Radiation dose-volume effects of the urinary bladder[END_REF][START_REF] Rosewall | The relationship between external beam radiotherapy dose and chronic urinary dysfunction a methodological critique[END_REF]Palorini et al., 2016b). As a result, it appears crucial to quantify uncertainties produced by inter-fraction bladder variations in order to determine dosimetric factors that affect late GU complications.

Several approaches have been developed in order to quantify and characterize the geometrical uncertainties produced by organ motion and deformation between fractions, including: serial imaging measurement of the organ during treatment course [START_REF] Roeske | Evaluation of changes in the size and A C C E P T E D M A N U S C R I P T location of the prostate, seminal vesicles, bladder, and rectum during a course of external beam radiation therapy[END_REF][START_REF] Van Herk | Quantification of organ motion during conformal radiotherapy of the prostate by three dimensional image registration[END_REF]Ten Haken et al., 1991), fiducial markers [START_REF] Crook | Prostate motion during standard radiotherapy as assessed by fiducial markers[END_REF][START_REF] Balter | Measurement of prostate movement over the course of routine radiotherapy using implanted markers[END_REF], margins of organs at risk (OAR) (Thor et al., 2010;[START_REF] Mageras | Computerized design of target margins for treatment uncertainties in conformal radiotherapy[END_REF][START_REF] Stroom | Inclusion of geometrical uncertainties in radiotherapy treatment planning by means of coverage probability[END_REF], rigid-body motion (Thor et al., 2013a;[START_REF] Killoran | A numerical simulation of organ motion and daily setup uncertainties: Implications for radiation therapy[END_REF][START_REF] Craig | Considerations for the implementation of target volume protocols in radiation therapy[END_REF], parametrization of the organ structure [START_REF] Hoogeman | A model to simulate day-to-day variations in rectum shape[END_REF][START_REF] Mageras | A method of incorporating organ motion uncertainties into three-dimensional conformal treatment plans[END_REF][START_REF] Pavel-Mititean | A geometric model for evaluating the effects of inter-fraction rectal motion during prostate radiotherapy[END_REF], as well as biomechanical [START_REF] Yan | A model to accumulate fractionated dose in a deforming organ[END_REF] and statistical models (Fontenla et al., 2001a,b;[START_REF] Söhn | Dosimetric treatment course simulation based on a statistical model of deformable organ motion[END_REF][START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF]Thörnqvist et al., 2013;[START_REF] Bondar | Statistical modeling of CTV motion anddeformation for IMRT of earlystage rectalcancer[END_REF][START_REF] Rios | Characterization of bladder motion and deformation in prostate cancer radiotherapy[END_REF]. Recently, a method based on weighted scenarios of the fundamental directions of the patient's geometric variability was implemented [START_REF] Söhn | Modelling individual geometric variation based on dominant eigenmodes of organ deformation: implementation and evaluation[END_REF][START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF]. These directions of variability, called motion/deformation modes, were obtained by applying principal component analysis (PCA) to a data set of (pre-)treatment organ geometries, usually parameterized by a set of corresponding surface points. This approach known as a point-distribution model (PDM) was introduced in statistical shape modelling and organ segmentation [START_REF] Cootes | Active Shape Models-Their Training and Application[END_REF][START_REF] Lorenz | Generation of point-based 3d statistical shape models for anatomical objects[END_REF][START_REF] Davatzikos | A framework for predictive modeling of anatomical deformations[END_REF][START_REF] Pekar | Automated model-based organ delineation for radiotherapy planning in prostatic region[END_REF].

PCA-based motion/deformation models were initially used by [START_REF] Söhn | Modelling individual geometric variation based on dominant eigenmodes of organ deformation: implementation and evaluation[END_REF] for modeling individual geometric variations between fractions of the rectum, bladder, and prostate. The main goal was to estimate the patient's variance along each mode of geometric variability. These modes were derived by diagonalization of a data matrix, called the covariance matrix, whose columns were made of patient's motion/deformation vectors (i.e., vectors that describe the observed geometric variation of the organ at the on-treatment CTs/cone-beam CTs (CBCTs) relative to the mean shape). New organ geometries were then generated by adding a weighted sum of a few dominating modes to the mean organ shape, where each weight obeyed a Gaussian distribution with mean of zero and the corresponding eigenvalue as variance. PCA has also been applied to generate individual 4D statistical models for organs undergoing respiratory motion like the lung, liver, and heart (Zhang et al., 2007a;[START_REF] Badawi | Optimizing principal component models for representing interfraction variation in lung cancer radiotherapy[END_REF][START_REF] Mcclelland | Respiratory motion models: A review[END_REF]. The same method was in turn applied to describe the geometric variations of organs like prostate, seminal vesicles, and pelvic lymph nodes (Thörnqvist et al., 2013).

A C C E P T E D M A N U S C R I P T

However, this individual-based methodology was limited by the required number of CT scans per patient for training the model. In clinical practice, only one CT scan is available at the planning treatment for a typical patient, and on-treatment CTs/CBCTs are seldom available during treatment. [START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF] thus proposed a population-based method that learned the geometric variations from a population database seeking to infer the possible organ variations for a typical patient despite his small number of CTs/CBCTs. [START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF] applied his method to describe clinical target volume (CTV) motion and deformation between fraction in prostate cancer using a longitudinal database. It is worth mentioning that longitudinal data set consists of repeated observations of a set of homologous objects or variables at several time points (Durrleman et al., 2013a;[START_REF] Pinheiro | Mixed-Effects Models in S and S-PLUS, 1st Edition[END_REF]. In our context, it corresponds to repeated observations of the organ at different fractions for a patient population. This methodology was also applied to model motion and deformation of the CTVs and bladder in rectal and prostate cancer radiotherapy, respectively [START_REF] Bondar | Statistical modeling of CTV motion anddeformation for IMRT of earlystage rectalcancer[END_REF][START_REF] Rios | Characterization of bladder motion and deformation in prostate cancer radiotherapy[END_REF]. In addition, [START_REF] Hu | Population-based prediction of subject-specific prostate deformation for MR-to-ultrasound image registration[END_REF] proposed a population-PCA model to describe prostate deformation in magnetic resonance (MR)-tumor-targeted biopsies using a longitudinal database of prostates obtained from MR images and biomechanical models. Similarly, several researchers have introduced so-called 4D populationbased or cross-population models for organs undergoing respiratory motion [START_REF] He | Online 4-D CT estimation for patient-specific respiratory motion based on real-time breathing signals[END_REF][START_REF] Mcclelland | Respiratory motion models: A review[END_REF][START_REF] Preiswerk | Model-guided respiratory organ motion prediction of the liver from 2d ultrasound[END_REF]. In population-based models, modes were derived from the population covariance matrix made by the patients' motion/deformation vectors, i.e., all the rigid and non-rigid organ displacements observed for each patient around his mean bladder were stacked as columns. For a given patient, new organ samples were also generated by a weighted sum of a few dominating modes to the mean shape, whereby each weight also obeyed a Gaussian distribution with mean of zero and the corresponding eigenvalue as variance (i.e., the total population variance). However, this may lead to a suboptimal estimation of individual-specific variability, given that that the mode variance is considered to be shared by all the individuals. This is an assumption that is not necessarily verified as the patients may exhibit heterogeneous organ variability.

In this paper, we sought to properly characterize patient' s variance along each mode of geometric variability from the population variance. Next, we have proposed a longitudinal analysis using mixed-effects (ME) models to separate intraand inter-patient variability along each mode in order to control confounding effects. We have also proposed an early step of dimensionality reduction using a SPHARM parameterization followed by PCA. We have thus found the modes of

A C C E P T E D M A N U S C R I P T
geometric variation in the SPHARM coordinates rather than in the surface points as PDM-based shape models. This step enabled us to reduce the number of parameters required to represent the bladder surface, further encoded with the dominant modes derived by PCA. This step also allowed us to further reduce data storage and processing time, among others aspects.

Our main contribution regarding groundwork was to introduce ME models in order to estimate individual-specific variations between fractions from a reducedorder model (SPHARM/PCA based model). Using the ME model the following aspects were addressed for each mode: i) the estimation and separation of individual-specific variance from that of the population, i.e., the inherent correlation of intra-individual observations; ii) the reduction of the observed population variance by grouping repeated observations per patient by means of a "patient factor"; iii) the characterization of the individual-specific directions of motion and deformation from the population; iv) the study of longitudinal observations ordered in time, which may reflect temporal trajectories in the latent space. This study is also the first explicitly addressing the issue of predicting a region of probable inter-fraction bladder motion/deformation using solely the planning CT. A quantitative comparison of predicted vs observed probability maps was proposed using our model and other previous motion/deformation models.

Our proposed model is based on the assumptions stated in the methodology developed by [START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF]: firstly, that the bladder moves and deforms in a limited number of directions, constrained by the body's anatomy; secondly, that the inter-fraction geometric organ variations are randomly distributed along the set of variability directions; thirdly, that despite the potential variations in organ size and shape across the population, principal directions of geometric variation prove to be similar in all patients; fourthly, that geometric variations between fractions are independent random variables, where the time sequence of the observations does not matter. Hence, each bladder within the population can be characterized by its variation along each direction.

The paper is organized as follows. Section 1 describes the two databases (training and testing) included in this study. Section 2 provides a detailed description of the proposed method together with the validation framework. Section 3 presents the results and comparison to groundwork methods. Lastly, Section 4 discusses the experimental findings, with a conclusion given in Section 5 in addition to final considerations and future perspectives.

A C C E P T E D M

A N U S C R I P T

Data

We included two population databases of patients treated for prostate cancer with external beam radiotherapy (RT). We trained the model using data from a database of 20 patients, which included a planning CT and several on-treatment CTs (5-8) for each patient. We also used an independent validation database of 28 patients described as follows: 25 patients had a planning CT as well as 6-9 ontreatment CTs, and 3 patients with a planning CT as well as 35-39 on-treatment CBCTs for each. The contours of the bladder, rectum, prostate and seminal vesicles (SV) were manually delineated following the same protocol. The prescribed dose was computed in a standard treatment planning system (TPS) step then resampled into the CT native space. We used 135 × 215 × 55 CTs and CBCTs with 1 mm × 1 mm × 2 mm voxels. From the delineated CT and CBCT scans, observedprobability maps called PM obs were computed for comparison with the prediction obtained using the proposed method, as explained in the validation section.

Methods

Population model training

The training steps are depicted in Fig. 1. Firstly, rigid-spatial normalization was performed in order to anatomically align the pelvic structures (prostate and bladder); secondly, bladder surface parameterization was conducted using SPHARM, followed by dimensionality reduction using PCA, thereby yielding a limited set of directions of maximum geometric variability. These directions, or modes, were validated by leave-one-out cross validation, aiming to determine whether patients exhibited common directions of bladder motion and deformation; finally, a ME model was fit to each direction to model the inter and intra-individual geometric variations obtained.

Prostate-based rigid registration

This was done to spatially align all the bladders in the same spatial referential. After defining the prostate barycenter of a patient as common coordinate system, henceforth called the template, all the patient structures were thus registered by aligning the prostate barycenters, as in [START_REF] Rios | Characterization of bladder motion and deformation in prostate cancer radiotherapy[END_REF], thus yielding a transformation T . The template was a patient from the training database whose prostate exhibited the lowest volume variation and the most similar mean volume to the population mean. The individual bladders were then propagated to the common template using the obtained transformation T . This enabled bladder motion and deformation to be quantified with respect to the template prostate barycenter. Fig. 2 depicts an example of the performed rigid registration/propagation step.

Spherical harmonics

An essential objective in shape modeling has been the reduction of parameters needed to describe a shape instance. As a first step towards this goal, PDM method was proposed to parameterize 3D shapes by means of a set of boundary points [START_REF] Cootes | Active Shape Models-Their Training and Application[END_REF][START_REF] Lorenz | Generation of point-based 3d statistical shape models for anatomical objects[END_REF]. However, this approach required thousands of points to properly describe shape surfaces and a one-to-one correspondence between these boundary points [START_REF] Yan | A model to accumulate fractionated dose in a deforming organ[END_REF][START_REF] Lorenz | Generation of point-based 3d statistical shape models for anatomical objects[END_REF][START_REF] Söhn | Modelling individual geometric variation based on dominant eigenmodes of organ deformation: implementation and evaluation[END_REF][START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF][START_REF] Bondar | Statistical modeling of CTV motion anddeformation for IMRT of earlystage rectalcancer[END_REF]. SPHARM, the extension of 2D Fourier techniques to three dimensions, was proposed to obtain smooth, accurate and fine parametrization of 3D shapes with spherical topology with a fewer number of parameters [START_REF] Ballard | Computer Vision[END_REF][START_REF] Brechbühler | Parametrization of Closed Surfaces for 3-D Shape Description[END_REF]Székely et al., 1996;[START_REF] Kelemen | Elastic model-based segmentation of 3-D neuroradiological data sets[END_REF]. This method consisted in expressing a surface as a linear combination of a reduced number of spherical harmonics basis functions, where different levels of approximation errors can be obtained by truncating the spherical harmonic series. As opposed to other methods like PDM, it offers the advantage of delivering a shape parametrization with an implicit correspondence between shapes on the boundary. Besides, the spherical harmonics define an orthogonal basis functions whose spanned space has a reduced dimension that enables us to describe and quantify shape variations. This implies that the number of parameters (i.e., SPHARM coefficients) needed to describe a shape instance is lower than other shape modeling methods like PDM.

A C C E P T E D M A N U S C
In medical image, SPHARM has played an crucial role in 3D-shape representation and modeling of anatomic structures in several applications, including computer assisted diagnosis [START_REF] Shen | A surface-based approach for classification of 3d neuroanatomical structures[END_REF][START_REF] Styner | Boundary and medial shape analysis of the hippocampus in schizophrenia[END_REF]Styner et al., , 2006;;Shen et al., 2009b,a), rigid registration [START_REF] Dillenseger | Spherical harmonics based intrasubject 3-D kidney modeling/registration technique applied on partial information[END_REF][START_REF] Shen | Fourier method for large scale surface modeling and registration[END_REF] and organ segmentation [START_REF] Brechbühler | Parametrization of Closed Surfaces for 3-D Shape Description[END_REF]Székely et al., 1996;[START_REF] Kelemen | Elastic model-based segmentation of 3-D neuroradiological data sets[END_REF][START_REF] Gerig | Shape analysis of brain ventricles using SPHARM[END_REF]Tateyama et al., 2012), among others. In most of these applications, statistical shape models were developed by applying PCA to a training population of organ surfaces parameterized with SPHARM coefficients, where eigenmodes of maximum variation around the mean shape were calculated. In these applications, the surfaces of these binary organs were initially parameterized by surface boundary points, and subsequently expanded into series

A C C E P T E D M A N U S C R I P T of spherical harmonics.
In our study, SPHARM was used to parameterize the surface of each bladder into a vector of SPHARM coordinates. Each bladder surface had the same radial surface sampling in spherical coordinates, maintaining the correspondences across all the points. Following mesh parameterization, each bladder represented by function f ∈ R 3 was encoded as a reduced number of spherical harmonic functions, as follows:

f (θ, φ) = ∞ l=0 l m=-l c m l Y m l (θ, φ), (1)
≈ L l=0 l m=-l c m l Y m l (θ, φ), 0 < L < ∞
where θ and φ denoted the spherical coordinates, and with

Y m l (θ, φ) = (-1) m 2l + 1 4π (l -m)! (l + m)! P m l (cos(θ))e imφ (2)
defining a spherical harmonic function of degree l and order m, where P m l was an associated Legendre polynomial; L was the number of spherical harmonic functions used to approximate f ; and coefficients c m l corresponded to the coordinates in the space spanned by the spherical harmonic basis. The c m l coefficients were obtained as the projection of the mesh function f onto the spherical harmonic basis, i.e, c m l =< f (θ, φ), Y m l (θ, φ) >. By fixing a harmonic degree to L, the bladder surface at the j-th CT scan of the i-th patient, denoted by x i , j ∈ R d 1 ×d 2 ×d 3 , was uniquely represented in spherical harmonic coordinates as follows:

c i j = c 0 i j ,0 c 0 i j ,1 c 1 i j ,1 c 0 i j ,2 • • • c L i j ,L c -1 i j ,1 c -1 i j ,2 • • • c -L i j ,L T (3)
where

c i j ∈ R p ∧ p = 2 (L + 1) 2 - (L + 1) 2 -(L + 1) 2 << d = d 1 × d 2 × d 3 (4)
Observe that we have taken only the SPHARM coefficients that were non-zero. The value of L was optimally chosen by analyzing its impact on the approximation

A C C E P T E D M A N U S C R I P T
error in the training cohort. The approximation error was assessed as a function of L in terms of dice score (DS) and the Hausdorff distance (HD).

Principal component analysis

This step was aimed to represent, in a reduced space, all possible rigid and non-rigid bladder displacements. After obtaining SPHARM coefficients, an empirical covariance matrix C of local patient motion-deformation vectors was thus defined in order to calculate the directions of maximum variability by using matrix diagonalization. Here, a local patient motion-deformation vector defined the difference between the observed bladder and patient mean bladder. Thus, denoting the observed bladder of the i-th patient at the j-th CT as c i j , the empirical covariance matrix C was defined as follows:

C empirical = 1 n n i=1 1 j i -1 j i j=1 (c i j -ci )(c i j -ci ) T (5) = 1 n n i=1 1 j i -1 j i j=1 d i j d T i j = UDU T
where d i j = c i j -ci , n was the number of patients, j i was the number of observations available for the i-th patient, and ci the average shape of the i-th patient. Similarly, U was a matrix composed of eigenvectors ϕ k of C empirical , and D was the diagonal matrix constructed from its corresponding eigenvalues λ k , i.e, C empirical ϕ k = λ k ϕ k . Eigenvectors ϕ k defined the directions of geometric variability of the bladder, called modes. As a result, each bladder in the spherical harmonic space was expressed as a linear combination of the modes, as follows:

c i j = ci + z i j 1 ϕ 1 + • • • + z i j r ϕ r (6)
where z i jk = d T i j ϕ k , and r was the rank of covariance matrix C empirical , and equal to the number of available images in the population database. Modes {ϕ k } k defined a new coordinate system, where each vector d i j was represented by new r coordinates z i j 1 , . . . , z i j r . The eigenvalues of matrix D defined a significance ranking for the data variability represented by each mode, attesting var(ϕ

1 ) = λ 1 ≥ var(ϕ 2 ) = λ 2 ≥ • • • ≥ var(ϕ r ) = λ r ,
and showing that data variability was dominated by the first q modes with q << r.

A C C E P T E D M A N U S C R I P T

In order to test whether the modes described common directions of bladder motion and deformation among all patients, we performed a leave-one out cross validation on the training database. Thus, the error assessing the difference between the original and reconstructed motion-deformation vector was defined as follows:

e i = d i -di d i (7)
where d i and di were the i-th bladder motion-deformation vector and its approximation, respectively. If the reconstruction of d i was exact, then e i was equal to 0; otherwise, the value of e i was between 0 and 1. The mean of the reconstruction error e i in all the motion-deformation vectors was also calculated, i.e., M = (1/N) • N i=1 e i , where e i denoted the approximation error of an i-th motiondeformation vector d and N the number of the available images, either in the PCA model or left-out patients.

Linear mixed-effects model

Several longitudinal studies have been designed to investigate anatomical or functional changes over time in objects or variables that are observed repeatedly at several time points (Durrleman et al., 2013b;[START_REF] Laird | Random-effects models for longitudinal data[END_REF][START_REF] Hedeker | Longitudinal Data Analysis[END_REF]Fitzmaurice et al., 2011;[START_REF] Frees | Longitudinal and Panel Data: Analysis and Applications in the Social Sciences, 1st Edition[END_REF]. For instance, these repeated observations may be brain volume, response to clinical treatments, disease progression, blood pressure, heart motion or tumor evolution, among others [START_REF] Prabhu | Effect of the addition of chemotherapy to radiotherapy on cognitive function in patients with low-grade glioma: Secondary analysis of RTOG 98-02[END_REF][START_REF] Hanlon | Modeling postradiation prostate specific antigen level kinetics[END_REF][START_REF] Bastogne | Phenomenological modeling of tumor diameter growth based on a mixed effects model[END_REF][START_REF] Bernal-Rusiel | Spatiotemporal linear mixed effects modeling for the mass-univariate analysis of longitudinal neuroimage data[END_REF][START_REF] Ziegler | Estimating anatomical trajectories with Bayesian mixed-effects modeling[END_REF]. In these longitudinal studies, observations may be obtained under changing experimental conditions that are not possible to fully control, leading to considerable variations among subjects in the number and time points of the observations. The resulting data sets, which are known as unbalanced data, have often analyzed using ME models that offer to simultaneously model both the random effects that differentiate one subject from a population and the inherent correlation of repeated measurements in the same subject [START_REF] Liu | Chapter 3 -linear mixed-effects models[END_REF][START_REF] Frees | Longitudinal and Panel Data: Analysis and Applications in the Social Sciences, 1st Edition[END_REF][START_REF] Pinheiro | Mixed-Effects Models in S and S-PLUS, 1st Edition[END_REF]Fitzmaurice et al., 2011;[START_REF] Hedeker | Longitudinal Data Analysis[END_REF][START_REF] Liang | Longitudinal data analysis using generalized linear models[END_REF][START_REF] Laird | Random-effects models for longitudinal data[END_REF]. A linear ME model is then a statistical model that incorporates both "fixed effects", which are parameters related to the entire population or a certain level of classification factor, as well as "random effects", which are parameters associated with differences among individual subjects drawn at random from a population.

A C C E P T E D M A N U S C R I P T

The use of basis functions and ME have already proposed in the literature to account for random effects in longitudinal data. These basis function methods include: smoothing splines, wavelet basis, polynomial basis, and functional principal components [START_REF] Guo | Functional Mixed Effects Models[END_REF][START_REF] Morris | Wavelet-based functional mixed models[END_REF][START_REF] Wu | Local Polynomial Mixed-Effects Models for Longitudinal Data[END_REF][START_REF] Aston | Linguistic pitch analysis using functional principal component mixed effect models[END_REF]. In our study, the longitudinal data was encoded in the PCA modes, where the observed motion-deformation vectors in the training database were projected over each mode while a linear ME model was fit to each. The score resulting from each projection could be interpreted as a measure of how much a direction of geometric variability was presented in a motion-deformation vector. The implemented ME model thus aimed the characterization of population and individual bladder variations along each mode by introducing a "patient" factor modeled as a random effect.

Let z be the measure of motion/deformation of an organ along a mode. Considering the patient as the unit of observation state (subject), we used an index i, ranging from 1 to n, to differentiate among patients, and an index j to differentiate between observation times in a patient. The encoded longitudinal data was thus defined as follows:

{z 11k , . . . , z 1 j 1 k }, {z 21k , . . . , z 2 j 2 k }, • • • , {z n1k , . . . , z n j n k }, k = 1, . . . , q (8)
with the linear ME model proposed for each mode described as follows:

z i jk = µ k + b ik + ε i j k , i = 1, . . . , n, j = 1, . . . , j i , (9) b ik ∼ N(0, σ 2 bk ), i jk ∼ N(0, σ 2 k ),
where µ k was the mean projection along the k-th mode for the whole population, b i k was a random variable representing the deviation of the i-th patient mean from the population mean, and ε i j, k was a random variable representing the deviation for j-th projection of the i-th patient from the i-th patient mean. In addition, variances σ 2 bk and σ 2 k denote inter and intra-patient variability, respectively. The population based model is thus described as follows:

c = ci + Wz (10) c = c 0 0 c 0 1 • • • c -L L , W = ϕ 1 ϕ 2 • • • ϕ q , z = z 1 z 2 • • • z q where ci was the patient mean shape, z k ∼ N(µ k , σ 2 bk + σ 2 k ), for k = 1, . . . , q. A C C E P T E D M A N U S C R I P T 2.

Out-of-sample problem

We sought to predict a likely motion-deformation region of the bladder for a new patient using only his planning CT information. We described regions of motion/deformation in terms of 3D probability maps (PM) that express the probability of a voxel being occupied by the bladder during treatment. As shown in Fig. 3, the first step was to spatially normalize the out-of-sample patient P l to match the template patient; secondly, the SPHARM and PCA representations of his delineated bladder c l,0 and z l,0 were obtained; and thirdly, each ME model was adapted to the new patient by estimating his patient-specific mean and variance μlk and σl jk .

In ME models, random subject-specific effects are summarized to predict the response of a new subject. The predictor is a linear combination of both fixed and random effects, and it can be interpreted as a shrinkage estimator. Let us suppose that the l-th patient was included in the data set and we want to estimate the following linear model:

z l jk = µ k + b l k + l j k ∼ N( μl k , σ2 l k ), j = 1, . . . , j l , k = 1, . . . , q (11)
A shrinkage estimator of μl k = µ k + b l k is then given by:

μl k = µ k + ζ k (µ l k -µ k ), ζ k = j l j l + σ 2 k /σ 2 bk , k = 1, . . . , q (12)
where µ l k is the mean of the l-th patient at the k-th mode, and j l is the number of patient's observations. Term ζ is known as the credibility factor. Similarly, the variance of the prediction error σ2 l k = var( l j k ) is given by:

σ2 l j k = 1 - σ 2 bk σ 2 k + σ 2 bk n σ 2 k + σ 2 bk 1 - σ 2 bk σ 2 k + σ 2 bk - (σ 2 bk ) 2 σ 2 k + σ 2 bk + σ 2 k + σ 2 bk ( 13
)
In our study, cl = c l,0 , µ l = W T (c l -c l,0 ) = z l = 0 and j l = 1 due to we only had the information obtained from the planning CT. In this way, probable structures of the bladder during treatment were thus inferred after sampling S times the Normal distributed vector z. Subsequently, each sample of z was mapped to a vector c in the spherical harmonic space using matrix W. Then, binary images representing probable bladder configurations were obtained from the parameterized surfaces that were coded in vectors c. Following this, an estimation of the motion/deformation region was calculated by the ratio between the number of times that a voxel was occupied by the bladder and the S generated samples.

A C C E P T E D

Evaluation of prediction performance

Two additional population PCA-based models were used to compare the model's prediction performance. Considering motion-deformation vectors computed between the observed bladders and population mean bladder, the PCA model implemented by [START_REF] Söhn | Modelling individual geometric variation based on dominant eigenmodes of organ deformation: implementation and evaluation[END_REF] was extended to the population data as follows:

c = cpop + q k=1 z k ψ k , (14a)
Global PCA model := z k ∼ N(0, α k ) ∧ cpop = 1 N n i=1 j i j=1 c i j (14b) 1 N -1 n i=1 j i j=1 (c i j -cpop )(c i j -cpop ) T ψ k = α k ψ k , (14c)
where N = n i=1 j i was the total number of observations in the training database, cpop was the population mean bladder, and vectors ψ k were the eigenvectors of the empirical covariance matrix centered in terms of the population mean bladder. Observe that in this model it was thus considered that all the observed bladders in

A C C E P T E D M A N U S C R I P T
the training set came from a single individual. The second model proposed was taken from [START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF], defined as follows:

c = ci + q k=1 z k ϕ k (15a) Local PCA model := z k ∼ N(0, λ k ) ∧ ci = 1 j i j i j=1 c i j (15b) C empirical ϕ k = λ k ϕ k , (15c)
where ci was the mean bladder of the i-th patient, and vectors ϕ i were the eigenvectors of empirical covariance C empirical . Then, given a new patient with planned delineated bladder z l , both PCA models were adapted as follows:

Global PCA model : c = c l,0 + q k=1 z k ψ k , ∧ z k ∼ N(z l , α k ) (16) Local PCA model : c = c l,0 + q k=1 z k ϕ k , ∧ z k ∼ N(z l , λ k ) (17)
Both PCA based models used the modes to transform the delineated bladder on the planning CT and predict new bladder structures. It is worth mentioning that both sets of modes {ψ k } and {ϕ k } expressed different directions of organ motion and deformation. For example, modes {ψ k } described motion/deformation deltas regarding the population mean bladder whereas modes {ϕ k } described local patient motion/deformation deltas.

The models' prediction performance was evaluated by following leave-one-out cross validation in the training database and using the independent database. We thus estimated probability maps with each model and compared with the observed PM in both databases. Estimated PMs were derived for each model and patient, as follows: firstly, 40 bladder structures were generated using the bladder observed at the planning CT (i.e. S = 40, see Fig. 3); and secondly, a distribution of 15 PMs were derived from the 40 estimated structures, where each PM was obtained from a sub-set of estimated bladders (without replacement) whose cardinality was equal to the number of available CT/CBCTs.

Metrics of similarity between PMs

Two metrics were proposed for assessing the similarity between two PMs X and Y.

Mutual information-based metric

dist intensity (X, Y) = 1 - I(X, Y) H(X, Y) (18) H(X, Y) = - x,y p XY (x, y) log p XY (x, y) I(X, Y) = H(X) + H(Y) -H(X, Y) A C C E P T E D M A N U S C R I P T
where p(X, Y) is the joint probability distribution, and I(X, Y) and H(X, Y) denoted mutual information (MI) and joint entropy between X and Y, respectively. Similarly, x and y corresponded to intensity values in the voxels of X and Y, respectively. H(X) and H(y) denoted marginal entropies of each image. MI measures the amount of information that both images share (redundant information) based on the intensities of corresponding voxels, which is assumed to be maximal if both images are geometrically aligned. Similarly, joint entropy measures the uncertainty or dispersion between images. dist intensity also provided a distance based on the ratio of dependence and dispersion between two PMs, attesting dist intensity (X, Y) = 0 iff X = Y, and dist intensity (X, Y) ≤ 1 for all pairs (X, Y) [START_REF] Kraskov | MIC: Mutual information based hierarchical clustering[END_REF][START_REF] Russakoff | Image similarity using mutual information of regions[END_REF][START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF]. It offers the advantage of not making any prior assumptions regarding the functional relationship between images. Since we were dealing with binary images, each voxel was considered as a discrete random variable. Joint and marginal entropies were calculated by normalizing the joint, and marginal histograms, i.e a (joint) probability distribution of intensity values were estimated by counting the number of times each value occurred in the images and dividing this value by the total number of voxels.

Metric of misestimated voxel

We also proposed an error metric to assess whether a voxel was overestimated or underestimated. Denoting v as a voxel, let X b = {v/X(v) > 0} be the binary image obtained by those voxels with a probability greater than 0. The proposed metric is defined as follows:

Accuracy(X, Y) = |X p ∩ Y p |+|X p ∩ Y p | d ( 19
)
where d denoted the total number of voxels (see Eq. 4), |•| denoted the cardinality, and X p and Y p denoted the complement of X p and Y p , respectively. As a result, this metric provided the proportion of the properly estimated voxels (i.e. both true positives and negatives voxels in the motion/deformation region) among the total number of voxels.

Results

The analysis was carried out in MATLAB (R2014a). The eigenvectors ϕ k and ψ k were obtained using the function eigs while the ME model of each score z i jk was fit using the function fitlme with restricted maximum likelihood estimation (REML) as estimation method (Pinheiro andBates, 1996, 2000).

SPHARM and PCA representation

In the first phase of dimensionality reduction, an L of 15 spherical harmonics was selected to approximate bladder surfaces in the population database. With L = 15 in Eq. 1, the average DS and HD values in the training database were 0.954 ± 0.0522 and 6.112 ± 4.997 mm, respectively. Conversely, validation database equivalents were 0.97 ± 0.0099 and 4.32 ± 2.27 mm, respectively. Any bladder in the databases can thus be represented by p = 272 coefficients instead of d = 1596375 voxel variables. largest variance (approximately 60% of the variability); similarly, the second, third, and fourth modes contributed, on average, 9%, 7.8% and 5% of the variability, respectively. The first eight modes contributed 90% of the cumulative variance, on average. Meanwhile, Fig. 5 shows that both mean and standard deviation of the approximation error M decreased as the number of modes was increased in both the PCA model and left-out patients. However, mean and standard deviation of the reconstruction error M did not reach the value of zero in the left-out patients, and they exhibited a lower decreasing rate after the first eight modes. Considering only the first 60 modes, the mean reconstruction error M was 10% and 20% in the PCA model and left-out patient, respectively. When truncating the first 40, it still obtained an average accumulated variability of 98%

A C C E P T E D M A N U S C R I P T

and an average reconstruction error of 15% and 23% for the PCA model and left-out patients, respectively. Therefore, after performing two phases of dimensionality reduction, any bladder was represented as a vector of q = 40 scores (q = 40 << p = 272 << d = 1596375). Fig. 5 also shows bands (dotted lines) around the mean approximation error M, which represent the mean plus and minus one standard deviation of e i (see Eq. 7). It can be observed that the variance also decreased as the number of modes was increased in both the PCA model and left-out patients. dorsal regions. Fig. 7 depicts the projection of the motion/deformation vectors of the training patients on the space generated by the first two modes, where each point represents the observed geometric variations along these two modes. This figure also presents the observed motion/deformation vectors of Patient 6, indicating that the geometric variation followed a random trajectory. In addition, we can observe that Patient 6 and 17 had variance that differed both between them and from the mode's population variance. For example, Patient 6's bladder exhibited greater changes in volume during treatment than that of Patient 17, as shown in the PMs of both patients, proving that Patient 6's specific-variance along the first mode was greater than that of Patient 17. In addition, we also observed that Patient 17's variance along the second mode was very slight due to the geometric variations of his bladder being greater in the inferior-medial and dorsal regions.

Evaluation of performance prediction with groundwork models

Fig. 8 and 9 demonstrates the PMs obtained by all three models and the available observations of two patients as out-sample patients (one for each database). In both patients, the Global and Local PCA model overestimated more voxels than our proposed model. Tables 1 and2 present the averages obtained by these metrics for these patients. lues. The Table 3 shows the number of patients that had significant difference in both databases. Besides, we also obtained four variance distributions by deriving the variance of the 15 values for each metric and patient in both databases. Table 4 presents the p-values of the significant tests that were made between these variance distributions of the PCA models using the Wilcoxon rank sum test.

Discussion

We have hereby proposed a new methodology for predicting bladder motion and deformation between fractions in prostate cancer radiotherapy via population analysis. We calculated bladder motion and deformation regions based solely on the planning CT scan. In order to predict which regions were likely to exhibit motion/deformation, we first performed a prostate-based rigid-registration ing database. We used leave-one-out cross validation on the training database and an external database to validate the model. The proposed model was then compared to two additional population PCA based models. Motion/deformation regions were characterized by means of probability maps (PMs) and metrics were provided in order to measure similarities between observed and estimated PMs. Previous studies have matched CT scans for modeling geometric variations of the pelvic organs based on body anatomy, soft-tissue registration, or fiducial markers. For example, [START_REF] Söhn | Modelling individual geometric variation based on dominant eigenmodes of organ deformation: implementation and evaluation[END_REF] used bone anatomy to match for each patient the repeat CT scans to the planning CT in order to obtain an intra-patient model for the prostate/bladder/rectum. Thörnqvist et al. (2013) used fiducial markers to

A C C E P T E D M A N U S C R I P T
A C C E P T E D M A N U S C R I P T
model individual geometric variations of the prostate. In contrast, [START_REF] Hu | Population-based prediction of subject-specific prostate deformation for MR-to-ultrasound image registration[END_REF] proposed anatomical apex and base points of the prostate. Meanwhile, [START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF] defined the prostate barycenter as the origin for each CT scan in order to derive a population-based model of the CTV for prostate cancer. However, we believe that using the barycenter as a reference point is not feasible for all pelvic organs to obtain population-based models. For this reason, [START_REF] Bondar | Statistical modeling of CTV motion anddeformation for IMRT of earlystage rectalcancer[END_REF] used a setup based on bone anatomy for each patient, to calculate intra-patient motion-deformation vectors and develop a population model of CTVs for rectum cancer, as did [START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF]. Nevertheless, this method may not be the most suitable setup for analyzing in the same framework all the different sizes and shapes that an organ can manifest in patients [START_REF] Rios | Characterization of bladder motion and deformation in prostate cancer radiotherapy[END_REF]. In our study, we did not use bone anatomy as it can significantly differ from one patient to another, affecting the rigid registration and distribution of the pelvic organs. On the other hand, the bladder is attached to the prostate across the urethra, and the surface in contact with the prostate is generally stable, or is at least considered the least mobile part, while the top and anterior part can expand. We thereby sought to align this fixed bladder section in order to model both bladder motion and deformation. We believe that the ideal setup could be to match CTs using the boundary between the bladder and prostate, yet that type of rigid-registration is not possible. For intra-patient analysis, we think that the prostate barycenter offers a better means to perform this alignment due to the prostate volume being near constant during treatment, despite its rigid displacements [START_REF] Roeske | Evaluation of changes in the size and A C C E P T E D M A N U S C R I P T location of the prostate, seminal vesicles, bladder, and rectum during a course of external beam radiation therapy[END_REF][START_REF] Melian | Variation in prostate position quantitation and implications for three-dimensional conformal treatment planning[END_REF]

(see Fig. A.14a in appendix A).

For inter-patient analysis, we consider that differences in prostate size and shape in the population can be reduced by defining the template as a patient presenting a prostate with the lowest volume variation and the most similar mean volume to the population mean (see Fig. A.14a in appendix A). In our training database, this patient prostate volume was 42.14 ± 3.33cm 3 , whereas the population volume was 42.65 ± 17.43cm 3 . For example, Fig. A.14b shows 3D sagittal views of some segmented patient prostates that were plotted as outside the first and third quartiles of the boxplot of the average prostate volumes. In Fig . A.14a, it can be observed that the size and shape differences were not significant for patients 17, 7, and 20, particularly in the upper anterior part of the prostate (the surface in contact with the bladder). The highest dissimilarity was obtained for Patient 8, an outlier like Patient 17. We thus consider that size and shape differences among prostates in the population are not sufficiently greater than those observed in bladders to introduce disturbances, at least in the first eigenmodes. For example, it can be observed in Fig. 6 that the first three modes did not reveal geometric variations

A C C E P T E D M A N U S C R I P T
in the lower posterior bladder section, i.e., the region that presented the geometric variations in the bladder produced by the barycenter alignment (see Figs. A.13 and A.14). Also,Fig. A.13 demonstrates that for intra-patient analysis, the first mode was also dominated by the bladder volume, meaning that volume changes in the prostate do not affect at least the most important mode in the bladder. It is worth mentioning that we did not introduce any rotation, nor conduct affine transformation, as the bladder does not present rotations as the prostate does, nor do volume changes in the prostate correlate with volume changes in the bladder. Besides, if all the bladders in the population are transformed using affine or nonrigid registration, inter-patient geometric variability can disappear, which is the main advantage of using ME models, where individual effects are separated from the mean. We therefore considered that the prostate barycenter also provided a suitable reference point in the pelvic zone to model bladder motion and deformation for inter-patient analysis, despite the potential variations of prostate in size and shape across the population.

Similar to other SPHARM-based shape models, our study was significantly aided by encoding the information of the bladder shape into vectors of SPHARM coefficients with orders of 10 2 , instead of vectors of boundary points with orders of 10 3 obtained from a meshing process, as used in [START_REF] Söhn | Modelling individual geometric variation based on dominant eigenmodes of organ deformation: implementation and evaluation[END_REF][START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF]. The SPHARM parameterization thus enabled us to reduce not only computation time but also numerical errors, by affording finite arithmetic operations, and data storage while preserving shape information. However, it is important to highlight that not all the bladder shapes were properly parameterized using the radial sampling in spherical coordinates. In our training database, some bladders with non-convex dilations in the inferior-ventral region resulted in higher reconstruction errors. Future work should thus be dedicated to improve bladder-shape modeling. One way would be to extend the SPHARM methodology to model a much larger class of simply-connected closed shapes using an area-preserving, distortion minimizing spherical mapping, as expressed in [START_REF] Brechbühler | Parametrization of Closed Surfaces for 3-D Shape Description[END_REF]Shen et al., 2009a). Alternatively, a second possibility is to use a different computational technique called isometric analysis, which generalizes and improves on the standard finite element method (Zhang et al., 2007b;[START_REF] Nguyen | Isogeometric analysis: An overview and computer implementation aspects[END_REF]. In population data, this technique offers the advantage of obtaining a surface parameterization with corresponding points in all organ instances without requiring a mesh relaxation procedure. This methodology defines a transformation to a canonical domain, where parametrization is carried out and extrapolated to any organ instance.

When considering population-based models developed for pelvic organs, we

A C C E P T E D M A N U S C R I P T
found that 15 modes were required to describe 90% of CTV motion and deformation for prostate cancer [START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF]. For rectum cancer, 28 modes were found to be required to describe 90% of the accumulated variability in CTVs [START_REF] Bondar | Statistical modeling of CTV motion anddeformation for IMRT of earlystage rectalcancer[END_REF]. In contrast, eight modes were sufficient in our study to describe the same percentage of accumulated variability. On analysis of the three methods' results, we noted that there were less fundamental directions of geometric variability in the bladder compared to other pelvic organs, despite its great volume and shape variations. Nevertheless, we were interested to note that more directions were observed in the prostate/SV composition, potentially due to changes in the bladder and rectum volumes, or due to patient-specific variations and delineation errors. Regarding the rectum, its greater number of directions may result from its flexible and vertical structure, enabling it to deform in any direction along its major axis [START_REF] Hoogeman | A model to simulate day-to-day variations in rectum shape[END_REF], in addition to patient-specific variations and delineation errors. Future studies using other population databases would be required to confirm or contradict the number of principal directions of motion and deformation reported in these organs.

On the other hand, as in [START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF], the reconstruction error in the PCA model and left-out patients in our study tended to decrease towards 0 as the number of modes was increased. However, the error did not completely disappear in the left-out patients, meaning that not all individual directions were accounted for by the PCA modes. In addition, as observed in Fig. 5, the rate in the left-out patient curve decreased after the first eight modes, implying that the order of the modes was not preserved in all patients. As a result, both of these issues may significantly affect the performance of any population PCA model. For example, the lowest measures of similarity were observed in Patients 6, 12, and 20 in the training database, whose bladders exhibited geometric variations in the superior-ventral region (see Patient 6's PM in Fig. 7). These variations were not represented by at least the first four modes (81% of the accumulated variability). It would therefore be interesting for future studies to first make a classification of bladders with similar shapes [START_REF] Rios | Characterization of bladder motion and deformation in prostate cancer radiotherapy[END_REF], then design the proposed model locally for each subgroup. For example, [START_REF] Bishop | A hierarchical latent variable model for data visualization[END_REF] introduced a methodology for developing hierarchical PCA models where each observation is plotted in the latent space to determine the number of sub-PCA models in the inferior levels.

It is also crucial to emphasize that we applied a quadratic relative error in the SPHARM space to assess the modes in our study, unlike [START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF], who used a quadratic error normalized by the number of points. Nevertheless, we considered it possible to compare the behavior of both errors. For example,

A C C E P T E D M

A N U S C R I P T [START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF] required 37 modes to reduce the initial error by 62% in the left-out patients, whereas 60 modes were needed in our study to reduce it by 80%. Both results may indicate there was significant intra-patient variability that was not sufficiently detected by the first modes obtained from the population. This also means that intra-patient variability may be higher in the bladder, potentially due to anatomical differences among patients and delineation errors. In this way, a trade-off was made between accumulated variability and the reconstruction error for determining the number of modes (40) in our PCA model, where an additional 20 modes only reduced the reconstruction error by 3%.

As established in the model assumptions, each mode corresponded to a direction of geometric variability that dilated or contracted certain regions of the bladder. Unlike the first mode, not all the rest were found to define an affine transformation that may be initially applied in the setup, such as rotation along an axis or a uniform dilation/contraction. In addition, we showed that the bladder motion and deformation between fractions follows random trajectory, as the data projections showed in the latent space (see Fig. 7). It was also possible, in this space, to characterize the geometric variations of each patient, where some directions were clearly found to be more dominant than others.

Although several studies have previously addressed the estimation of motion/deformation regions of the pelvic organs by means of the coverage matrix concept [START_REF] Stroom | Inclusion of geometrical uncertainties in radiotherapy treatment planning by means of coverage probability[END_REF][START_REF] Mageras | Computerized design of target margins for treatment uncertainties in conformal radiotherapy[END_REF][START_REF] Craig | Considerations for the implementation of target volume protocols in radiation therapy[END_REF][START_REF] Price | A method to calculate coverage probability from uncertainties in radiotherapy via a statistical shape model[END_REF], this is the first, to our knowledge, to introduce a quantitative and qualitative comparison between estimated and observed regions using metrics of similarity and joint histograms. As observed in Figs. 10, 11, and 12, our proposed model obtained the lowest measures of dissimilarity and misestimated voxels. Previous population PCA models have thus overestimated or underestimated more voxels, primarily those not occupied by the bladder. As a result, our proposed model reduced uncertainties in estimations of the probable region of motion and deformation. This was achievable as the ME models enable us to reduce the observed population variance along each mode by grouping repeated observations per patient. For example, we obtained more patients presenting significant differences between the Global PCA model and our proposed model in both databases than the Local PCA model; moreover, we significantly reduced uncertainties in estimations of motion/deformation region regarding the Global PCA model (see Table 4). We therefore concluded that the Global, Local, and Local PCA-ME models provided progressive decrease in uncertainties for estimating likely motion/deformation regions, while our proposed model led to the lowest reduction in all the patients (see Fig. 8 and9). When more CTs/CBCTs scans become

A C C E P T E D M A N U S C R I P T
available for a patient during treatment, this information could be used to improve estimation of the motion/deformation region by personalizing modes, variance, and mean bladder shape. It is worth mentioning that one limitation to use PCA models for modeling organ motion and deformation is their strong dependence on the organ's mean shape, as planning CT scans may not provide a good estimation. Future studies should be focused on analyzing the methodology's robustness for handling outlier patients whose variance along a mode is almost similar to the population variance. This kind of patient may lead to a misestimation of the intrapatient variability σ 2 affecting the performance of the ME model (see Patient 6 in the training database in Fig. 7). A robust formulation of PCA, the so-called robust PCA, could help overcome the problems associated with outlier patients.

PCA-based shape models are typically efficient for organs whose shape variations can be captured by a reasonable number of modes [START_REF] Heimann | Statistical shape models for 3d medical image segmentation: A review[END_REF]. For instance, two studies have shown that 15 and 16 modes were required to describe most of the inter-and intra-patient motion and deformation of CTVs in prostate and rectum cancer, respectively [START_REF] Budiarto | A population-based model to describe geometrical uncertainties in radiotherapy: applied to prostate cases[END_REF][START_REF] Bondar | Statistical modeling of CTV motion anddeformation for IMRT of earlystage rectalcancer[END_REF]. However, in population databases concerning highly-varying soft-tissue volumes, a large number of modes may be required to properly capture the complex shape variations observed in the training database. In our study, for instance, we required 40 modes to generate an inter-patient PCA model for the bladder. We believe that the large number of modes may render interpretation of higher modes difficult, as they may not describe directions of geometric variability observed in the population rather than patient-specific directions. In addition, the large number of modes may also hinder inference and analysis of the results for each ME model. However, as in [START_REF] Aston | Linguistic pitch analysis using functional principal component mixed effect models[END_REF], the modes define an orthogonal basis that enables us to consider each score z i jk as an independent random variable across k that can be easily modeled as a response variable [START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF][START_REF] Jolliffe | Principal Component Analysis[END_REF].

We have previously mentioned that the scores can be interpreted as measures of geometric variability along the modes. Future studies should therefore focus on determining covariates that may help to describe patient-specific bladder deformation. The results displayed in Fig. 7 and published by [START_REF] Rios | Characterization of bladder motion and deformation in prostate cancer radiotherapy[END_REF] appear to indicate that there are some covariates that may underlie the motion/deformation region in the latent space or mode's variance for any arbitrary patient. For example, bladder size or volume could be considered as fixed parameters that may determine changes in volume and directions of geometric variations. [START_REF] Rios | Characterization of bladder motion and deformation in prostate cancer radiotherapy[END_REF] and (Palorini et al., 2016a) reported, in line with Fig. 7, that bladders with large volumes exhibited geometric variations in the superior-ventral region and A C C E P T E D M A N U S C R I P T higher variance along the first mode. In an adaptive radiotherapy scheme, a time covariate can also be added, as a random or fixed variable, to describe the daily volume changes of the bladder.

As previously stated, [START_REF] Hu | Population-based prediction of subject-specific prostate deformation for MR-to-ultrasound image registration[END_REF] proposed a population-PCA model to describe prostate deformation in MR-tumor-targeted biopsies. Concerning this methodology, we found two shared characteristics: firstly, that the scores were also considered as parameters that describe intra-subject organ motion; secondly, the added objective of characterizing subject-specific probability density functions (SSPDFs) of motion and deformation (parameterized by means of multivariate Gaussian distributions). However, [START_REF] Hu | Population-based prediction of subject-specific prostate deformation for MR-to-ultrasound image registration[END_REF] proposed a fixed-effects model to describe SSPDF of the scores across the population, while such a model consists of a multivariate non-linear regression involving the scores of the reference shapes defined as covariates. [START_REF] Hu | Population-based prediction of subject-specific prostate deformation for MR-to-ultrasound image registration[END_REF] also decided to model the SSPDF parameters instead of the subject-specific scores, as proposed in our study. Nevertheless, we believe that it could be interesting for future studies to add the scores of the bladder segmented at the planning CT as fixed variables, as this may help to correlate the subject-specific variance and directions of geometric variability. For future works, we wish to underline that both studies have suggested the possibility of including other parameters in the learning framework, such as organ size or any temporal information.

For future developments, it is worth noting that our main goal has been to predict late GU toxicity based only on planning treatment information, and we thus only used the planning CT scan to predict bladder motion/deformation regions. Nonetheless, we have also performed some simulations using two and three observations of the out-of-sample patient, with the aim of simulating an adaptive radiotherapy treatment, and the results achieved better estimation of the region. We think that the number of observations for the out-of-sample patient helps to decrease the variance of the estimate for the shrinkage predictor. Thus, the more observations we have for the out-of-sample patient, the more certain we can be about how the patient differs from the population mean. In these simulations, it is also worth mentioning that we only used these observations to improve the estimation of the patient's average bladder and shrinkage estimator of each mode, i.e. to improve the model adaptation yet not retrain the model's parameters (modes ϕ k , and inter-intra-patient variances). We considered that more observations should be available to fit the model's parameters again, as five CTs/CBCTs has been established as the minimum number of observations required for intra-patient models [START_REF] Söhn | Modelling individual geometric variation based on dominant eigenmodes of organ deformation: implementation and evaluation[END_REF].

In this study, we described regions of motion/deformation in terms of 3D pro-

A C C E P T E D M

A N U S C R I P T bability maps (PM) that were obtained by sampling the model's distribution or using patient images. We have thus considered motion/deformation regions as histograms that can be compared using the proposed metrics. We believe that the next step should be focused on estimating the uncertainty of the motion/deformation region using the model's joint distribution. However, we considered that the large number of modes may hinder inference of the associated joint distribution, as the number of combinations of scores at the boundaries of the confidence intervals will significantly grow when we add modes to the joint distribution. Nonetheless, we consider that the uncertainty region can also be estimated using only those modes that have a strong influence on variables, like the accumulated delivered dose, i.e., that the uncertainty region can be restricted to those modes that have a strong influence on the mean accumulated dose. Finally, it may also be interesting to add linearly modes that describe geometric variations between the mean bladders of the patients as the modes obtained from the covariance matrix C empirical only described directions observed on the intra-patient level. These models are known as multi-level PCA models [START_REF] Di | Multilevel functional principal component analysis[END_REF][START_REF] Greven | Longitudinal functional principal component analysis[END_REF].

Conclusion

We proposed a population based model to predict bladder motion and deformation between fractions using solely the planning CT. In comparison with previous studies, our proposed model was able to decrease uncertainty in the estimation of the region where the bladder will likely move and deform. We also demonstrated that, by following a longitudinal study with ME models, it was possible to separate the patient-specific variance from population variance and thus reduce the total variance. The potential applications of this model include margin evaluation, delivered dose estimation, toxicity prediction, and the design of robust treatment plans, among others. Future studies should now be conducted to validate the proposed model by means of a large cohort of patients who have undergone prostate cancer radiotherapy. It could also be valuable to focus on the relationship between inter-fraction bladder motion/deformation, delivered dose, and late GU toxicity in the prostate cancer radiotherapy context.

Acknowledgements

This research was supported by the Departamento Administrativo de Ciencia, Tecnología e Innovación (COLCIENCIAS) -Colombia with the "Becas de Doctorado en Colombia 2010 -Convocatoria 511" grant, and also by the Institut National du Cancer(INCa) -France with the "STIC IGRT-P" project.

Figure 1 :

 1 Figure 1: Workflow of the method used to train the population-based model. 0) Prostate-based registration, 1) spherical harmonics (SPHARM) representation, 2) principal component analysis (PCA) reduction, 3) ME model.

Figure 2 :

 2 Figure 2: Overlay of a typical individual bladder registered to the template.

Figure 3 :

 3 Figure 3: Schematic representation for estimating new structures based on the planning CT for a new patient.

  Four PMs were thus calculated for each out-sample patient: one PM estimated by means of the proposed model, labeled PM Local PCA-ME ; one PM obtained A C C E P T E D M A N U S C R I P T from the patient's observed available images, labeled PM obs ; one PM estimated by means of the global PCA model, labeled PM Global PCA ; and one PM estimated by using the local PCA model, labeled PM Local PCA .

Figure 4 :

 4 Figure 4: Relative mode variances obtained from the PCA in the training database (relative values: sum of all eigenvalues normalized to 100%).

Figure 5 :

 5 Figure 5: Mean approximation error M as a function of the number of modes in the PCA model and left-out patients. There are bands (dotted-lines) that represent the mean approximation error M (solid-lines) plus and minus one standard deviation of e i (see Eq. 7).

Fig. 4 Figure 6 :

 46 Figure6: 3D sagittal views of the first three motion/deformation modes applied on the population mean bladder (center column). The left and right columns correspond to geometric variations along the respective mode ϕ i , defined as population mean bladder ± 2σ i ϕ i , where σ 2 i is the mode variance. The prostate (red color) and rectum (green color) of the template patient were also included as spatial reference.

Figure 7 :

 7 Figure 7: Scatterplot of the patients from the training database projected in the first two modes (black points). Color circles and squares correspond to the vectors of two individuals. d 6, j denotes the observed motion/deformation vectors of Patient 6's bladder.

Fig. 6

 6 Fig.6presents sagittal views of the first three modes over the mean bladder of the training database. The first mode was associated with bladder volume. Similarly, modes with greater indices, like the second and third, represented directions of dilation and contraction in some regions. For example, the second mode indicated dilation/contraction in the left-dorsal and right-ventral regions of the bladder whereas the third mode indicated dilation/contraction in the inferior-ventral and

  Figure 9: 3D sagittal views of PMs obtained for one patient in the validation database as an outof-sample patient. Segmented prostate (blue) and rectum (green) at the planning CT scan are overlaid. Top left figure also depicts the segmented bladder (cyan) at the planning CT.

Figure 10 :Figure 11 :

 1011 Figure 10: Average dist intensity values of observed and estimated PMs for patients in both databases.

  Fig. 12 also presents the average joint histograms of the voxel intensities of PM obs vs PM Local PCA-ME -PM Global PCA and PM obs vs PM Local PCA-ME -PM Local PCA in both database. This presents the number of times that a joint intensity value was repeated between the observed and estimated PMs. For example, Figs. 12a and 12c reveal that the global PCA overestimated more voxels than our proposed model in both databases in probabilities p ranging between 0 and 0.6 instead of A C C E P T E D M A N U S C R I P T

Figure 12 :

 12 Figure 12: Average joint histograms of observed and estimated PMs for patients in the training and validation database. Left column shows the average joint histograms between PM obs vs PM Local PCA-ME -PM Global PCA while right column shows the average joint histograms between PM obs vs PM Local PCA-ME -PM Local PCA .

  R I P T

			Registered and template
			patient
	Unregistered
	patient	rigid Prostate based	bladder Template	Registered bladder
	Template	registration	Registered prostate	Template prostate
	patient

Table 1 :

 1 Averaged values from the metric dist intensity for two patients between estimated and observed PMs.

		PM obs	PM obs	PM obs
		vs	vs	vs
	Patient	PM Global PCA PM Local PCA PM Local PCA-ME
	Training database	0.691	0.676	0.66
	Validation database	0.74	0.66	0.65

Table 2 :

 2 Averaged values from the metric Accuracy for two patients between estimated and observed PMs.

		PM obs	PM obs	PM obs
		vs	vs	vs
	Patient	PM Global PCA PM Local PCA PM Local PCA-ME
	Training database	0.968	0.972	0.98
	Validation database	0.946	0.97	0.975
	the observed value of p = 0. Meanwhile, Figs. 12b and 12d demonstrate that
	our proposed model was also more capable of estimating voxels observed with
	probability p = 0.

Table 3 :

 3 Number of patients exhibiting significance difference using the Wilcoxon rank sum test in both databases.

		Local PCA-ME model Local PCA model Local PCA-ME model
		vs	vs	vs
		Global PCA model	Global PCA model	Local PCA model
	Database	dist intensity /Accuracy	dist intensity /Accuracy	dist intensity /Accuracy
	Training	17/18	17/14	0/0
	Validation	20/27	20/26	5/3

Table 4 :

 4 p-values of the Wilcoxon rank sum tests for the variance distributions between PCA models in both databases. dimensionality reduction steps, and fitting regression. The first dimensionality reduction step consisted of surface parameterization of the bladder in the spherical harmonic space. The second step of dimensionality reduction was to determine the directions of bladder geometric variability by means of PCA. These directions, called modes, were validated by leave-one-out cross validation.

	Local PCA-ME model Local PCA model Local PCA-ME model
		vs	vs	vs
		Global PCA model	Global PCA model	Local PCA model
	Database	dist intensity /Accuracy	dist intensity /Accuracy	dist intensity /Accuracy
	Training	0.0411/0.0001	0.0720/0.003	0.6/0.32
	Validation	0.037/3 * 10 -5	0.10/0.0001	0.51/0.44
	alignment, two

A Székely, G., Kelemen, A., Brechbühler, C., Gerig, G., Mar. 1996