Calculating Curved Needles Deformation in Automated High- Dose Rate Brachytherapy Operations

Corresponding Author: Galina S. Kireeva Russian State Scientific Center for Robotics and Technical Cybernetics, St. Petersburg, Russia Email: galinakireyeva@mail.ru Abstract: The article outlines the approach to applying the mechatronic technologies and robotic systems in carrying out the brachytherapy operations. To calculate the elastic state of the flexible needles of different geometry in the prostate model, a mathematical model of the flexible needle deformation has been applied which took into account the medium reaction and was based on the theory of flexible and elastic rods. The numeric calculations for different types of needles have been introduced for the first time. Also for the first time the advantages of applying flexible needles for brachytherapy procedure have been mathematically justified.


Introduction
Over a span of many decades the basic and, as it was, the only one available method of treatment for those suffering from prostate cancer used to be the surgical approach, implicating complex traumatizing operations often resulting in various side-effects: In 19-45% of cases leading to serious urination problems, in 60-93% of cases resulting in evident sexual dysfunction (D'Amico et al., 1998;Kupelian et al., 1997;Oesterling et al., 1995;Gleason and VACURG, 1977;Lee et al., 1995). The range of the available therapeutic solutions has been diversified considerably up to now (Kirschner et al., 2014;Skowronek, 2013;Sylvester et al., 1997;Stock et al., 1998;Tejwani et al., 2012). Notably, for the patients with localized-type prostate cancer the radiation therapy is considered to be either a good alternative or an important supplement to the surgical treatment. The task of delivering high tumoricidal doses into the prostate and to the tumor, which is located inside it, with considerably less radiation load on the surrounding tissues, is now successfully solved by means of modern methods of interstitial brachytherapy, a version of which is represented by High-Dose Brachytherapy (HDB) (Fig. 1) (Kanayev and Novikov, 2013;Khmelyovskiy et al., 2008;Bracarda et al., 2005;Jereczek-Fossa and Orecchia, 2007;Kovacs et al., 2005;Martinez et al., 2011). An important advantage of this method is that when performing high-dose brachytherapy there are favorable conditions for the intrastat needles introduction beyond the prostate capsule, which facilitates radiating the seminal vesicle and the periprostatic area. However, a serious disadvantage to the existing technique of implanting the radioactive seeds is the direct trajectory of the needle targeting, which is done manually.

Fig. 2.
A picture of the prostate taken after implanting the needles superimposed on the initial radiation plan: aprostate contour; b-urethra contour; c-rectum wall contour; d-urinary bladder contour; e-isodose curve (100% of the prescribed dose); f-intrastat needle contour Due to this fact, in the course of an operation a considerable alteration of the radiation plan occurs as compared to the pre-operation plans (Fig. 2), resulting in an inadequate radiation dose distribution on the tumor and, as a consequence, in the poor results of treatment in general. Mathematical modeling of the physical principles constituting the basis of the medical soft/hardware complexes is a modern approach to solving complicated procedural tasks, including surgery. This study indicates the approach to mathematical modeling of the flexible needle behavior in the model of the prostate.

Materials and Methods
To describe the elastic state of a medical needle, the thin flexible rod curvature equations shall be applied (Birger and Panovko, 1988;Turner and Ford, 1960;Salzmann, 1946;Marçal and Turner, 1961).
Consider a flexible curved rod in xOy plane, shown in Fig. 3.
The axis of the rod is given in parametric representation: The formula for calculating Lame coefficient (Birger and Panovko, 1988;Turner and Ford, 1960;Salzmann, 1946;Marçal and Turner, 1961) is as follows: Rod axis curvature shall be computed using the following formula: Rod equilibrium equation should be represented as follows: Here T, N are the tangential and the shear strains; M is the moment of deflection; p T (α)и p N (α) are the tangential and the normal components of the distributed load; A(α)-Lame parameter; R (α) = 1/k(α)-rod curvature radius. Positive directions of forces and moments are shown in Fig. 4.
For the rod center line deformation the following formulae are applied: Here, ε and κ are tensile deformation and flexural deformation accordingly; θ is the cross section flexion angle; u and w are the tangential and the normal displacement components.
The connection between the forces in the rod section (4) and the deformations (5) is set forth by the following correlations: , where, B and G are tensile stiffness and flexural stiffness. For continuous annular section rods now we have the following: Here, D is the section diameter; E is Young's modulus.
For a annular section rod the following will take place: Here, d is the inner diameter of the section. The specific character of the task under consideration stipulates some extra preconditions for the distributed load, which is the response of the medium to the impact of the needle. Obviously, the tangential load p T from the part of the medium on the lateral surface of the needle is a negligibly small value for calculating the strained state, as compared to the normal component. Therefore, further we assume the following: For the normal component of load at the stretch of the needle 0< a < α 0 , which does not get in contact with the medium, we also have p N (α) ≡ 0.
At α 0 < a < L, that is at the stretch of the needle, which gets into contact with the medium, the normal component of load is stipulated by the response of the medium to the displacement component. As a hypothesis, we assume that the response of the medium is in direct proportion to the value of this displacement. The proportionality factor k is a characteristic of the medium. Thus, we have the following: For modeling the behavior of the medical needle in a patient's body with the help of the Equation 4-7 we have to supplement these with the boundary conditions corresponding to the specific conditions of this task.
The left end of the rod (α = 0) is assumed to be rigidly fixed: At the right end of the rod (α = L) the boundary conditions of forces will be set as follows: The equation system (4-7) is of the sixth order. Together with the boundary conditions (8-9) it makes up a boundary value problem, by solving which we shall find the shape of the deformed needle.

Straight Needle
With a straight needle, which axis is located along x axis, the Equation 4 become more simple: α = x, R = ∞, A = 1. The equilibrium equations will be as follows: The formula for deformations will be put in the following way: Thus, the equation of strain and the equation of flexion can be solved independently and the corresponding systems are of the second and of the forth order accordingly. Of course, the correlations (10-11) will coincide with the well-recognized equations of strain and flexion of the straight rods (Birger and Panovko, 1988;Turner and Ford, 1960;Salzmann, 1946;Marçal and Turner, 1961).

Annular Needle
In case of a annular needle with radius R, a point on the center plane will be identified using the angular coordinate ϕ. Then A = R and the equilibrium equations will be as follows: For deformations we have the following formula: Thus, the annular needle equations are linear differential equations with constant factors. The solution for the later system is expressible in an essentially closed form in the same way as for the straight rod.

Elliptical Needle
The parametric equations of ellipse are as follows (0≤ α≤2π): Here, a and b are the long and the short half-axes of ellipse. Then, according to the formulae (2-3), we have the following: Later equations, when put into the formulae (4-5), will result in the linear equation system with constant factors and this system will have to be solved by methods of approximation.

Boundary Value Problem
Let us write down the boundary value problem of the curved needle deformation in the medium, which models a patient's body. For this purpose we shall enable the correlations (4-5) pertaining the derivatives, use the correlations (6) and take into account the specific features of the distributed load, characteristic for the task, using (7). As parameter α we shall choose the curve length s of the center line of the needle. Then A ≡ 1. This will result in the following: Let us supplement the system (12) with the boundary conditions: The system (12-13) can be solved by the sweep method, i.e., the solution to the boundary value problem shall be found as a linear combination of the solutions for three auxiliary initial tasks, satisfying the conditions at the left end of the interval of integration (13a).
To make the numeric solution more convenient, we shall write down the system (12) as a matrix: Where the following indications are introduced: Now, U 1 (s), U 2 (s), U 3 (s) are three solutions to the system (14), meeting the initial conditions: The boundary value problem will be solved as follows: Obviously, the solution (15) meets the conditions (13a) at the left end of the interval of integration.
To meet the boundary conditions at the right end of the interval of integration the following correlation has to be valid: To solve the later system, we shall determine the values of the variables c i in the Equation 15. The values of displacement and the angular deflection at the right end of the interval of integration will be expressed by the equations as follows:

Straight Needle
Consider the straight needle with the length L. We shall study the deformation of such needle, which takes place when the needle is pointed in the medium modeling a patient's body. At the first stage, we shall consider the problem occurring when the straight rod loses stability at meeting a rigid obstacle. The equilibrium Equation 4 have to be related to the deformed axis of the rod and they will, as a result, contain the terms, accounting for the curvatures, thus resulting in a non-linear problem. Neglecting the details, we shall assume that the specified computations have been carried out and that the linearization of the resulting equations has been implemented. The boundary conditions will be taken from the form (8-9). Let us find the critical force T L . The first equilibrium Equation 10 shows (provided that no distributed tangential load is present) a continuous axial force in the rod: T(x) = T L . The second Equation 4 will be as follows: Given the third equilibrium equation, we shall get the following: The later equation coincides with the equation for determining the formula of the rod stability loss (Birger and Panovko, 1988), but it takes into account the supporting force. Characteristic equation has two realvalued and two pure imaginary solutions: The solution to the Equation 16, meeting the conditions (13a) at the left end, will be put as follows:  Therefore:   1  2  2  1  2  1  2  2  3  3  3  2  2  1  1  2  2  1  1   sin  cos  0  cos sin Equation 18 together with the correlations (17) make it possible to find the value of the critical force T L , at which the needle loses stability. Consider an important special situation when there is no medium response, i.e., when the needle perforates a surface. Thus,

Annular Needle
The computation of the tension state for the annular needle has been carried out for the following input data: Material of the needle is steel, E = 2.0⋅1011 N/m 2 ; radius of the needle, R = 10 cm; outer diameter of the needle, D = 1 mm; inner diameter of the needle, d = 0.8 mm; value of the force at the tip of the needle, T L = 0.01 N; medium response, k = 10 N/m 2 . Length of the needle L = 78.5 mm, which corresponds to 1/8 of the circumference: 0≤ϕ≤π/4. The results of the computation are shown in Fig. 6.
The maximum value of a tangential displacement is reached at the tip of the needle and it amounts to app. 0.01 mm. The shearing force and the moment of deflection are turned to zero at the tip of the needle. Such boundary conditions secure the minimum traumatic impact on the tissues.
The second graph shows that the maximum value of normal displacement assumes the value of app. 0.07 mm.

Elliptical Needle
The computation of the tension state for the elliptical needle with wide eccentricity have been carried out for the input data as follows: Material of the needle is steel, E = 2.0⋅10 11 N/m; the short and the long half-axis of the ellipse, a = 1 mm and d = 100 mm; outer diameter of the needle, D = 1 mm; the inner diameter of the needle, d = 0.8 mm; value of force at the tip of the needle, T L = 0.01 N; medium response, k = 10 N/m 2 . 0≤ϕ≤π/2.
The results of the calculations are shown in Fig. 7. It could be observed, that close to the right end of the needle the displacement and the force graphs have the appearance of a characteristic boundary effect, at the same time the angle of the section deflection and the moment of deflection are altered slowly over the length of the needle.

Discussion
As of today, at the stage of theoretical and practical research, a scientific and technical basis is formed in the area of developing the flexible needles and the methods for controlling their movement in the robotic complexes for executing brachytherapy operations. With mathematical modeling and under the phantom model experimental conditions the achieved precision of needle implanting is quite high with minimum mistakes at targeting. However, under the clinical research conditions the loss of targeting precision should be expected, as there appear other influencing factors, such as physiological movement of the tissue, flows of biological liquids and the tissue heterogeneity. The functions of the control system can be enhanced up to detecting the patient's movements during the needle targeting process (breath, flows of biological liquids).

Conclusion and Further Researches
The numerical calculations can be used for solving the problem of minimum tissue traumatizing at implanting the radioactive seeds, as well as for calculating the controls of the needle operating manipulator.
It has been found, that in the practically important case with the needle of annular shape inserted along the annular trajectory, the shearing force and the moment of deflection are turned to zero at the tip of the needle, thus ensuring the minimum traumatizing effect on the tissues and, as a consequence, improving the safety of the procedure. With minor deviations this has been also confirmed for the elliptical needle.
In the course of performing the calculations for the straight needle a dependency of the critical force on the length of the needle has been established, which, provided that the input data on the tissue are available, makes it possible to determine the value of this force to ensure the perforation.
This article considered the issue of calculating the curved needle deformation as a problem of interaction between the tissue and the controlled needle. In future, additional calculations are planned to be performed to formulate the curved trajectory lines with different number of inflections. This in turn will bring up the question of extra investigations to be performed on traumatic effects on the tissue.