Abstract
To reduce the side effects and to improve the efficiency of radiation therapy in lung cancer, a pinpoint radiation therapy system is under development. In the system, the movement of lung tumor during breathing could be estimated by employing a suitable numerical modeling technique. This paper presents a gridless numerical technique called Moving Particle Semi-implicit (MPS) method to simulate the lung deformation during breathing. The potential of the proposed method to employ in the future pinpoint radiation therapy system has been explored. Deformation of lung during breathing was dynamically tracked and compared against the experimental results at two different locations (upper lobe and lower lobe). Numerical simulations showed that the deformation of lung surface ranged from less than 4 mm to over 20 mm depending on the location at the surface of lung. The simulation showed that the lower section of lung exhibited comparatively large displacement than the upper section. Comparing with the experimental data, the lung surface displacement during inspiration process was predicted reasonably well. Comparison of numerical prediction with experimental observations showed that the root mean squared error was about 2 mm at lower lobe and less than 1 mm at upper lobe at lung surface.
1. Introduction
Death related to cancer is increasing annually in Japan. In 2005, thenumber of deaths related to cancer was approximately 326000. In terms with cancer site, lung was leading with 23% for male. To treat lung cancer, radiation therapy is widely used. In radiation therapy cancer cells are irradiated with X-ray beam while limiting injury to surrounding healthy tissues. However, the issue of limiting injury to healthy neighboring tissues, in the treatment of lung cancer, is more complicated as lung undergoes volumetric expansion and contraction during breathing. Moreover, displacement of lung surface is not uniform throughout and depends on its location. Hence, a proper estimation of lung displacement at the specific location is very crucial for planning radiation therapy to treat lung cancer.
To reduce the side effects and to improve the efficiency in radiation therapy, our group is developing a pinpoint radiation therapy system [1] for lung cancer (Figure 1). In this system, a numerical method could be applied to estimate the deformation of lung and tumor during breathing. The displacement data obtained from the simulation would be provided to the robotic arm holding an accelerator which produces a very thin X-ray beam. The accelerator producing a very thin X-ray beam, hence, would be able to move along with the movement of tumor in the lung, consequently, reducing the chance of irradiating healthy cells neighboring the tumor.
Soft tissues like lung, liver, kidney are a complex structure and the major constraint of soft tissue modeling is its correct representation of deformability. Generally, lung is assumed to be elastic material whose mechanical properties are described with Young modulus and Poisson’s ratio.
Mori et al. did study on lung tumor displacement by using four-dimensional Computed Tomography (CT) techniques on 14 patients. They noticed that the average respiratory cycle for all patients was 3.8 seconds with a standard deviation of 0.6 seconds. Their results showed that the displacement in each lobe differed in Left-Right (LR), Anterior-Posterior (AP), and Superior-Inferior (SI) direction. Moreover, the tumor in lower lobes (segment 6–segment 10) showed greater SI displacement than that in other lobes. The total displacement in the lower lobe was around twice the middle or upper lobe [2].
Nakao et al. treated lung as an elastic object and analyzed the deformation of the lower part of the lung, during breathing, based on linear finite element method [3]. Although their numerical simulation matches quite well with the experimental results, their proposed model was not able to simulate the variation in lung deformation depending on the location at lung’s surface.
Didier et al. considered lung tissue as homogeneous, isotropic and its constitutive equation defined by a linear elastic behavior and simulated lung deformation by using finite element method [4]. The maximum error with this technique was around 6 mm. However, for clinical purpose error less than 4 mm is considered to be acceptable. Moreover, the ability of their proposed model to dynamically track the deformation during breathing was not demonstrated.
Villard et al. used finite element method to simulate lung deformation during breathing and measured the influence and relevance of mechanical parameters [5]. However, the potential of the proposed model to employ in clinical application was not cleared as the simulated results still needed to be compared against clinical data. They suggested that Young modulus can be chosen arbitrarily; however, Poisson’s ratio plays a significant role and may not be chosen arbitrarily for correct prediction. In literatures, Poisson’s ratio ranging from 0.25 to 0.47 has been chosen to model lung deformation during breathing [5, 6].
Moreover, in the above approaches of modeling lung deformation, grid generation was necessary. Since soft biological tissues, like lung, exhibit the property of high deformation under application of external pressure, grid technique like finite element method might suffer from grid distortion at some point. Hence, a gridless method called Moving Particle Semi-implicit (MPS) method is proposed in this paper as an alternative method to simulate lung deformation during breathing. In gridless method, like MPS, grid generation is not required. The elastic solids are represented by a collection of finite number of particles. During the simulation, the changes in physical properties of particles while they undergo displacement are calculated. Hence the problem of grid distortion while modeling high deformation is not of the concern.
In 1996, Koshizuka and Oka [7] introduced the Moving Particle Semi-implicit (MPS) method for fluid dynamics. In 2001, Chikazawa et al. presented the particle method for elastic and visco-plastic structures and fluid-structure interaction [8]. In their study particle method was proposed for thick elastic and visco-plastic structures based on the concept of MPS which provided particle interaction model for differential operators. Dynamic simulation of elastic solids was performed by applying MPS method [9, 10]. They applied the modeling technique to simulate large deformation and fracture of elastic solids of different material.
Three-dimensional elastic analysis and large deformation of soft material were performed by [11]. Kondo et al. applied symplectic scheme to three-dimensional elastic analysis using MPS method [12]. They found that the scheme conserved energy and momentum. Recently, Suzuki and Koshizuka introduced Hamiltonian particle method for nonlinear elastodynamics [13]. The particle method inherits the symplectic structure possessed by their governing equation for both compressible and incompressible materials.
2. Materials and Methods
2.1. Governing Equations
As in fluid mechanics, the governing equation of motion in elastic solid can be written as [9, 10]
Here,
where is density of material, an isotropic pressure which is obtained in particle location , unisotropic components of stress tensor which is obtained between the particles, and the Lame constant, external force acting on the particles, Young modulus, and Poisson’s ratio.
In MPS method, the differential operators on the right-hand side of (1) is replaced by particles interaction and calculated by using Lagrangian approach. The external force can be applied directly at particles. The procedure is explained briefly in the Numerical Model section.
2.2. Numerical Model
In MPS method, as shown in Figure 2, elastic solid is represented as a collection of finite number of particles. Particles behave as if they are connected by normal and shear springs which are governed by the parallel (normal) and perpendicular (shear) displacement of the particles.
As an initial condition, coordinate, velocity, angle, and angular velocity are attributed to each particle as degrees of freedom. In order to simulate large deformation, displacement is not given as degrees of freedom but it is calculated from the coordinates. Hence at each time step, each particle has two systems of coordinates; one of initial nondeformed stage and the other one is the coordinate of later or deformed stage.
The positions and the velocities are vectors. The particle number density and the pressure are scalars.
As show in Figure 2, particles interact with its neighboring particles covered with a weight function as.
Since the area that is covered by this kernel function is bounded, a particle interacts with a finite number of neighboring particles (Figure 2). The interaction radius is determined by .
The global accuracy in particle methods depends on the initial distribution of particles. It is recommended that a finer particle distribution model and approximately large particle number to be used according to the available computer resources to obtain better computational results [14, 15].
2.2.1. Displacement Calculation
The initial position vector is represented as . As the particles move to a new position, the new position vector is represented as
Here, and represent particle and its neighbor . and represent position of particles and .
The rotated vector of the initial inter particle location vector is evaluated as in
where orientation matrix is represented by using unit quaternion [10]. The displacement vector () is calculated as
2.2.2. Calculation of Isotropic Pressure
Isotropic pressure is calculated as (2) in the particle method [10]. Strain term in (2) is calculated as
Here, is the number of spatial dimension, is the initial weight function, and is the initial particle number density. Particle number density at coordinate where particle is located is defined as
2.2.3. Calculation of Unisotropic Stress
Unisotropic stress (3) is calculated by using
2.2.4. Calculation of Divergence Terms and Time Integration
The divergence terms in the right-hand side of (1) is calculated by using Lagrangian approach
The velocity and the corresponding coordinates at a new time step are then updated by using a fourth-order Runge Kutta scheme.
Rotation of particle is calculated by updating the angular velocity and rotation angle [10].
Time step is controlled in the computation to satisfy the following Courant condition:
where is the original distance between two particles and the maximum velocity among all particles. The time step is set at the beginning of computation and kept constant throughout the computation.
2.3. Lung Deformation and Numerical Experiments
Lungs are divided into lobes and each lobe is composed of segments. Right lung, which is larger than left lung, consists of 3 lobes which are subdivided in 10 segments (S1–S10). Left lung consists of 2 lobes which are subdivided in 8 segments (S1–S8).
Respiration process can be both voluntary and involuntary. At the beginning of inspiration process, diaphragm moves downward and intercoastal muscles around the ribcage move forward. Hence, a negative pressure develops around the lung surface. To compensate to that negative pressure, air rushes into the lungs.
Alveoli are the final branch of the respiratory tree where gas exchange takes place and fill the lung volume. During the inspiration phase of breathing, alveolar expansion takes place. Generally, alveoli are better approximated by spherical shape. Hence during inspiration, the direction of alveoli expansion is radial outward. Hence, the direction of expansion of lung surface is also radial outward. However, at the lower lobe, the displacement can be affected by the diaphragm motion as well.
In the numerical modeling, to simulate the behavior of lung expansion due to alveoli expansion, a 3D lung model represented by particles was generated. To create a 3D lung model, CT images of a lung was used. The 2D CT images were then converted to a 3D model (Figure 3) by an image processing and measurement software called 3D-Doctor manufactured by Able Software Corporation. The 3D model of lung were then converted to a particle model by considering each voxel of the 3D model as a particle. In the numerical simulation, the lung was represented by a set of 14 298 uniformly distributed particles or mass points along lung’s contour. The particles were distributed uniformly with the spacing of 1 mm.
(a)
(b)
A radial outward force was applied to each mass point at the surface of the lung. Forces acting on the points generated displacements.
In this study, Poisson’s ratio () was selected as 0.45. Young Modulus () and the density of the material were selected as 0.01 kPa and 700 kg/m3, respectively.
3. Results
Displacement of lung surface during breathing depends largely on its location. To simulate the deformation of lung during breathing, a radial outward force at each particle was applied and the displacement was tracked until the final stage was reached. In the numerical simulation, the movement of lung in LR and AP direction was restricted to 2 mm and 4 mm to assimilate the influence of rib cage around lung [2]. However, in the SI direction, deformation was allowed without any restrictions. Numerical results of lung deformation along SI direction have been compared against experimental results. The numerical results demonstrated that during breathing, lower lobe exhibited higher displacement along SI direction compared to the upper lobe.
The overall deformation of lung surface by the end of inspiration is shown in Figure 4. The simulation showed that the deformation at the lung surface, in SI direction, was higher in lower part compared to the upper part. Deformation was ranged from lower than 4 mm to over 20 mm depending on the location at the surface of lung.
The numerical prediction of lung deformation, at positions S6 (lower lobe) and S3 (upper lobe), along SI direction has been compared against experimental results obtained by tracking the movement of gold marker placed at those locations.
During inspiration, the deformation rate of lung surface is slower at first and then it picks up after certain time interval. This phenomenon occurs through out the lung surface. However, deformation rate depends on the location. This process has been well predicted by the numerical method.
The position S6 at the lung surface in the numerical simulation is depicted in Figure 5. In the figure, the top figure shows the simulated results of lung at the beginning of inspiration and the bottom figure shows the deformed lung at the end of inspiration.
(a)
(b)
The graph of displacement versus time at location S6 is depicted in Figure 6. The experimental result showed that during the first 0.2 seconds of inspiration, the displacement was less than 1 mm. After that the deformation rate increased gradually and by 1.3 seconds, or the end of inspiration, the S6 location was displaced, along SI direction, by a little over 14 mm. The graph in Figure 6 shows that the numerical prediction generally followed the experimental result with reasonable accuracy.
The numerical result showed that by the end of inspiration process the position at S6 was displaced by a little over than 14 mm. Though the maximum displacement by the end of inspiration was accurately predicted, by the numerical simulation, during the other phases of inspiration, the displacement was slightly under predicted. The maximum error obtained from the simulation was mm which occurred at 0.7 seconds after the beginning of inspiration. Comparing the numerical prediction against the observed data, the root mean squared error was 2 mm, at position S6.
The position S3 at the lung surface in the numerical simulation is depicted in Figure 7. The top figure shows the lung at the beginning of inspiration and the bottom one shows the deformed lung at the end of inspiration obtained by simulation.
(a)
(b)
The graph of displacement versus time at location S3 is depicted in Figure 8. At S3, the experimental results (Figure 8) show that during the first 0.2 seconds of inspiration, the displacement of lung surface was almost zero. After that the displacement picked up and by the end of inspiration, or at 1.3 seconds, the S3 position was distended by 6 mm along SI direction. This trend has been well predicted by the simulation (Figure 8). The simulation shows that until 0.4 seconds of inspiration the deformation of lung surface at S3 position was almost zero. The numerical result showed that by the end of inspiration process, the position at the surface of the lung was distended by 6.7 mm in SI direction. Comparing the numerical results with the observed data, the displacement during inspiration process has been well predicted. Comparing the numerical prediction against the observed data, the root mean squared error was 0.72 mm, at position S3.
4. Concluding Remarks
Experimental methods, like 4D CT or gold marker tracking, to estimate or to dynamically track the movement of lung or tumor during breathing are not a very comfortable process for patients. Moreover those techniques are also considered quite risky as the patients have to be exposed in radiation for a long time. Hence, a numerical method that accurately predicts the lung or tumor deformation during breathing could serve as an alternative tool in radiation treatment. So far, most of the researches on modeling lung deformation during breathing employ a grid-based technique like finite element method. The grid-based technique sometimes suffers from grid distortion while modeling very high deformation. To respite from the grid distortion problem, a gridless technique could be an alternative. Hence, a gridless technique called MPS method is presented, in this paper, to dynamically track the deformation of lung surface during breathing. The potential of the proposed method to employ in the future pinpoint radiation therapy system has been explored. Variation of lung deformation at the lower and upper lobes was demonstrated, and the simulated results were compared against the experimental results. The comparison showed that the deformation was predicted with reasonable accuracy by the proposed method. Comparison of numerical prediction with experimental observations showed that the root mean squared error was about 2 mm at location S6 (lower lobe) and less than 1 mm at location S3 (upper lobe) at lung surface.
The application of this technique for the future pinpoint radiation system is still primal as the rigorous quantitative validation of numerical results is required to clear the limitation of the model before applying it in clinical application. In future, along with the quantitative validation test, the potential of the proposed model will be explored to dynamically track the displacement of cancer cell inside a lung during breathing.
Acknowledgments
The authors would like to heartily thanks the Associate Professor Dr. Masayori Ishikawa at the Department of Medical Physics and Engineering, Hokkaido University for providing the experimental data on lung deformation. The authors would also like to express sincere appreciation to the anonymous reviewers for their valuable comments and suggestion to improve our work.