On the flexible needle insertion into the human liver.

In the present research, the navigation of a flexible needle into the human liver in the context of the robotic-assisted intraoperative treatment of the liver tumors, is reported. Cosserat (micropolar) elasticity is applied to describe the interaction between the needle and the human liver. The theory incorporates the local rotation of points and the couple stress (a torque per unit area) as well as the force stress (force per unit area) representing the chiral features of the human liver. To predict the deformation of the needle and the liver, the elastic properties of the human liver have been evaluated. Outcomes reveal that considering smaller deformations of the needle and the liver results in better needle navigation mechanism. The needle geometry can enhance the penetration.

Introduction

By suggesting the flexible bee needles as useful tools to transport drugs into the liver tumors, the needle navigation performance enhances[1,2]. We refer to the insertion trajectory of the needle which should avoid the ribs, blood vessels, and other organs to protect the liver[3-6] (Fig. 1a). The bee needle has the advantage to reduce the insertion forces and to ensure small tissue deformations. The literature reports a number of interesting papers on the surgical needle navigation into the liver[7-9]. The bee needle is displayed in Fig. 1b. The front angle has 157°, the back angle, 110°, the height is 0.5 mm, and the tip thickness is 0.15 mm.

Figure 1

(a) The suggested trajectory towards the liver tumor; (b) Honeybee barbed needle inspired from[1,2].

Another topic refers to the collision free trajectory of the needle to the target. This topic requires experience in imaging the tumor location, in the liver structure and in the microstructural interaction between the needle and the liver. Elastic properties of the liver, minimum execution time, minimum energy of the needle navigation, load carrying capacity are some topics of interest.

The needle flexibility is essential for a good precision in the handling. Important concentration in strain and stress and the topological changes of the liver are not to be neglected during the needle navigation towards the tumor[10-13]. Details of the forces during needle insertion into the liver are find in[14], the real time collision detection for virtual surgery in[15] and the minimal hierarchical collision detection in[16]. Optimization is required to modify the needle trajectory in order to protect the liver[17,18], to manage the tumor risk[19], and to change the robot architecture[20-22]. The inverse sonification problem for capturing hardly detectable details in a medical image is treated in[23], and the control in[24-27].

Microscopic investigation of the human liver offers details of its microanatomy with emphases to the granular, fibrillar components and irregular solid–fluid interfaces[28-30]. The basic functional unit of the liver is the hepatic lobule which comprises a hexagonal and a portal triad-portal vein, hepatic artery, bile duct[31,32]. Lobuli form a two layers membrane with internal space of 100 A and the cellular elements with twisted, spiraling fibers braided into the helical and screw-shaped gaps (pores) of 40–100 µm in size[33-36] (Fig. 2).

Figure 2

Representation of the hepatic lobule—basic functional unit of the liver.

A straight line for the needle trajectory is typically used in the treatment of the liver tumors. But avoiding the obstacles and reaching the regions currently inaccessible using straight line trajectory require the motion planning algorithm for flexible needle insertion with integrative models of human liver[37] with knowledge of the interaction forces raised during needle insertion[38]. The feedback force from the needle can be formulated as the gradient of the potential energy of the soft tissue based on particle constraint[39]. The trajectories inside the liver must be recorded by a camera to compare with the simulation trajectories in order to reduce the errors between the experimental and simulation trajectories less than 0.8 mm[40].

Development of a robotic system and control algorithms is discussed in[41] with different topics related to needle steering. Mechanics-based models are adapted from beam theories[42-44]. In these models, the fact that the needle deflection and tissue deformation are coupled is considered. 3D-needle shape reconstruction with an array of fiber brag grating sensors is developed in[45] and inventions were reported for guiding devices for needle placement and performing the percutaneous computed tomography-guided surgical activity[46-48].

Two challenges persist until now: what kind of material the liver is made of, and which is the interaction between the needle and the liver. The present paper puts together results on the elastic properties of the human liver, deformation of the needle and the liver, respectively.

Investigations on the human liver confirm a chiral (noncentrosymmetric) behavior providing evidence of isotropy with respect to the coordinate rotations but not with respect to inversions. The Cosserat (micropolar) elasticity is the appropriate theory that recognizes and describe the rotation of the cellular components as well as the translation, the couple per unit area as well as force per unit area of the hepatic membranes, the size effect in tension and bending, and the stress concentration near discontinuities[49-54]. It should be added that these features cannot be described by conventional elastic theories. The chirality leads to a vibrational amplification of the displacements and stresses which can be explained by a suitable adaptation of the liver dynamics to the attractive and repulsive forces. The chirality-triggered oscillations suggest that the linearity is preserved at the microscopic level, while becoming strongly nonlinear to the macroscopic scale.

Deformation of the needle

In this research, the common focus of a serial surgical robot composed of a revolute joint and a flexible needle, was set to a reference Lagrange frame of base vectors and origin in the entry point of the skin (Fig. 3). The Euler frame with origin in the joint and the base vectors is attached to the needle. The angle between the flexible arm and axis is . Bending and torsion of the needle are described by the strain functions . The robot has degrees of freedom , where 1 is the generalized coordinate of the rigid system and  = 3 are the degrees of freedom of the flexible needle.

Figure 3

Schematic of the needle with the Lagrange coordinate system and the Euler coordinate system attached.

According to the previous studies of Munteanu and Donescu[55,56], the orientation of the Euler axes relative to the Lagrange axes are expressed by Euler angles and

Strain functions measure the bending and torsion of the needle aswhere means the partial differentiation with respect to s—the coordinate along the central line of the needle.

Especially to our problem, and measure the bending of the needle, and the function measures the torsion of the needle. Therefore, and are components of the curvature of the central line denoted by corresponding to the planes and while is the torsion of the needle denoted by

In this way, the needle is rigid along the tangential direction and the total length of the needle is invariant, the ends being fixed by external forces.

The link between the position vector and unit tangential vector is obtained as , or

To write the equations which describe the needle deformation, we introduce the inertia of the needle characterized by the functionswhere is the mass density per unit volume, A0 the area of the cross section, are geometrical moments of inertia around the axis, which is perpendicular to the central axis and respectively around the central axis.

The exact set of equations of the needle with the ends fixed by the force with becomeswhere A and C are the bending stiffness and respectively the torsional stiffness of the needle, related to the Lame constants , by where is the Young's elastic modulus, and a is the radius of the cross section of the needle, and

The system of Eqs. (7)–(11) is exactly solved by the cnoidal method. As a result, the closed form solutions of the Euler angles and are derived[55]with and where is the normal elliptic integral of the third kind. Functions are solutions of the equationwith

For the present study, our objective is to determine the functions which measure the bending of the needle ( and ), and the function which measures the torsion of the needle (). This can be done by (2). The strain profile of the needle is computed for two needle routes as considered in Fig. 4. For the first route, the tumor is red and the entry point is A and for the second route the tumor is blue with entry at point B.

Figure 4

Two needle routes: the first route with a red tumor and the entry point A, and the second route longer with a blue tumor and the entry point B.

Figure 5 represents the bending functions and of the needle which depend on the coordinate along the central line of the needle , and the time , for the route (a) and for the route (b), respectively. In a similar matter, Fig. 6 represents the torsion functions of the needle.

Figure 5

Functions and for (a) the first route and (b) the second route.

Figure 6

Function for (a) the first route and (b) the second route.

It is noteworthy that Figs. 5 and 6 show that the deformation of the needle depends on the needle trajectory. For both cases, the deformation of the needle is small and finite without any tendency to grow towards chaos. The strains operate in a solitonic regime in which the localized waves propagate for a long time without changes. The soliton is a localized wave with an infinite number of degrees of freedom that conserve their properties even after interaction among them, and then act somewhat like particles[55]. The basic idea is that the Eqs. (7)–(11) have unique properties that are locally preserved such as an infinite number of exact solutions expressed in terms of the Jacobi elliptic functions or the hyperbolic functions, and the simple formulae for nonlinear superposition of explicit solutions.

Degree of deformation of the honeybee barbed needle during insertion into the liver is caused by the flexibility of the needle. The main parameters that determine the flexibility of the needle are the height and the tip thickness of the needle (see Fig. 1b).

For study the effect of the needle parameters and on the stress change during the insertion, simulations are carried out under different values for these parameters, as 0.5, 0.55, 0.6 and 0.15, 0.2, 0.25. The simulation results are shown in Fig. 7 for ; a) 0.5 and 0.15; b) 0.55 and 0.2; c) 0.6 and 0.25.

Figure 7

Insertion simulation results for different values of and ; (a) 0.5 and 0.15; (b) 0.55 and 0.2; (c) 0.6 and 0.25.

We can see from the Fig. 7 that if the parameters and of the needle increase, the deformation of the needle in the tissue also increases. The stress order of magnitude for 0.5 and 0.15 is 0.01, for 0.55 and 0.2 is about 0.1, while for 0.6 and 0.25 is about 0.17. We conclude that the larger parameters and of the needle, the greater is the stress distribution.

Elastic constants of the liver

Cosserat theory makes our work distinct from previous studies describing the elasticity of the human liver. In spite of the fact that the liver was usually modeled as an elastic medium, we noted a rich dynamic behavior of the human liver involving nonlinearity and chirality. Our findings show that the effect of the chirality leads to a vibrational amplification of the elastic response of the liver due to the needle penetration, explained by the adaptation of the liver dynamics to the attractive and repulsive forces (inter-atomic bonds). The human liver is viewed as a spatial overlap of mesoscopic subsystems in which interferences overlap or are lost, as such it naturally vibrates through overlapping the oscillations—and this result is a natural effect of chirality. Therefore, the human liver can be regarded as an elastic chiral Cosserat material (noncentrosymmetric) characterized by the fourth rank elastic constant[49-56]where is the potential energy of deformation per unit volume (elastic potential) of the hepatic lobule, is the total energy of the liver, is the volume of the hepatic lobule, and is the Lagrangian strain tensor, andwith the material or Lagrangian coordinates and the Eulerian coordinates.

After that, the potential theory[45] is defined by four terms: the free-electron energy , then electrostatic Coulomb energy which is often called the Mandelung energy, the band-structure energy and the Born–Mayer ion-core repulsive energy [57-62]

Jankowski and Tsakalakos[59] advanced that is the predominant term for the evaluation of the elastic constants. The sum is extended to all nearest neighbors points of coordinates located at distances , with respect to the green atom from Fig. 7, , is the parameter of the repulsive energy and is the repulsive range parameter[57]. The is measured in Ryd (Rydberg) 1 Ryd = 13.6 eV = 2.092 10 J, and in atomic units [ua]. These parameters are evaluated from a genetic algorithm that use the experimental value of the elastic modulus of the liver 5.9359 Pa and Poisson’s ratio 0.49 available in[58].

The aim of the inverse problem is to use the difference between the experimental and theoretical values of the Young’s modulus and the Poisson’ratio to provide a procedure that leads to the least disaccord between predictions and experimental observations. We consider that and are approximated by polynomials of five degree , , characterized by coefficients , . To calculate , from experimental data, an objective function is chosen to measure the agreement between theoretical and experimental datawhere , are the predicted values of and , and , , are the experimental values of and . The functions , , are estimated by a genetic algorithm. This algorithm assures an iteration scheme that guarantees a closer correspondence of predicted and experimental values of and at each iteration. A binary vector with 42 genes representing the real values of the parameters , , is used. The length of the vector is six places after the decimal point. The domain of parameters with length is divided into a least 15,000 equal size ranges. That means that each parameter , , is represented by a gene (string) of 22 bits . One individual consists of the row of 42 genes, that is, a binary vector with 22 42 components given by . From one generation to the next the genetic algorithm usually decreases the objective function of the best model. The starting population is usually randomly generated. Then, new descendant populations are iteratively created with the goal of an overall objective function decrease from generation to generation. Each new generation is created from the current one by the main operations: selection, crossover and reproduction, mutation and fluctuation. The alternation of generations is stopped when convergence is detected. If no convergence the iteration process continues until the specified maximum number of generations is reached.

We report here the results of the genetic algorithm. The unknown and are obtained after 31 iterations. The results of the genetic algorithm are Ryd and 10.22 ua.

Describe the constitutive equations for the isotropic centrosymmetric Cosserat solids as follows[49-53]where is the Cauchy stress tensor, is the couple stress tensor, is the macrostrain strain tensor, is the displacement vector and is the Levi–Civita permutation tensor. Further, the microrotation vector is distinct from the macrorotation vector , i.e. refers to the rotation of points themselves, while refers to the rotation associated with movement of nearby points. The comma denotes differentiation with respect to spatial coordinates and a superposed dot indicates the time rate.

Involved elastic constants are the Lame elastic constants , and , the Cosserat rotation modulus , and the first, second and third microrotation constants .

For a healthy liver without tumors, the estimates of the elastic constants and are shown in Table 1.

Table 1

Estimates for elastic constants in the () coordinates as a function of strain.

Strain	λ[Pa]	μ[Pa]	κ[Pa] ×10^2	α [N] ×10^4	β[N] ×10^4	γ[N] ×10^4
− 0.040	97.7806	2.1033	459	2.0976	12.0833	2.0865
− 0.036	97.7703	2.0876	432	2.0878	12.0775	2.0814
− 0.032	97.7657	2.0456	401	2.0502	12.0656	2.0800
− 0.028	97.7587	1.9981	396	2.0481	11.9961	1.9976
− 0.024	97.7498	1.9978	388	1.9934	11.9908	1.9974
− 0.020	97.7480	1.9964	387	1.9914	11.9904	1.9969
− 0.016	97.7398	1.9956	381	1.9897	11.9895	1.9958
− 0.012	97,7341	1.9933	374	1.9883	11.9891	1.9913
− 0.008	97.7034	1.9928	371	1.9878	11.9792	1.9900
− 0.004	97.6351	1.9925	366	1.9865	11.9725	1.9895
0.000	97.6037	1.9919	354	1.9849	11.9719	1.9819
0.004	97.5939	1.9909	351	1.9839	11.9709	1.9809
0.008	97.5744	1.9874	342	1.9814	11.9704	1.9776
0.012	97.3901	1.9832	340	1.9802	11.9695	1.9772
0.016	96.8372	1.9802	339	1.9782	11.9692	1.9707
0.020	96.6333	1.9745	336	1.9735	11.9655	1.9705
0.024	96.4403	1.9698	332	1.9718	11.9598	1.9688
0.028	96.3012	1.9453	329	1.9653	11.9553	1.9552
0.032	95.7854	1.9245	324	1.9625	11.9545	1.9265
0.036	95.7840	1.9203	322	1.9603	11.9409	1.9213
0.040	97.7737	1.9189	321	1.9518	1.9408	1.9177

It is interesting to see from the Table 1 that the elastic constants of liver reduces all the way when the strain increases from negative (compression) to zero and then to positive(extension) values. This is explained by the dependency of the liver elasticity on the size of the inhomogeneities and the surface/interface stress at the micro-scale.

Finally, the results show that even small strains can affect the values of the elastic constants. The increase in the strain reduces the value of the elastic constants.

Deformation of the liver during the needle insertion

The deformation of the needle is coupled with the deformation of the liver. To evaluate the liver deformation, the connection between the needle and the liver is modeled as a spring layer with a very small thickness[63-65]. The tractions are continuous but displacements can be discontinuous across the layer.

The evaluation of the liver deformation is simplified by discretization of the interface between the needle and the liver, into tiny homogeneous cells[65]. Across this grid of points, a network of springs is introduced to ensure that the behaviour inside each component is elastic and, in the case of perfect contact interfaces, a perfect contact among different components.

Figure 9 represents a 2D spring model with a generic point the tip of the needle. The nearest neighbours in the liver are labelled from 1 to 8. It is thus possible to obtain the iteration equations for the deformation of the liver starting from the deformation of the needle.

Figure 9

2D spring model for the interface between the tip of the needle and the liver[65].

A Cartesian coordinates system is attached to the human liver. The needle moves in the direction of the -axis, in the origin of the cartesian coordinate. The concentrated needle force has a magnitude . For simplicity without loss of generality we consider that . Also, we assume the particular case in which all quantities depend only on and .

To this end, the constitutive relations (21) and (22) are completed with the kinematic description which includes the microrotations , , as independent degrees of freedom in addition to the usual displacements , . The transfer of loading between points in the liver is achieved through the couple stresses , and the classical Cauchy stresses , .

The equilibrium equations are in the absence of body forces and body couples

Define the kinematic variables as the strain tensorand the strain gradient tensor

As we said before, the basic functional unit of the liver is the hexagonal hepatic lobule which comprises the portal triad-portal vein, hepatic artery, bile duct.

The solutions in displacements and of the Eqs. (23)–(26) are taken under the form

where the first term represents a linear superposition of cnoidal functions and the second term , a nonlinear superposition of the cnoidal functions[55]

The analytic expressions of , and are available once the displacements and are specified. Variation of displacements and with respect to for the human liver is presented in Fig. 8. Variations of the normal stresses ,, and with respect to are presented in Fig. 9

Figure 8

Representation of the twelve nearest neighbors of the green atom.

Graphs from Figs. 10 and 11 illustrate an oscillatory behavior. After 1.5 s the graphs do not change significantly, this means that the soliton functions with which these solutions are expressed are stabilized as form and identity. The oscillatory phenomenon is a result of the emission of acoustic signals rich in ultrasound components through which the liver atoms communicate between them to balance the liver state. The role of vibrations in the strain and the stress fields depends on the ability of the liver to respond to loading conditions by making use of the vibratory communication between tissue and the needle. This communication is made by means of excitation signals rather that the response signals alone. To our knowledge, this vibratory communication of the human liver with the needle has not previously been recorded.

Figure 10

Plots of and with respect to .

Figure 11

Plots of and with respect to .

Snapshots at various time of the amplitude of the stress in the tip of the needle are shown in Fig. 12.

Figure 12

Snapshots at various time of the amplitude of the stress in the tip of the needle.

The deformation of the needle at 400, 600, 800 and 1000 s. are represented in Fig. 13.

Figure 13

The deformation of the needle at 400, 600, 800 and 1000 s.

Conclusions

As mentioned in the introduction, the purpose of the present study is to investigate the navigation of a flexible needle into the human liver. In our analysis, we adopted the Cosserat elasticity to describe the interaction between the needle and the human liver. The theory incorporates the local rotation of points and the couple stress as well as the force stress representing the chiral properties of the human liver. However, since the liver is a deformable body that needs to be mechanically characterized, its elastic properties and deformation are evaluated. Estimates of the elastic constants have been made for a healthy liver in the absence of tumors. Computations result from the pseudopotential energy by retaining only the predominant terms in the evaluation of the elastic constants.