PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Trans Biomed Eng. Author manuscript; available in PMC 2010 April 23.
Published in final edited form as:
PMCID: PMC2859043
NIHMSID: NIHMS192959

Modeling and Control of Needles with Torsional Friction

Abstract

A flexible needle can be accurately steered by robotically controlling the bevel tip orientation as the needle is inserted into tissue. Friction between the long, flexible needle shaft and the tissue can cause a significant discrepancy between the orientation of the needle tip and the orientation of the base where the needle angle is controlled. Our experiments show that several common phantom tissues used in needle steering experiments impart substantial friction forces to the needle shaft, resulting in a lag of over 45° for a 10 cm insertion depth in some phantoms; clinical studies report torques large enough to cause similar errors during needle insertions. Such angle discrepancies will result in poor performance or failure of path planners and image-guided controllers, since the needles used in percutaneous procedures are too small for state-of-the-art imaging to accurately measure the tip angle. To compensate for the angle discrepancy, we develop an estimator using a mechanics-based model of the rotational dynamics of a needle being inserted into tissue. Compared to controllers that assume a rigid needle in a frictionless environment, our estimator-based controller improves the tip angle convergence time by nearly 50% and reduces the path deviation of the needle by 70%.

I. Introduction

Accurate needle tip placement is essential for many needle-based diagnostic and therapeutic procedures, such as biopsies and prostate brachytherapy [1], [2]. Guiding a maneuverable needle in real time can improve placement accuracy, enable access to multiple targets for a single insertion, and expand the range of procedures performed with needle interventions by providing the ability to avoid obstacles and reach previously inaccessible subsurface targets.

One mechanism for steering needles involves harnessing the unbalanced reaction forces imparted by tissue on an asymmetric-tip needle as it is inserted into tissue. For sufficiently flexible needles, these forces cause the needle to follow a nearly circular are of constant curvature [3]. The instantaneous direction of curvature can be controlled by rotating the needle about its shaft at the base of the needle (outside the tissue), thereby reorienting the needle tip. The needle tip can then be guided to a desired location inside the body while avoiding obstacles. See [4] for a survey that includes robotically steered needles and insertion modeling.

Specifically for prostate brachytherapy, the steerable needle could be precisely driven to a target where a radioactive seed is delivered. Without removing the needle tip from the tissue, the tip can then be steered to subsequent targets for further seed placement. Biopsy procedures typically employ large needles that cannot maneuver inside tissue. In this case, a small steerable needle could accurately reach the target and the larger biopsy needle could reach the target by sliding over the steerable needle. The steerable needle would then be removed and the tissue excised using standard methods.

When rotating (steering) the needle inside tissue, friction between the tissue and the inserted length of the needle results in the angle of the needle tip “lagging” the angle of the needle base as illustrated in Figure 1. We show that several tissues impart sufficient frictional forces to the needle shaft to cause significant discrepancies between the base and tip angles. During clinical needle insertions for prostate brachytherapy, Podder et al. [5] found torques large enough to cause a flexible needle to similarly twist in human tissues.

Fig. 1
Torsional friction

When steering a needle using an asymmetric tip, any inconsistency between base and tip orientations or misalignment of the needle tip can result in poor performance or failure of controllers and path planners. Reed et al. [6] experimentally demonstrate that using a pre-bent needle tip instead of a bevel-tip needle causes deviations in the tip location when the needle is rotated. The small discrepancy causes the robotically controlled needle to puncture the obstacles and miss the targets. Pre-bent needles are often used since they afford greater bending due to the larger tip asymmetry.

Current medical imaging technology does not enable direct measurement of the bevel angle around the needle shaft due to the small diameter of the needles used in clinical procedures. Ultrasound and other imaging techniques can achieve a resolution of 0.8 mm [7], which is similar to the diameter of needles used in many percutaneous procedures. At these resolutions, imaging alone is unable to determine the roll of the bevel tip. Kallem and Cowan [8] detail an image-based state estimation approach to recover the needle shaft orientation in real time, but their controller assumes that changes in the base angle directly relate to changes in the tip angle. They demonstrate that the real system converges slower than simulations, which they speculate is due to unmodeled torsion dynamics. Our torsion controller is designed to work with controllers such as their planar controller to decrease the convergence time.

Several papers [8], [9], [10], [11], [12], [13] mention the possibility that torsional friction along the inserted portion of the needle could be a problem, but the authors are not aware of any work quantifying or modeling this effect during continuous insertion. Abolhassani et al. [9] and Reed [10] describe a method of adjusting the needle rotation based on the torque measured at the base of the needle while the needle is not moving in translation. Several studies examine the friction required to insert a needle into tissue [14], [15], but the authors are not aware of any previous studies that analyze or model the torsional friction during continuous insertions.

This paper develops a fundamental model of a long thin object rotating inside a continuous medium. We demonstrate the practical significance of torque generated by the needle–tissue interface and formulate a model of the rotational dynamics that can be incorporated into path planners and position controllers for steerable needles like the system described in [6].

II. Quantifying Needle Lag

A. Materials and Methods

We use the device shown in Figure 2 for our measurements and experimental validation. The setup is capable of rotating and inserting a needle into tissue using two DC servo motors. A 6-axis sensor (ATI Nano 17) measures the forces and torques at the base of the needle.

Fig. 2
Experimental setup

For the experiments described throughout this paper, we use a solid 0.59 mm diameter Nitinol wire (Nitinol Devices and Components, Fremont, CA) with a bevel tip, unless stated otherwise. The needle’s polar moment of inertia is defined as J = πd4/32, so J = 0.012 mm4. We determine the shear modulus, G, by rigidly fixing the needle tip at a fixed angle while rotating the needle base. The remainder of the needle is unrestrained. The angular change (in radians) over the length l, with torque τ, is given by

θ=0lτ(x)dxJ(x)G(x).
(1)

A 90° rotation of the 41 cm needle resulted in a torque of 1.24 N-mm, so G is 27.2 GPa. The density of the Nitinol needle is 6.45g/cm3.

To quantify the effect of torsion, we inserted the needle through the five tissue types described in Table I. All tissues, except for the chicken, were homogenous. We inserted the needle in a straight path by rotating the needle 180° every 1 cm. A needle inserted in a highly curved path under large torsion would be more likely to buckle [16], but we have not seen this effect during our testing. After inserting, we instrumented the protruding needle tip with two magnets. A hall effect sensor (Honeywell Sensing and Control, Golden Valley, MN, USA) measured the angle of the needle tip with a resolution of at least 0.6° over the range used during these experiments. A platform maintained the proper distance between the hall effect sensor and the magnets on the needle during insertion and retraction. The platform and needle translated together. Extending the needle through the tissue eliminates the possible effect of the bevel tip on rotation, thus isolating the effect of shaft friction.

TABLE I
Tissue and Phantom Descriptions

We performed two types of needle rotation experiments: (1) the needle was not translating and (2) the needle was translating. In the non-translating experiments, the needle was inserted such that the distance between the needle base and tissue was 23 cm. Each material was tested with a tissue depth of 10 cm and the porcine gelatin and plastisol were additionally tested at 2.5, 5, and 7.5 cm. One test consisted of rotating the needle base by 90° and another consisted of rotating the needle base clockwise and counterclockwise by 90°. Both variations ended by retracting the needle 0.5 cm because we suspected that retracting the needle slightly would break the stiction, allowing the needle to release much of its torsional energy.

For the translating experiment, the needle base started 5 cm from the tissue, as shown in Figure 2, and retracted at 0.25 cm/sec. To reduce the effects of the needle buckling inside the support sheath, we retracted the needle instead of inserting. Buckling would have caused the needle to press against the support sheath and cause a further lag between the tip and the base. The needle tip is not cutting since it is outside the tissue, so the force should be the same during insertion and retraction. Since there is a possibility of added forces from the support sheath during insertion, we chose retraction to focus on the interaction forces. Future needle steering systems could be designed without the problems inherent in this support sheath. The plastisol was tested at depths of 5, 10, and 20 cm.

B. Tip Lag without Inserting the Needle

1) Tissue Effects

Torsional friction is a significant problem for all the phantoms listed in Table I. Figure 3 shows the lag associated with slowly rotating the needle base 90°. The plastisols have a lag of about 10° and show considerable creep once the base has stopped rotating. SimTest gel shows a similar, yet more significant, creep effect during the same time period. Porcine gelatin exhibits the most dramatic “stick-slip” behavior of the tissue samples tested, indicative of substantial stiction; for this material the needle almost immediately reaches steady state when the base rotation halts. Chicken has tissue structures relatively similar to human tissue, so it may provide results similar to human tissue. However, the chicken tested had been packaged in a fluid that may have significantly reduced its coefficient of friction. A lag of 3° is not likely to cause significant placement errors, but, as we will discuss in Section II-D, there is likely to be a larger and problematic angle lag in human tissues.

Fig. 3
Torsion behavior in multiple materials

2) Depth Effects

Torsional friction is a function of depth since more tissue interacts with the needle at deeper insertions [15]. Figure 4 shows the lag associated with rotating the needle base at four different depths in porcine gelatin. As expected, the depth of insertion affects the torsion lag.

Fig. 4
Torsion behavior for multiple insertion depths

The torque measured at the base also increases with the depth during both the constant angular velocity and also once the tip reaches a steady offset. Table II shows the average torque while the needle is not rotating in porcine gelatin (t = 20 to 30 sec) and the damping coefficient during the constant rotation in the plastisol (t = 10 to 20 sec). The static and damping forces increase linearly with depth. Compared to a least squares fit line, the R2 values are 0.98 and 0.97 for porcine and plasticol, respectively.

TABLE II
Rotational friction forces

Figure 4 shows other complications associated with torsion windup. At t = 30 sec, the base begins to rotate back to 0° while the tip also starts to move in the same direction as the base. The tip is briefly leading the base. At a depth of 10 cm, the needle tip has rotated by almost 20° before the base overtakes the tip. This lead effect is likely due to tissue elasticity: the tissue briefly pulls the needle in the same direction as the base. The tissue, rather than the base actuator, is causing the tip to rotate during this early reversal period.

3) Overcoming Torsional Stiction

The final motion in Figures 3 and and44 is the needle base retracting 0.5 cm. This motion breaks the stiction between the needle and tissue and allows the needle to snap near to the final base angle. The required retraction distance is a function of the elasticity of the tissue – more elastic materials will require a larger retraction.

One method to overcome the lag is, thus, to retract the needle a small distance after rotating and then reinsert the needle the same distance. Retracting the needle leaves a precut path. Upon reinsertion, the needle follows this path and returns to the previous location with the bevel tip at the desired orientation. Additional insertion causes the needle to follow a new path defined by the new tip angle. However, this method is not recommended when using a pre-bent needle in real tissue since the needle may cause additional tissue damage and may follow a different path upon reinsertion.

Another method to reduce the effect of stiction is to use dithering, which adds a high frequency periodic signal to the insertion [17]. Dithering would keep the controlled surface moving and thus prevent the controlled surface from encountering stiction. High frequency motions are likely unsuitable for needle steering due to the possible additional trauma associated with these vibrations. Some clinicians use a similar strategy as they insert needles by slowly twisting the needle back and forth during insertion.

C. Tip Lag with Needle Motion

Even when the needle is continuously moving through tissue, the tip still lags the base when rotated. Figure 5 shows the lag associated with rotating a needle in soft plastisol while the needle is being moved through tissue at 0.25 cm/sec. The base motion is tuned to provide the fastest response without oscillating. Similar to the non-translating case, more tissue contact causes a larger tip lag. Whereas the non-translating case remains at a constant offset angle due to static friction, continuously inserting the needle prevents any static friction effects since the needle is constantly sliding through the tissue, so the tip eventually reaches the desired angle.

Fig. 5
Torsion behavior during continuous needle motion

The insertion velocity used in this experiment is slower than a prostate brachytherapy performed by surgeons. However, Podder et al. [5] note that the slower speeds typically used in robotic needle insertions do not significantly affect the procedure compared to the higher velocities used by surgeons. The needle tip will eventually reach the base angle as long as the needle is being inserted fast enough to prevent stiction effects.

Rotating the needle 180° while inserting causes the needle tip to deviate from its current plane of motion. As the needle rotates about its axis, the tip angle points in an undesired direction, thus the needle moves away from the desired plane of motion. During the 2 seconds it takes for the needle tip to converge to the base angle, the needle will have moved forward and followed a path perpendicular to the desired plane. We demonstrate that this slow convergence causes a needle to deviate by 1.2 cm over a 7cm insertion and develop a model and controller to improve the performance of the tip convergence.

D. Clinical Implications

Accurate needle placement is vital to many medical diagnoses and treatments. Even small deviations in the placement of a needle tip during a biopsy can result in misdiagnosis or ineffective placement of radiation during a prostate brachytherapy. Many percutaneous procedures are performed by hand and the forces/torques are not recorded. However, Podder et al. [5] measured the forces and torques in vivo during several prostate brachytherapy needle interventions. Compared to the torque in the other two directions, the torque around the insertion axis appears minuscule, but the 7N-mm they measured is actually higher than the torque we measured in SimTest gel (1.0 N-mm), which had a lag of nearly 50°.

Using the torque measured at the base of a needle, Reed [10] demonstrates that the lag between the base and the tip of the needle is calculated as

θbaseθtip=τbaseJG(lin+lout2),
(2)

where lin is the distance from the needle base to the entry point in the tissue, lout is the distance from the needle base to the tip of the needle or the point where the needle exits the tissue, τbase is the torque measured at the base of the needle, G is the shear modulus of the needle, and J is the polar moment of inertia of the needle.

Based on (2), the 7 N-mm of torque reported by Podder et al. for a human intervention would cause a lag of about 11° in a 1.27 mm diameter Nitinol needle of the same 20 cm length inserted 5 cm into tissue. Podder et al. used larger 1.27 mm diameter stainless steel needles, which reduces the amount of lag due to increased diameter and thus, torsional stiffness. However, bevel and bent-tip steerable needles require the superelastic properties of materials such as Nitinol. Currently, there is no steering during a prostate brachytherapy since the large stainless steel needles are much more rigid than a steerable needle. When a flexible needle is used to allow the needle to steer, the torques will be large enough to cause significant errors and it will be necessary to compensate for the developed inaccuracies. Larger steerability will cause larger angle errors, thus increasing the need for torsion compensation.

III. Dynamic Needle Torsion Model

Currently, it is not possible to measure the angle of the needle tip inside tissue. The only angle measurement is at the base of the needle. At certain instances, the torque at the needle base directly relates to the steady-state lag at the needle tip as shown in (2), but this equation says nothing about the trajectory of the needle tip. Therefore, in order to estimate the motion of the needle tip for feedback control, we formulate a mechanics-based model of the rotational dynamics of a needle rotating inside tissue.

The experiments in Section II-C demonstrate that the needle tip eventually reaches the base angle when rotated while translating. Stiction has a minimal effect because the surface of the needle is continuously sliding past the tissue, but the slow angle convergence will cause the needle to move in an unspecified direction. Our model assumes that each part of the needle is moving relative to the local tissue such that the interaction with the tissue is viscous only and the effects of stiction will not be present.

A. Lumped Mass Model

As a brute force approach to modeling the needle inserting and rotating inside a tissue, we use a lumped mass model with viscous friction. This model requires many unnecessary states, but allows us to show that a more compact modal coordinate representation can be used.

Each lumped mass of the needle is connected to neighboring sections by a torsional spring and is connected to the ground through a damper. The dynamics of the lumped mass system are given by a vector second-order differential equation:

(ηlsIm×m)θ¨+(dIm×m)θ.+Kθ=[1,0,0,,0]Tu(t)
(3)

where θ is a vector of the angles of the m masses, ls is the length between masses, d is the damping exerted on each mass from the needle-tissue interaction, and u(t) is an input exerted only at the first mass. The following terms frequently appear together, so we define η = ρJ and κ = JG where ρ is the density of the needle and J and G are previously defined needle parameters. The stiffness matrix K is given by

K=[κκ0κ2κκκ2κκ0κκ].
(4)

For many types of damping, the modal coordinates are coupled through the damping terms and cannot be decoupled into independent states. However, for a proportionally (viscously) damped system, the matrix of eigenvectors decouples the system into modal coordinates [18]. Caughey and O’Kelly [19] show that systems have the same mode shapes (i.e., eigenvectors) as the same system without damping if they satisfy

KΛ1D=DΛ1K,
(5)

where K is the stiffness matrix, Λ is the inertia matrix, and D is the damping matrix.

The system described by (3) is viscously damped and satisfies (5), thus it can be described conveniently in modal coordinates as linear independent combinations of the normal modes, regardless of the amount of friction. Describing a system in modal coordinates reduces the number of states necessary to adequately describe the system. In Section III-E, we highlight the benefits of the modal model as we compare it with the lumped parameter model.

B. Continuous Model

Figure 6 shows our model of the needle inside tissue and the torque exerted on a small element of the needle. Writing Newton’s Second Law for this small element results in

Fig. 6
Needle model

τxdxβθtdx=η2θt2dx,
(6)

where τ is the torque exerted on this element from neighbor elements (internal torque), β is viscous damping exerted on a small element, t is time, and θ is the angle.

As shown in Section II-B2, the total lag increases with insertion depth since additional tissue interacts with the needle. Assuming homogeneous tissue, the total torsional friction increases linearly since each additional piece of tissue exerts an additional resistive force. Thus we model the friction at each small element as

β=bl,
(7)

where b is the effective damping coefficient along the entire needle-tissue interface of length l.

Mechanics principles show that the torque required to twist a shaft of length l by an angle θ is determined by (1). The needle (shaft) has constant material properties, so J(x) = J and G(x) = G. Taking the second partial derivative of (1) with respect to x results in

τx=κ2θx2.
(8)

Substituting (8) into (6) and integrating results in a homogenous partial differential equation that defines the motion of the needle through time and space:

κ2θx2=η2θt2+βθt
(9)

C. Forced Modal Analysis

Modal analysis is often used to analyze many structures such as cantilever beams [18], [20], [21]. The authors are not, however, aware of an analysis on long slender rods being controlled on one edge with distributed damping, such as needles inserted inside tissue. In fact, modal analysis of members in torsion is not studied nearly as extensively as cantilevered beams.

To control the orientation of the needle, a torque is exerted at the base. We represent the torque as a spatial impulse function:

v(x,t)=δ(x)u(t),
(10)

where δ is the Dirac Impulse function and u(t) is the input torque. The forced partial differential equation is

η2θt2+βθtκ2θx2=v(x,t).
(11)

We assume the solution is separable in time and space, so we can write the solution as a standard Fourier Series expansion:

θ(x,t)=12ψ0(x)q0(t)+k=1ψk(x)qk(t),
(12)

where ψk(x) is typically taken as the kth mode shape, and qk(t) is the associated time dependent factor for each mode.1 The zeroth mode corresponds to a constant angle along the spatial coordinate, which is directly measured at the base of the needle. The general solution for each mode shape is

ψk(x)=cos(ωkx),
(13)

where ωk is the frequency of the kth mode and

ωk=kπlk=0,1,2,3,
(14)

Plugging (12) into (11) results in

12(ηψ0q¨0+βψ0q.0)+k=1{ηψkq¨k+βψkq.kκψkqk}=v(x,t).
(15)

We apply the orthogonality principle by multiplying (15) by an arbitrary mode shape (ψs) and integrating over the length. Since each mode is orthogonal to every other mode, we have

0lψk(x)ψs(x)dx={l2,k=s>0l,k=s=00,otherwise
(16)

and

0lψk(x)ψs(x)dx={l2ωk2,k=s0,otherwise.
(17)

We use the sifting property of the Dirac function at x = 0:

0lv(x,t)ψs(x)dx=u(t)0lδ(x)ψs(x)dx=u(t),
(18)

which means that any input at the base excites all cosine modes; the sinusoidal modes are not excited with this base input and similarly do not have any effect on the tip angle. The final differential equation for the kth mode is then

ηq¨k(t)+βq.k+κωk2qk(t)=2lu(t).
(19)

Truncating to n modes, we have

Mq¨+Dq.+Kq=Pu(t)
(20)

where

M=diag(η,η,,η)
(21a)
D=diag(β,β,,β)
(21b)
K=diag(0,κω12,κω22,,,κωn12)
(21c)
P=2l[1,1,,1]T
(21d)
q=[q0(t),q1(t),q2(t),,qn1(t)]T
(21e)

Each state (q) in the modal system represents a mode shape, so observability and controllability correspond to the ability to observe and control the modes of the system. The system can be shown to be fully controllable and observable for any number of modes.

D. Inertialess Modal Model in State Space Form

In the case of needle steering, we observe that the system is over-damped and the inertia of the needle may be insignificant compared to the stiffness and damping. To test the effect of needle inertia, we formulate an inertialess system following the same procedure as above, where η is considered zero:

x.=Ax+Bu
(22a)

y=Cq
(22b)

where

A=D1K
(23a)

B=D1P
(23b)

C=[ψ0(l),ψ1(l),ψ2(l),,ψk1(l)]
(23c)

where D, K, and P are defined in (21).

E. Model Comparison

Table III shows the eigenvalues of the lumped mass model with both 64 and eight masses and of the fourth-order modal models with and without inertia for a 20 cm insertion. The slowest poles of the modal model with mass are within 0.001 % of the poles of the inertialess modal model. The eigenvalues of the 64 mass lumped system are within 0.5% of the modal models. For an eight mass lumped system, the eigenvalues deviate by more than 10%. The system shows similar results at other depths. Using a modal model in this case allows the system to be more accurately described using fewer states than the lumped parameter model. Without sacrificing accuracy, the inertialess modal model further simplifies the system representation, which is beneficial for realtime control of the needle tip. Thus, we use the inertialess modal formulation.

TABLE III
Comparison of Models for 20 cm Insertion

The Hankel singular values of a system quantify the contribution of each state to the system behavior. The ratio of the minimum and maximum values provide a good measure of where to truncate the system. Table III shows that a fifth-order system has a min/max ratio on the order of 10−3, which means the fifth mode only accounts for 0.1 % of the system dynamic response. The contribution of the fifth mode decreases for shorter needles, so a fourth-order model is sufficient to adequately describe and control the system. Thus, we use a fourth-order inertialess modal model to control the angle of the needle tip as described in Section V.

The vital parameter for torsion compensation is the angle at the tip of the needle, so our analysis focuses on the settling time for the tip angle. Settling time is a reasonable comparison since we are mainly focused on the rotations necessary to steer the needle. Slow settling time causes the needle to deviate from the desired path.

IV. Needle-Tissue Interaction Model

The true behavior of a needle during insertion is dependent not only on the needle dynamics, but also on the dynamics of the tissue. We use the finite element model shown in Figure 7 to simulate the tissue properties and the needle-tissue interaction so we can determine if tissue effects are significant. We assume the elements are small enough that tissue interaction torques along each element occur discretely at the associated nodes. The angle between each node is based on (1). The properties are constant throughout the needle, so the angle between two neighboring nodes, i and i + 1, is

Fig. 7
Finite element model

θi+1θi=τlJG
(24)

where τ is the torque exerted from neighboring nodes and the damping along the element, l is the length between nodes, and J and G are the needle’s material properties determined in Section II-A. The input torque is exerted only to the base node and tissue friction is exerted along each element as described in Section IV-A1.

The simulated needle properties and lengths match the experimental setup for the 20 cm insertion described in Section II-C. The FE model was constructed of 200 disk elements inside the tissue, each of length 1mm, and one element of length 50 mm outside the tissue. The element outside the tissue was simply treated as a long torsional member since no external forces acted along its length. The simulation was run with a time step of 10−7 s. Decreasing the time step and/or elements does not significantly change the results.

A. Material Characteristics

1) Needle-Tissue Interaction

In previous studies [10], [15], the interaction between the needle and the tissue was modeled as Karnopp friction, which is composed of dynamic friction everywhere except for a constant static friction within a small velocity near zero. As demonstrated in Figure 5, the needle tip eventually reaches the desired angle. Static friction causes the needle to maintain a constant offset. Since the needle is constantly inserted through the tissue, the static friction region is never reached, so we model the dynamic friction as viscous damping.

We experimentally determine the viscous damping. During steady state needle rotation, all applied torque from the tissue is from damping, so we can use

τb=bθ.
(25)

to estimate the damping coefficient, b. Equation (25) corresponds to the zeroth mode from the modal model (Section III-C). As shown in Table II, the total rotational damping force increases linearly with depth, much like how the insertion force increases linearly with depth [15]. For the 20 cm insertion modeled here, we use a total damping coefficient of b = 3.45 N-mm-s. The total damping is divided equally between each of the elements inside the tissue, so the torque exerted due to friction at each needle element interacting with the tissue, τfi, is

τfi=bn(θ.iØ.i)
(26)

where Ø̇i is the angular velocity of the tissue touching the ith element and n is the number of nodes inside the tissue.

2) Tissue Model

To simulate the tissue properties, we use the Kelvin-Voigt model, which uses a spring and damper in parallel to model the creep characteristics of tissue [22]. We use the Rheometrics Solids Analyzer II (RSA) experiments performed by Misra et al. [23] for a similar plastisol gel to determine the tissue stiffness and damping parameters. The RSA test measures the force as a displacement is applied to a known surface area of the material. This method determines the elasticity and creep of the tissue. Misra et al. show that the tissue elasticity (Et) is 45.2 kPa and Poisson’s ratio is 0.45. Mechanics principles show that the shear modulus of the tissue (Gt) is calculated as

Gt=Et2(1+ν)=15.6kPa.
(27)

As the needle rotates inside the tissue, some of the soft tissue also moves. For this simulation, we assume the needle rotation affects tissue up to 0.5 mm away, roughly the size of the needle diameter. We then use the shear modulus to approximate the tissue as a spring surrounding each element of the needle. For a force applied perpendicular to the tissue (such as the needle rotating), mechanics principles show that

FΔx=GtAh
(28)

where F is the applied force from the tissue, Δx is the distance the surface of the tissue moves, h is the assumed height of the affected tissue, and A is the surface area of the needle-tissue interface. FΔx is an equivalent spring for the tissue over each element. Using the surface area of each needle element inside tissue (A = 1.9 mm2), the interaction height (h = 0.5 mm), and the radius of the needle (r = 0.295 mm), the tissue surrounding each node of the needle is approximated by a rotational spring (kt) where

kt=Fr2Δx=5.2×103Nmmrad.
(29)

The RSA test linearly increases the strain on the tissue for 5 seconds and then maintains a constant position for 5 seconds. During the 5 seconds at a constant position, the measured force decreases as the material creeps with a time constant of 0.96. Since the Kelvin-Voigt model is a spring and damper in series, the time constant, a, is

a=ktbt
(30)

where bt is the tissue damping. The damping coefficient (bt) for each element is then solved to be 5.4 × 10−3N-mm-s/rad.

B. FE Simulation

We ran two simulations of the needle with all parameters held constant except for the stiffness and damping of the tissue. We ran the FEM with an input similar to that used in the experimental tests. Figure 8 shows the tip angles for the two simulated tissue types: one hard and one soft. The soft tissue is modeled after the soft plastisol gel described in Section IV-A and the hard tissue has a stiffness and damping coefficient 100 times larger than the modeled plastisol, which is much stiffer than any tissue a steerable needle would be inserted through.

Fig. 8
Finite element simulation

The needle tip rotates slightly differently in the soft and hard tissues. The needle in soft tissue initially rotates faster than the hard tissue because the soft tissue moves slightly with the needle, thus providing less impedance. The simulated soft tissue rotates to within 3° in 1.05 sec and the hard tissue in 1.22 sec, a difference of 15%, but the needle in soft tissue does not reach the desired angle for several seconds.

The soft tissue starts to creep back toward its original location as the needle approaches the desired angle. As the tissue returns, it pulls the needle slightly away from the base angle, causing the needle to briefly maintain a small offset. The tissue reaches its resting state around 2.7 sec and the needle then reaches the desired angle. For multiple rotations in alternating directions, the restoring force of the tissue applies a torque in the direction of the second rotation and briefly aids in needle rotation, much like the multiple rotation experiments described in Section II-B2 and in [10]. Overshooting the desired angle can alleviate much of tissue effect since the final motion causes the tissue to return to its resting state as the needle reaches the desired angle.

Another metric to determine the impact of the angle lag is the integral of the tip error during a rotation. This integral indicates how much the needle diverges from the desired path. The integral of the error in the soft tissue is 43.8 sec° and the hard tissue is 48.6 sec°, a difference of 10%.

Although tissue elasticity causes a 10–15 % difference in rotation between a soft and hard tissue, it is unlikely to cause a significant deviation in human tissue. One of the target regions for needle steering is the prostate, which Krouskop et al. [24] shows has a Young’s Modulus around 60 kPa, larger than the tissue used in our simulations. The higher stiffness in in vivo tissue will reduce the effects of tissue elasticity when a needle is rotated. Since the controller we discuss next naturally develops a small overshoot and the tissues are likely to be stiffer than those tested here, we do not foresee tissue elasticity causing a significant angle difference. Any remaining effects from soft tissue can be eliminated by assuming a slightly more flexible needle to account for tissue elasticity. We will assume no tissue effects and will use the continuous model for the remainder of this paper.

V. Control of the Needle Tip

The fourth-order inertialess modal model without tissue effects described in Section III-C allows better control of the needle tip than using the base angle measurement alone. Since there is no feedback of the needle tip angle, we will use a state estimator to predict the current states (mode shapes) of the needle.

For comparison, we develop two estimators: (1) feedforward and (2) feedback from the tip angle sensor described in Section II-A. The accompanying controller is designed to control the modes of the system, as opposed to the typical method of controlling positions or velocities of discrete points. The block diagram for the estimator-controller is shown in Figure 9. The dotted line indicates the tip angle feedback that is not available during a needle insertion. In both observers, the zeroth mode is the needle angle directly measured at the needle base and the higher modes are estimated based on the derived dynamics. The control input is specified based on the the estimated position of the needle tip, which is the sum of each mode evaluated at the needle tip. Without feedback, the non-zero-mode system dynamics are estimated as

Fig. 9
Torsion controller

x^.=Ax^+Bu
(31a)

y^=Cx^
(31b)

where A, B, C, and x are defined in (23).

In the comparative case with tip feedback, we used a Luenberger observer:

x^.=Ax^+Bu+L(yy^)
(32)

where the state matrices are the same as (31) and the feedback gains (L) were chosen so that the eigenvalues of ALC are five times faster than the system (A) eigenvalues.

We used the same needle described in Section II-A with the damping values determined from Table II. Figure 10 shows the needle base and average tip trajectories resulting from three trials of a rotation with and without torsion compensation at two depths. The torsion compensator is able to rotate the needle tip to within 3° in nearly half the time of only rotating the base to the desired angle. The needle tip convergence times are shown in Table IV. The benefit of compensating for torsion increases at deeper insertions, since the total friction, and thus the tip lag, increases with distance inserted. The overshoot in the base angle alleviates much of the tissue effects since the final motion is in the same direction as the tissue relaxation.

Fig. 10
Estimator-controller with feedback and feedforward
TABLE IV
Needle Tip Settling Time

The feedback estimator is faster than the feedforward estimator. Although no tip feedback is available for the rotation angle during a needle insertion, comparing the two estimators shows that the feedforward mechanics-based model is capable of correcting a large portion of the tip angle error.

The integral of the tip angle error during the rotations indicates how slow needle tip convergence affects the path. Table V shows the integral of the tip angle error over two seconds for the needle with and without torsion compensation. Note that Figure 10 shows the absolute value of the tip error, but this integral incorporates negative angles because a negative angle will counteract the out-of-plane motion caused by a positive angle. Adding torsion compensation decreases the the integral error by about 40%. A control scheme could be to drive this error to zero, ensuring minimal motion out of the plane.

TABLE V
Integral of Tip Angle Error

The faster convergence time and decreased tip error with torsion compensation allow more accurate control of the needle position in the tissue. We tested the feedforward torsion compensation by inserting a needle 7 cm at 0.5 cm/sec. We used a pre-bent needle with a 45° bevel and a 15° pre-bend 12 mm from the tip in a complimentary direction. Reed et al. [6] experimentally determined the radius of curvature to be 6.1 cm. Figure 11 shows the final positions from the camera for each case. The baseline case consists of an insertion without the needle rotating. The baseline was run twice: (1) with the needle curved to the right and (2) with the needle rotated 180° before insertion so the needle curved to the left. In each of these baseline trials the needle did not drastically deviate from the starting plane. In the second case, the needle base rotated 180° after each of the first four centimeters without torsion compensation. After a 7 cm insertion, the needle had deviated by 12 mm out of the plane. The same four rotations with torsion compensation only caused a 4 mm deviation over a 7 cm insertion. Figure 12 shows the height of the needle for each case. Each point represents the triangulated and filtered deviation from the desired plane recorded at 7.5 Hz. All cases started at the same point, but the cameras were unable to accurately track the tip during the first cm due to edge occlusions in the tissue.

Fig. 11
Insertions with and without torsion compensation
Fig. 12
Needle height deviation

Although the torsion compensation does not take into account the added torques caused by the pre-bent tip, the controller was able to reduce the height deviation by nearly 70%. A bevel-tip needle has less curvature, so the deviations from a plane would be smaller, but would show a similar proportion of deviation with and without torsion compensation. The rotations in this example could have alternated directions between the tip pointing up and the tip pointing down. Alternating the rotation direction would reduce the height deviation, but this example demonstrates the reduced deviation from the plane using torsion compensation.

VI. Conclusions

Friction is detrimental to accurate control of steerable needles, both in the phantom tissues tested here and likely in human tissues, particularly when using superelastic needles to enhance steering. Torsional friction causes a discrepancy of over 45° between the base and tip angles in certain phantom materials. Based on previous experiments with steel needles during a prostate brachytherapy [5], we estimate that torsion will cause a 10–15° discrepancy in human tissues. It is possible that torsional forces could be amplified by lateral forces on the needle during insertion into human tissue. Such significant errors will likely imperil image-guided controllers and path planners designed for flexible needle steering. Unfortunately, there is a trade off between the flexibility of needles for steering and the stiffness of needles to rotate as a rigid body inside tissue. Our model and controller for the torsional dynamics can alleviate a large portion of the angle lag due to torsional friction, thus allowing the use of more flexible needles with enhanced steerability.

We developed an estimator-controller based on a modal model of the needle dynamics to predict the tip angle for a needle rotating inside tissue. The controller quickly drives the needle tip to the desired angle. The faster convergence allows the needle to maintain motion in a prescribed plane significantly better than without torsion compensation. Our torsion model is designed to work in conjunction with the planar controller derived by Kallem and Cowan in [8]. They demonstrate that the real system takes twice as long to converge than simulations, which they speculate is due to unmodeled torsion dynamics. Incorporating the torsion dynamics developed here is likely to significantly increase the convergence speed. Future experiments using in vivo tissues with multiple tissue types and membrane layers under medical imaging will be necessary for complete validation.

If a steerable needle or catheter is controlled by a robot or a clinician via teleoperation, estimating the torque may be important. Clinicians steering a needle will likely occasionally stop and rotate before further insertion. When they subsequently insert further, the inaccuracies following a rotation can cause path deviations. It may be necessary to automatically compensate or provide torque feedback to the doctor. One feedback mode might be an amplified version of the torsional stiffness. This would be particularly important for pre-bent needles that are not being inserted at that instant.

Acknowledgments

This work is supported in part by the National Institutes of Health under Grant R01-EB006435. Portions of this work were previously presented at the IEEE/RAS-EMBS International Conference on Biomedical Robotics and Biomechatronics, 2008 and the IEEE International Conference on Robotics and Automation, 2009.

The authors thank Sarthak Misra, John Swensen, and Tom Wedlick for help with this research.

Biographies

An external file that holds a picture, illustration, etc.
Object name is nihms192959b1.gif

Kyle B. Reed received the B.S. degree from the University of Tennessee, Knoxville, TN in 2001, and M.S. and Ph.D. degrees from Northwestern University, Evanston, IL, in 2004 and 2007, respectively, all in Mechanical Engineering. From 2007 to 2009, he was a Postdoctoral Fellow in the Laboratory for Computational Sensing and Robotics at The Johns Hopkins University, Baltimore, MD. In 2009, he joined the faculty of the University of South Florida, Tampa, FL, as an Assistant Professor. He was the recipient of the 2001 NSF Graduate Fellowship. Dr. Reed’s interests include haptics, human-machine interaction, rehabilitation engineering, medical robotics, and engineering education.

An external file that holds a picture, illustration, etc.
Object name is nihms192959b2.gif

Allison M. Okamura received the B.S. degree from the University of California at Berkeley, and the M.S. and Ph.D. degrees from Stanford University, in 1994, 1996, and 2000, respectively, all in mechanical engineering. She is a Professor of Mechanical Engineering and the Decker Faculty Scholar at The Johns Hopkins University, Baltimore, MD. Her research interests include haptics, teleoperation, robot-assisted surgery, tissue modeling and simulation, rehabilitation robotics, and prosthetics. Dr. Okamura received the 2004 National Science Foundation CAREER Award, the 2005 IEEE Robotics and Automation Society Early Academic Career Award, and the 2009 IEEE Technical Committee on Haptics Early Career Award. She is an associate editor of the IEEE Transactions on Haptics.

An external file that holds a picture, illustration, etc.
Object name is nihms192959b3.gif

Noah J. Cowan received the B.S. degree in electrical engineering from the Ohio State University, Columbus, in 1995, and the M.S. and Ph.D. degrees in electrical engineering and computer science from the University of Michigan, Ann Arbor, in 1997 and 2001, respectively.

From 2001 to 2003, he was a Postdoctoral Fellow at the University of California, Berkeley. In 2003, he joined the faculty of The Johns Hopkins University, Baltimore, MD, as an Assistant Professor. Prof. Cowan’s research interests include multisensory control in animals and machines.

Prof. Cowan received the NSF CAREER award in 2009, the William H. Huggins Award for excellence in teaching in 2004, and a Rackham Doctoral Fellowship from the University of Michigan in 2000.

Footnotes

Personal use of this material is permitted. However, permission to use this material for any other purposes must be obtained from the IEEE by sending an email to gro.eeei@snoissimrep-sbup.

1Of course, one does not require that these be mode shapes — any convenient set of orthogonal basis functions will work; we drop the sine terms for reasons that will be clear shortly.

References

1. Nath S, Chen Z, Yue N, Trumpore S, Peschel R. Dosimetric effects of needle divergence in prostate seed implant using I and Pd radioactive seeds. Med Physics. 2000 May;27:1058–1066. [PubMed]
2. Youk J, Kim E, Kim M, Lee J, Oh K. Missed breast cancers at US-guided core needle biopsy: How to reduce them. Radiographics. 2007 Jan–Feb;27:79–94. [PubMed]
3. Webster R, III, Kim JS, Cowan NJ, Chirikjian GS, Okamura AM. Nonholonomic modeling of needle steering. Int J Robot Res. 2006;25(5–6):509–525.
4. Abolhassani N, Patel R, Moallem M. Needle insertion into soft tissue: A survey. Med Engr and Physics. 2007;29:413–431. [PubMed]
5. Podder T, Clark D, Sherman J, Fuller D, Messing E, Rubens D, Strang J, Liao L, Ng W, Yu Y. In vivo motion and force measurement of surgical needle intervention during prostate brachytherapy. Med Physics. 2006;33(8):2915–2922. [PubMed]
6. Reed KB, Kallem V, Alterovitz R, Goldberg K, Okamura AM, Cowan NJ. Integrated planning and image-guided control for planar needle-steering. Proc IEEE Conf Biomed Robotics. 2008 October;:819–824. [PMC free article] [PubMed]
7. Ding M, Cardinal H. Automatic needle segmentation in three-dimensional ultrasound image using two orthogonal two-dimensional image projections. Med Physics. 2003;30(2):222–234. [PubMed]
8. Kallem V, Cowan NJ. Image guidance of flexible tip-steerable needles. IEEE Trans Robotics. 2009;25:191–196. [PMC free article] [PubMed]
9. Abolhassani N, Patel R, Ayazi F. Minimization of needle deflection in robot-assisted percutaneous therapy. Int Journal of Medical Robotics and Computer Assisted Surgery. 2007;3:140–148. [PubMed]
10. Reed KB. Compensating for torsion windup in steerable needles. Proc. IEEE Conf. Biomed. Robotics; 2008. pp. 936–941. [PMC free article] [PubMed]
11. Alterovitz R, Lim A, Goldberg K, Chirikjian GS, Okamura AM. Steering flexible needles under markov motion uncertainty. IEEE/RSJ Int. Conf. on Intelligent Robots and Systems (IROS); 2005. pp. 120–125.
12. Webster R, III, Memisevic J, Okamura AM. Design considerations for robotic needle steering. IEEE Int. Conf. on Robotics and Automation (ICRA); 2005. pp. 3588–3594.
13. Kallem V, Chang D, Cowan NJ. Task-induced symmetry and reduction in kinematic systems with application to needle steering. IEEE/RSJ Int. Conf. Intelligent Robots and Systems (IROS); 2007. pp. 3302–3308. [PMC free article] [PubMed]
14. Kataoka H, Washio T, Chinzei K, Mizuhara K, Simone C, Okamura AM. Measurement of the tip and friction force acting on a needle during penetration. Proc. Medical Image Computing and Computer Assisted Intervention (MICCAI); 2002. pp. 216–223.
15. Okamura AM, Simone C, O’Leary M. Force modeling for needle insertion into soft tissue. IEEE Trans Biomed Engr. 2004;110(51):1707–1716. [PubMed]
16. Thompson JMT, Champneys AR. From helix to localized writhing in the torsional post-buckling of elastic rods. Proc: Mathematical, Physical and Engineering Sciences. 1996;452:117–138.
17. Olsson H, Åström K, Canudas De Wit C, Gäfvert M, Lichinsky P. Friction models and friction compensation. European Journal of Control. 1998;4(3):176–195.
18. Inman D. Vibration with Control. West Sussex, England: John Wiley & Sons, Ltd; 2006.
19. Caughey T, O’Kelly M. Classical normal modes in damped linear dynamic systems. J Appl Mech. 1965;32:583–588.
20. Salehi-Khojin A, Bashash S, Jalili N. Modeling and experimental vibration analysis of nanomechanical cantilever active probes. J Micromech Microeng. 2008;18
21. Bashash S, Salehi-Khojin A, Jalili N. Forced vibration analysis of flexible Euler-Bernoulli beams with geometrical discontinuities. American Control Conference; June 2008; pp. 4029–4034.
22. Misra S, Ramesh KT, Okamura AM. Modeling of tool-tissue interactions for computer-based surgical simulation: A literature review. Presence: Teleoperators and Virtual Environ. 2008;17:463–491. [PMC free article] [PubMed]
23. Misra S, Reed KB, Schafer B, Ramesh KT, Okamura AM. Mechanics of robotically steered flexible needles through soft tissue. Int J Robot Res. submitted. [PMC free article] [PubMed]
24. Krouskop T, Wheeler T, Kallel F, Garra B, Hall T. Elastic moduli of breast and prostate tisues under compression. Ultrasonic Imaging. 1998;20:260–274. [PubMed]