Mathematical and computer modeling of endovenous laser treatment

Abstract

An initial-boundary value problem for quasilinear equations of radiative-conductive heat transfer, simulating the process of endovenous laser ablation, is considered. The unique solvability of the problem is proved and an algorithm for finding its solution is proposed. The eficiency of the algorithm is illustrated by numerical experiments.

Keywords

• eol>numerical simulation
• endovenous laser ablation
• quasilinear equations of heat transfer

1. Introduction

The endovenous laser ablation (EVLA) is an efective and popular treatment of varicose veins. During EVLA a laser optical fiber is introduced into a damaged vein. Afterwards, the optical fiber transmitting laser radiation is drawn out of the vein. The fiber tip emitting radiation is usually coated with a carbonized layer, which splits the laser energy, such that a part of the energy is absorbed by the carbonized layer with the release of heat, the other part is spent on radiation. The heat from the carbonized layer is transmitted into blood by the conductive heat transfer. Heat transfer is significantly intensified by the flow of bubbles formed on the hot fiber tip. The radiation entering into the blood and the surrounding tissue is partially absorbed with the release of heat. As a result, the heat energy generated by diferent mechanisms causes essential heating the vein, that leads to its obliteration (closure).

Optimization of radiation parameters during EVLA ensures successful vein obliteration with minimum frequency and severity of complications. In other words, optimal radiation must provide a rather high temperature inside the vein for its obliteration and simultaneously the generated temperature ifeld must be relatively safe for the tissue surrounding the vein. The principal parameters influencing the eficiency and safety of the laser ablation procedure are the laser power, radiation wavelength, pullback speed of the optical fiber, and the ratio between the laser powers spent on radiation and ifber tip (carbonized layer) heating. As a rule, the laser ablation is performed by a radiation with a wavelength from 810 to 1950 nm. The suficiently widely used ranges of the pullback speed of the ifber and laser radiation power are 1–3 mm/s and 5–15 W, respectively.

The classic model describing the interaction between optical radiation and a biological tissue is a reaction-difusion model including the radiation transfer and heat conduction equations [1, 2]. Correspondingly, the calculation of radiation and temperature fields involves finding the bulk density of absorbed radiation energy from the radiation transfer equation and subsequent solving the heat conductivity equation, in which the absorbed radiation energy form heat sources. The theoretical and numerical analysis of reaction-difusion models describing various phenomena in physics, biology, and ecology can be found in [3, 4, 5, 6]. The reaction-difusion models of heat transfer accounting the radiation efects are studied in [7, 8]. In [9, 10, 11, 12], to simulate EVLA, the authors apply the reaction-difusion heat transfer equation, where the bulk density of the absorbed energy of radiation is calculated using an explicit formula for intensity of radiation induced by a point source in infinite homogeneous medium.

The principal efects, which are usually taken into account in EVLA modeling, are the conductive heat transfer, radiation transfer and absorption with the release of heat, and heat transfer by the flow of bubbles formed at the hot fiber tip. In the works [10, 11, 12], on the base of an estimation of experimental data, heat transfer by the flow of bubbles is modelled by the use of a piecewise constant heat conductivity coeficient depending on the temperature as follows: when the temperature at a certain point reaches 95◦C and higher, the heat conductivity coeficient increases by 200 times. We will use this approach in the model of the reaction-difusion describing process of EVLA. The corresponding quasilinear initial-boundary value problem of EVLA is studied in the current work. The unique solvability of the problem is proved and an algorithm for finding its solution is proposed. The eficiency of the algorithm is illustrated by numerical examples. Theoretical and numerical analysis of the considered model are also conducted in [13, 14], where extremal problems for the quasilinear model of EVLA are studied. 2. Mathematical model of endovenous laser ablation Let us formulate an initial-boundary value problem for the following equations modeling process of EVLA in a bounded three-dimensional domain Ω with the boundary Γ = Ω: / − div( ( )∇ ) − = 1 ,

−div( ∇ ) + = 2 , ∈ Ω, 0 < < , ( ) + ( − ) = 0, + 0.5 = 0 on Γ, | =0 = 0. ( 1 ) ( 2 ) ( 3 ) Here, is the temperature, is the radiation intensity averaged over all directions, is the difusion coeficient for optical radiation, is the absorption coeficient, ( ) is the thermal conductivity coeficient, (, ) is the product of the specific heat capacity and the density of the medium, is the heat transfer coeficient, 1 describes the source power used on heating the fiber tip, and 2 is the source power going to the radiation, is the characteristic function of the fiber tip area divided by the volume of the fiber tip. The functions , 0 define the boundary and initial temperature distributions, respectively. By we denote the derivative in the direction of outward normal to the boundary Γ. ( ) 0 < 0 ≤ ≤ 1, | / | ≤ 2, = . ( ) 0 < 0 ≤ ( ) ≤ 1, | ′( )| ≤ 2, ∈ ℝ, = .

0, ∈ ∞(Σ).

We define a nonlinear operator ∶ → ′, using the equality valid for any , ∈ : ( ( ), ) = ( ( )∇, ∇ ) + ∫ Γ = (∇ℎ( ), ∇ ) + ∫

Γ. Here, ℎ( ) = ∫ ( ) .

Let ∈ ∞(0, ; ′), (, ) = ∫ Γ.

In what follows, as the inner product in , we will use the bilinear form (, ) = (∇, ∇ )+∫

The following weak formulation is valid for problem (4)–(5).

The corresponding norm is equivalent to the standard norm of the space .

Definition. Function ∈ 2(0, ; ) is a weak solution of the problem ( 4 )–( 5 ) if ′ ∈ 2(0, ; ′) 3. Solvability of the initial-boundary value problem 3.1. Quasilinear equation of heat transfer

Let us consider the problem:

In what follows, we assume that Ω is a Lipschitz bounded domain, Γ = Ω = Ω × (0, ), Σ = Γ × (0, ). , By , 1 ≤ ≤ ∞

, we denote Lebesgue space, through 1 the Sobolev space the Lebesgue space of functions of the class , defined on (0, ), with values in the Banach space . = 2(Ω), = 1(Ω), by ′ we denote the dual space of . The space is identified with the 21, and by (0, ; ) space ′, so that ⊂ of the functional = ′ ⊂ ′. We denote by ‖ ⋅ ‖ the standard norm in , and by (, ) the value ∈ ′ at the element ∈ , which coincides with the inner product of if ∈ . ( 4 ) ( 5 ) Hereinafter, ′ = / , = + .

Note that ( )′ = ′ +

/ Therefore, the initial condition (0) = 0 ∈ 2(0, ; ′). Then

makes sense.

′ + ( ) = , (0) = 0.

Theorem 1. Let conditions (i)–(iii) hold. Then the problem ( 4 )–( 5 ) is solvable.

Proof. Let us define the Galerkin approximations of a solution of the problem ( 4 )–( 5 ) and of the following Cauchy problem: derive a priori estimates which are necessary for proving the solvability. In the space , we consider an orthonormal in basis 1, 2, … , = span{ 1, ..., }. Let ( ) ∈ , ∈ (0, ), be a solution

( ′ + ( ) − , ) = 0 ∀ ∈ , | =0 = 0 .

Here, 0 is the orthogonal projection in of the function 0 onto the subspace .

∈ 2(0, ; ′) and hence ∈ ([0, ]; ′).

Let = min{ 0, 0} > 0. Since and integrate with respect to time from 0 to . Then

Let us derive a priori estimates which are necessary for proving the solvability. Put = in ( 7 ) ‖ ‖2 + ∫ ( 0 ( ( ))∇ ( ), ∇ ( )) + ∫ 2 Γ =

) 0 1 2 then, accounting for the conditions ( ), ( ), ( ), and the definition of the norm in , we derive the 0‖ ( )‖2 + ∫ ‖ ( )‖2 ≤ + 2 ∫ ‖ ( )‖2.

Here, = 1‖ 0‖2 + −1‖ ‖2 ′. By the Gronwall’s inequality, for ‖ ( )‖2, we obtain the estimate The resulting inequality allows one to estimate the right-hand side in ( 8 ) and therefore 0 0 exp 2 a.e. on (0, ).

0 2 ∫ ‖ ( )‖ ≤ exp

Here and below, in the proof of the theorem, by we denote constants independent of . The obtained estimates ( 11 ) allow us to assert that passing, if necessary, to a subsequence, there is a function such → weakly in 2(0, ; ), ∗ –weakly in ∞(0, ; ), ℎ( ) → weakly in 2(0, ; ). ( 12 ) that the limit function ∈ 2(0, ; ) is such that ′ ∈ 2(0, ; ′), and the equality Convergence results ( 12 ) are suficient for passing to the limit as → ∞ in the system ( 7 ) and proving ( ′, ) + (∇ , ∇ ) + ∫

Γ = (, ) ∀ ∈ and the initial condition hold. To prove that is a weak solution of the problem ( 4 )–( 5 ), it is suficient to check that the equality = ℎ( ) holds. We obtain an estimate guaranteeing the compactness of the where Accounting for the non-negativity of (∇ℎ( ( )), ∇ ( )) + ∫ 2 ( ) Γ, we obtain the inequality (, ) = (∇ℎ( ( )), ∇( ( ) − ( ))) + ∫ ( )( ( ) − ( )) Γ− ( ( ) + ( ( ) − ( )), ( ) − ( )).

1 2 sequence in 2( ). In the system ( 7 ), we put = ( ) − ( ) and integrate with respect to on the interval (, + ) and with respect to on (0, − ), assuming that > 0 is small enough. Then 1

√

), 2 ‖ ( )‖2 ≤ 0 exp 2 a.e. on (0, ), ∫ ‖ ( )‖2 ≤ exp (14) ∞( ), then there are no other solutions of the problem.

Let us show that the solution is unique in the class of functions with bounded gradient.

Theorem 2. Let conditions (i)–(i ) hold. If ∗ is a weak solution to problem ( 4 )–( 5 ) such that ∇ ∗ ∈ (, ) ≤ 1|( ( ), ( )) | + ( ( ), ( ) − ( )) + 2 2‖ ( ) − ( ))‖2. we have Here, 1 = max{ 1, ‖ ‖ ∞(Γ)}. Therefore, taking into account the continuity of the embedding ⊂

To estimate the integrals of the terms depending on , it is enough to change the order of integration. Using the boundedness of the sequence in 2(0, ; ), we obtain the estimate of equicontinuity: ∫ ‖ ( + ) − ( )‖ ≤ . (13) obtained in the proof of Theorem 1.

Indeed, by virtue of equality

2( ) and hence = ℎ( ). The theorem is proved.

From the estimate (13) (passing, if necessary, to a subsequence), it follows that

Therefore, by virtue of the inequality |ℎ( ) − ℎ( )| ≤ 1| − |, we conclude that ℎ( ) → ℎ( ) in → in 2( ).

Remark. If is an arbitrary weak solution of the problem ( 4 )–( 5 ), then it satisfies the estimates ( , ) = ( ′, ) + 2 ( , ), 1 of Theorem 1, we derive where = / , the function → ( ,

)/ is integrable on (0, ). Multiplying the first equation in ( 6 ), in the sense of the inner product of , by and integrating with respect to time, as in the proof Proof. Let 1 be another solution of the problem ( 4 )–( 5 ), = 1 − ∗. Then

′ + ( 1) − ( ∗) = 0, (0) = 0. to time. Then We multiply the first equation, in the sense of the inner product of , by and integrate with respect Using the constraints on the functions ,

and their derivatives, we obtain the inequality 1

0 0 2 ∫ ‖ ‖2 + 2‖∇ ∗‖ ∞( ) ∫ ‖ ‖‖∇ ‖.

) 0

0 0 1 4 2

=

0 2‖∇ ∗‖ ∞( ) .

We take into account that ‖ ‖‖∇ ‖ ≤ ‖∇ ‖2 +

‖ ‖2 and estimate the last term assuming Then it follows from Gronwall’s inequality that = 0 and the solution 1 coincides with ∗. 3.2. Initial-boundary value problem modeling EVLA Let us reduce the problem ( 1 )–( 3 ) to the problem ( 4 )–( 5 ). For this purpose, we first consider the boundary value problem − div( ∇ ) + = , ∈ Ω, + 0.5 = 0 on Γ. (15)

Let the following conditions hold:

( ) 0 < 0 ≤ ( ) ≤ 1, 0 < 0 ≤ ( ) ≤ 1, ∈ Ω.

Let us define the operator ∶

→ ′, (, ) = ( ∇, ∇ ) + ( , ) + 0.5 ∫ Γ ∀ ∈ .

The problem (1)–(3) is now reduced to problem (4)–(5) if we put

→ is continuous. equation operator −1 ∶ It follows from the Lax-Milgram lemma that for any function ∈ there is a unique solution of the = , which is a weak solution of the boundary value problem (15). Moreover, the inverse = 1 + 2 −1 , 1,2 ∈ ℝ.

Theorem 3. Let conditions (i)–(i ), (r) hold. Then there exists a solution of the problem ( 1 )–( 3 ), {, } ∈ 2(0, ; ) ∩ ∞(0, ; ) × ∞(0, ; ). If, in addition, ∇ ∈ ∞( ), then the solution is unique.

4. Numerical simulation

Let us present an iterative algorithm for solving the initial-boundary value problem ( 1 )–( 3 ). At each iteration of this algorithm, a linear initial-boundary value problem is solved: Here, −1, = 1, 2, ... is the temperature field obtained at the previous iteration. The solvability of the linear problem (16) is well known, and similarly to the proof of Theorem 1, the following estimates are derived: ‖ ( )‖2 ≤ exp 2

This allows us to go to the limit in (16), concluding that is a solution of the problem ( 1 )–( 3 ).

The eficiency of laser ablation can be estimated by the behavior of temperature profiles at diferent points of the computational domain [10, 11, 12]. The analysis of temperature profiles provides the possibility to estimate whether the temperature level inside the vein is suficient for obliteration and, at the same time, whether the thermal efect outside the vein will be safe for living tissue. The calculations based on Eqs. ( 1 )–( 3 ) and for the computational domain presented in Fig. 1 are performed with the optical and thermophysical problem parameters presented in [10, 11]. The parameters and 0 are equal to 37, and

= 1. In all the calculations, the initial position of the optical fiber tip corresponds to = 5, and its motion velocity is 2 mm/s. at four observation points (1.5, 10) (inner vein wall), ( 2, 10 ), (2.5, 10) (outer vein wall), and (3.5, 10) (perivenous tissue). The source power is set as ( 1, 2) = ( 3, 7 ) (hereinafter, the power is given in watts). As can be seen from these plots, the perivenous tissue temperature remains quite safe despite a high temperature inside the vein. wavelengths: 810 nm, 1064 nm, 1470 nm, and 1950 nm. The source power is set as ( 1, 2) = ( 3, 7 ) in all cases. As can be seen from the figure, a change of the radiation wavelength has a significant efect on the behavior of the temperature profile. It is possible to ensure the same boiling duration (when the temperature is more than 95◦C) for temperature profiles corresponding to diferent wavelengths by changing the radiation power as it is shown in Fig. 4. At the same time, the temperature level in the perivenous tissue is quite safe. Here, in all cases, the ratio 1/ 2 remains unchanged as 3/7.

As can be seen from the experiments, the use of computer simulation is a promising way to determine the optimal radiation parameters that ensure an eficient and safe procedure of EVLA.

Acknowledgments

This work was supported by the Russian Foundation for Basic Research (project no. 20-01-00113 (a)) and by the Ministry of Science and Higher Education of the Russian Federation (contract no. 075-022020-1482-1, 21.04.2020) outcomes, and issues for debate, Lasers in Medical Science 29 (2014) 393–403. doi:10.1007/ s10103-013-1480-5. [13] A. E. Kovtanyuk, A. Y. Chebotarev, A. Astrakhantseva, A. Sushchenko, Optimal control of endovenous laser ablation, Optics and Spectroscopy 128 (2020) 1508–1516. doi:10.1134/ S0030400X20090131. [14] A. E. Kovtanyuk, A. Y. Chebotarev, A. Astrakhantseva, Inverse extremum problem for a model of endovenous laser ablation, Journal of Inverse and Ill-Posed Problems (2021). doi:10.1515/ jiip-2020-0118. 60

References

[1] I. V. Krasnikov , A. Y. Seteikin , A. P. Popov , Simulation of the efect of photoprotective titanium dioxide (TiO2) and zinc oxide (ZnO) nanoparticles on the thermal response and optical characteristics of skin , Optics and Spectroscopy 118 ( 2015 ) 668 - 673 . doi: 10 .1134/S0030400X15040116.

[2] A. E. Kovtanyuk , G. V. Grenkin , A. Y. Chebotarev , The use of the difusion approximation for simulating radiation and thermal processes in the skin , Optics and Spectroscopy 123 ( 2017 ) 205 - 210 . doi: 10 .1134/S0030400X17080136.

[3] A. G. Maslovskaya , Physical and mathematical modeling of the electron-beam-induced charging of ferroelectrics during the process of domain-structure switching , Journal of Surface Investigation 7 ( 2013 ) 680 - 684 . doi: 10 .1134/S1027451013040125.

[4] R. V. Brizitskii , Z. Y. Saritskaya , Optimization analysis of the inverse coeficient problem for the nonlinear convection-difusion-reaction equation , Journal of Inverse and Ill-Posed Problems 26 ( 2018 ) 821 - 833 . doi: 10 .1515/jiip-2017-0011.

[5] R. V. Brizitskii , Z. Y. Saritskaya , Inverse coeficient problems for a non-linear convectiondifusion-reaction equation , Izvestiya Mathematics 82 ( 2018 ) 14 - 30 . doi: 10 .1070/IM8517.

[6] A. G. Maslovskaya , L. I. Moroz , A. Y. Chebotarev , A. E. Kovtanyuk , Theoretical and numerical analysis of the landau-khalatnikov model of ferroelectric hysteresis , Communications in Nonlinear Science and Numerical Simulation 93 ( 2021 ) 105524 . doi: 10 .1016/j.cnsns. 2020 . 105524 .

[7] A. Y. Chebotarev , A. E. Kovtanyuk , G. V. Grenkin , N. D. Botkin , K.-H. Hofmann , Nondegeneracy of optimality conditions in control problems for a radiative-conductive heat transfer model , Applied Mathematics and Computation 289 ( 2016 ) 371 - 380 . doi: 10 .1016/j.amc. 2016 . 05 .036.

[8] A. Y. Chebotarev , G. V. Grenkin , A. E. Kovtanyuk , N. D. Botkin , K.-H. Hofmann , Inverse problem with finite overdetermination for steady-state equations of radiative heat exchange , Journal of Mathematical Analysis and Applications 460 ( 2018 ) 737 - 744 . doi: 10 .1016/j.jmaa. 2017 . 12 . 015.

[9] S. R. Mordon , B. Wassmer , J. Zemmouri, Mathematical modeling of endovenous laser treatment (ELT ), BioMedical Engineering OnLine 5 ( 2006 ) 26 . doi: 10 .1186/ 1475 -925X- 5 -26.

[10] P. W. M. van Ruijven , A. A. Poluektova , M. J. C. van Gemert , H. A. M. Neumann , T. Nijsten , C. W. M. van der Geld , Optical-thermal mathematical model for endovenous laser ablation of varicose veins , Lasers in Medical Science 29 ( 2014 ) 431 - 439 . doi: 10 .1007/ s10103-013-1451-x.

[11] A. A. Poluektova , W. Malskat , M. J. C. van Gemert , M. Vuylsteke , C. Bruijninckx , H. A. M. Neumann , C. W. M. van der Geld , Some controversies in endovenous laser ablation of varicose veins addressed by optical-thermal mathematical modeling , Lasers in Medical Science 29 ( 2014 ) 441 - 452 . doi: 10 .1007/s10103-013-1450-y.

[12] W. Malskat , A. A. Poluektova , C. W. M. van der Geld , H. A. M. Neumann , R. Weiss , C. Bruijninckx , M. J. C. van Gemert , Endovenous laser ablation (evla): A review of mechanisms, modeling