KR102441143B1 - Simulation method of human respiratory system - Google Patents

Abstract

According to the present invention, a method for simulating the human respiratory system comprises: a first step in which a computer system extracts a plurality of slices of an image of the human respiratory system at maximum inspiration and maximum expiration; a second step of generating 3D models, respectively, at maximum inspiration and maximum expiration by using the plurality of slices extracted in the first step; a third step of performing a hole filling operation on the 3D model at maximum inspiration, generated in the second step, by using a Gaussian technique; a fourth step of performing a smoothing operation on the 3D model that has undergone the third step by using the Gaussian technique; a fifth step of determining the density and size of elements for the 3D model that has undergone the fourth step; a sixth step of generating a FE model by performing a meshing operation on the 3D model that has undergone the fifth step; a seventh step of applying an Ogden N6 model, which is a superelastic material model, to the lung parenchyma in the FE model generated in the sixth step, and applying elasticity properties to the bronchi; an eighth step of setting numerical analysis conditions to the FE model that has undergone the seventh step; and a ninth step of simulating the FE model to which the numerical analysis conditions are set, for a first reference time set in advance. Therefore, the method can realize fast and accurate minimally invasive surgery.

Description

Simulation method of human respiratory system

The present invention relates to a method for simulating the respiratory organs of a human body, and more particularly, to a method for simulating the respiratory organs of a human body through finite element analysis based on a hyperelastic configuration model.

Recently, various lung diseases are occurring due to the rapid increase of the elderly population and the deterioration of the environment, and the early diagnosis of lung diseases is being actively performed according to the development of imaging technology and lung cancer screening tests.

According to the study results, lesions such as indeterminate lung nodules and ground glass lung nodules were identified in 7.2% of the lung cancer screening population. Typically, a pulmonary nodule refers to a localized lung shadow of 3 cm or less in size. Recently, pulmonary nodules are increasingly found in the form of a semi-solid (part-solid).

In particular, when the semi-solid nodule is a tumor causing lung cancer, a biopsy through surgery is inevitably performed for pathological confirmation and treatment. However, if the lung nodule is not solid or the size of the lung nodule is small, so that the lung nodule cannot be exposed through minimally invasive surgery, the conversion to thoracotomy rather than minimally invasive surgery is estimated to reach 50%. is known

For successful minimally invasive surgery, accurate exposure of the tumor must be preceded, and the following methods have been applied, but it is reported that each method has fatal limitations.

The existing percutaneous medical device insertion method has a problem in that it does not have an anesthetic process, causing psychological anxiety and pain in the patient, and has a limitation in insertion depth, making it impossible to apply to central tumors of the lungs.

In addition, the radioactive chemical injection method has a problem in that excessive radiation exposure to patients and medical staff for material tracking occurs.

In addition, the water-soluble chemical injection method has a problem in that localization of the lung tumor location is impossible due to the high diffusion rate of the chemical itself, as well as increasing the number of operations.

This makes it difficult to target lung nodules in minimally invasive procedures.

To overcome this, a method for predicting three-dimensional lung deformation during regular and irregular breathing cycles is needed.

It is difficult to perfectly estimate the three-dimensional lung strain. Currently, three-dimensional lung deformation is approximated by a one-dimensional or multi-dimensional signal of a device such as a spirometer, an abdominal pressure belt, an external marker or an image.

Therefore, in order to minimize the incision of only the lung lesion area with minimally invasive surgery, it is necessary to evaluate/predict the lung behavior according to the patient's respiration and the accurate motion of the lesion site.

US 10-2010-0096470 A

The present invention has been devised to solve the above problems, and an object of the present invention is to provide a simulation method of a human respiratory system capable of simulating the respiratory system of the human body through finite element analysis based on a hyperelastic configuration model. have.

In addition, it is an object of the present invention to simulate a human respiratory system that can predict the physical behavior of the respiratory organs of the human body through finite element analysis based on a mathematical model according to the volume and pressure changes of the respiratory organs of the human body without opening the chest or inserting a material is to provide a way.

In order to solve the above technical problems, the method for simulating the human respiratory tract according to the present invention includes a first step of extracting a plurality of slices for an image of a human respiratory tract by a computer system during maximum inspiration and maximum expiration, the In the second step of generating 3D models of maximum inspiration and maximum expiration, respectively, with the plurality of slices extracted in the first step, hole filling is performed using the Gaussian technique on the 3D models of maximum inspiration generated in the second step a third step of performing a smoothing operation using a Gaussian technique on the 3D model that has undergone the third step, a fifth step of determining the density and size of elements with respect to the 3D model that has undergone the fourth step, In the sixth step of generating an FE model by performing meshing on the 3D model that has undergone the fifth step, the Ogden N6 model, which is a superelastic material model, is applied to the lung parenchyma in the FE model generated in the sixth step, and , a seventh step of applying elastic properties to the bronchus, an eighth step of setting numerical analysis conditions for the FE model that has passed through the seventh step, and simulating the FE model in which the numerical analysis conditions are set for a preset first reference time It is characterized in that it includes a ninth step.

In addition, in the fifth step, a tetrahedral element is applied as a mesh element when working with the mesh, and in the sixth step, the elastic modulus and Poisson's ratio of the bronchus are 5 kPa and 0.44, respectively.

In addition, in the eighth step, the bronchus of the FE model that has undergone the 8-1 step and the 8-1 step of restricting the rotational and translational motions in the three directions of the X, Y, and Z axes with respect to the upper end of the bronchi of the FE model, respectively. Step 8-2, which allows rotation and translation of the Z-axis with respect to the center of Step 3, Step 8-4, applying the transpulmonary pressure load condition to the FE model that has undergone Step 8-3, and Full Newtonian static analysis to the type of numerical analysis of the FE model that has gone through Step 8-4. Step 8-5 of applying, Step 8-6 of applying the nonlinear large deformation model to the lungs of the FE model that has undergone steps 8-5, Contact between the lungs and bronchi of the FE model that has undergone Steps 8-6 Steps 8-7 of applying a fixed face-to-face contact condition to the site and steps 8-8 of applying the master surface and slave surface to the lungs and bronchi of the FE model that have undergone the steps 8-7, respectively characterized.

In addition, the load condition of the transpulmonary pressure is a condition in which the transpulmonary pressure increases from 0 to a preset maximum pressure during the first reference time, and decreases from the maximum pressure to 0 during a second preset reference time. do it with

In addition, the first reference time is a time interval for exhalation in which the respiratory tract is converted from maximum inspiration to maximum expiration, and the second reference time is a time interval for inspiration in which the respiratory tract is converted from maximum inspiration to maximum inspiration, The second reference time is set to be greater than the first reference time.

In addition, the simulation method of the human respiratory system further comprises a tenth step of verifying the CT model at the time of maximum expiration and the FEM model at the time of the maximum expiration corresponding to the CT model, wherein the tenth step is the CT model at the time of the maximum expiration. Step 10-1 of overlapping the FEM model at the time of maximum expiration corresponding to the CT model on the model, and comparing the CT and FEM models overlapped in step 10-1 with a plurality of preset criteria from the upper part to the lower part of the respiratory tract. Step 10-2 of obtaining a plurality of reference cross-sections by cutting with a line segment, Step 10-3 of calculating the cross-sectional areas of the plurality of reference cross-sections obtained in Step 10-2 as landmark plane slices, respectively, the 10th It is characterized in that it includes a step 10-4 of calculating and comparing relative errors for each left lung, right lung, and bronchus from the cross-sectional areas of a plurality of reference sections calculated in step -3.

In addition, the CT model is a 3D model of the subject's respiratory organs generated in the second step, the FEM model is an FE model simulated under the load condition of the transpulmonary pressure, and the plurality of reference line segments are in the horizontal direction. It is characterized in that it is a line segment in which a plurality of line segments are offset to be spaced apart by a preset reference distance.

In addition, the relative error is a relative error between the cross-sectional area of the reference cross-section of the CT model (Area_CT) and the cross-sectional area of the reference cross-section of the FEM model (Area_FEM) corresponding to the cross-sectional area of the reference cross-section of the CT model, (Equation) (Area_FEM- It is characterized in that it is calculated by Area_CT)/(Area_CT)×100.

In addition, in the simulation method of the human respiratory system according to the present invention, the physical behavior of the lungs of the human body is predicted through finite element analysis based on a hydraulic model according to the volume and pressure changes of the lungs of the human body without opening the chest or inserting a material, and Clinically verified, it is possible to track the location of lung tumors after collapse, which has the effect of realizing rapid and accurate minimally invasive surgery.

In addition, the simulation method of the human respiratory system according to the present invention has the effect of accurately predicting the location of a tumor in the lung even in an extreme volume change such as collapse of the lung.

1 is a schematic diagram of the inspiratory and exhalation procedures in the human respiratory tract.
2 is a flowchart of a method for simulating a human respiratory system according to the present invention.
3 is a view showing 3D models of (a) bronchi, (b) lungs, and (c) respiratory tract, in which hole filling and smoothing operations are completed by a 3D slicer.
4 is a diagram illustrating an FE model of a respiratory organ generated by HyperMesh.
5 is a graph comparing the strain energies of the test data and the four models.
6 is a graph comparing the strain energies of test data and six models.
7 is a graph comparing the strain energies of test data and five models.
8 is a graph comparing the strain energy of the test data and the Ogden N6 model.
9 is a flowchart for step S800 of setting numerical analysis conditions in the FE model.
10 is a diagram in which loads, boundaries, and contact conditions are defined in the FE model.
11 is a graph showing the differential pressure between the alveoli and the pleura applied during respiration.
12 is a flowchart for step S900 of verifying the CT and FEM models during maximum expiration.
13A is a diagram showing the displacement results in the FE model for the respiratory tract in Case 1 during maximum expiration.
13B is a diagram showing the displacement results in the FE model for the respiratory tract in Case 2 during maximum expiration.
13C is a diagram showing the displacement results in the FE model for the respiratory tract in Case 3 during maximum expiration.
14A is a diagram showing the result values of the vector field in the FE model for the respiratory tract of Case 1 during maximum expiration.
14B is a diagram illustrating the result values of the vector field in the FE model for the respiratory tract of Case 2 during maximum expiration.
14C is a diagram showing the result values of the vector field in the FE model for the respiratory tract in Case 3 at the time of maximum exhalation.
FIG. 15 is a diagram showing overlapping CT and FEM models of Cases 1 to 3 during maximum expiration.
16 is a cross-sectional view of the CT and FEM models of Cases 1 to 3 taken along the line C1 to C7, respectively, during maximum exhalation.
17 is a graph showing cross-sectional areas and relative errors of C1 to C7 in the CT and FEM models of Case 1. FIG.
18 is a graph showing the cross-sectional areas and relative errors of C1 to C5 in the CT and FEM models of Case 2;
19 is a graph showing cross-sectional areas and relative errors of C1 to C5 in the CT and FEM models of Case 3;

Hereinafter, an embodiment of the present invention will be described with reference to the accompanying drawings in order to describe in detail enough that a person of ordinary skill in the art to which the present invention pertains can easily implement the technical idea of the present invention.

However, the following examples are merely examples to help the understanding of the present invention, thereby not reducing or limiting the scope of the present invention. In addition, the present invention may be embodied in several different forms and is not limited to the embodiments described herein.

First, the anatomy and physiology of pulmonary respiration will be described.

The mechanical movement of pulmonary respiration is controlled by two muscles: the diaphragm and the intercostal muscles attached to the ribs.

1 is a schematic diagram of the inspiratory and exhalation procedures in the human respiratory tract.

Referring to FIG. 1 , the diaphragm refers to a dome-shaped muscle fiber membrane that separates the chest from the abdominal cavity, and is composed of a central tendon and peripheral muscle fibers.

First, on inspiration, the diaphragm moves downward, creating negative pressure around the chest cavity and reducing the pressure in the lungs and abdominal organs. Then, the diaphragm moves upward during exhalation, increasing the pressure in the lungs and abdominal organs.

And, the lungs are connected to the ribs through the diaphragm and pleura. And, the pleura is composed of the visceral pleura and the body wall side pleura. Here, the visceral pleura refers to the pleura attached to the lung, and the body wall side pleura refers to the diaphragm and the pleura in contact with the chest wall.

At this time, the incompressible fluid flows between the visceral pleura and the body wall side pleura while the lung moves along with the movement of the diaphragm and the ribs.

Next, a method for simulating a human respiratory system according to the present invention will be described.

2 is a flowchart of a method for simulating a human respiratory system according to the present invention.

Referring to FIG. 2 , the computer system extracts a plurality of slices for images of human respiratory organs during maximum inspiration and maximum expiration. (S100)

After that, a 3D model of maximum inspiration and maximum expiration is generated using the plurality of slices extracted in step S100, respectively. (S200)

In this case, the 3D model generation program is preferably a 3D slicer. For example, a 3D model of maximum inspiration and maximum expiration can be extracted from 120 slices.

3 is a view showing 3D models of (a) bronchi, (b) lungs, and (c) respiratory tract, in which hole filling and smoothing operations are completed by a 3D slicer.

Referring to FIG. 3 , a hole filling operation is performed using a Gaussian technique so that a plurality of holes formed on the surface of the 3D model can be filled with respect to the 3D model at the time of maximum intake generated in step S200 ( S300 ). )

After that, a smoothing operation is performed using a Gaussian technique to smooth the surface of the 3D model for the 3D model that has undergone step S300 (S400).

In this case, the Gaussian hole filling and smoothing operation program is preferably a 3D slicer.

In the surface of the 3D model created in step S200, not only abnormal holes are sporadically formed, but also the surface is very rough. In order to compensate for the above problems, hole filling and smoothing are performed.

After that, the density and size of the elements are determined for the 3D model that has gone through the step S400 (S500).

In this case, it is preferable that the program for determining the density and size of the element is Design X or Hypermesh.

4 is a diagram illustrating an FE model of a respiratory organ generated by HyperMesh.

Referring to FIG. 4 , a mesh operation is performed on the 3D model that has undergone step S500 to generate a finite element (FE) model (S600).

In the mesh work of step S600, the mesh element is preferably a C3D4 element that is a tetrahedral element. Because a number of curves are formed on the surface of the respiratory tract, it is determined that a tetrahedral element is more suitable than a hexahedral element in order to express the respiratory tract similarly to the real thing.

And, the FE model generation program is preferably a mesh mixer (Meshmixer) or a hypermesh (Hypermesh).

After that, in the FE model generated in step S600, the Ogden N6 model, which is a superelastic material model, is applied to the lung parenchyma, and elastic properties are applied to the bronchi. (S700)

In the S700 series, the bronchial modulus and Poisson's ratio are preferably 5 kPa and 0.44, respectively.

After that, numerical analysis conditions are set for the FE model that has undergone step S700 (S800).

After that, the FE model in which the numerical analysis conditions are set is simulated for the first reference time set in advance. (S900)

Next, in the FE model generated in step S700, the reason for applying the Ogden N6 model, which is a superelastic material model, to the lung parenchyma will be described.

Most soft tissues under load exhibit a nonlinear stress-strain relationship with finite elasticity applied. Carter, Werner et al. modeled lung deformation by applying linear Hooke's law to explain a simple linear stress-strain relationship.

However, according to a study by Pinart et al., rat lungs exhibit nonlinear behavior under certain physiological conditions, such as pulmonary fibrosis, and Freed et al. is said to be suitable for

And, according to a study of 14 lung cancer patients, in the case of patients who exercise diaphragm during inspiration, most of the deformation is concentrated in the lower lobe of the lung and rapidly deforms within a short distance from the contact part of the lung and the diaphragm, transformation proceeds to the upper lobe of

Similar behavior was also observed in a comprehensive study of 152 lung cancer patients.

In the present invention, the simulation of the respiration of the lungs is performed from the time of maximum inspiration to the time of maximum expiration, that is, until the lungs expand to the maximum and contract to the maximum. When breathing in the lungs changes from inspiration to exhalation, and the lungs contract, enormous strain occurs in the lungs.

In this case, experimental data measured by Zeng, Gao, etc. are used for the physical properties of lung tissue. Lung tissue was subjected to tensile-compression experiments in biaxial mode, and experiments in uniaxial mode were also performed if necessary.

In general, the behavior of a superelastic material can be described as a potential of strain energy. In this case, strain energy means energy stored in a unit volume of the material.

5 is a graph comparing the strain energies of the test data and 4 models, FIG. 6 is a graph comparing the strain energies of the test data and 6 models, and FIG. 7 is a comparison of the strain energies of the test data and 5 models. It is one graph.

5 to 7 , in order to select a superelastic material model suitable for the lung parenchyma, experimental data and multiple models were evaluated with strain energy.

In this case, the strain energy evaluation program is preferably ABAQUS (ABAQUS 2020, Dassault System, USA).

At this time, the plurality of models are 1st polynomial (Polynomial N1), 2nd polynomial, Ogden N1, Ogden N2, Ogden N3, Ogden N4, Ogden N5, Ogden N6, Reduced Polynomial N1 , Neo Hooke), second-order reduced polynomial, third-order reduced polynomial (Yeoh), fourth-order reduced polynomial, fifth-order reduced polynomial, sixth-order reduced polynomial, and van der Waals.

When curve fitting multiple models to experimental data in the strain energy curve, first-order polynomial, Ogden N1, N2, and N6 models are the most stable.

On the other hand, in a strain energy curve with a nominal strain of 0.5 or more, nonlinearity is largely generated. When curve fitting multiple models to experimental data on curves with nominal stress of 0.5 or more, Ogden N3, N4, N5, and N6 models match the experimental data best.

And, when the least squares method is applied to multiple models, the sum of the squares of the error between the approximate solution and the actual solution in the Ogden N6 model becomes the minimum.

8 is a graph comparing the strain energy of the test data and the Ogden N6 model.

Referring to FIG. 8 , it can be seen that it is most preferable to apply the Ogden N6 model to the lung parenchyma in consideration of stability and least squares method.

Next, the numerical analysis model theory for explaining the respiratory system will be described.

The strain energy density of a superelastic material is the energy stored per unit volume and can be explained by the stress-strain behavior. The strain energy function based on the three invariants (I ₁ , I ₂ , I ₃ ) of the strain tensor is shown in Equation (1) below.

where W is the strain energy density, and I ₁ , I ₂ , I ₃ are the invariants of the strain tensor given below.

??

λ ₁ , λ ₂ , and λ ₃ in Equation (2) represent three principal stretch ratios, and in incompressible materials such as respiratory organs, λ ₂ λ ₃ =1, λ ₁ =λ, and λ ₂ =λ ₃ =λ ^-0.5 .

Therefore, if I ₃ =1 and k=0 are substituted into Equation (1), the last term of Equation (1) disappears, resulting in Equation (3).

In the present invention, 15 hyperelastic models were used to analyze the nonlinear stress-strain behavior of the lung according to the pressure load, and among them, the Ogden N6 model was used for numerical analysis.

In the Ogden model, it is assumed that λ ₁ λ ₂ λ ₃ =1, and the superelastic behavior of the material is expressed through Equation (4).

Here, μ _i and a _i are material constants that describe the shear force.

In the present invention, the Ogden N6 model to which the strain energy density potential (W) of N=6 is applied was used.

Next, the step S800 of setting numerical analysis conditions in the FE model will be described.

9 is a flowchart for step S800 of setting numerical analysis conditions in the FE model, and FIG. 10 is a view in which the load, boundary and contact conditions are defined in the FE model.

First, as shown in FIG. 10, rotation and translational motions in three directions of X, Y, and Z axes are constrained, respectively, with respect to the upper end (BC-1) of the bronchi of the FE model. (S801)

Here, the X, Y, and Z axes are formed in the left, forward, and opposite directions of gravity, respectively.

After that, rotation and translation of the Z-axis are allowed with respect to the center (BC-2) of the bronchi of the FE model, respectively. (S802)

After that, the rotation and translational motions of the X and Y axes are restricted to the center (BC-2) of the bronchi of the FE model, respectively. (S803)

The reason that steps S802 and S803 are applied is that the center of the bronchi is observed on CT to contract and relax like a bellows during respiration.

After that, the load condition of transpulmonary pressure is applied to the FE model. (S804)

11 is a graph showing the differential pressure between the alveoli and the pleura applied during respiration.

As shown in FIGS. 10 and 11 , the load condition of transpulmonary pressure applied to the FE model is that the transpulmonary pressure increases from 0 to the preset maximum pressure during the preset first reference time, and the preset second It means a condition that decreases from the maximum pressure to 0 during the reference time.

In this case, the transpulmonary pressure means a pressure value obtained by subtracting the internal pressure of the pleura from the internal pressure of the alveoli.

In addition, the first reference time refers to a time interval for exhalation in which the respiratory tract is converted from maximum inspiration (maximum expansion) to maximum expiration (maximum contraction). In addition, the second reference time refers to a time interval for inspiration during which the respiratory tract is switched from maximum expiration (maximum contraction) to maximum inspiration (maximum expansion).

And, the second reference time is set to be greater than the first reference time. At this time, it is preferable that the first reference time and the second reference time are 2 seconds and 3 seconds, respectively, and the maximum pressure is 10 cmH 2 O (=0.98 [kPa]).

Referring to FIG. 9, full Newtonian static analysis is applied to the type of numerical analysis of the FE model to obtain a solution of the nonlinear equilibrium equation (S805).

After that, a nonlinear large deformation model is applied to the lungs of the FE model to reflect the characteristics of the lungs with very large deformation (S806).

After that, fixed tie-to-surface contact conditions are applied to the contact areas of the lungs and bronchi of the FE model (S807).

After that, since the quantity of bronchial elements is smaller than the quantity of lung elements in the lungs and bronchi, a master surface and a slave surface are applied respectively (S808).

Next, a method for verifying the CT model, which is clinical data, and the FEM model, which is the result of simulation, will be described.

The simulation method of the human respiratory system according to the present invention is configured to further include a step S1000 of verifying the CT model at the time of the maximum expiration and the FEM model at the time of the maximum expiration corresponding to the CT model.

Here, the CT model is a 3D model of the subject's respiratory organs generated in step S200, and the FEM model means an FE model simulated under the load condition of transpulmonary pressure.

12 is a flowchart for step S1000 of verifying the CT and FEM models during maximum expiration.

Referring to FIG. 12 , first, the CT model at the time of maximal expiration is overlapped with the FEM model at the time of maximal expiration corresponding to the CT model (S1001).

Specifically, the CT model at the time of maximum exhalation is loaded on the 3D slicer. (S1001-1)

After that, the FEM model at the time of maximum expiration, which is the result of the simulation, is extracted as an OBJ file (S1001-2).

In this case, the FEM model at the time of maximum expiration may be extracted using the Abacus program.

After that, the FEM model at the time of maximal expiration extracted in step S1001-2 may be loaded on the 3D slicer, and the CT and FEM models at the time of maximal exhalation may be overlapped (S1001-3).

After that, in step S1001, the overlapping CT and FEM models are cut into a plurality of preset reference line segments from the upper part to the lower part of the respiratory tract to obtain a plurality of reference cross-sections. (S1002)

Here, the plurality of reference line segments are line segments in which a plurality of horizontal line segments are offset to be spaced apart by a preset reference distance.

For example, the plurality of reference line segments may be named C1, C2, C3, C4, C5, C5, C7, and the like, and the reference distance may be set to 28.6 mm.

After that, the cross-sectional areas of the plurality of reference cross-sections obtained in step S1002 are respectively calculated as landmark plane slices. (S1003)

Specifically, the markups module included in the landmark plane slice uses a plurality of resampled points to generate a plurality of landmark points at the edge of the reference section.

After that, the markups module sequentially calculates the perimeter, maximum curvature, minimum curvature, and cross-sectional area of the reference cross-section.

After that, the relative errors are calculated for each left lung, right lung, and bronchus from the cross-sectional areas of the plurality of reference cross-sections calculated in step S1003, respectively, and compared (S1004).

Here, the relative error is a relative error between the cross-sectional area of the reference cross-section of the CT model (Area_CT) and the cross-sectional area of the reference cross-section of the FEM model (Area_FEM) corresponding to the cross-sectional area of the reference cross-section of the CT model, (Equation) (Area_FEM-Area_CT) )/(Area_CT)×100.

Meanwhile, the simulation result is the analysis result obtained from the finite element discretization equation and can be expressed as a displacement vector (U).

13A is a diagram showing the displacement result value in the FE model for the respiratory tract of Case 1 during maximum expiration, and FIG. 13B shows the displacement result value in the FE model for the respiratory tract of Case 2 during maximum expiration. It is a diagram, and FIG. 13C is a diagram showing the displacement result value in the FE model for the respiratory tract of Case 3 during maximum exhalation.

14A is a diagram showing the result values of the vector field in the FE model for the respiratory tract of Case 1 at the time of maximum expiration, and FIG. It is a diagram, and FIG. 14C is a diagram showing the result value of the vector field in the FE model for the respiratory tract of Case 3 at the time of maximum exhalation.

13A to 13C , the displacement field of the displacement vector can be quantitatively confirmed in the FE model for the respiratory tract. And, the arrows in FIGS. 14A to 14C mean the final positions and directions of all nodes of the finite elements when the simulation is completed, and the length of the arrow is proportional to the coefficient of the displacement vector.

13A to 13C , the lung can be largely classified into three parts: an upper lobe, a middle lobe, and a lower lobe.

13A to 13C , during maximum exhalation (complete contraction), little deformation occurs in the upper lobe of the lung, whereas the greatest deformation occurs in the lower lobe of the lung with the diaphragm. At this time, the lung displacement value is calculated as an average of displacement values in the U1 (X-axis), U2 (Y-axis), and U3 (Z-axis) directions.

13A to 13C , the displacement values of the upper lobe are, respectively, Case 1: 0 to 7.3 mm, Case 2: 0 to 16 mm, and Case 3: 0 mm to 11 mm, respectively, and almost no deformation occurs in the upper lobe. On the other hand, the displacement values of the lower lobe are Case 1: 40.4-55.1 mm, Case 2: 56-96 mm, and Case 3: 48-63 mm, respectively, and the largest deformation occurs in the lower lobe. Through this, it can be seen that the largest deformation occurs in the dorsal part of the human body, which is the end of the U2 (Y-axis) direction.

In order to verify the proposed technique in the same way as above, the CT and FEM models of Cases 1 to 3 were compared and verified during maximum exhalation.

At this time, the FEM model refers to the FE model at the time of maximum expiration to which the load condition of transpulmonary pressure is applied. And, the transpulmonary pressure means a pressure value obtained by subtracting the internal pressure of the pleura from the internal pressure of the alveoli.

15 is a diagram showing the overlapping CT and FEM models of Cases 1 to 3, respectively, during maximum expiration, and FIG. 16 is a view showing the CT and FEM models of Cases 1 to 3 during maximum expiration, respectively, cut along the baseline segments C1 to C7. It is a cross section.

In this case, the reference line segments C1 to C7 are line segments offset by 28.6 mm from the upper part to the lower part of the respiratory tract.

In FIG. 16 , the left and right cross-sections are the right and left lungs of the human body, respectively, and the small circle provided between the left and right cross-sections means the bronchus.

In case 1 of FIG. 16 , the section C1 to C3 is the upper lobe of the lung and includes the lungs and bronchi, and the section C4 to C7 is the middle and lower lobes of the lung without bronchus. And, in Cases 2 and 3, the section C1 to C2 is the upper lobe of the lung with lungs and bronchi, and the section C3 to C5 is the section without bronchus as the middle and lower lobes of the lung.

In this case, since cases 2 and 3 have a shorter vertical length than case 1, only the sections C1 to C5 are displayed.

After sequentially calculating the perimeter, maximum curvature, minimum curvature, and cross-sectional area of the reference section from the reference sections of C1 to C7 in Case 1 and C1 to C5 in baseline segments C1 to C5 in Cases 2 and 3, the relative error between clinical data and simulation results ( relative 1error) was calculated.

Here, the relative error is a relative error between the cross-sectional area of the reference cross-section of the CT model (Area_CT) and the cross-sectional area of the reference cross-section of the FEM model (Area_FEM) corresponding to the cross-sectional area of the reference cross-section of the CT model, (Equation) (Area_FEM-Area_CT) )/(Area_CT)×100.

17 is a graph showing the cross-sectional areas and relative errors of C1 to C7 in the CT and FEM models of Case 1. FIG.

Referring to FIG. 17 , it can be seen that the total volumes of the CT and FEM models of Case 1 are 3,072.93 cm ³ and 3,318.05 cm ₃ , respectively, showing a slight difference of 245.12 cm ³ .

Referring to FIG. 17 , the relative errors of the upper lobe of the left lung, C1~C3, are 46.2%, 6.8%, and 7.1%, respectively, and the relative errors of the middle lobe, C4~C5, are 1.9%, 0.2%, respectively, and the lower lobe C6~ The relative errors of C7 are 29.2% and 98.9%, respectively.

And, the relative errors of the upper lobe, C1~C3, of the right lung, are 12.6%, 14.8, and 17.6%, respectively, the relative errors of C4~C5, the middle lobe, are 20.5% and 16.8%, respectively, and the relative errors of the lower lobe, C6~C7, are 108.9% and 66.8%, respectively.

And, the relative errors of bronchial C1 to C3 are 13%, 3%, and 53%, respectively. And, the total mean relative errors of the right lung and the left lung are 36.8% and 27.1%, respectively.

18 is a graph showing the cross-sectional areas and relative errors of C1 to C5 in the CT and FEM models of Case 2;

Referring to FIG. 18, the relative errors of the upper lobe C1 to C2 of the left lung are 46.2% and 18.3%, respectively, the relative error of C3, the middle lobe, is 7.4%, the relative errors of the lower lobe C4 to C5 are 6.8%, respectively, 5.9%.

And, the relative errors of the upper lobe C1~C2 of the right lung are 44.9% and 8.9%, respectively, the relative errors of the middle lobe C3 are 8.3%, and the relative errors of the lower lobe C4~C5 are 20.5% and 34.9%, respectively.

And, the total mean relative errors of the right lung and the left lung are 23.5% and 16.9%, respectively.

19 is a graph showing cross-sectional areas and relative errors of C1 to C5 in the CT and FEM models of Case 3;

Referring to FIG. 19 , the relative errors of the upper lobe C1 to C2 of the left lung are 68% and 17.6%, respectively, the relative errors of the middle lobe C3 are 1.2%%, and the relative errors of the lower lobe C4 to C5 are 7.2%, respectively. , 21.2%.

And, the relative errors of the upper lobe C1 to C2 of the right lung are 38.8% and 5.1%, respectively, the relative errors of the middle lobe C3 are 9.6%, and the relative errors of the lower lobe C4 to C5 are 16.1% and 80%, respectively.

And, the total mean relative errors of the right lung and the left lung are 29.9% and 23%, respectively.

17 to 19 , it can be seen that the right lung has a larger relative error than the left lung. This is because, when the lung contracts, the right lung is bounded and supported by the heart and liver, so that the right lung deforms more than the left lung.

And, it can be seen that the relative error is large in the order of the lower lobe, the upper lobe, and the middle lobe of the lung.

In addition, the shape of the bronchi is deformed according to the respiration of the lungs, and the phenomenon that the bronchi expand and contract like a bellows can also be confirmed in CT and FEM models.

In the present invention, a finite element analysis was performed by extracting a 3D model from the CT model and applying the physical properties of the human lung to the 3D model. And, in the present invention, based on the result of the finite element analysis, the relative error of the cross-sectional area between the CT and FEM models was calculated, and the CT and FEM models were compared and verified.

In addition, in the simulation method of the human respiratory system according to the present invention, the physical behavior of the lungs of the human body is predicted through finite element analysis based on a hydraulic model according to the volume and pressure changes of the lungs of the human body without opening the chest or inserting a material, and Clinically verified, it is possible to track the location of lung tumors after collapse, which has the effect of realizing rapid and accurate minimally invasive surgery.

In addition, the simulation method of the human respiratory system according to the present invention has the effect of accurately predicting the location of a tumor in the lung even in an extreme volume change such as collapse of the lung.

As described above, the present invention has a main technical idea to provide a simulation method of the human respiratory system, and the embodiment described above with reference to the drawings is only one embodiment, and the true scope of the present invention is the patent. However, based on the claims, it will also extend to equivalent embodiments that may exist in various ways.

10: Simulation device of the human respiratory system according to the present invention

Claims (8)

A method for simulating a human respiratory system by a computer system, the method comprising:
The simulation method of the human respiratory system is
a first step of extracting, by the computer system, a plurality of slices for images of human respiratory organs during maximum inspiration and maximum expiration;
a second step of generating 3D models of maximum inspiration and maximum expiration, respectively, from the plurality of slices extracted in the first step;
a third step of performing a hole filling operation using a Gaussian technique on the 3D model at the time of maximum intake generated in the second step;
a fourth step of performing a smoothing operation using a Gaussian technique on the 3D model that has undergone the third step;
a fifth step of determining the density and size of elements with respect to the 3D model that has undergone the fourth step;
a sixth step of generating an FE model by performing a meshing operation on the 3D model that has undergone the fifth step;
a seventh step of applying the Ogden N6 model, which is a superelastic material model, to the lung parenchyma in the FE model generated in the sixth step, and applying elastic properties to the bronchi;
an eighth step of setting numerical analysis conditions for the FE model that has undergone the seventh step; and
A ninth step of simulating the FE model in which the numerical analysis conditions are set for a preset first reference time; includes;
The eighth step
Step 8-1 of constraining the rotational and translational motions in three directions of the X, Y, and Z axes with respect to the upper end of the bronchus of the FE model, respectively;
an 8-2 step of allowing rotation and translation of the Z-axis with respect to the center of the bronchus of the FE model that has undergone the step 8-1;
a step 8-3 of constraining the rotation and translation movements of the X and Y axes with respect to the center of the bronchus of the FE model that has undergone the step 8-2;
a step 8-4 of applying a load condition of transpulmonary pressure to the FE model that has undergone step 8-3;
an 8-5 step of applying a Full Newtonian static analysis to the type of numerical analysis of the FE model that has undergone the 8-4 step;
a step 8-6 of applying a nonlinear large deformation model to the lungs of the FE model that has undergone steps 8-5;
an 8-7 step of applying a fixed face-to-face contact condition to the contact area of the lung and bronchi of the FE model that has undergone steps 8-6; and
and an 8-8th step of applying a master surface and a slave surface to the lungs and bronchi of the FE model that have undergone the 8th-7th steps, respectively. The method of claim 1,
In the sixth step,
When working with mesh, a tetrahedral element is applied as a mesh element,
In the seventh step,
A simulation method of the human respiratory system, characterized in that the elastic modulus and Poisson's ratio of the bronchi are 5 kPa and 0.44, respectively. delete The method of claim 1,
The load condition of the transpulmonary pressure is
The human respiratory system simulation method, characterized in that the transpulmonary pressure increases from 0 to a preset maximum pressure during the first reference time and decreases from the maximum pressure to 0 during a second preset reference time. 5. The method of claim 4,
The first reference time is
The time interval between exhalation during which the respiratory tract transitions from maximal inspiratory to maximal exhalation,
The second reference time is
It is the time interval for inspiration during which the respiratory tract transitions from maximum expiration to maximum inspiration,
The second reference time is
A simulation method of a human respiratory system, characterized in that it is set to be greater than the first reference time. The method of claim 1,
The simulation method of the human respiratory system is
A 10th step of verifying the CT model during the maximum expiration and the FEM model during the maximum expiration corresponding to the CT model; further comprising,
The tenth step is
a step 10-1 of overlapping the CT model at the time of maximal expiration with the FEM model at the time of maximal expiration corresponding to the CT model;
a step 10-2 of obtaining a plurality of reference sections by cutting the CT and FEM models overlapped in step 10-1 into a plurality of preset reference line segments from the upper part to the lower part of the respiratory tract;
a 10-3 step of calculating the cross-sectional areas of the plurality of reference cross-sections obtained in the 10-2 step as landmark plane slices, respectively;
a step 10-4 of calculating and comparing the relative errors for each left lung, right lung, and bronchi from the cross-sectional areas of the plurality of reference cross-sections calculated in step 10-3; and
The relative error is
As a relative error between the cross-sectional area of the reference cross-section of the CT model (Area_CT) and the cross-sectional area of the reference cross-section of the FEM model (Area_FEM) corresponding to the cross-sectional area of the reference cross-section of the CT model,
(Equation) (Area_FEM-Area_CT)/(Area_CT)×100 A simulation method of a human respiratory system, characterized in that it is calculated. 7. The method of claim 6,
The CT model is
It is a 3D model of the subject's respiratory organs generated in the second step,
The FEM model is
It is an FE model simulated under the load condition of transpulmonary pressure,
The plurality of baseline segments are
A simulation method of a human respiratory system, characterized in that the line segments in the horizontal direction are offset by a plurality of lines that are spaced apart by a preset reference distance.
delete

Images