H Ladjal
email: hamid.ladjal@liris.cnrs.fr

M Beuve

P Giraud

Towards Non-invasive Lung Tumor Tracking Based on Patient Specific Model of Respiratory System

Hamid Ladjal, Michael Beuve, Philippe Giraud and Behzad Shariat (memberIEEE)

Abstract-The goal of this paper is to calculate a complex internal respiratory and tumoral movements by measuring respiratory air flows and thorax movements. In this context, we present a new lung tumor tracking approach based on a patient-specific biomechanical model of the respiratory system, which takes into account the physiology of respiratory motion to simulate the real non-reproducible motion. The behavior of the lungs, is directly driven by the simulated actions of the breathing muscles, i.e. the diaphragm and the intercostal muscles (the rib cage). In this paper, the lung model is monitored and controlled by a personalized lung pressure/volume relationship during a whole respiratory cycle. The lung pressure and rib kinematics are patient specific and obtained by surrogate measurement. The rib displacement corresponding to the transformation which was computed by finite helical axis method from the end of exhalation (EE) to the end of inhalation (EI). The lung pressure is calculated by an optimization framework based on inverse finite element analysis, by minimizing the lung volume errors, between the respiratory volume (respiratory airflow exchange) and the simulated volume (calculated by biomechanical simulation). We have evaluated the model accuracy on five public datasets. We have also evaluated the lung tumor motion identified in 4D CT scan images and compared it with the trajectory that was obtained by finite element simulation. The effects of rib kinematics on lung tumor trajectory were investigated. Over all phases of respiration, our developed model is able to predict the lung tumor motion with an average landmark error of 2.0 ± 1.3mm. The results demonstrate the effectiveness of our physics-based model. We believe that this model can be potentially used in 4D dose computation, removal of breathing motion artifacts in positron emission tomography (PET) or gamma prompt image reconstruction.

Index Terms-Respiratory motion -Biomechanics modeling -Finite Element Method, Radiation therapy

I. INTRODUCTION

One of the major difficulties in radiation therapy is the treatment of moving tumours. The modern external beam radiotherapy techniques are available for cancer treatment, including Three-Dimensional Conformal Radiation Therapy (3D-CRT), Intensity Modulated Radiation Therapy (IMRT), Helical Tomotherapy (HT) and Hadron THerapy (HTH). The HTH has been the subject of growing interest, primarily due to the favorable ballistic properties of ion-matter interactions, which allow for a high degree of dose conformality in the tumor while minimizing the dose in the healthy tissue. Unfortunately, in many cases the advantage of the increased dose conformity is severely compromised due to organ motion. The internal anatomy of the patient can undergo intrafractional and inter-fractional deformation and displacement [START_REF] Mori | Effects of Interfractional Anatomical Changes on Water-Equivalent Pathlength in Charged-Particle Radiotherapy of Lung Cancer[END_REF], [START_REF] Ehrhart | 4D Modeling and Estimation of Respiratory Motion for Radiation Therapy[END_REF]. Moreover, in the case of Hadron THerapy (HTH), apart the ballistic problem caused by the target motion, the organs density variations due to respiratory motion affect the particle range which can lead to unwanted dose distributions [START_REF] Grzinger | Volume conformal irradiation of moving target volumes with scanned ion beams[END_REF], [START_REF] Ac | Range Verification Methods in Particle Therapy: Underlying Physics and Monte Carlo Modeling[END_REF], overshoot to critical organs or beam undershoot to tumour. The error due to density changes depends the position of the tumour and the beam direction and may exceed 1 or 2 cm [START_REF] Moskvin | A semiempirical model for the therapeutic range shift estimation caused by inhomogeneities in proton beam therapy[END_REF]. The effects of density variations that occur along the beam path on dose distribution are more important in the case of HTH than in the case of conventional radiotherapy [START_REF] Manescu | Four-dimensional radiotherapeutic dose calculation using biomechanical respiratory motion description[END_REF]. In this paper we are going to focus on breathing-induced intrafractional motion. Respiration is the most important source of intra-fractional motion in the patient's body. Lung tumour displacement can range up to more than 2 cm, while lung volume was being varied by about 25% during respiration [START_REF] Langen | Organ motion and its management[END_REF], [START_REF] Shirato | Speed and amplitude of lung tumor motion precisely detected in four-dimensional setup and in real-time tumor-tracking radiotherapy[END_REF]. In this context, various intra-fractional motion mitigation techniques have been developed. One solution is to reduce the organ motion of the patient using abdominal compression or partial/complete stop of the patient respiration [START_REF] Mohatt | A Quantitative Assessment of Lung Tumor Motion under Abdominal Compression using 4D-CT Imaging[END_REF]. However, these techniques would substantially increase the patient setup time as well as the associated clinical workload. Moreover, the abdominal compression methods can be hard to reproduce. Breath hold level and breathing patterns can drift between time of simulation and treatment or even within a treatment session as the patient becomes more relaxed. Other existing solutions consist on beam gating or beam tracking. The beam gating refer to the irradiation during a selected period of the breathing cycle known as the gating window [START_REF] Modiri | Inversed-Planned Respiratory Phase Gating in Lung Conformal Radiation Therapy[END_REF], [START_REF] Mizuhata | Respiratory-gated Proton Beam Therapy for Hepatocellular Carcinoma Adjacent to the Gastrointestinal Tract without Fiducial Markers[END_REF]. This technique has the major disadvantage that the treatment time is prolonged due to the frequent interruptions of the beam delivery. Moreover, the beam gating, as well as mitigation techniques, are based on previously acquired images and therefore cannot take into consideration the variations of the patient's breathing pattern. An example of a beam tracking system used in the clinical environment is the robotic Cyberknife Synchrony system applied to xray radiosurgery [START_REF] Kilby | The CyberKnife Robotic Radiosurgery System in 2010[END_REF], [START_REF] Yang | Target margin design for realtime lung tumor tracking stereotactic body radiation therapy using CyberKnife Xsight Lung Tracking System[END_REF]. This technique increases dose to the patient from imaging and is an invasive procedure [START_REF] Shirato | Speed and amplitude of lung tumor motion precisely detected in four-dimensional setup and in real-time tumor-tracking radiotherapy[END_REF]. Furthermore, the density changes along the beam path cannot be properly determined which is essential in the case of HTH.

Other existing solutions consist of tracking by implantation of fiducial markers. These solutions are either intrusive or invasive and there is a risk of infection and pneumothorax [START_REF] Shirato | Speed and amplitude of lung tumor motion precisely detected in four-dimensional setup and in real-time tumor-tracking radiotherapy[END_REF], [START_REF] Giroux | Patient-Specific Biomechanical Modeling of the Lung Tumor for Radiation Therapy[END_REF].

An ideal solution would continuously adapt the beam to the variation of the tumour position. An alternative solution is to use a motion model that is able to model the correspondence between the internal displacement and deformations of the organs and respiratory surrogate signals such as spirometry and deformation of the skin surface. Moreover, in the case of four dimensional radiotherapy treatment planning that aims to estimate the evolution of the dose distribution during the respiratory cycle on one hand the position of the organs, on the other hand the displacements and the deformations of the tissues throughout the respiratory cycle need to be estimated. Breathing is a dynamic and a complex process in which the respiratory motion is non-reproducible [START_REF] Shirato | Speed and amplitude of lung tumor motion precisely detected in four-dimensional setup and in real-time tumor-tracking radiotherapy[END_REF].

As shown by many studies, the breathing periodicity and amplitude of patients can vary during a single imaging or treatment session [START_REF] Vedam | Acquiring a four dimensional computed tomography dataset using an external respiratory signal[END_REF], [START_REF] George | The application of the sinusoidal model to lung cancer patient respiratory motion[END_REF]. Furthermore, lung tumors can even present hysteresis in their trajectories [START_REF] Seppenwoolde | Precise and real time measurement of 3D tumor motion in lung due to breathing and heartbeat, measured during radiotherapy[END_REF], making them more difficult to locate with precision. Methods to estimate respiratory organ motions can be divided into 4 major categories; deformable image registration [START_REF] Han | A hybrid patient-specific biomechanical model based image registration method for the motion estimation of lungs[END_REF], [START_REF] Chen | internal motion estimation by internal-external motion modeling for lung cancer radiotherapy[END_REF], biomechanical models [START_REF] Behr | A three-dimensional human trunk model for the analysis of respiratory mechanics[END_REF], [START_REF] Li | Biomechanical model-based 4DCT simulation[END_REF], [START_REF] Fuerst | Patient-Specific Biomechanical Model for the Prediction of Lung Motion from 4D CT Images[END_REF], [START_REF] Behnaz | A Biomechanical Model of Human Lung Deformation Utilizing Patient-Specific Elastic Property[END_REF], [START_REF] Coelho | A biomechanical approach for in vivo diaphragm muscle motion prediction during normal respiration[END_REF], [START_REF] Chen | internal motion estimation by internal-external motion modeling for lung cancer radiotherapy[END_REF], hybrid models with biomechanics and deformable image registration [START_REF] Kk | Accuracy and sensitivity of finite element model-based deformable registration of the prostate[END_REF], [START_REF] Al-Mayah | Toward efficient biomechanical-based deformable image registration of lungs for imageguided radiotherapy[END_REF] and statistical models including or not the machine learning techniques [START_REF] Rostampour | Markerless respiratory tumor motion prediction using an adaptive neuro-fuzzy approach[END_REF], [START_REF] Terunuma | Novel real-time tumor-contouring method using deep learning to prevent mistracking in X-ray fluoroscopy[END_REF], [START_REF] Li | Prediction of Lung Motion From Four-Dimensional Computer Tomography (4DCT) Images Using Bayesian Registration and Trajectory Modelling[END_REF]. Currently, the existing methods do not explicitly take into account the information related to breathing physiology and physical properties of organ tissues, and all those methods assume a reproducible motion of the respiratory system and cannot fully take into account the variability of the respiratory motion. Recently, the authors in [START_REF] Fuerst | Patient-Specific Biomechanical Model for the Prediction of Lung Motion from 4D CT Images[END_REF] have proposed a generative lung biome-chanical model for half respiration cycle, where the motion is not constrained by any fixed boundary condition. The authors have used 4 and 16 pressure zones on the sub-diaphragm and thoracic cavity, respectively. However, this method, could hardly be correlated with respiratory surrogate signals during treatment. In order to overcome these limitations, the aim of this study is to propose a non-invasive approach for lung tumor tracking based on a patient-specific model of the respiratory system, which can be used during External Beam Radiation Therapy (EBRT). Fig. 1 illustrates our global methodology to predict and track the internal motion (lung tumor motion), from external measurable parameters, such as spirometry and Real-time Position Management (RPM), during treatment. In this study, we have developed a lung tumor tracking approach based on a patient-specific biomechanical model of the respiratory system, which takes into account the physiology of respiratory motion to simulate the real non-reproducible breathing motion. The behavior of the lungs, is driven directly by simulated actions of the breathing muscles, i.e. the diaphragm and the intercostal muscles (the rib cage). The lung model is monitored and controlled by a personalized lung pressure/volume relationship (lung compliance) during a whole respiratory cycle. The lung pressure and diaphragmatic forces are patient-specific, and calculated for different breathing states. However, as these parameters are not directly measurable, we have used the volume of the lungs, which is a measurable parameter, and through a personalized compliance pressure/volume curve, to estimate the lung pressure and to calculate diaphragmatic forces at different respiratory states. In this order, we have developed an optimization framework based on inverse finite element analysis, by minimizing the lung volume errors, between the respiratory volume (measured respiratory airflow exchange) and the simulated volume (calculated by biomechanical simulation) [START_REF] Giroux | Patient-Specific Biomechanical Modeling of the Lung Tumor for Radiation Therapy[END_REF]. Finally, our model is controlled by two physiological parameters: the respiratory volume and the ribs displacement. The respiratory volume can be replaced by surrogate signal such as spirometry during treatment and the ribs displacement by Real-time Position Management (RPM) system. In the context of HTH, our model takes into account the density variations due to breathing. However, in this paper, we do not address the real time challenging issues but we focus on modeling aspects. First, in section 2, we present the patient-specific anatomical and biomechanical model, as well as the lung-pressure/diaphragm-force optimization algorithms based on inverse finite element. Then, in section 3, a qualitative and quantitative analysis and an experimental validation of our approach are presented. Finally, we give some concluding remarks and directions for future work.

II. MATERIALS AND METHODS

A. Anatomy and physiology of the respiratory system

The lung is a passive organ which is divided into two halves, the right and left lung. It is situated in the thorax on either side of the heart. The pleural cavity is surrounded by the chest wall on the sides, and the diaphragm on the bottom. This space contains pleural fluid which facilitates near frictionless sliding at this boundary. The diaphragm is a dome-shaped musculofibrous membrane concave toward the lungs which separates the thorax from the abdominal cavity. It is composed of a peripheral part (muscular fibre) and a central part (tendon). Lungs are linked to the diaphragm and to the ribs through the pleura. The mechanics of human breathing involves two steps that alternate with each other: inhalation (inspiration) and exhalation (expiration). Negative pressure in the pleural cavity (natural breathing) initiates when the diaphragm and chest wall move away from the lung. The negative pressure expands lung volume, dropping the internal lung pressure, allowing air to enter passively in the lung. The ability of the lungs to expand is expressed by using a measure known as the lung compliance (Fig. 2). Lung compliance (C L ) is the relationship between how much pressure is required to produce a degree of volume change of the lungs. It is usually expressed in ml/cmH2O. The ratio between the volume variation △V and the pressure △P represents the compliance:

C L = △V △P (1)
The lung compliance is affected by the elastic properties of the lung. The pulmonary compliance therefore reflects the lungs ability to develop in response to an increase in pressure.

The lung compliance can be written as:

K = E 3(1 -2ν) = -V 0 dP dV = - V 0 C L (2)
where K, E, ν, V 0 are the bulk modulus, the Young's modulus, Poissons ratio and the initial lung volume respectively.

B. Patient-specific anatomical model of respiratory system

One of the major difficulties of the finite element method applied to computational biomechanics is the complexity associated with the building of patient-specific geometrical models. In this order, we have developed previously a global strategy to generate patient-specific anatomical 3D models of the respiratory system from medical images (CT scan data) for finite element simulation [START_REF]Physiological and Biomechanical Model of Patient Specific Lung Motion Based on 4D CT Images[END_REF], [START_REF] Ladjal | Biomechanical Modeling of the Respiratory System: Human Diaphragm and Thorax[END_REF], [START_REF] Giroux | Patient-Specific Biomechanical Modeling of the Lung Tumor for Radiation Therapy[END_REF]. For this, after segmentation, a 3D surface mesh is created for each volume, using the marching cubes algorithm. Due to the excessive number of nodes and large number of bad quality elements, which are common features in mesh-based models, a CAD-based approach has been developed: the meshes are rebuilt as a solid using a procedure of semi-automatic surface creation with Non-Uniform Rational Basis Splines (NURBS). Using the resulting smooth surface, a quality mesh with four-node tetrahedral elements is generated using ABAQUS packages. In this paper, we did not focus on the anatomical geometry reconstruction that do not constitute the core of this contribution.

C. Patient-specific biomechanical model of the respiratory system

First, we describe the biomechanics of the respiratory system including the mechanical behavior of the lung tissue, the diaphragm muscles and rib kinematics, as well as the boundary conditions which are used in our model. Then, we present the developed automatic tuning and optimization algorithm, to calculate the personalized physiological parameters diaphragm forces and the lungs pressure during the whole respiratory cycle. 1) Respiratory system and boundary conditions : Our anatomical model of the respiratory system comprises the lungs, tumor, thorax (rib kinematics), diaphragm behavior, skin and mediastinum. The boundary conditions (BC) are inferred from the anatomy and identified by medical experts. For the diaphragm, we have considered the direction of muscle forces to the direction of muscle fibers (Fig. 5(a)). The force is applied to the muscular part of the diaphragm and simple homogeneous Dirichlet boundary condition is applied to the lower part of the diaphragm and the Lagrange multiplier method is used for the contact model. To simulate the sliding of the lungs, a surface-to-surface contact model is applied to the lung-chest cavity (Fig. 5(b)). The frictionless contact surfaces are used to simulate the pleural fluid behavior. The rib cage kinematics have been applied as displacement boundary conditions during the whole respiratory cycle. Then, we have used the personalized compliance to apply pressure boundary conditions to the surface of the lung model. The lung model is monitored directly by simulated actions of the breathing muscles: the diaphragm and the intercostal muscles (the rib cage). In our work the amplitude of the lung pressure and diaphragm force are patient-specific, determined at different respiratory states by an optimization framework based on inverse FE analysis methodology, using lung volume variations. We have segmented the respiratory system at end expiration (EE, the reference state). Also, we have segmented lungs at all states for a full cycle (10 states from DIR-Lab Dataset [START_REF] Castillo | Four-dimensional deformable image registration using trajectory modeling[END_REF]). Based on personalized pressurevolume curve (Fig 5(c)), we can estimate the lung pressure for each patient as:

P Lung (i) = -K V i -V 0i V 0i = - E 3(1 -2ν) V i -V 0i V 0i (3)
where P Lung (i) presents the lung elastic pressure at each intermediate respiratory state. The K, E and ν are bulk modulus, Young Modulus and Poisson ratio respectively. V 0i is the initial lung volume at each intermediate respiratory state and V i represents the current lung volume.

2) Rib kinematics: In our study, the thoracic cage model includes all skeletal structures: ribs and associated costal cartilages, the sternum and the thoracic vertebrae (Fig. 5(a)). The ribs are considered as articulated rigid bodies. The motion of a rigid body between two positions can be characterized by an Euclidean transformation. Several authors have investigated the ribs kinematics. In [START_REF] Wilson | Respiratory effects of the external and internal intercostal muscles in humans[END_REF], human ribs displacement were studied and the transformation parameters have been defined using planes attached to the ribs. In our previous works [START_REF] Didier | A chest wall model based on rib kinematics[END_REF], [START_REF] Ladjal | Appropriate Biomechanics and kinematics Modeling of the respiratory System: Human Diaphragm and Thorax[END_REF], we have presented a methodology to study rib kinematics, using the finite helical axis method. The idea is to predict, from the transformation, the evolution the ribs position and orientation parameters with time. Each rib transformation parameters are automatically computed between the initial and final state (Fig. 5(b)). As the rotation angle is small, a linear interpolation between the two states is applied to predict the rib motion at any intermediate breathing states. For more details about finite helical axis method, one can refer to [START_REF] Didier | Mechanical role of pleura on lung motion during breathing[END_REF], [START_REF] Didier | A chest wall model based on rib kinematics[END_REF].

3) Lung and diaphragm Behaviors: Based on the results published [START_REF] Giroux | Patient-Specific Biomechanical Modeling of the Lung Tumor for Radiation Therapy[END_REF], the Saint-Venant Kirchhoff (SVK) hyperelastic model can describe the mechanical behavior of the lung tissue. Then, we have used the personalized pressure-volume curve to apply pressure boundary conditions to the surface of the lung model based on our computational framework presented in section II-C4. The muscles of the diaphragm are skeletal muscle type, and their action has the effect of shortening or lengthening the muscles along their fiber direction. The muscles of the diaphragm are considered as compressible isotropic hyperelastic materials (SVK) [START_REF] Ladjal | Appropriate Biomechanics and kinematics Modeling of the respiratory System: Human Diaphragm and Thorax[END_REF], [START_REF] Ladjal | Biomechanical Modeling of the Respiratory System: Human Diaphragm and Thorax[END_REF] and the tendon part is considered as linear elastic behavior (Fig. 5(a)). For an isotropic elastic material, the elastic energy, noted W , can be written as:

W (E) = λ 2 ( trE ) 2 + µ tr ( E 2 ) ( 4
)
where E is the Green-Lagrange strain tensor, λ and µ are the Lame's coefficients. For small deformations, the Green-Lagrange strain tensor is linearized into the infinitesimal strain tensor (ϵ ϵ ϵ):

ϵ ϵ ϵ = 1 2 (grad U + grad T U ) ( 5
)
where U is the displacement field. The Hooke's law or Hookean materials is the relation between the stress tensor (σ) and the strain tensor (ϵ ϵ ϵ), can be written as:

ϵ ϵ ϵ = 1 + ν E σ - ν E tr(σ)I (6)
With E: Young's modulus, ν: Poisson's ratio, and I is the identity matrix. Other expression can be written, introduce Lame's constants, where µ:shearing coefficient and λ: compression coefficient.

µ = E 2(1 + ν) λ = ν E (1 -2ν)(1 + ν) (7)
For the non-linear hyper elastic behavior, the Saint-Venant Kirchhoff law extends the Hooke's law for large displacement. The relation between the second Piola-Kirchhoff stress tensor, the Green-Lagrange strain tensor and the rate effect leads to:

S = λ (tr E) I + 2 µ E (8)
S is the 2nd Piola-Kirchhoff stress tensor, can be derived from the relation: S = ∂W ∂E . For the dynamic simulation of the respiratory system, the equation of motion of a vertex v of the model can be written:

M v { üv } + γ v { uv } + ∑ τ ∈Vv ( F int v ) = F v ext ( 9
)
Where M v , γ v are respectively the mass and damping coefficients of each vertex. The V v is the neighborhood of vertex v (i.e. the tetrahedra τ containing node v). The F int v are the internal forces calculated by FE method and the F v ext are the imposed forces calculated by our developed algorithm based on inverse FE. To solve the dynamic system, and for more stability, the implicit finite scheme has been chosen. In our work, the mass density of each tissue is patientspecific, calculated and determined directly from CT scan images, based on density mapping algorithm defined and developed in our last work [START_REF] Manescu | Four-dimensional radiotherapeutic dose calculation using biomechanical respiratory motion description[END_REF]. First, organs tetrahedral meshes are generated from segmented CT scanner images. Next, the Hounsfield values issued from CT scanner images are converted into density values that are mapped to these mesh nodes, respecting the principles of mass conservation (Fig. 3). For more information related to density mapping algorithm, one may refer to [START_REF] Manescu | Human liver multiphysics modeling for 4d dosimetry during hadrontherapy[END_REF], [START_REF] Manescu | Four-dimensional radiotherapeutic dose calculation using biomechanical respiratory motion description[END_REF].

4) Automatic tuning and optimization:

Physiological parameters of the respiratory system are patient-specific and vary from one patient to another. Based on the planning 4D CT, we calculate firstly the initial lung elastic pressure at each intermediate respiratory state, using the equation (3) as a function of the volume and mechanical properties (Young modulus and Poisson ratio) of the lungs. Then, based on our developed model using the inverse finite element analysis methodology, we recalculate and optimize the personalized lung pressure/volume curve by minimizing the lung volume errors at each intermediate respiratory state, between the segmented volume (from 4D CT) and the simulated volume (calculated by biomechanical simulation). The schematic of the computational framework developed for this study is illustrated in (Fig. 4). This framework is designed to calculate automatically the personalized diaphragm forces and the lungs pressure during the whole respiratory cycle and for each intermediate respiratory state based on inverse finite element optimization. The optimization loops were achieved with a program coded in Python coupled within finite element Abaqus solver. Besides data acquisition (4D scan images), there are three successive steps involved in this framework. The first step is the estimation of the initial parameters:

• the initial geometry is segmented and the organs tetrahedral meshes are generated from segmented CT images at EE, • through a personalized compliance (pressure/volume curve), we have calculated and applied the negative lung pressure (not directly measurable) at EE and at different respiratory states, • the material density values calculated from CT data are mapped on the vertices of the tetrahedral meshes in such a way that the mass of each organ is preserved. • we apply the diaphragm forces, these forces are applied on each node of the muscular part of the diaphragm, parallel to muscle orientation. • the kinematic model of the rib cage based on the finite helical axis method (FHAM), which play a significant role in the pleura outer-surface motion and therefore in the lung motion. The second step, the biomechanical simulation (simulation loop) based on finite element analysis of the respiratory system is generated using the current parameters. The outputs in this step are:

• the simulated lung volume,

• the diaphragm forces,

• the displacement field of all organs,

• the 3D lung tumor position Then, in the third step, based on an iterative algorithm (the optimization loop), in which the diaphragm forces and lung pressures are adapted and corrected for each phase of respiration. The optimization algorithm starts with initial diaphragm forces and updates them iteratively. Using the estimated diaphragm force, the simulation loop calculates the deformations of the respiratory system (lung motion) using a finite element model. The simulation loop, then returns the estimated parameters of the diaphragm forces and the simulated lungs volume. The mean absolute error (MAE) was used as cost function to minimize the the errors between the simulated lung volume (V sim ) mesh and the segmented lung volume (V exp ) mesh extracted from the CT scan images of the same phase of respiration. If the cost function does not converge, an estimation of the parameters are updated, and the new parameters are fed back into the simulation loop.

III. RESULTS AND EXPERIMENTAL VALIDATION

In order to demonstrate the validity of our patient specific biomechanical model, a quantitative and qualitative analysis of simulations were conducted. We have compared the results of a simulated motion with the experimental data provided directly from the 4DCT scan images DirLab database [START_REF] Castillo | Four-dimensional deformable image registration using trajectory modeling[END_REF]. Fig. 6 illustrates the different specific compliances (pressurevolume curve) for each patient, calculated and identified at each state directly from 4D CT scan images. Then, these compliances are used an input in our biomechanical patient specific model to simulate the full breathing cycle.

In our FE simulation, we define the simulation time for the inspiration phase is 2 seconds and for the expiration phase is 3 seconds. The mechanical properties and behaviors of the different organs used in our simulations are settled in the Table .I. Fig. 7 shows the total deformation and the maximum displacement components of the lungs and diaphragm during breathing. We can observe the maximum displacement of the diaphragm on the right-posterior (RP) and left-posterior (LP) sides. It is also possible to notice a slightly larger (RP) side motion than (LP) side motion, in concordance with the physiological anatomy. For the lungs deformation, the maximum displacement occurs in the posterior region along the superior-inferior (SI) direction (diaphragm direction).

A. Anatomical landmarks evaluation at end of inhalation (EI) and the end of exhalation (EE)

To measure the precision of our biomechanical simulation, we have evaluated the model accuracy on five selected lung cancer patients, from DIRLab Dataset [START_REF] Castillo | Four-dimensional deformable image registration using trajectory modeling[END_REF], where the entire thorax was visible and with small and large breathing amplitudes(Table.II). We have compared the finite element simulation results on 300 landmarks, at end inspiration (EI), end expiration (EE) states, and 75 landmarks at each intermediate respiratory state, obtained by manual delineations [START_REF] Castillo | Four-dimensional deformable image registration using trajectory modeling[END_REF]. Table .IV illustrates the comparative study between our FE simulation results and the anatomical landmarks displacement vectors for five patients. The average error for anatomical landmarks inside the lung at end inspiration (EI) and end expiration (EE) states for five patients is less than 2.5 ± 1.5mm. In addition to the need for accuracy, we were also evaluated the % errors distribution on the 300 landmarks for the same five lung cancer patients, and the mean errors for anatomical landmarks in the 2 lobes of the left lung and 3 lobes of the right lung, along the AP (Anterior/Posterior), RL (Right/Left), and SI (Superior/Inferior) directions. Fig. 8 and Table .V summarize the histogram of landmark errors distribution and the mean landmark errors in the 2 lobes of the left lung and the 3 lobes of the right lung in the LR, AP and SI directions between EE and EI phases, respectively. The results show that more than 95% of mean landmark error values are less than 3mm for each lobe in the LR, AP and SI directions.

B. Anatomical landmarks evaluation at intermediate states between EI and EE

The performance of the proposed biomechanical model has been evaluated by comparing the simulation results with ground truth (CT images) on 75 landmarks at five intermediate states between EI and EE. displacement vectors for five patients. In our simulation, we have obtained an average mean error less than 2.0 ± 1.3mm. These results show that the developed physiological model coupled with the personalized lung-pressure/diaphragm-force optimization algorithm of the respiratory system is in a good agreement with the experimental data, produces more accurate predictions with lower errors compared to other works ( [START_REF] Fuerst | Patient-Specific Biomechanical Model for the Prediction of Lung Motion from 4D CT Images[END_REF], [START_REF] Li | Biomechanical model-based 4DCT simulation[END_REF]) that used the same datasets.We have also evaluated the influence of the variation of lung tissue Young's modulus and Poisson's ratio commonly used in the literature [START_REF] Nakao | Simulating lung tumor motion for dynamic tumor-tracking irradiation[END_REF], [START_REF] Kk | Accuracy and sensitivity of finite element model-based deformable registration of the prostate[END_REF], [START_REF] Villard | Interactive simulation of diaphragm motion through muscle and rib kinematics[END_REF], [START_REF] Werner | Patient-specific finite element modeling of respiratory lung motion using 4D CT image data[END_REF]. Young modulus can vary between 0.1kPa to 10kPa and Poisson's ratio between 0.2 to 0. ratio=0.49, is three times larger than the lung landmarks error obtained by Poisson's ratio=0.33.

C. Tumor tracking and impact of the rib kinematics on lung tumor motion

In order to evaluate the impact of the rib kinematics on lung tumor motion, we have evaluated the tumor motion and 3D trajectory identified in 4D CT scan images on two selected patients (6 and 10) from the same Dataset DIR-Lab [START_REF] Castillo | Four-dimensional deformable image registration using trajectory modeling[END_REF], where the tumor location is visible. These trajectory are compared and evaluated with the trajectory obtained by FE simulation, during the whole cycle of breathing (10 phases between the EI and EE). To track the lung tumor movement in 4D CT images, we have used the affine registration applied to the segmented lung tumor volume at different respiratory states (Fig. 9). Then, the accuracy of the proposed tumor tracking method is evaluated by comparing and calculating the Hausdorff distance between the 3D mesh surface of the segmented tumor and predicted FE lung tumor, including or not the rib kinematics. Fig. 11 and Fig. 10 show the mean errors ± standard deviation of lung tumor position and comparison study between the 3D hysteresis trajectories of the lung tumor during the whole cycle of the breathing without and within the rib kinematics compared to the trajectory calculated directly from 4D CT images. The figures demonstrate that our patient specific biomechanical model for lung tumor position estimation is accurate (less than 3 mm). However, it is important to note that the results for these two patients are slightly better with the rib kinematics but the difference is not significant during the whole cycle (for patient 6 and 10). This is because the respiration for these patients is mainly diaphragmatic. In addition, another investigation conducted to study the impact of rib kinematics on the lung tumor motion prediction. This investigation was conducted by varying the ribs displacement by 10%, 25% and 50% during the whole breathing cycle. Fig. 12 illustrates the impact of the rib kinematics on lung tumor trajectory for the patient 6 and patient 10, and the robustness of the proposed approach based on biomechanical patient specific model of the respiratory system.

IV. CONCLUSION AND FUTURE WORK

In this paper, a patient specific biomechanical model of the respiratory system was proposed for lung tumor tracking for the whole respiratory cycle. Our results are quite realistic compared to the 4D CT scan images. We can observe that the proposed physically-based FE model is able to predict correctly at high precision the lung tumor displacements including the personalized and physiological lung-pressure/diaphragm-force and boundary conditions of the organs. This could be a potential tool to provide valuable location-specific tumor motion information for medical physicians to reduce the margins between clinical target volume (CTV) and planning target volume (PTV). One of the limitations of our work is the 3D reconstruction organs delineation and treatment of the multiple organs mesh, requiring time consuming, manual operations for each patient. In order to avoid manual contouring and 3D geometry segmentation for different organs, and to reduce the computational costs without lowering the quality, we plan to develop and use a realistic atlas-based 3D shape reconstruction of the respiratory system based on statistical training or machine learning, to get a fast and automatic patient-specific model adaptable taking into account complexities such as geometrical irregularities and organs contacts. Currently, we are working on real time computation of our patient-specific model. The model should be controlled by external non-invasive surrogates, to predict patient's lung tu- mor motion, during treatment. In this context, it is important to investigate the accuracy and reproducibility of RPM versus surface imaging to find an appropriate surrogate for 3D rib motion. In the future, the model should also be validated pre-clinically on a realistic deformable anthropomorphic lung phantom (LuCa). The results of our model can greatly contribute to motion compensation algorithms for diagnosis, beam therapy, dose calculation or online imaging systems.

Fig. 1 .

 1 Photon or hadron beam

Fig. 2 .

 2 Fig. 2. Theoretical lung compliance (pressure/volume curve) for the whole respiratory cycle

Fig. 3 .Fig. 4 .

 34 Fig. 3. Density mapping from regular voxel grid to our irregular tetrahedral elements: the mass of a tetrahedral element equals the sum of the masses of volumes of intersection between the tetrahedron and the grid of voxels: m(T k ) = m(I 1 k ) + m(I 2 k ) + m(I 4 k ) + m(I 5 k ) + m(I 6 k ) + m(I 8 k ) + m(I 9 k ).

Fig. 5 .

 5 Fig. 5.The boundary conditions (BC) of our patient specific biomechanical model of the respiratory system, including (a) diaphragm behavior: heterogeneous tissue with the muscles in peripheral part and the tendon in central part, (b) rib kinematics using the finite helical axis method, and (c) personalized compliance calculated for each patient.

Fig. 6 .

 6 Fig. 6. Patient specific compliances for five patients from DirLab database, calculated and identified at each respiratory state directly from 4D CT scan images.

Fig. 7 .

 7 Fig. 7. Qualitative analysis of patient specific biomechanical simulation; lungs and diaphragm deformations during the whole cycle of breathing for two patient P1 and P6.

Fig. 8 .

 8 Fig. 8. The histogram of landmark errors distribution for five patients from DirLab database [30] along the AP (Anterior/Posterior), RL (Right/Left), and SI (Superior/Inferior) directions, at end inspiration (EI).

Fig. 9 .Fig. 10 .

 910 Fig. 9. Evaluate lung tumor motion: rigid registration applied to the segmented lung tumor volume at different respiratory states

Fig. 11 .

 11 Fig. 11.Mean errors ± standard deviation of lung tumor position during the whole cycle of breathing (10 phases between the EI and EE) between the trajectory issued from 4D CT images compared to the trajectory calculated by biomechanical finite element simulation coupled with the lungpressure/diaphragm-force optimization for two patients P6 and P10. The error bars are always higher or same for the without ribs results.

10 Fig. 12 .

 1012 Fig. 12. 3D lung tumor trajectory (in mm) issued from 4D CT images compared to the trajectory calculated by biomechanical finite element model including rib kinematics displacement changed by 10%, 25% and 50% of the rib displacement calculated from CT scan images.

TABLE I MECHANICAL

 I PROPERTIES OF BREATHING SYSTEM: E YOUNGS MODULUS AND ν POISSON RATIO (FROM[START_REF] Kk | Accuracy and sensitivity of finite element model-based deformable registration of the prostate[END_REF],[START_REF] Villard | Interactive simulation of diaphragm motion through muscle and rib kinematics[END_REF],[START_REF] Ladjal | Biomechanical Modeling of the Respiratory System: Human Diaphragm and Thorax[END_REF],[START_REF] Giroux | Patient-Specific Biomechanical Modeling of the Lung Tumor for Radiation Therapy[END_REF]), ρ VOLUMETRIC DENSITY (CALCULATED FROM CT IMAGES) .

	Tissues	Mechanical E	ν	ρ
		behavior	(MPa)		(t/mm 3 )
	Lungs	HSVK	3.74x10 -3 0.3	3.0x10 -10
	Lung tumor	LE	49	0.4	1.5x10 -9
	Mediastinum LE	5.87x10 -3 0.4	1.0x10 -10
	Diaphragm	HSVK	5.32	0.33 1.0x10 -9
	muscle
	Diaphragm	LE	33	0.33 1.0x10 -9
	tendon
	Ribs	LE	5000	0.4	1.5x10 -9
	Cartilage	LE	49	0.3	1.0x10 -9
	Sternum	LE	11500	0.3	1.5x10 -9
	Vertebra	LE	9860	0.4	1.5x10 -9
	Skin	LE	5.32	0.4	1.0x10 -6

TABLE III SURFACE

 III ERRORS FOR THE DIAPHRAGM AND LUNGS SURFACES AT END INSPIRATION EI BETWEEN FEM SIMULATION SURFACE AND SEGMENTED SURFACES EXTRACTED FROM CT IMAGESFOR FIVE PATIENTS FROM DIRLAB DATABASE[START_REF] Castillo | Four-dimensional deformable image registration using trajectory modeling[END_REF].

	Patients	Diaphragm	Lungs
	End Inspiration (EI) Mean ± SD (mm)	Mean ± SD(mm)
	Patient 1	1.8 ± 1.5	2.0± 1.5
	Patient 4	2.0 ± 1.6	2.2 ± 1.7
	Patient 6	1.9 ± 1.6	2.1 ± 1.5
	Patient 9	2.1 ± 1.5	2.3 ± 1.6
	Patient 10	2.2 ± 1.1	2.5 ± 1.8

TABLE IV AVERAGE

 IV LANDMARKS LUNG ERROR (MM) DURING EXHALATION AT DIFFERENT RESPIRATORY STATES: THE FIRST STATE T00, THE END INSPIRATION (T50), THE END EXPIRATION (T10).

	Patients	Volume			Mean ± SD (mm)			Mean
		between EE and EI (%)	T10	T20	T30	T40	T50	All states
	Patient 1	12.8	2.1 ± 1.5 2.2 ± 1.2 2.1 ± 1.6 1.6 ± 1.4 1.1 ± 0.8	1.8± 1.3
	Patient 4	14.6	2.3 ± 1.2 2.5 ± 1.3 2.1 ± 1.2 1.8 ± 1.2 1.5 ± 1.2	2.0±1.2
	Patient 6	25.4	2.4 ± 1.5 2.3 ± 1.2 2.0 ± 1.6 1.9 ± 1.6 1.4 ± 1.1	2.0±1.3
	Patient 9	13.4	2.3 ± 1.4 2.2 ± 1.1 2.1 ± 1.3 1.8 ± 1.4 1.3 ± 0.9	1.9± 1.2
	Patient 10	16.7	2.1 ± 1.5 2.2 ± 1.2 2.1 ± 1.6 1.6 ± 1.5 1.1 ± 0.8	1.8±1.3

TABLE V MEAN

 V LANDMARK ERRORS IN THE 2 LOBES (INFERIOR (INF), SUPERIOR (SUP)) OF THE LEFT LUNG AND 3 LOBES ( (INFERIOR (INF), MIDDLE AND SUPERIOR (SUP)) OF THE RIGHT LUNG, IN THE AP (ANTERIOR/POSTERIOR), RL (RIGHT/LEFT), AND SI (SUPERIOR/INFERIOR) DIRECTIONS BETWEEN EE AND EI PHASES.

	Patients					Patient 1
	Lobes		SUP			Middle			INF
	Directions	L/R	S/I	A/P	L/R	S/I	A/P	L/R	S/I	A/P
	Landmarks error (right lung)(mm)	1.3	1.3	1.72	1.5	1.6	1.8	1.8	2.2	1.9
	Landmarks error (left lung)(mm)	1.1	1.1	1.5	-	-	-	1.3	1.6	2.3
	Patients					Patient 4
	Lobes		SUP			Middle			INF
	Directions	L/R	S/I	A/P	L/R	S/I	A/P	L/R	S/I	A/P
	Landmarks error (right lung)(mm)	1.4	1.5	1.9	1.4	1.2	2.3	1.6	2.1	2.5
	Landmarks error (left lung)(mm)	1.1	0.9	1.6	-	-	-	1.3	1.6	1.8
	Patients					Patient 6
	Lobes		SUP			Middle			INF
	Directions	L/R	S/I	A/P	L/R	S/I	A/P	L/R	S/I	A/P
	Landmarks error (right lung)(mm)	0.77 0.11 1.72	0.23 0.55 1.04	1.35 0.60 1.15
	Landmarks error (left lung)(mm)	0.84 0.95 1.38	-	-	-	0.69 0.65 1.03
	Patients					Patient 9
	Lobes		SUP			Middle			INF
	Directions	L/R	S/I	A/P	L/R	S/I	A/P	L/R	S/I	A/P
	Landmarks error (right lung)(mm)	1.3	1.3	1.7	1.5	1.8	2.0	1.7	2.1	2.1
	Landmarks error (left lung)(mm)	1.2	1.1	1.6	-	-	-	1.5	1.8	1.9
	Patients					Patient 10
	Lobes		SUP			Middle			INF
	Directions	L/R	S/I	A/P	L/R	S/I	A/P	L/R	S/I	A/P
	Landmarks error (right lung)(mm)	1.9	1.8	1.7	1.9	2.1	2.4	2.1	2.3	2.6
	Landmarks error (left lung)(mm)	1.6	1.8	1.8	-	-	-	1.9	1.9	2.3

  Table.VI shows the

TABLE VI COMPARISON

 VI BETWEEN OUR BIOMECHANICAL PATIENT SPECIFIC MODEL RESULTS AND THE RESULTS FROM FUERST ET AL.[START_REF] Fuerst | Patient-Specific Biomechanical Model for the Prediction of Lung Motion from 4D CT Images[END_REF] AND LI ET AL.[START_REF] Li | Biomechanical model-based 4DCT simulation[END_REF], ON PATIENT 6, PATIENT 9 AND PATIENT 10 ISSUED FROM DIRLAB DATASET[START_REF] Castillo | Four-dimensional deformable image registration using trajectory modeling[END_REF].

	References Cycle	Organs	Boundary	P6	P9	P10
			Modeled	conditions	(mm)	(mm)	(mm)
	[22]	half		4 pressure zones (sub diaphragm)
					3.27 ±1.0	2.97 ±1.0 2.83 ±1.0
			No ribs behavior	16 pressure zones (thoracic)
	[21]	half	Only lungs organ	One CT and uniform pressure	-	3.2±1.4	4.4±2.9
					-
			All respiratory system personalized pressure-volume curve
	Our model	Full	Thorax	Rib kinematics	2.0±1.3	1.9± 1.

This paper was presented in part at the MICCAI 2017. This research is supported by the LABEX PRIMES (ANR-11-LABX-0063), within the program Investissements dAvenir(ANR-11-IDEX-0007) operated by the French National Research Agency (ANR) and by France Hadron. H. Ladjal and B. Shariat are with LIRIS