Modeling of the Flexible Needle Insertion into the Human Liver

Abstract

The insertion of the needle is difficult because the deformation and displacement of the organs are the key elements in the surgical act. Liver and tumor modeling are essential in the development of the needle insertion model. The role of the needle is to deliver into the tumor an active chemotherapeutic agent. We describe in this chapter the deformation of the needle during its insertion into the human liver in the context of surgery simulation of the high- robotic-assisted intraoperative treatment of liver tumors based on the integrated imaging-molecular diagnosis. The needle is a bee barbed type modeled as a flexible thread within the framework of the Cosserat (micropolar) elasticity theory.

1. Introduction

The flexible bee needles are useful tools to transport drugs into the liver tumors [1, 2]. The insertion trajectory of the needle must avoid the ribs, blood vessels, and other organs to protect the liver [3, 4, 5, 6] (Figure 1a). The bee needle assures reduced insertion forces and small tissue deformations because of the tip deflections. The furthermost current publications on the surgical needle navigation into the liver can be demonstrated in [7, 8, 9]. The bee needle is shown in Figure 1b. The front angle has 157 deg., the back angle, 110 deg., the height is 0.5 mm, and the tip thickness 0.15 mm.

Figure 1.

a) Trajectory towards the liver tumor; b) honeybee barbed needle [1, 2].

A number of scientific researches have been carried out on the collision free trajectory of the needle to the target. The surgical event requires experience in imaging the tumor location based on the liver structure and the microstructural interaction between the needle and the liver. Several studies have revealed that the needle flexibility is essential to achieve a good precision in the handling.

The strain and stress fields 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 on its microanatomy with emphases to the granular, fibrillar components and irregular solid–fluid interfaces [28, 29, 30]. The basic unit of the liver is the hepatic lobule which is a hexagonal element with comprised the portal triad -portal vein, hepatic artery and the bile duct [31, 32]. Lobuli form two layers membranes with internal space of 100A and the cellular elements with twisted, spiraling fibers braided into the helical and screw-shaped gaps (pores) of 40–100 $\mathrm{\mu m}$ in size (Figure 2) [33, 34, 35, 36].

Figure 2.

Hepatic lobule - basic unit of the liver.

In this chapter we try to answer a few questions such as how is the deformation of the needle and how the free-collision trajectories are determined.

2. Deformation of the needle

The needle is a bee barbed needle and it is modeled as a flexible thread within the framework of the Cosserat (micropolar) elasticity theory [37, 38, 39, 40, 41, 42]. The Cosserat elasticity is applied to describe the interaction between the needle and the human liver. Let us consider a serial surgical robot composed of a revolute joint and a flexible needle. A Lagrange frame $\left(X,Y,Z\right)$ of base vectors $\left(e_1,e_2,e_3\right)$ and origin $O$ in the entry point of the skin is attached to the robot (Figure 3). The Euler frame $K\left(x,y,z\right)$ with origin in the joint and the base vectors $\left(d_1,d_2,d_3\right)$ is attached to the needle. The angle between the flexible arm and axis $x$ is $\mathrm{\theta }$. Bending and torsion of the needle are described by the strain functions $\left(u_1,u_2,u_3\right)$. The robot has $f$ degrees of freedom $f=f_r+f_e$, where $f_r=1$ is the generalized coordinate of the rigid system and $f_e=3$ are the degrees of freedom of the flexible needle [43, 44].

Figure 3.

Lagrange coordinate system $\mathit{\text{OXYZ}}$ and the Euler coordinate system $\mathit{\text{oxyz}}$ attached to the needle.

The Euler axes are oriented with respect to the Lagrange axes by Euler angles $\mathrm{\upsilon },\mathrm{\psi }$ and $\mathrm{\phi }$ [43].

$$\begin{array}{c}{d}_1=\left(-sin\mathrm{\psi }sin\mathrm{\phi }+cos\mathrm{\psi }cos\mathrm{\phi }cos\mathrm{\upsilon }\right)e_1+\\ +\left(cos\mathrm{\psi }sin\mathrm{\phi }+sin\mathrm{\psi }cos\mathrm{\phi }cos\mathrm{\upsilon }\right)e_2-sin\mathrm{\upsilon }cos\mathrm{\phi }{e}_3,\end{array}$$

$$\begin{array}{c}{d}_2=\left(-sin\mathrm{\psi }cos\mathrm{\phi }-cos\mathrm{\psi }sin\mathrm{\phi }cos\mathrm{\upsilon }\right)e_1+\\ +\left(cos\mathrm{\psi }cos\mathrm{\phi }-sin\mathrm{\psi }sin\mathrm{\phi }cos\mathrm{\upsilon }\right)e_2+sin\mathrm{\upsilon }sin\mathrm{\phi }{e}_3,\end{array}$$E1

$$d_3=sin\mathrm{\upsilon }cos\mathrm{\psi }{e}_1+sin\mathrm{\upsilon }sin\mathrm{\psi }{e}_2+cos\mathrm{\upsilon }{e}_3.$$

The functions $\left(u_1,u_2,u_3\right)$ measure the bending and torsion of the needle as

$$\begin{array}{l}{u}_1={\mathrm{\upsilon }}^{\prime }sin\mathrm{\phi }-{\mathrm{\psi }}^{\prime }sin\mathrm{\upsilon }cos\mathrm{\phi },\\ u_2={\mathrm{\upsilon }}^{\prime }cos\mathrm{\phi }+{\mathrm{\psi }}^{\prime }sin\mathrm{\upsilon }sin\mathrm{\phi },\end{array}$$E2

$$u_3={\mathrm{\phi }}^{\prime }+{\mathrm{\psi }}^{\prime }cos\mathrm{\upsilon },$$

where $\left(\text{'}\right)$ means the partial differentiation with respect to s which is the coordinate along the central line of the needle. The functions $u_1$ and $u_2$ describe the bending of the needle, and the function $u_3$ the torsion of the needle. In addition, $u_1$ and $u_2$ are components of the curvature $\mathrm{\kappa }$ of the central line corresponding to the planes $\left(\mathit{yz}\right)$ and $\left(\mathit{xz}\right)$

$${\mathrm{\kappa }}^2=u_1^2+u_2^2={{\mathrm{\upsilon }}^{\prime }}^2+{{\mathrm{\psi }}^{\prime }}^2{sin}^2\mathrm{\upsilon }.$$E3

The function $u_3$ measures the torsion $\mathrm{\tau }$ of the needle

$$u_3=\mathrm{\tau }={\mathrm{\phi }}^{\prime }+{\mathrm{\psi }}^{\prime }cos\mathrm{\upsilon }.$$E4

So, the needle is rigid along the tangential direction and the total length of the needle $l$ is invariant, the ends being fixed by the force $F=-f$ with $f=\left(f_1,f_2,f_3\right)$. This force describes the contact between the needle and the tissue $f=p_{c}n=p_c{g}^{\prime }$, where $g$ is a gap function with respect to $s$.

The link between the position vector $r=\left(x,y,z\right)$ and unit tangential vector $d_3$ is $r=\underset{0}{\overset{s}{\int }}{d}_3\mathrm{d}s$, or

$$x\left(s\right)=\underset{0}{\overset{s}{\int }}cos\mathrm{\psi }sin\mathrm{\upsilon }\mathrm{d}s,y\left(s\right)=\underset{0}{\overset{s}{\int }}sin\mathrm{\psi }sin\mathrm{\upsilon }\mathrm{d}s,z\left(s\right)=\underset{0}{\overset{s}{\int }}cos\mathrm{\upsilon }\mathrm{d}s.$$E5

We introduce the inertia of the needle characterized by

$$\left({\mathrm{\rho }}_0{A}_0\right)\left(s\right),\left({\mathrm{\rho }}_0{I}_{\mathbf{1}}\right)\left(s\right),\left({\mathrm{\rho }}_0{I}_2\left(s\right)\right),$$E6

where ${\mathrm{\rho }}_0$ is the mass density per unit volume, A₀ 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 equations which describe the deformation are

$$-\mathrm{\rho }\ddot{r}-{\mathrm{\lambda }}^{\prime }=0,$$E7

$$\begin{array}{l}{k}_1\left({\dot{\mathrm{\psi }}}^{2}sin\mathrm{\upsilon }cos\mathrm{\upsilon }-\ddot{\mathrm{\upsilon }}\right)-k_2\left(\dot{\mathrm{\phi }}+\dot{\mathrm{\psi }}cos\mathrm{\upsilon }\right)\dot{\mathrm{\psi }}sin\mathrm{\upsilon }-\\ -A\left({{\mathrm{\psi }}^{\prime }}^{2}sin\mathrm{\upsilon }cos\mathrm{\upsilon }-{\mathrm{\upsilon }}^{″}\right)+C\left({\mathrm{\phi }}^{\prime }+{\mathrm{\psi }}^{\prime }cos\mathrm{\upsilon }\right){\mathrm{\psi }}^{\prime }sin\mathrm{\upsilon }-\\ -{\mathrm{\lambda }}_{1}cos\mathrm{\upsilon }cos\mathrm{\psi }-{\mathrm{\lambda }}_{2}cos\mathrm{\upsilon }sin\mathrm{\psi }+{\mathrm{\lambda }}_{3}sin\mathrm{\upsilon }=0,\end{array}$$E8

$$\begin{array}{l}-\frac{\partial }{\partial t}\left\{k_1\dot{\mathrm{\psi }}{sin}^2\mathrm{\upsilon }+k_2\left(\dot{\mathrm{\phi }}+\dot{\mathrm{\psi }}cos\mathrm{\upsilon }\right)cos\mathrm{\upsilon }\right\}+\\ +\frac{\partial }{\partial s}\left\{A{{\mathrm{\psi }}^{\prime }}^2{sin}^2\mathrm{\upsilon }+C\left({\mathrm{\phi }}^{\prime }+{\mathrm{\psi }}^{\prime }cos\mathrm{\upsilon }\right)cos\mathrm{\upsilon }\right\}+\\ +{\mathrm{\lambda }}_{1}sin\mathrm{\upsilon }sin\mathrm{\psi }-{\mathrm{\lambda }}_{2}sin\mathrm{\upsilon }cos\mathrm{\psi }=0,\end{array}$$E9

$$-k_2\frac{\partial }{\partial t}\left(\dot{\mathrm{\phi }}+\dot{\mathrm{\psi }}cos\mathrm{\upsilon }\right)+C\frac{\partial }{\partial s}\left({\mathrm{\phi }}^{\prime }+{\mathrm{\psi }}^{\prime }cos\mathrm{\upsilon }\right)=0.$$E10

where A and C are the bending stiffness and respectively the torsional stiffness of the needle, related to the Lame constants $\mathrm{\lambda }$, $\mathrm{\mu }$ by $A=\frac{1}{4}\mathrm{\pi }{a}^{4}E,C=\frac{1}{2}\mathrm{\pi }{a}^4\mathrm{\mu }$, $A=\frac{1}{4}\mathrm{\pi }{a}^{4}E,C=\frac{1}{2}\mathrm{\pi }{a}^4\mathrm{\mu }$, where $E=\frac{\mathrm{\mu }\left(3\mathrm{\lambda }+2\mathrm{\mu }\right)}{\mathrm{\lambda }+\mathrm{\mu }}$ is the Young’s elastic modulus, and a is the radius of the cross section of the needle, and

$$\mathrm{\rho }=A_0{\mathrm{\rho }}_0=\mathrm{\pi }{a}^2{\mathrm{\rho }}_0,k_1=I_1{\mathrm{\rho }}_0=\frac{\mathrm{\pi }{a}^4}{4}{\mathrm{\rho }}_0,k_2=I_2{\mathrm{\rho }}_0=\frac{\mathrm{\pi }{a}^4}{2}{\mathrm{\rho }}_0.$$E11

The Eqs. (7–11) are solved by using the cnoidal method [43].

In short, this method is reducible to a generalization of the Fourier series with the cnoidal functions as the fundamental basis function. This is because the cnoidal functions are much richer than the trigonometric or hyperbolic functions, that is, the modulus $m$ of the cnoidal function, $0\le m\le 1$, can be varied to obtain a sine or cosine function $\left(m\cong 0\right)$, a Stokes function $\left(m\cong 0.5\right)$ or a solitonic function, sech or tanh.

To understand the cnoidal method, consider now a nonlinear system of equations that govern the motion of a dynamical system

$$\frac{{\mathrm{d\theta }}_i}{\mathrm{d}t}=F_i\left({\mathrm{\theta }}_1,{\mathrm{\theta }}_2,\dots ,{\mathrm{\theta }}_n\right),i=1,\dots ,n,n\ge 3,$$E12

with $x\in {\mathrm{R}}^n$, $t\in \left[0,T\right]$, $T\in R$, where $F$ may be of the form

$$\begin{array}{l}{F}_i=\sum _{p=1}^n{a}_{\mathit{ip}}{\theta }_p+\sum _{p,q=1}^n{b}_{\mathit{ipq}}{\theta }_p{\theta }_q+\sum _{p,q,r=1}^n{c}_{\mathit{\text{ipqr}}}{\theta }_p{\theta }_q{\theta }_r+\\ +\sum _{p,q,r,l=1}^n{d}_{\mathit{\text{ipqrl}}}{\theta }_p{\theta }_q{\theta }_r{\theta }_l+\sum _{p,q,r,l,m=1}^n{e}_{\mathit{\text{ipqrlm}}}{\theta }_p{\theta }_q{\theta }_r{\theta }_l{\theta }_m+\dots ,\end{array}$$E13

where $i=1,2,\dots ,n$, and $a,b,c\dots$ constants.

This system of equations can be reduced to Weierstrass equations of the type

$${\dot{\mathrm{\theta }}}_2=P_n\left(\mathrm{\theta }\right),$$E14

We introduce the function transformation

$$\mathrm{\theta }=2\frac{{\mathrm{d}}^2}{\mathrm{d}t}log{\mathrm{\Theta }}_n\left(t\right),$$E15

where the theta function ${\mathrm{\Theta }}_n\left(t\right)$ are defined as

$${\mathrm{\Theta }}_1=1+exp\left({\mathrm{i\omega }}_{1}t+B_{11}\right),$$

$${\mathrm{\Theta }}_2=1+exp\left({\mathrm{i\omega }}_{1}t+B_{11}\right)+exp\left({\mathrm{i\omega }}_{2}t+B_{22}\right)+exp\left({\mathrm{\omega }}_1+{\mathrm{\omega }}_2+B_{12}\right),$$

$$\begin{array}{c}{\mathrm{\Theta }}_3=1+exp\left({\mathrm{i\omega }}_{1}t+B_{11}\right)+exp\left({\mathrm{i\omega }}_{2}t+B_{22}\right)+\\ +exp\left({\mathrm{i\omega }}_{3}t+B_{33}\right)+exp\left({\mathrm{\omega }}_1+{\mathrm{\omega }}_2+B_{12}\right)+\\ +exp\left({\mathrm{\omega }}_1+{\mathrm{\omega }}_3+B_{13}\right)+exp\left({\mathrm{\omega }}_2+{\mathrm{\omega }}_3+B_{23}\right)+\\ +exp\left({\mathrm{\omega }}_1+{\mathrm{\omega }}_2+{\mathrm{\omega }}_3+B_{12}+B_{13}+B_{23}\right),\end{array}$$E16

and

$${\mathrm{\Theta }}_n=\sum _{M\in \left(-\infty ,\infty \right)}exp\left(\mathrm{i}\sum _{i=1}^n{M}_i{\mathrm{\omega }}_{i}t+\frac{1}{2}\sum _{i<j}^n{B}_{\mathit{ij}}{M}_i{M}_j\right),$$E17

$$exp{B}_{\mathit{ij}}={\left(\frac{{\mathrm{\omega }}_i-{\mathrm{\omega }}_j}{{\mathrm{\omega }}_i+{\mathrm{\omega }}_j}\right)}^2,exp{B}_{\mathit{ii}}={\mathrm{\omega }}_i^2.$$E18

Further, we write the solution (15) under the form

$$\mathrm{\theta }\left(t\right)=2\frac{{\partial }^2}{\partial t^2}log{\mathrm{\Theta }}_n\left(\mathrm{\eta }\right)={\mathrm{\theta }}_{\mathit{lin}}\left(\mathrm{\eta }\right)+{\mathrm{\theta }}_{\mathit{int}}\left(\mathrm{\eta }\right),$$E19

for $\mathrm{\eta }=-\mathrm{\omega }t+\mathrm{\varphi }$. The first term ${\mathrm{\theta }}_{\mathit{lin}}$ represents a linear superposition of cnoidal waves. Indeed, after a little manipulation and algebraic calculus, obtain

$${\mathrm{\theta }}_{\mathit{lin}}=\sum _{l=1}^n{\mathrm{\alpha }}_l\left[\frac{2\mathrm{\pi }}{K_l\sqrt{m_l}}\sum _{k=0}^{\infty }{\left[\frac{q_l^{k+1/2}}{1+q_l^{2k+1}}cos\left(2k+1\right)\frac{{\mathrm{\pi \omega }}_{l}t}{2K_l}\right]}^2\right].$$E20

In (20) we recognize the expression [43].

$${\mathrm{\theta }}_{\mathit{lin}}=\sum _{l=1}^n{\mathrm{\alpha }}_l{\mathrm{cn}}^2\left[{\mathrm{\omega }}_{l}t;m_l\right],$$E21

with

$$q=exp\left(-\mathrm{\pi }\frac{K^{\prime }}{K}\right),$$

$$K=K\left(m\right)+\underset{0}{\overset{\mathrm{\pi }/2}{\int }}\frac{\mathrm{d}u}{\sqrt{1-m{sin}^{2}u}},$$

$$K^{\prime }\left(m_1\right)=K\left(m\right),m+m_1=1.$$

The second term ${\mathrm{\theta }}_{\mathrm{int}}$ represents a nonlinear superposition or interaction among cnoidal waves. We write this term as

$$2\frac{{\mathrm{d}}^2}{\mathrm{d}t}log\left[1+\frac{F\left(t\right)}{G\left(t\right)}\right]\approx \frac{{\mathrm{\beta }}_k{\mathrm{cn}}^2\left(\mathrm{\omega }t,m_k\right)}{1+{\mathrm{\gamma }}_k{\mathrm{cn}}^2\left(\mathrm{\omega }t,m_k\right)}.$$E22

If $m_k$ take the values $0$ or $1$, the relation (22) is directly verified. For $0\le m_k\le 1$, the relation is numerically verified with an error of $\left|e\right|\le 5\times {10}^{-7}$. Consequently, we have

$${\mathrm{\theta }}_{\mathit{\text{nonlin}}}=\frac{{\sum }_{k=0}^n{\mathrm{\beta }}_k{\mathrm{cn}}^2\left[\mathrm{\eta };m_k\right]}{1+{\sum }_{k=0}^n{\mathrm{\lambda }}_k{\mathrm{cn}}^2\left[\mathrm{\eta };m_k\right]}.$$E23

As a result, the cnoidal method yields to solutions consisting of a linear superposition and a nonlinear superposition of cnoidal waves.

Therefore, by applying the cnoidal method, the closed form solutions of the Euler angles $\mathrm{\theta },\mathrm{\psi }$ and $\mathrm{\phi }$ are obtained [43].

$$\begin{array}{c}cos\mathrm{\upsilon }=\mathrm{\zeta }={\mathrm{\zeta }}_2-\left({\mathrm{\zeta }}_2-{\mathrm{\zeta }}_3\right){\mathrm{cn}}^2\left(\sqrt{\frac{\left|{\mathrm{\lambda }}_3\right|}{2A}\left({\mathrm{\zeta }}_1-{\mathrm{\zeta }}_3\right)}\left(\mathrm{\xi }-{\mathrm{\xi }}_3\right),m\right)=\\ ={\mathrm{\zeta }}_2-\left({\mathrm{\zeta }}_2-{\mathrm{\zeta }}_3\right){\mathrm{cn}}^2\left[w\left(\mathrm{\xi }-{\mathrm{\xi }}_3\right),m\right],\end{array}$$E24

where $m=\frac{{\mathrm{\zeta }}_2-{\mathrm{\zeta }}_3}{{\mathrm{\zeta }}_1-{\mathrm{\zeta }}_3}$ and $w=\sqrt{\frac{\left|{\mathrm{\lambda }}_3\right|}{2A}\left({\mathrm{\zeta }}_1-{\mathrm{\zeta }}_3\right)}$,

$$\begin{array}{l}\mathrm{\psi }=\frac{1}{4{\left(A-k_1{v}^2\right)}^2{w}^2}\left\{-\frac{\beta +\left(C-k_2{v}^2\right)\tau }{1-{\zeta }_3}\prod \left[w\left(\mathrm{\xi }-{\mathrm{\xi }}_3\right),\frac{{\mathrm{\zeta }}_2-{\mathrm{\zeta }}_3}{1-u_3},m\right]-\right.\\ \left.-\frac{\mathrm{\beta }-\left(C-k_2{v}^2\right)\mathrm{\tau }}{1+\mathrm{\zeta }{u}_3}\prod \left[w\left(\mathrm{\xi }-{\mathrm{\xi }}_3\right),\frac{{\mathrm{\zeta }}_2-{\mathrm{\zeta }}_3}{1+u_3},m\right]\right\},\end{array}$$E25

$$\begin{array}{c}\mathrm{\phi }=-\frac{\mathrm{\tau }\left[C-A-\left(k_2+k_1\right)v^2\right]}{A-k_1{v}^2}\mathrm{\xi }+\frac{1}{4{\left(A-k_1{v}^2\right)}^2{w}^2}\left\{\frac{\mathrm{\beta }+\left(C-k_2{v}^2\right)\mathrm{\tau }}{1-{\mathrm{\zeta }}_3}\times \right.\\ \left.\times \prod \left[w\left(\mathrm{\xi }-{\mathrm{\xi }}_3\right),\frac{{\mathrm{\zeta }}_2-{\mathrm{\zeta }}_3}{1-{\mathrm{\zeta }}_3},m\right]-\frac{\mathrm{\beta }-\left(C-k_2{v}^2\right)\mathrm{\tau }}{1+{\mathrm{\zeta }}_3}\prod \left(w\left(\mathrm{\xi }-{\mathrm{\xi }}_3\right),\frac{{\mathrm{\zeta }}_2-{\mathrm{\zeta }}_3}{1+{\mathrm{\zeta }}_3},m\right)\right\},\end{array}$$E26

with $\prod \left(x,z,m\right)=\underset{0}{\overset{x}{\int }}\frac{\mathrm{d}y}{1-z{\mathrm{sn}}^2\left(y,m\right)}$ the normal elliptic integral of the third kind. Functions ${\mathrm{\zeta }}_1,{\mathrm{\zeta }}_2,{\mathrm{\zeta }}_3$ are solutions of the equation

$$\frac{1}{2}{{\mathrm{\zeta }}^{\prime }}^2=a{\mathrm{\zeta }}^3+b{\mathrm{\zeta }}^2-a\mathrm{\zeta }+c,$$E27

$$a=-\frac{{\mathrm{\lambda }}_3}{A}\ne 0,b=\frac{1}{2A}\left(\mathrm{\gamma }-\frac{C^2{\mathrm{\tau }}^2}{A}\right),c=-\frac{1}{2A}\left(\mathrm{\gamma }-\frac{{\mathrm{\beta }}^2}{A}\right).$$E28

Our objective is to determine the functions which measure the bending of the needle ($u_1$ and $u_2$), and the torsion $\left(u_3\right)$. To visualize the strain profile of the needle, we chose two routes (Figure 4). For the first route the tumor is red and the entry point is A. The second route is restricted by the presence of blood vessels that should not be touched and has the tumor (blue) with entry at point B.

Figure 4.

Two needle trajectories: For first route the tumor is red and the entry point A, and for the second route tumor is blue and the entry point B.

Figures 5 and 6 show that the deformation of the needle for both routes. We see that the deformation is small with no tendency to chaos. The strains are described by localized solitons which propagate for a long time without changes. The soliton is a localized wave with an infinite number of degrees of freedom. This wave conserves its properties even after interaction with another wave. In short, this wave acts somewhat like particles [43]. The system of Eqs. (7–11) has unique properties. These properties 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.

Figure 5.

Functions $u_1$ and $u_2$ for a) first route and b) second route.

Figure 6.

Function $u_3$ for a) first route and b) the second route.

3. Determination of the free-collision trajectories

Let us present a model for determining the collision-free trajectories by using the Fibonacci sequence. The trajectories are determined from the restrictions of avoiding the collisions with blood vessels, ribs and surrounding tissues, and also the interference of needles with each other. We consider that more needles are planned to be inserted into the liver [45].

The trajectory of each needle $j=1,\dots ,n$, is defined as a set of segments connecting the insertion point with the tumor. Two binary control parameters are introduced on each needle. The first parameter is the length of the $k^{\mathrm{th}}$ segment of the $j^{\mathrm{th}}$ needle’s trajectory, $l_{\mathit{kj}}=f_k/r^{\mathit{kj}}$, where $f_k$ is the $k^{\mathrm{th}}$ Fibonacci number and $r_{\mathit{kj}}$ a scaling number. The second control parameter is the angle $\mathrm{\omega }$ between the current needle and the previous one, $\mathrm{\omega }\in \left(0,\mathrm{\pi }\right)$ [46, 47].

The sequence terms in the Fibonacci problem are 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144. It is clear that each term is a sum of two preceding sequence term. In fact, the sequence can be recursively defined in the form $f_n=f_{n-1}+f_{n-2}$, $f_0=0,f_1=1$. The limit of the ratio of two consecutive terms in the Fibonacci string tends to the gold ratio $\mathrm{\phi }=\left(1+\sqrt{5}\right)/2$ [48, 49].

The Fibonacci sequence was highlighted in nature for example in the arrangement of the flower petals, in seeds and in the spiral arrangement of pine cones and pineapple. The keys on a piano are divided into Fibonacci numbers, and numerous classical compositions implement the golden section. Such an example is found in the Alleluia Choir in Handel’s Messiah and in many of Chopin’s preludes [50].

Development of a target drug delivery technique usually consisted of three steps. In step 1, a 2D ultrasound image of the tumor is obtained. The size and position of the tumor are analyzed by the surgeon who decides the number of required needles. In step 2, the insertion points on the skin of the needles and the positions of the targets in the tumor are established. In step 3, once the needles are placed at the target points on the skin, they are guided according to a precise surgical planning based on the optimization of the collision-free trajectories to avoid the ribs, blood vessels and tissues in the abdominal area.

The base coordinate frame is established on the first needle with the vertical axis $\mathit{Oz}$ and the origin in the insertion point (Figure 7). The workspace boundary for the needle is a 1D curve in each $\mathit{xy}$- plane means the $z_h$-plane for $z=h$, $h\in \left[0,h_{\max }\right]$ with a clinical value $h=200\mathrm{mm}$. The $S$ is the plan of the needle insertion trajectory, $\left({\mathrm{\theta }}_k,{\mathrm{\phi }}_k\right)$ are the rotation angles with respect to $y$ and $z$ axes ${\mathrm{\theta }}_{\min }\le {\mathrm{\theta }}_k\le {\mathrm{\theta }}_{\max }$, ${\mathrm{\varphi }}_{\min }\le {\mathrm{\phi }}_k\le {\mathrm{\phi }}_{\max }$ [45].

Figure 7.

The base coordinate frame.

The first control parameter is the length of the $k^{\mathrm{th}}$ segment of the $j^{\mathrm{th}}$ needle, $l_{\mathit{kj}}=f_k/r^{\mathit{kj}}$, $k=1,\dots ,m_{\mathit{lk}}$, $j=1,\dots ,n$ where $f_k$ is the $k^{\mathrm{th}}$ Fibonacci number

$$f_0=f_1=1,f_{k+2}=f_{k+1}+f_k,k\ge 0.$$E29

and $r_{\mathit{kj}}$ a scaling number. The second control parameter is the angle $\mathrm{\phi }$ between the current needle and the previous one, $\mathrm{\phi }\in \left(0,\mathrm{\pi }\right)$. The kinematic constraint of the $j^{\mathrm{th}}$ needle, $j=1,\dots ,n$, is given by

$$\mathrm{\varphi }=\left[\begin{array}{c}{x}_j-x_{0j}-{\mathrm{\delta }}_{\mathit{xj}}-\left(h-z_{0j}-{\mathrm{\delta }}_{\mathit{zj}}\right)tan{\mathrm{\theta }}_{j}cos{\mathrm{\phi }}_j\\ y_j-y_{0j}-y_{\mathit{yj}}-\left(h-z_{0j}-{\mathrm{\delta }}_{\mathit{zj}}\right)tan{\mathrm{\theta }}_{j}cos{\mathrm{\phi }}_j\end{array}\right]=0,$$E30

where ${\left[x_{0j}+{\mathrm{\delta }}_{\mathit{xj}},y_{0j}+{\mathrm{\delta }}_{\mathit{xj}},y_{0j}+{\mathrm{\delta }}_{\mathit{xj}}\right]}^T$ is the actual target of the tip of needle $\left(j-1\right)$, $\left({\mathrm{\theta }}_j,{\mathrm{\phi }}_j\right)$ are the rotation angles with respect to $y$ and $z$ axes ${\mathrm{\theta }}_{\min }\le {\mathrm{\theta }}_j\le {\mathrm{\theta }}_{\max }$, ${\mathrm{\varphi }}_{\min }\le {\mathrm{\phi }}_j\le {\mathrm{\phi }}_{\max }$ and $\left({\mathrm{\delta }}_{\mathit{xj}},{\mathrm{\delta }}_{\mathit{yj}},{\mathrm{\delta }}_{\mathit{zj}}\right)$ denote the deformation of the liver.

The choice of the scaling number $r_{\mathit{kj}}$ is done by a binary control

$$l_{\mathit{kj}}=‖x_k-x_{k-1}‖=\frac{u_{\mathit{kj}}{f}_k}{r^{\mathit{kj}}}.$$E31

The possible collision point between the needle and the tissue is analyzed by an identifier to check the minimum distance between needle and the surrounding tissue [45]. The minimum distance is expressed as

$$min\left(\frac{1}{2}{\left(r_1-r_2\right)}^T\left(r_1-r_2\right)\right),$$E32

with ${\mathrm{g}}_1\left(r_1\right)\le 0$, ${\mathrm{g}}_2\left(r_2\right)\le 0$, $r_1$, $r_2$, the position vectors of two points belonging to the needle and the tissue, respectively, and $g_1$, $g_2$, the surfaces to the needle and the tissue, respectively. The interference distance or penetration is defined as

$$min\left(-d\right),{\mathrm{g}}_1\left(r_1\right)\le -\frac{d}{2},{\mathrm{g}}_2\left(r_2\right)\le -\frac{d}{2},$$E33

where $d$ is the penetration. The configuration of the collision-free trajectories of each needle is defined as a sequence of trajectories corresponding to a particular choice for the control $\left(u,v\right)$. The set of all collision-free trajectories are computed based on the kinematic constraint (30) as

$$R=\left\{\sum _{k=0}^n\frac{u_k{f}_k}{r^k}exp\left(-\mathrm{i\omega }\sum _{j=0}^k{v}_j|\left(v_j\right)\in \left\{0,1\right\}\right.\right\}.$$E34

As application, the case of a tumor with a difficult location in the vicinity of the portal tree of the vascular territory in the liver, is considered (Figure 8a). The tumor image seen on the microscope is shown in Figure 8b. White and gray denote forbidden areas while the shade of purple are safe regions. The tumor is drawn in red. The Fibonacci algorithm is applied to three-needles with restrictions to avoid the collision with the tissues, blood vessels, ribs and previously inserted needles.

Figure 8.

a) Location of the tumor; b) tumor image seen on the microscope.

The task of our simulation is to determine the boundaries of each needle as a collision-free surface which represents the feasible insertion area based on given constraints. Then, the optimal trajectory of each needle can be chosen in this surface automatically.

Once all needles are placed at the predetermined epidermis and the ordering of entry is chosen to be 1, 2 and 3, following the first needle’s insertion, the operation is repeated for the needles 2 and 3.

Simulation of the collision-free trajectories for three needles is presented in Figure 9. The insertion scheme is determined by the Fibonacci spirals. A set of free-collision trajectories (red) in the immediate vicinity of the epidermis, is suggested.

Figure 9.

Simulation of the collision-free trajectories for three-needles.

From these possible collision-free trajectories (red) the green paths corresponding to the Fibonacci spirals (black) are chosen. These trajectories avoid the blood vessels (purple) and the coasts (brown) in all directions until the tumor. Optimal solution with 3 collision-free trajectories to the target is displayed in Figure 10. The Fibonacci spirals with the centers in ribs (brown) are displayed for needles 1 and 3. For the needle 2, the Fibonacci.

Figure 10.

Optimal solution with 3 collision-free trajectories to the target.

Three locally optimal collision-free trajectories for the surgical needle corresponding to three different entry points into the skin A, B and C are displayed in Figure 11.

Figure 11.

Three optimal collision-free trajectories for the needle robot.

4. Conclusions

The study investigates the navigation of a flexible needle into the human liver. The role of the needle is to deliver into the tumor an active chemotherapeutic agent. The deformation of the needle during its insertion into the human liver is describe in this chapter in the context of intraoperative treatment of liver tumors based on the integrated imaging-molecular diagnosis. The needle is a bee barbed type modeled as a flexible thread within the framework of the Cosserat (micropolar) elasticity theory. The Cosserat elasticity describes the interaction between the needle and the human liver by incorporating the local rotation of points and the couple stress as well as the force stress representing the chiral properties of the human liver.

Acknowledgments

This work was supported by a grant of the Romanian Ministry of Research and Innovation project PN-IIIP2-2.1-PED-2019-0085 CONTRACT 447PED/2020 (Acronim POSEIDON).

Conflict of interest

The authors of this paper certify that they have no affiliations with or involvement in any organization or entity with any financial or nonfinancial interest in the subject matter or materials discussed in this manuscript.