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.

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.
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][4][5][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][8][9] . The bee needle is displayed in Fig. 1b. The front angle has 157°, the back angle, 110°, the height h is 0.5 mm, and the tip thickness b is 0.15 mm.
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][11][12][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][21][22] . The inverse sonification problem for capturing hardly detectable details in a medical image is treated in 23 , and the control in [24][25][26][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][29][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][34][35][36] (Fig. 2).
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][43][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][47][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. www.nature.com/scientificreports/ 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][50][51][52][53][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 (X, Y , Z) of base vectors (e 1 , e 2 , e 3 ) and origin O in the entry point of the skin (Fig. 3). The Euler frame K(x, y, z) with origin in the joint and the base vectors   Especially to our problem, u 1 and u 2 measure the bending of the needle, and the function u 3 measures the torsion of the needle. Therefore, u 1 and u 2 are components of the curvature of the central line denoted by κ corresponding to the planes (yz) and (xz) while u 3 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 l is invariant, the ends being fixed by external forces.
The link between the position vector r = (x, y, z) and unit tangential vector d 3 is obtained as r = s 0 d 3 ds , or To write the equations which describe the needle deformation, we introduce the inertia of the needle characterized by the functions where ρ 0 is the mass density per unit volume, A 0 the area of the cross section, I 1 , I 2 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 F = − with = ( 1 , 2 , 3 ) becomes (1) d 1 = (− sin ψ sin ϕ + cos ψ cos ϕ cos υ)e 1 + + (cos ψ sin ϕ + sin ψ cos ϕ cos υ)e 2 − sin υ cos ϕ e 3 , d 2 = (− sin ψ cos ϕ − cos ψ sin ϕ cos υ)e 1 + +(cos ψ cos ϕ − sin ψ sin ϕ cos υ)e 2 + sin υ sin ϕe 3 , d 3 = sin υ cos ψe 1 + sin υ sin ψe 2 + cos υe 3 . (2) − ∂ ∂t {k 1ψ sin 2 υ+k 2 (φ+ψ cos υ) cos υ}+ ∂ ∂s {Aψ ′ 2 sin 2 υ+C(ϕ ′ +ψ ′ cos υ) cos υ}+ 1 sin υ sin ψ− 2 sin υ cos ψ = 0, 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 where �(x, z, m) = For the present study, our objective is to determine the functions which measure the bending of the needle ( u 1 and u 2 ), and the function which measures the torsion of the needle ( u 3 ). 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. 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 h and the tip thickness of the needle b (see Fig. 1b).
For study the effect of the needle parameters h and b on the stress change during the insertion, simulations are carried out under different values for these parameters, as h = 0.5, 0.55, 0.6 and b = 0.15, 0.2, 0.25. The simulation results are shown in Fig. 7 for b ; a)

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][50][51][52][53][54][55][56] where V = U is the potential energy of deformation per unit volume (elastic potential) of the hepatic lobule, U is the total energy of the liver, is the volume of the hepatic lobule, and ε ij is the Lagrangian strain tensor, and with X i the material or Lagrangian coordinates and x i the Eulerian coordinates. After that, the potential theory 45 is defined by four terms: the free-electron energy U fe , then electrostatic Coulomb energy U es which is often called the Mandelung energy, the band-structure energy U bs and the Born-Mayer ion-core repulsive energy U r [57][58][59][60][61][62] Jankowski and Tsakalakos 59 advanced that U r is the predominant term for the evaluation of the elastic constants. The sum is extended to all nearest neighbors points of coordinates (X, Y , Z) located at distances R (n) , R = X 2 + Y 2 + Z 2 with respect to the green atom from Fig. 7, R (n) , α 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 −21 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 E = 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 P i (b 6i−5 , . . . , b 6i ) , i = 1, 2, . . . , 7 , characterized by coefficients b j , j = 1, 2, . . . , 42 . To calculate P i , i = 1, 2, ...7 from experimental data, an objective function ℑ is chosen to measure the agreement between theoretical and experimental data (17)  (1) (42) 0 . 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 α × 10 6 = 0.19 Ryd and β = 10.22 ua.
Describe the constitutive equations for the isotropic centrosymmetric Cosserat solids as follows [49][50][51][52][53] where σ kl is the Cauchy stress tensor, m kl is the couple stress tensor, e kl = 1 2 (u k.l + u l,k ) is the macrostrain strain tensor, u is the displacement vector and ε klm is the Levi-Civita permutation tensor. Further, the microrotation vector ϕ k is distinct from the macrorotation vector r k = 1 2 ε klm u m,l , i.e. ϕ k refers to the rotation of points themselves, while r k 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.
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 m kl = αϕ r,r δ kl + βϕ k,l + γ ϕ l,k , Table 1. Estimates for elastic constants in the ( X, Y , Z ) coordinates as a function of strain. www.nature.com/scientificreports/ 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][64][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 O 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. A Cartesian coordinates system (x, y, z) is attached to the human liver. The needle moves in the direction of the z-axis, in the origin of the cartesian coordinate. The concentrated needle force has a magnitude F = F 0 δ(x)δ(y)δ(t) . For simplicity without loss of generality we consider that τ 1 = τ 2 = · · · = τ 8 = F 0 . Also, we assume the particular case in which all quantities depend only on x and z.
To this end, the constitutive relations (21) and (22) are completed with the kinematic description which includes the microrotations ϕ k , k = 1, 2, 3 , as independent degrees of freedom in addition to the usual displacements u i , i = 1, 2, 3 . The transfer of loading between points in the liver is achieved through the couple stresses m ij , i, j = 1, 2, 3 and the classical Cauchy stresses σ ij , i, j = 1, 2, 3.
The equilibrium equations are in the absence of body forces and body couples Define the kinematic variables as the strain tensor and 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 u 1 , u 2 and u 3 of the Eqs. (23)- (26) are taken under the form where the first term z lin represents a linear superposition of cnoidal functions and the second term z nonlin , a nonlinear superposition of the cnoidal functions 55 The analytic expressions of U 1 = u 1 F 0 , U 2 = u 2 F 0 and U 3 = u 3 F 0 are available once the displacements u 1 , u 2 and u 3 are specified. Variation of displacements U 1 , U 2 and U 3 with respect to x for the human liver is presented in Fig. 8. Variations of the normal stresses T 13 = σ 13 F 0 ,T 31 = σ 31 F 0 ,T 32 = σ 32 F 0 and T 33 = σ 33 F 0 with respect to x are presented in Fig. 9 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 (23) σ ji,j = 0, (24) m ji,j + ε jk σ jk = 0, i, j, k = 1, 2, 3.

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.