Research Article | Open Access

Subas Chhatkuli, Seiichi Koshizuka, Mitsuru Uesaka, "Dynamic Tracking of Lung Deformation during Breathing by Using Particle Method", *Modelling and Simulation in Engineering*, vol. 2009, Article ID 190307, 7 pages, 2009. https://doi.org/10.1155/2009/190307

# Dynamic Tracking of Lung Deformation during Breathing by Using Particle Method

**Academic Editor:**Ewa Pietka

#### 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/m^{3}, 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.

#### References

- M. Uesaka, K. Mizuno, A. Sakumi et al., โPinpoint keV/MeV X-ray sources for X-ray drug delivery system,โ in
*Proceedings of the IEEE Particle Accelerator Conference (PAC '07)*, pp. 2793โ2795, Albuquerque, NM, USA, June 2007. View at: Publisher Site | Google Scholar - S. Mori, M. Endo, S. Komatsu, T. Yashiro, S. Kandatsu, and M. Baba, โFour-dimensional measurement of lung tumor displacement using 256-multi-slice CT-scanner,โ
*Lung Cancer*, vol. 56, no. 1, pp. 59โ67, 2007. View at: Publisher Site | Google Scholar - M. Nakao, A. Kawashima, M. Kokubo, and K. Minato, โSimulating lung tumor motion for dynamic tumor-tracking irradiation,โ in
*Proceedings of the IEEE Nuclear Science Symposium Conference Record*, vol. 6, pp. 4549โ4551, Honolulu, Hawaii, USA, October-November 2007. View at: Publisher Site | Google Scholar - A.-L. Didier, P.-F. Villard, J.-Y. Bayle, M. Beuve, and B. Shariat, โBreathing thorax simulation based on pleura physiology and rib kinematics,โ in
*Proceedings of the 4th International Conference Medical Information Visualisation: BioMedical Visualisation (MediViz '07)*, pp. 35โ40, Honolulu, Hawaii, USA, October 2007. View at: Publisher Site | Google Scholar - P.-F. Villard, M. Beuve, B. Shariat, V. Baudet, and F. Jaillet, โSimulation of lung behaviour with finite elements: influence of bio-mechanical parameters,โ in
*Proceedings of the 3rd International Conference on Medical Information Visualisation-BioMedical Visualisation (MediVis '05)*, pp. 9โ14, July 2005. View at: Publisher Site | Google Scholar - A. Al-Mayah, J. Moseley, and K.K. Brock, โContact surface and material nonlinearity modeling of human lungs,โ
*Physics in Medicine and Biology*, vol. 53, no. 1, pp. 305โ317, 2008. View at: Publisher Site | Google Scholar - S. Koshizuka and Y. Oka, โMoving-particle semi-implicit method for fragmentation of incompressible fluid,โ
*Nuclear Science and Engineering*, vol. 123, no. 3, pp. 421โ434, 1996. View at: Google Scholar - Y. Chikazawa, S. Koshizuka, and Y. Oka, โA particle method for elastic and visco-plastic structures and fluid-structure interactions,โ
*Computational Mechanics*, vol. 27, no. 2, pp. 97โ106, 2001. View at: Publisher Site | Google Scholar - M. Song, S. Koshizuka, and Y. Oka, โA particle method for dynamic simulation of elastic solids,โ in
*Proceedings of the 6th WCCM in Conjunction with APCOM*, Beijing, China, 2004. View at: Google Scholar - M. S. Song, S. Koshizuka, and Y. Oka, โMoving particle semi-implicit (MPS) method for solid simulation,โ in
*Proceedings of the American Nuclear Society of International Congress on Advances in Nuclear Power Plants (ICAPP '05)*, pp. 1125โ1132, Seoul, Korea, May 2005. View at: Google Scholar - S. Koshizuka, M. Song, and Y. Oka, โA particle method for three- dimensional elastic analysis,โ in
*Proceedings of the 6th WCCM in Conjunction with APCOM*, Beijing, China, 2004. View at: Google Scholar - M. Kondo, S. Koshizuka, and Y. Suzuki, โApplication of symplectic scheme to three-dimensional elastic analysis using MPS method,โ
*Transactions of the Japan Society of Mechanical Engineers, Part A*, vol. 72, no. 4, pp. 425โ431, 2006. View at: Google Scholar - Y. Suzuki and S. Koshizuka, โA Hamiltonian particle method for non-linear elastodynamics,โ
*International Journal for Numerical Methods in Engineering*, vol. 74, no. 8, pp. 1344โ1373, 2008. View at: Publisher Site | Google Scholar - J. J. Monaghan, H. E. Huppert, and M. G. Worster, โSolidification using smoothed particle hydrodynamics,โ
*Journal of Computational Physics*, vol. 206, no. 2, pp. 684โ705, 2005. View at: Publisher Site | Google Scholar - R. Speith and H. Riffert, โThe viscous gas ring as an astrophysical test problem for a viscous SPH-code,โ
*Journal of Computational and Applied Mathematics*, vol. 109, no. 1-2, pp. 231โ242, 1999. View at: Google Scholar

#### Copyright

Copyright © 2009 Subas Chhatkuli et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.