Loïc Chaumont
email: loic.chaumont@univ-angers.fr

Jacek Małecki
email: jacek.malecki@pwr.edu.pl

THE ENTRANCE LAW OF THE EXCURSION MEASURE OF THE REFLECTED PROCESS FOR SOME CLASSES OF LÉVY PROCESSES

Keywords: 2010 Mathematics Subject Classification. 60G51 Lévy process, supremum process, reflected process, Itô measure, entrance law, stable process, subordinate Brownian motion, integral representation

We provide integral formulae for the Laplace transform of the entrance law of the reflected excursions for symmetric Lévy processes in terms of their characteristic exponent. For subordinate Brownian motions and stable processes we express the density of the entrance law in terms of the generalized eigenfunctions for the semigroup of the process killed when exiting the positive half-line. We use the formulae to study indepth properties of the density of the entrance law such as asymptotic behavior of its derivatives in time variable.

Introduction

It follows from excursion theory that the trajectories of a Lévy process can be decomposed using the excursions of the process reflected in its past infimum. This result justifies the importance of knowing the excursion measure of the reflected process and more particularly, the entrance law of this measure. There are also several interesting applications of this entrance law. First it is directly related to the potential measure of the time space ladder height process, see Lemma 1 in [START_REF] Chaumont | On the law of the supremum of Lévy processes[END_REF]. Moreover it provides a useful expression of the distribution density of the supremum of the Lévy process before fixed times, [START_REF] Chaumont | On the law of the supremum of Lévy processes[END_REF], [START_REF] Chaumont | On the asymptotic behavior of the density of the supremum of Lévy processes[END_REF]. More recently it has been involved in the study of the probability of creeping through curves of Lévy processes, [START_REF] Chaumont | Creeping of Lévy processes through deterministic functions[END_REF].

In this article we obtain integral representations of the densities and the Laplace transforms of the entrance laws of the reflected excursions for two classes of real valued Lévy processes. The first class consists of symmetric Lévy processes, with a particular emphasis on subordinate Brownian motions, when the Lévy measure of the underlying subordinator has a completely monotone density. The other class is that of stable processes. The presented formulae for symmetric processes are given in terms of the corresponding Lévy-Kchintchin exponent Ψ(ξ) and the related generalized eigenfunctions introduced by M. Kwaśnicki in [START_REF] Kwaśnicki | Spectral analysis of subordinate Brownian motions in half-line[END_REF]. In the stable case, we based the calculations on the generalized eigenfunctions studied recently by A. Kuznetsov and M. Kwaśnicki in [START_REF] Kuznetsov | Spectral analysis of stable processes on the positive half-line Electron[END_REF]. Then we used the formulae obtained for the entrance law densities to derive corresponding integral representations for supremum densities. Although the theory of Lévy processes is very rich and abounds in numerous general relationships, as those coming from the Wiener-Hopf factorizations, there are few examples where the explicit representations of the related densities are available. Apart from Brownian motion and Cauchy process, some series representations were recently found in [START_REF] Hubalek | A convergent series representation for the density of the supremum of a stable process[END_REF], [START_REF] Kuznetsov | On the density of the supremum of a stable process[END_REF], [START_REF] Hackmann | A note on the series representation for the density of the supremum of a stable process[END_REF] in the case of stable processes. A different approach was presented in [START_REF] Kwaśnicki | First passage times for subordinate Brownian motions[END_REF], where the theory of Kwaśnicki's generalized eigenfunctions were used to described the first passage time density through a barrier for subordinate Brownian motions with regular Lévy measures. This concept was generalized to non-symmetric stable processes in [START_REF] Kuznetsov | Spectral analysis of stable processes on the positive half-line Electron[END_REF]. In the present paper we stay in this framework and show that a similar approach leads to integral representations of the entrance law density, the supremum density or the density of joint distribution of the process itself and its supremum. Then we apply the obtained formulae to study the asymptotic behavior of the derivatives in time variable of the entrance law densities of the reflected excursions. Let us finally mention that these formulae can be used to perform numerical simulations and study in-depth properties of the process coupled with its past supremum.

Preliminaries

Let (X, P) be the real valued Lévy process whose characteristic exponent Ψ(ξ) is characterized in terms of the Lévy triplet (a, σ 2 , Π) by the Lévy-Kchintchin formula

Ψ(ξ) = -aiξ + 1 2 σ 2 ξ 2 - R\{0} e iξx -1 -iξ1 {|x|<1} Π(dx), ξ ∈ R.
We write P x for the law of the process starting from x ∈ R. We denote by X * = -X the dual process and P * x stands for its law with respect to P x . The past supremum and past infimum of X before a deterministic time t ≥ 0 are

X t = sup{X s ; 0 ≤ s ≤ t}, X t = inf{X s ; 0 ≤ s ≤ t}.
For given t > 0 we write f t (dx) = P(X t ∈ dx) for the corresponding distribution and f t (x) stands for its density with respect to the Lebesgue measure on (0, ∞) whenever it exists. Recalling that the reflected processes X -X and X -X are Markovian, we write L t and L * t for their local times at 0 respectively, where these are normalized in the following way

E ∞ 0 e -t dL t = E ∞ 0 e -t dL * t = 1.
We write n (and n * ) for the Itô measure of the excursions away from 0 of the reflected process X -X (resp. X -X). Our main objects of studies are the corresponding entrance laws defined by

q t (dx) = n(X t ∈ dx, t < ζ), q * t (dx) = n * (X t ∈ dx, t < ζ), t > 0,
where ζ is the life time of the generic excursion and q t (x), q * t (x) denotes the densities on (0, ∞) of q t (dx) and q * t (dx), whenever they exist. In this paper, it will always be assumed that 0 is regular for both half-lines (-∞, 0) and (0, ∞). In this case, the double Laplace transform of q t (dx) is given by

∞ 0 ∞ 0 e -ξx e -zs q s (dx)ds = 1 κ(z, ξ) , (2.1)
where κ(z, ξ) is the Laplace exponent of the ladder process (L -1 t , H t ), t < L(∞). Here L -1 t denotes the right continuous inverse of L t (ladder time process) and the ladder-height process is defined by H t = X L -1 t . Analogous relations hold for q * t (dx) and the Laplace exponent κ * (z, ξ) for the ladder process ((L * t ) -1 , H * t ). Formula (2.1) actually shows that q s (dx)ds is the potential measure of (L -1 , H). We denote by h the renewal function of the ladder height process H, that is

h(x) = ∞ 0 P(H t ≤ x) dt, x ≥ 0.
In the light of Theorem 6 in [START_REF] Chaumont | On the law of the supremum of Lévy processes[END_REF], the entrance laws q t (dx) and q * t (dx) seem to be basic objects in the study of the supremum distributions. More precisely, under our assumption that 0 is regular for both negative and positive half-lines, the representation (4.4) from [START_REF] Chaumont | On the law of the supremum of Lévy processes[END_REF] reads as

P(X t ∈ dx, X t -X t ∈ dy) = t 0 q * s (dx)q t-s (dy)ds, (2.2)
which, in particular, implies

f t (dx) = t 0 n(t -s < ζ) q * s (dx). (2.3)
Finally, for x > 0, we denote by Q * x the law of the processes killed when exiting the positive half-line, i.e.

Q * x (Λ, t < ζ) = P x (Λ, t < τ - 0 ), Λ ∈ F t , where τ - 0 = inf{t > 0 : X t < 0}.
The law Q x is defined in the same way, but with respect to the dual process. The corresponding semigroups are defined as

Q * t f (x) = Q * x f (X t ), Q t f (x) = Q x f (X t )
, for non-negative Borel functions f . We also write q * t (x, dy), q t (x, dy) and q * t (x, y), q t (x, y) for the corresponding transition probability measures and their densities whenever they exist. Recall that whenever q * t (x, •) and q t (x, •) are absolutely continuous, the duality relation holds q * t (x, y) = q t (y, x).

Symmetric Lévy processes and subordinated Brownian motions

This section is devoted to symmetric Lévy processes with some addition regularity assumptions on the Lévy-Kchintchin exponent Ψ(ξ) presented in details below. We also exclude compound Poisson processes from our considerations. Note that the symmetry assumptions simplify the general exposure presented in Preliminaries, where, roughly speaking, we can remove the notation with * . Moreover, the ladder time process is the 1/2-stable subordinator for every symmetric Lévy process, which implies that

n(t < ζ) = t -1/2 √ π , t > 0. (3.1)
Finally, we recall the integral representation of the Laplace exponent of the ladder process

κ(z, ξ) = √ z exp 1 π ξ log(1 + Ψ(ζ) z ) ξ 2 + ζ 2 dζ , z, ξ ≥ 0, (3.2)
where, in the symmetric case, Ψ(ξ) is a real-valued function.

Our first result gives the expression for the Laplace transform of q t (dx) (for fixed t > 0) in the case of symmetric Lévy processes with increasing Lévy-Khintchin exponent. This is an analogue of Theorem 4.1 in [START_REF] Kwaśnicki | Suprema of Levy processes[END_REF], where the corresponding formula for X t was derived. Note that even though the formulae for the Laplace transforms of q t (dx) and P(X t ∈ dx) seem to be similar, passing from one to the other by using (2.3) and (3.1) is not straightforward. Theorem 1. Let (X, P) be a symmetric Lévy process that is not a compound Poisson process. Assume that the Lévy-Khintchin exponent Ψ(ξ) of (X, P) is increasing in ξ > 0.

Then ∞ 0 e -ξx q t (dx) = 1 π ∞ 0 λΨ (λ) λ 2 + ξ 2 exp 1 π ∞ 0 ξ log λ 2 -u 2 Ψ(λ)-Ψ(u) ξ 2 + u 2 du e -tΨ(λ) dλ. (3.3)
Proof. The proof is based on the same idea as in the proof of Theorem 4.1 in [START_REF] Kwaśnicki | Suprema of Levy processes[END_REF] with a slight modification of the arguments. For the completeness of the exposure and the convenience of the reader we present it below. We put ψ(ξ) = Ψ( √ ξ) for ξ > 0 and define

ϕ(ξ, z) = exp - 1 π ∞ 0 ξ log(1 + ψ(ζ 2 ) z ) ξ 2 + ζ 2 dζ , which is √ z(κ(z, ξ)) -1 by (3.2).
Obviously, for fixed ξ > 0, the function ϕ(ξ, z) is a holomorphic function of z, which is positive for z > 0. Note also that lim z→0+ ϕ(ξ, z) = 0 (by monotone convergence) and lim z→∞ ϕ(ξ, z) = 1 (by dominated convergence). Moreover, as it was shown in [START_REF] Kwaśnicki | Suprema of Levy processes[END_REF], that for Im z > 0 we have

Arg ϕ(ξ, z) = - 1 π ∞ 0 ξ Arg(1 + ψ(ζ 2 )/z) ξ 2 + ζ 2 dζ ∈ (0, π/2). Thus Arg( √ zϕ(ξ, z)) ∈ (0, π) for Im z > 0. This is equivalent to h ξ (z) = ϕ(ξ, z)/ √ z being a Stieltjes function (for fixed ξ). In general, a function g(z) is said to be a Stielties function if g(z) = b 1 z + b 2 + 1 π ∞ 0 1 z + ζ ν(dζ), z ∈ C \ (-∞, 0], (3.4)
where

b 1 , b 2 ≥ 0 and ν(dζ) is a Radon measure on (0, ∞) such that min(1, ζ -1 )µ(dζ) < ∞.
The constants and a measure appearing in the definitions of Stielties functions are given by

b 1 = lim z→0+ zg(z), b 2 = lim z→∞ g(z), ν(dζ) = lim ε→0+ Im(-g(-ζ + iε)dζ). (3.5)
Note that the last limit is understood in the sense of weak limit of measures. Since

lim z→0+ zh ξ (z) = lim z→0+ √ zϕ(ξ, z) = 0, lim z→0+ h ξ (z) = lim z→0+ ϕ(ξ, z)/ √ z = 0
the constants appearing in the representation (3.4) for Stielties function h ξ (z) are zero. Moreover, for z = ψ(λ 2 ) we get

h + ξ (-z) = lim ε→0 + h ξ (-z + iε) = i ψ(λ 2 ) λ(λ + ξi) λ 2 + ξ 2 exp   1 π ∞ 0 ξ log ψ(λ 2 ) λ 2 λ 2 -u 2 ψ(λ 2 )-ψ(u 2 ) ξ 2 + u 2 du   = λi -ξ λ 2 + ξ 2 exp 1 π ∞ 0 ξ log λ 2 -u 2 ψ(λ 2 )-ψ(u 2 ) ξ 2 + u 2 du .
Therefore, by (3.5), for every z > 0 we have

ϕ(ξ, z) √ z = 1 π ∞ 0 Im h + ξ (-ζ) 1 z + ζ dζ = 1 π ∞ 0 2λψ (λ 2 ) Im h + ξ (-ψ(λ 2 )) 1 z + ψ(λ 2 ) dλ = 2 π ∞ 0 λψ (λ 2 ) λ λ 2 + ξ 2 1 z + ψ(λ 2 ) exp 1 π ∞ 0 ξ log λ 2 -u 2 ψ(λ 2 )-ψ(u 2 ) ξ 2 + u 2 du dλ = ∞ 0 e -tz 1 π ∞ 0 λΨ (λ) λ 2 + ξ 2 exp 1 π ∞ 0 ξ log λ 2 -u 2 Ψ(λ)-Ψ(u) ξ 2 + u 2 du e -tΨ(λ) dλ dt.
Thus, the Laplace transform of the right-hand side of (3.3) is equal to 1/κ(z, ξ) and the theorem follows from uniqueness of the Laplace transform.

From now on, for the rest of the section, we will follow the approach presented in [START_REF] Kwaśnicki | First passage times for subordinate Brownian motions[END_REF] and restrict our consideration to the case where (X, P) is a subordinate Brownian motion whose underlying subordinator has a complete monotone density. The process (X, P) has the latter form if and only if its characteristic exponent Ψ(ξ) can be written as Ψ(ξ) = ψ(ξ 2 ) for a complete Bernstein function ψ (see Proposition 2.3 in [START_REF] Kwaśnicki | Spectral analysis of subordinate Brownian motions in half-line[END_REF]). A function ψ(ξ) is called a complete Bernstein function (CBF) if

ψ(z) = a 1 + a 2 z + 1 π ∞ 0 z z + ζ µ(dζ) ζ , z ∈ C \ (-∞, 0], (3.6)
where a 1 ≥ 0, a 2 ≥ 0 and µ(dζ) is a Radon measure on positive half-line such that min(ζ -1 , ζ -2 )µ(dζ) is finite. As in the Stielties function representation, the above-given constants and the measure µ are determined by suitable limits as follows

a 1 = lim z→0+ ψ(z), a 2 = lim z→∞ ψ(z) z , µ(dζ) = lim ε→0+ Im(ψ(-ζ + iε)dζ). (3.7)
The spectral theory of subordinate Brownian motion on a half-line was developed by M. Kwaśnicki in [START_REF] Kwaśnicki | Spectral analysis of subordinate Brownian motions in half-line[END_REF], where the generalized eigenfunctions F λ (x) of the transition semigroup Q t of the process (X, P) killed upon leaving the half-line [0, ∞) were constructed. Some additional properties of F λ (x) were also studied in [START_REF] Kwaśnicki | First passage times for subordinate Brownian motions[END_REF]. For a fixed CBF ψ and λ > 0 the generalized eigenfunctions of Q t with eigenvalue e -tψ(λ 2 ) are given by

F λ (x) = sin(xλ + ϑ λ ) -G λ (x), x > 0 (3.8)
where the phase shift ϑ λ belongs to [0, π/2) and is given by

ϑ λ = - 1 π ∞ 0 λ λ 2 -u 2 log ψ (λ 2 )(λ 2 -u 2 ) ψ(λ 2 ) -ψ(u 2 )
du, λ > 0.

Recall the following upper-bounds (Proposition 4.3 and Proposition 4.5 in [START_REF] Kwaśnicki | First passage times for subordinate Brownian motions[END_REF] respectively)

ϑ λ ≤ sup ξ>0 ξ|ψ (ξ)| ψ (ξ) π 4 .
(3.9)

ϑ λ ≤ π 2 -arcsin λ 2 ψ (λ 2 ) ψ(λ 2 ) , λ > 0. (3.10)
The function G λ is the Laplace transform of the finite measure

γ λ (dξ) = 1 π Im λψ (λ 2 ) ψ(λ 2 ) -ψ + (ξ 2 ) exp - 1 π ∞ 0 ξ ξ 2 + u 2 log ψ (λ 2 )(λ 2 -u 2 ) ψ(λ 2 ) -ψ(u 2 )
du dξ.

Here ψ + denotes the holomorphic extension of ψ in the complex upper half-plane. The Laplace transform of F λ (x) is given by

LF λ (ξ) = λ λ 2 + ξ 2 exp 1 π ∞ 0 z z 2 + u 2 log ψ (λ 2 )(λ 2 -u 2 ) ψ(λ 2 ) -ψ(u 2 ) du . (3.11)
Recall also the following estimates

|LF λ (ξ)| ≤ |λ + ξ| |λ 2 + ξ 2 | , x > 0, Re ξ > 0.
(3.12) Proposition 5.4 in [START_REF] Kwaśnicki | First passage times for subordinate Brownian motions[END_REF] states that for unbounded ψ such that lim sup λ→0 + ϑ λ < π/2 we have the following limiting behavior

lim λ→0 + F λ (x) λ ψ (λ 2 ) = h(x), x ≥ 0, (3.13)
and the convergence is locally uniform in x ≥ 0.

The functions F λ (x) were used to find the integral representations for the density function of τ - 0 and its derivatives (see Theorem 1.5 in [START_REF] Kwaśnicki | First passage times for subordinate Brownian motions[END_REF]). In the next Theorem we show that an analogous representation can be obtained for the density of the entrance law.

Theorem 2. Let (X, P) be a symmetric Lévy process whose Lévy-Khintchin exponent

Ψ(ξ) satisfies Ψ(ξ) = ψ(ξ 2 ) for a complete Bernstein function ψ(ξ). If there exists t 0 > 0 such that ∞ 1 e -t 0 ψ(λ 2 ) λ ψ (λ 2 )dλ < ∞, (3.14)
then, for every t ≥ t 0 , q t (dx) has a density with respect to the Lebesgue measure given by the formula

q t (x) = 2 π ∞ 0 e -tψ(λ 2 ) F λ (x)λ ψ (λ 2 )dλ, x > 0. (3.15)
Proof. The definition (3.8) of F λ (x) and the fact that G λ (x) is a Laplace transform of finite measure and

G λ (0) = sin(ϑ λ ) entail that |F λ (x)| ≤ 2 for all x, λ > 0.
The assumption (3.14) together with the estimate

1 0 e -tψ(λ 2 ) λ ψ (λ 2 )dλ ≤ 1 ψ (1) 1 0 e -tψ(λ 2 ) λψ (λ 2 )dλ = 1 -e -tψ(1) 2t ψ (1)
give that the function e -ξx e -tψ(λ 2 ) |F λ (x)|λ ψ (λ 2 ) is jointly integrable on (λ, x) ∈ (0, ∞) 2 and consequently the Laplace transform of the integral appearing on the right-hand side of (3.15) is given by

∞ 0 e -tψ(λ 2 ) LF λ (ξ)λ ψ (λ 2 )dλ = ∞ 0 λ 2 ψ (λ 2 ) λ 2 + ξ 2 exp 1 π ∞ 0 ξ log λ 2 -u 2 ψ(λ 2 )-ψ(u 2 ) ξ 2 + u 2 du e -tψ(λ 2 ) dλ.
This is just the right-hand side of (3.3) with Ψ(ξ) = ψ(ξ 2 ) divided by 2/π. The uniqueness of the Laplace transform ends the proof.

It is very easy to see that if ψ is regularly varying at infinity with strictly positive order then the exponential factor in (3.14) makes the integral convergent for every t > 0. Thus we derive the following result.

Corollary 1. If ψ is CBF regularly varying at infinity with order α ∈ (0, 1], then the measure q t (dx) has a density for every t > 0 and the formula (3.15) holds for every t > 0 and x > 0.

(ξ) = ξ α/2 , α ∈ (0, 2) (symmetric stable), ψ(ξ) = ξ α/2 + ξ β/2 , α, β ∈ (0, 2) (sum o two independent stable), ψ(ξ) = (m 2 + ξ) α/2 -m, α ∈ (0, 2) (relativistic stable), ψ(ξ) = log(1 + ξ α/2 ), α ∈ (0, 2], t > 1/α (geometric stable)
Then we obtain the following straightforward consequence of the previous result. Theorem 3. Let (X, P) be a symmetric Lévy process whose Lévy-Khintchin exponent Ψ(ξ) satisfies Ψ(ξ) = ψ(ξ 2 ) for a complete Bernstein function ψ(ξ). If (3.14) holds, then for every t > t 0 the distribution of (X t , X t -X t ) is absolutely continuous with respect to the Lebesgue measure on (0, ∞) 2 with density

4 π 2 (0,∞) 2 e -tψ(λ 2 ) -e -tψ(u 2 ) ψ(λ 2 ) -ψ(u 2 ) F λ (x)F u (y)λ u ψ (λ 2 )ψ(u 2 ))du dλ.
Moreover, we have

f t (x) = 2 π 3/2 ∞ 0 e -tψ(λ 2 ) tψ(λ 2 ) 0 e u du √ u F λ (x) λ ψ (λ 2 ) dλ,
for every t > t 0 .

Proof. The proofs of both formulae are direct consequences of the integral representation (3.15), the relations (2.2), (2.3) and (3.1) together with the Fubini's theorem, which can be applied due to the integral condition (3.14).

The representation (3.15) enables to compute the derivatives of q t (x) and examine its behavior in two asymptotic regimes: as t goes to infinity and x goes to 0. It is described in the following theorem. Theorem 4. Let (X, P) be a symmetric Lévy process whose Lévy-Khintchin exponent Ψ(ξ) satisfies Ψ(ξ) = ψ(ξ 2 ) for a complete unbounded Bernstein function ψ(ξ). If there exists t 0 > 0 such that (3.14) holds, then

(-1) n d n dt n q t (x) = 2 π ∞ 0 e -tψ(λ 2 ) F λ (x)λ(ψ(λ 2 )) n ψ (λ 2 )dλ, x > 0, (3.16)
for every t > t 0 and n = 0, 1, 2, . . .. Moreover, if additionally (a) ψ is increasing, regularly varying of order α 0 ∈ (0, 1) at 0, then the following holds

lim t→∞ t n+1 ψ -1 (1/t) d n dt n q t (x) = (-1) n π Γ n + 1 2α 0 -1 h(x), x ≥ 0, (3.17)
where ψ -1 denotes the inverse of ψ, and the convergence is locally uniform in x. This also holds for α 0 = 1 with the additional assumption

sup ξ>0 |ψ (ξ)| ψ (ξ) < 2. (3.18)
(b) ψ is regularly varying at infinity with index α ∞ ∈ [0, 1] and (3.18) holds, then

lim x→0 + d n dt n q t (x) = 1 Γ(1 + α ∞ ) d n dt n p t (0) t , t > t 0 , (3.19)
where p t denotes the density of the transition semigroup of (X, P).

Proof. To justify (3.16) it is enough to show that we can interchange the derivative and the integral in (3.15). However, taking any t > t 0 , where t 0 is such that (3.14) holds, we can find t 1 ∈ (t 0 , t) and a constant

c 1 = c 1 (t 0 , t, n) such that e -tψ(λ 2 ) (ψ(λ 2 )) n λ ψ (λ 2 ) ≤ c 1 e -t 1 ψ(λ 2 ) λ ψ (λ 2 )
and the claim follows from dominated convergence. Assuming additionally, that ψ is increasing and regularly varying at 0 with index α 0 ∈ (0, 1], we get that ψ -1 is regularly varying (at 0) with index 1/α 0 (see [START_REF] Bingham | Regular Variation[END_REF]). Thus, there exists a constant

c 2 = c 2 (α 0 ) > 1 such that 1 c 2 u 1/(2α 0 ) ≤ ψ(u/t) ψ(1/t) ≤ c 2 u 2/α 0 , u > 0 (3.20)
and t ∈ (0, 1). Recall also (3.13), which asserts that under the assumptions from point (a) the function F λ (x)/ λ 2 ψ (λ 2 ) extends to a continuous function for λ ∈ [0, 1]. Here we use the upper-bound given in (3.10) for α 0 ∈ (0, 1) and (3.9) in the case α 0 = 1 to show that lim λ→0 + ϑ λ < π/2, which is required to claim (3.13). Consider the measure

µ t (dλ) = t n+1 ψ -1 (1/t) e -tψ(λ 2 ) λ 2 (ψ(λ 2 )) n ψ (λ 2 )1 [0,1] (λ)dλ
and note that it tends to a point-mass at 0 (as t → ∞), as its density function tends to 0 uniformly on [ε, 1], for every ε > 0. The mass of µ t can be calculated as follows

||µ t || = t n+1 ψ -1 (1/t) 1 0 e -tψ(λ 2 ) λ 2 (ψ(λ 2 )) n ψ (λ 2 )dλ = 1 2 ψ -1 (1/t)
tψ( 1)

0 e -u u n ψ -1 (u/t)du.
Since

lim t→∞ ψ -1 (u/t) ψ -1 (1/t) = u 1/(2α 0 ) , using (3
.20) and dominated convergence we obtain

lim t→∞ ||µ t || = 1 2 Γ n + 1 2α 0 -1 .
Finally, the expression

∞ 1 e -tψ(λ 2 ) F λ (x)λ(ψ(λ 2 )) n ψ (λ 2 )dλ
can by bounded for every t > t 1 by 2e (t-t 1 )ψ( 1)

∞ 1 e -t 1 ψ(λ 2 ) λ(ψ(λ 2 )) n ψ (λ 2 )dλ ≤ c 3 e (t-t 1 )ψ(1) ∞ 1 e -t 0 ψ(λ 2 ) λ ψ (λ 2 )dλ,
with some c 3 = c 3 (n, t 0 , t 1 ) > 0, which together with the regularity of ψ -1 at zero and estimates (3.20) implies that 1

ψ -1 (1/t) ∞ 1 e -tψ(λ 2 ) F λ (x)λ(ψ(λ 2 )) n ψ (λ 2 )dλ
vanishes uniformly in x, as t → ∞. Collecting all together we arrive at

lim t→∞ (-1) n ψ -1 (1/t) d n dt n g t (x) = 1 π Γ n + 1 2α 0 -1 h(x), x ≥ 0.
Because the justification of the fact that under assumption from point (b) we have

lim x→0 + ψ(1/x 2 ) d n dt n q t (x) = (-1) n πΓ(1 + α ∞ ) ∞ 0 e -tψ(λ 2 ) (ψ(λ 2 )) n λ 2 ψ (λ 2 )dλ,
follows in the same way as in the proof of Theorem 1.7 in [START_REF] Kwaśnicki | First passage times for subordinate Brownian motions[END_REF], we omit the proof. Note that using (3.14) we can rewrite the last integral as

(-1) n d n dt n ∞ 0 e -tψ(λ 2 ) λ 2 ψ (λ 2 )dλ = (-1) n d n dt n 1 t ∞ 0 e -tψ(λ 2 ) dλ ,
where the last equality follows simply by integration by parts. Finally, the regular behavior of ψ at infinity implies that e -tψ(λ 2 ) is in L 1 (R, dλ), which in particular means that the transition probability density is given by the inverse Fourier transform

p t (x) = 1 2π
∞ 0 e -tψ(λ 2 ) e -ixλ dλ.

Combining all together we get (3.19), which ends the proof.

In addition to numerical applications of our results, they can be used to obtain more transparent representations as in the following example related to the Cauchy process.

Proposition 1. For the symmetric Cauchy process, i.e. ψ(ξ) = √ ξ, we have

q t (x) = 1 √ π sin π 8 + 3 2 arctan x t (t 2 + x 2 ) 3/4 + 1 2π 3/2 ∞ 0 y (1 + y 2 )(xy + t) 3/2 exp - 1 π ∞ 0 log(y + s) 1 + s 2 dy,
where x, t > 0. Then the density f t (x) of the past supremum at time t of (X, P) can be derived from the above expression together with (2.3) and (3.1).

Proof. Since ψ(ξ) = ξ 1/2 , ψ (ξ) = 1/(2 √ ξ) the formula (3.15) reads as

q t (x) = √ 2 π ∞ 0 e -tλ F x (λ) √ λdλ,
where we used the scaling property

F λ (x) = F 1 (λx) = F x (λ)
. By the Plancherel's theorem we get, for fixed b ∈ (0, t), that

q t (x) = ∞ 0 e -bλ F x (λ) e -(t-b)λ √ 2λ π dλ = 1 2π ∞ -∞ LF x (b + is)Lf (t -b + is)ds = 1 2πi LF x (b + is) b+i∞ b-i∞ LF x (z)Lf (t -z)dz, (3.21)
where f (x) = √ 2x/π. The Laplace transform of f can easily be computed as follows

Lf (z) = √ 2 π ∞ 0 e -zx √ xdx = 1 √ 2π 1 z 3/2 , Re(z) > 0. Formula (3.11) gives LF x (z) = 1 √ 2 x x 2 + z 2 exp 1 π ∞ 0 z log(1 + u/x) z 2 + u 2 du .
Substituting u = z/s in the last integral we get

∞ 0 z z 2 + u 2 log(1 + u/x)du = ∞ 0 log 1 + z xs 1 + s 2 ds = ∞ 0 log(z/x + s) 1 + s 2 ds,
where the last equality follows from the fact that

∞ 0 log s ds 1 + s 2 = 1 0 + ∞ 1 log s ds 1 + s 2 = 1 0 log s ds 1 + s 2 + 1 0 log(1/s) ds 1 + s 2 = 0.
Finally, the function

B(z) = 1 π ∞ 0 log(z + u) 1 + u 2 du, z ∈ C \ (-∞, 0],
studied in details in [START_REF] Kulczycki | Spectral properties of the Cauchy process on halfline and interval[END_REF], is holomorphic in the region. We recall (see (3.13) in [START_REF] Kulczycki | Spectral properties of the Cauchy process on halfline and interval[END_REF]) that

B(i) = log 2 2 + i π 8 (3.22)
and (see (4.1) in [START_REF] Kulczycki | Spectral properties of the Cauchy process on halfline and interval[END_REF]) that

e B(z) = (1 -izσ(z))e -B(-z) , (3.23)
where σ(z) = 1 for Im(z) > 0 and σ(z) = -1 for Im(z) < 0. Consequently, defining (for fixed x) the function of complex variable z

G x (z) = x x 2 + z 2 1 (t -z) 3/2 exp 1 π ∞ 0 log(z/x + u) 1 + u 2 du = x x 2 + z 2 1 (t -z) 3/2 e B(z/x) , it is easy to see that G x (z) is a meromorphic function on {z ∈ C : Re(z) < t} \ (-∞, 0]
with single poles at ix and -ix. To evaluate the integral (3.21) we integrate G x over the (positively oriented) curve consisting of (see figure 3 below) (i) four horizontal segments:

γ 1 = {z : Im(z) = n, Re(z) ∈ [-n, b]}, γ 2 = {z : Im(z) = -n, Re(z) ∈ [-n, b]}, γ 3 = {z : Im(z) = 1/n, Re(z) ∈ [-n, 0]}, γ 4 = {z : Im(z) = -1/n, Re(z) ∈ [-n, 0]}; (ii) three vertical segments: γ 5 = {z : Re(z) = -n, Im(z) ∈ [1/n, n]}, γ 6 = {z : Re(z) = -n, Im(z) ∈ [-n, -1/n]}, Γ = {z : Im(z) = b, Re(z) ∈ [-n, n]} (iii) a semi-circle: γ 7 = {z : |z| = 1/n, (z) ≥ 0}. - 6 ¡ b -in in -n 6 γ 1 γ 2 γ 5 γ 6 γ 3 γ 4 γ 7 Γ - q q ix -ix q t Figure 1. The contour of integration.
First we compute the residua of G x at ix and -ix. By (3.22), we have

Res (G x , ix) = 1 2i √ 2e iπ/8 (t -ix) 3/2 , Res (G x , -ix) = - 1 2i √ 2e -iπ/8 (t + ix) 3/2 .
Since (t ± ix) 3/2 = (t 2 + x 2 ) 3/4 e ±3i/2 arctan(x/t) , we arrive at

Res (G x , ix) + Res (G x , -ix) = √ 2 sin π/8 + 3 2 arctan(x/t) (t 2 + x 2 ) 3/4 . (3.24)
Using the relation (3.23) we obtain

γ 3 + γ 4 G x (z)dz n→∞ -→ 0 -∞ x((1 -iu/x) -(1 + iu/x))e -B(-u/x) (x 2 + u 2 )(t -u) 3/2 du,
where the last integral, after substituting y = -xu, is equal to

-2i ∞ 0 y 1 + y 2 1 (t + yx) 3/2 exp - 1 π ∞ 0 log(y + u) 1 + u 2 du dy. (3.25)
Using the bounds (3.12) we can write

|G x (z)| = 2 √ π|LF x (z)Lf (t -z)| ≤ c 1 |x + z| |x 2 + z 2 | 1 |t -z| 3/2 ≤ c 2 (Im z) -5/2 .
It implies that the integrals of G x (z) over γ 1 , γ 2 , γ 5 and γ 6 vanish as n goes to infinity. Since G x (z) is bounded in the neighborhood of 0 (Re(z) > 0), the same holds for the integral over the semi-circle γ 7 . Now we can finish the computations by applying the residue theorem in order to get

1 2πi lim n→∞ n -n G x (z)dz = Res (h x , ix) + Res (G x , -ix) -lim n→∞ γ 3 + γ 4 G x (z)dz.
Taking into account (3.25) and (3.24) and dividing both sides by 2 √ π lead to the result.

Remark 2. It is also possible to find similar formula for the entrance law density of the symmetric α-stable process with index α ∈ (0, 1). Using the scaling property F λ (x) = F 1 (λx) and writing

e -tλ α = ∞ 0 e -uλ g (α) t (u)du, t > 0,
where g

(α)
t (u) is the density of the α-stable subordinator we obtain

q t (x) = √ 2α π ∞ 0 ∞ 0 e -uλ F x (λ) λ α/2 dλ g (α)
t (u)du.

The inner integral can be evaluated similarly as in Proposition 1.

Stable processes

For the rest of the paper we focus on stable processes and use the theory of the corresponding generalized eigenfunctions developed in [START_REF] Kuznetsov | Spectral analysis of stable processes on the positive half-line Electron[END_REF]. We assume that X is a stable process with characteristic exponent Ψ(ξ) = |ξ| α e πiα(1/2-ρ)sign(ξ) , ξ ∈ R.

We exclude spectrally one-sided processes from our considerations, i.e. we assume that α ∈ (0, 1] and ρ ∈ (0, 1) or α ∈ (1, 2], but then we assume that ρ ∈ (1 -1/α, 1/α). We write ρ * = 1 -ρ and define non-symmetric analogous of F 1 (x) defined in Section 3 for stable processes as follows

F (x) = e π cos(πρ) sin(x sin(πρ) + πρ(1 -ρ * )/2) + √ α 4π S 2 (-αρ * )G(x), (4.1)
where

G(x) = ∞ 0 e -zx z αρ/2-1/2 |S 2 (1 + α + αρ * /2 + iα ln(z)/(2π))| 2 dz. (4.2)
The function S 2 (z) = S 2 (z; α) is the double sine function uniquely determined by the following functional equations

S 2 (z + 1) = S 2 (z) 2 sin(πz/α) , S 2 (z + α) = S 2 (z) 2 sin(πz)
together with the normalizing condition S 2 ((1 + α)/2) = 1 (see [9], [START_REF] Koyama | Values of the double sine function[END_REF] and Appendix A in [START_REF] Kuznetsov | Spectral analysis of stable processes on the positive half-line Electron[END_REF] for equivalent definitions and further properties). We define F * (x) and G * (x) by the same formulae as in (4.1) and (4.2) but with ρ replaced by ρ * (and consequently ρ * replaced by (ρ * ) * = ρ). Note that whenever ρ > 1/2 the oscillations of F coming from the sine function are multiplied by the exponentially decreasing factor, but then F * oscillates exponentially, when x → ∞ and the situation is reversed for ρ < 1/2. The behaviour of F at zero is described by (see the proof of Lemma 2.8 in [START_REF] Kuznetsov | Spectral analysis of stable processes on the positive half-line Electron[END_REF])

F (x) = √ α 2 S 2 (αρ) Γ(1 + αρ * ) • x αρ * (1 + o(1)), x → 0 + . (4.3)
Although the constant

√ α 2 S 2 (αρ)
Γ(1+αρ * ) was not specified in [START_REF] Kuznetsov | Spectral analysis of stable processes on the positive half-line Electron[END_REF], using (1.10) and (1.19) from [START_REF] Kuznetsov | Spectral analysis of stable processes on the positive half-line Electron[END_REF], we obtain that

∞ 0 e -zx F (x)dx = √ α 2 S 2 (αρ)z -1-αρ * (1 + o(1)), z → ∞.
Consequently, using the Karamata's Tauberian theorem and the Monotone Density Theorem we obtain (4.3). Moreover, if ρ > 1/2 then

F (x) = G(x) = O(x -α-1 ), F * (x) = O(e x cos(πρ * ) ), x → ∞. (4.4)
Even though the functions F and F * do not simultaneously belong to L 2 (0, ∞) (for ρ = 1/2), they can be understood as the generalized eigenfunctions of the semigroups Q * t and Q t , respectively (see Theorem 1.3 in [START_REF] Kuznetsov | Spectral analysis of stable processes on the positive half-line Electron[END_REF]). Moreover, using Theorem 1.1 in [START_REF] Kuznetsov | Spectral analysis of stable processes on the positive half-line Electron[END_REF] the transition probability density of the process X killed when exiting the positive half-line is given by

q * t (x, y) = 2 π ∞ 0 e -tλ α F (λx)F * (λy)dλ, x, y, t > 0, (4.5)
whenever α > 1. Note that the restriction on α ensures that the exponential oscillations of F (or F * ) are suppressed by the factor e -tλ α , which makes the integral convergent. The formula (4.5) is the analogue of the integral representation for subordinate Brownian motions presented in [START_REF] Kwaśnicki | Spectral analysis of subordinate Brownian motions in half-line[END_REF]. Note also that assuming ρ = 1/2 we have F (x) = F * (x) = F 1 (x), where F λ (x) is the generalized eigenfunction defined in Section 3 for symmetric α-stable process.

In the next theorem we present a relation between the entrance laws densities q * t (x), q t (x) and the functions F , F * . Theorem 5. Let (X, P) be a stable process with parameter α > 1 and ρ ∈ (1 -1/α, 1/α) or ρ = 1/2. Then

q t (x) = √ α π S 2 (αρ * ) ∞ 0 e -tλ α F (λx)λ αρ dλ, q * t (x) = √ α π S 2 (αρ) ∞ 0 e -tλ α F * (λx)λ αρ * dλ,
for every x, t > 0.

Proof. We will exploit formula (4.5) together with the relation (see Proposition 1 in [START_REF] Chaumont | On the asymptotic behavior of the density of the supremum of Lévy processes[END_REF])

lim x→0 + q * t (x, y) h * (x) = q * t (y), y, t > 0, (4.6)
where the renewal function h * (x) of the ladder height process H * is defined, in general, by h * (x) = ∞ 0 P(H t ≤ x)dt. In the stable case H * t is the αρ * stable subordinator and

h * (x) = EH -αρ * 1 • x αρ * = x αρ * Γ(1 + αρ * ) , x ≥ 0.
Choosing c > 0 small enough and using (4.3) we can write

1 {λx<c} e -tλ α F (λx)

x αρ * F * (λy) ≤ c 1 e -tλ α λ αρ * F * (λy), where the latter function is integrable over (0, ∞) (for fixed t and y) by (4.4). Thus, by the Lebesgue dominated convergence theorem and (4. By duality, we have the corresponding integral representation for q t (x) with F * (x) and ρ * replaced by F (x) and ρ.

The analogue of Theorem 3 can now be proved. Theorem 6. Let (X, P) be a stable process with parameter α > 1 and ρ ∈ (1 -1/α, 1/α) or ρ = 1/2. The density of (X t , X t -X t ) with respect to the Lebesgue measure dxdy on (0, ∞) 2 is given by 2α sin(πρ * ) π 2 (0,∞) 2 e -tλ α -e -tu α λ α -u α F (uy)F * (λx)λ αρ u αρ * du dλ. Proof. As previously, the result follows from the integral representations for q t (x) and q * t (x), the relations (2.2), (2.3) and Fubini's theorem. However, since the ladder time process (L * t ) -1 is ρ * -stable subordinator and n * (t < ζ) = π * (t, ∞), where π * is the Lévy measure of (L * t ) -1 we have

n(t < ζ) = 1 Γ(1 -ρ * )t ρ * , t > 0,
which gives the representations for f t (x). To find the constant in the other formula we use the relations S 2 (z)S 2 (1 + α -z) = 1 and 2 sin(πz/2)S 2 (z + 1) = S 2 (z) (see (A.7) and (1.9) in [START_REF] Kuznetsov | Spectral analysis of stable processes on the positive half-line Electron[END_REF]).

Recall that spectrally one sided Lévy processes are excluded in Theorem 6. However let us note that in this case, some expressions of the law of (X t , X t -X t ) in terms of the density of X t are given in Theorem 3.1 in [START_REF] Chaumont | Short proofs in extrema of spectrally one sided Lévy processes[END_REF].

Remark 1 .

 1 The condition(3.14) is satisfied for a large class of CBFs like ψ

e 2 S 2 1 2x α ∞ c 1 eee

 2211 -tλ α F (λx) x αρ * F * (λy)dλ = √ α (αρ) Γ(1 + αρ * ) ∞ 0 e -tλ α F * (λy)λ αρ * dλ.Moreover, by (4.4), we can write for x < 1 that∞ c 1 /x e -tλ α |F (λx)F * (λy)|dλ ≤ ∞ c 1 /xe -tλ α e λ(x∨y) cos(π(ρ∨ρ * )) dλ ≤ exp -tc α -tλ α /2 e λ(1∨y) cos(π(ρ∨ρ * )) dλ,where the last integral is convergent according to our assumption α > 1. It shows thatlim x→0 + 1 x αρ * ∞ c 1 /xe -tλ α F (λx)F * (λy)dλ = 0 and consequently, by (4.6), we obtain -tλ α F (λx)F * (λy)dλ , -tλ α F * (λy)λ αρ * dλ, y, t > 0.

  u ρ * F * (λx) dλ, for every t, x > 0.

J. Małecki is supported by the Polish National Science Centre (NCN) grant no. 2015/19/B/ST1/01457.


Pablo Alvarez

Saramati Narasimhan

Simon Rouzé

Jean-Louis Dillenseger

Yohan Payan

Michael I Miga

Matthieu Chabanas

BIPHASIC MODEL OF LUNG DEFORMATIONS FOR VIDEO-ASSISTED THORACOSCOPIC SURGERY (VATS)

Keywords: Nonrigid Image Registration, Biomechanical Model, Lung, Video-Assisted Thoracoscopic Surgery

Intraoperative localization of small, low-density or deep lung nodules during Video-Assisted Thoracoscopic Surgery (VATS) is a challenging task. Localization techniques used in current practice require an additional preoperative procedure that adds complexity to the intervention and might yield to clinical complications. Therefore, clinical practice may benefit from alternative, intraoperative localization methods. We propose a nonrigid registration approach for nodule localization. Our method is based on a biomechanical model of the lung, where the lung parenchyma is represented as a biphasic medium. Preliminary results are promising, with target registration errors reduced from 28.39 mm to 9.86 mm in median, and to 3.68 mm for the nodule in particular.

INTRODUCTION

Surgical resection of lung nodules via Video-Assisted Thoracoscopic Surgery (VATS) is one of the treatments available for early stage lung cancer. In comparison to open thoracotomy, this minimally invasive procedure reduces the length of hospitalization and minimizes post-operative complications [START_REF] Falcoz | Video-assisted thoracoscopic surgery versus open lobectomy for primary non-small-cell lung cancer: a propensity-matched analysis of outcome from the European Society of Thoracic Surgeon database[END_REF]. However, at the beginning of the procedure, the insertion of surgical ports and the artificial ventilation applied only on the contralateral lung, allow air to flow into the intrapleural space. This abnormal air inflow, known as a pneumothorax, induces a collapse of the lung towards the hilum, and, therefore large anatomical deformations. As a result, the intraoperative localization of small, deep or low-density nodules becomes considerably difficult [START_REF] Gould | Evaluation of Individuals With Pulmonary Nodules: When Is It Lung Cancer?[END_REF].

In current practice, localization techniques rely on preoperative procedures for nodule marking, generally involving the placement of hookwires or dyes under CT fluoroscopy [START_REF] Keating | Novel Methods of Intraoperative Localization and Margin Assessment of Pulmonary Nodules[END_REF]. However, these procedures not only increase the cost and complexity of care, but also occasionally lead to serious complications [START_REF] Park | Comparative Effectiveness and Safety of Preoperative Lung Localization for Pulmonary Nodules[END_REF]. As an alternative, several groups have recently proposed localization methods based on intraoperative imaging only, for example, Cone Beam CT (CBCT) [START_REF] Rouzé | Small pulmonary nodule localization with cone beam computed tomography during video-assisted thoracic surgery: a feasibility study[END_REF] or Ultrasound [START_REF] Wada | Thoracoscopic ultrasonography for localization of subcentimetre lung nodules[END_REF]. While promising, these methods rely entirely on nodule visibility in the images, which may be insufficient. For instance, the fuzzy borders and low-density of ground glass opacity nodules could make them indistinguishable, especially considering the increased density of the collapsed lung. A potential solution to address this limitation could be nonrigid image registration, to transfer information from the preoperative CT to the deformed intraoperative configuration.

Many studies have addressed the registration of lung images, with methods based on image intensity [START_REF] Murphy | Evaluation of Registration Methods on Thoracic CT: The EMPIRE10 Challenge[END_REF], biomechanical models [START_REF] Seyfi | Effect of gravity on subject-specific human lung deformation[END_REF][START_REF] Al-Mayah | Deformable image registration of heterogeneous human lung incorporating the bronchial tree: Deformable image registration of heterogeneous human lung[END_REF], or a combination of both [START_REF] Han | A hybrid patient-specific biomechanical model based image registration method for the motion estimation of lungs[END_REF]. The purpose of these methods is to compensate for the lung deformations during normal breathing, mainly for conformational radiation therapy. However, lung deformations are considerably larger during VATS [START_REF] Alvarez | Lung deformation between preoperative CT and intraoperative CBCT for thoracoscopic surgery: a case study[END_REF], which increases the difficulty of the registration problem. To the best of our knowledge, only one study has proposed the use of nonrigid image registration in the context of VATS using CBCT images of a pig lung [START_REF] Uneri | Deformable registration of the inflated and deflated lung in cone-beam CT-guided thoracic surgery: Initial investigation of a combined modeland image-driven approach[END_REF].

In this preliminary work, we propose to compensate for lung deformation using a nonrigid image registration approach guided by a biomechanical model of the lung. Using clinical data with deformations consistent to those of VATS, we modeled the lung as a biphasic medium and estimated its deformation via finite element simulation. The main objective was to evaluate the capacity of the model to estimate lung deformations during VATS. We intend to use this approach as an intraoperative nodule localization strategy for VATS.

MATERIALS AND METHODS

Data

This retrospective study was conducted on a single clinical case. After a CT image was acquired, the patient had a nee- dle biopsy of a solitary nodule located in his left lung. During this procedure, the patient developed a pneumothorax that was later detected with a post-interventional Low Dose CT (LDCT) image. Although this clinical case is not a VATS intervention, the resulting lung deformation is very similar. Specifically, there is a change of patient pose from supine during the CT imaging to lateral decubitus during the LDCT imaging, as well as a pneumothorax (Fig. 1).

Several structures were manually segmented from these images: the airway tree, lung parenchyma and nodule from the CT image; and the deflated lung parenchyma, thoracic cage, and nodule from the LDCT image.

Finite element model

In this work, the biphasic model proposed in [START_REF] Miga | Model-updated image guidance: initial clinical experiences with gravity-induced brain deformation[END_REF] was used to represent the lung. This model is based on Biot's theory of soil consolidation [START_REF] Biot | General Theory of Three Dimensional Consolidation[END_REF], and has been previously used for brain shift compensation [START_REF] Sun | Near real-time computer assisted surgery for brain shift correction using biomechanical models[END_REF]. The governing equations are:

∇ • (G∇u) + ∇ G (1 -2ν) (∇ • u) -∇p = -ρg (1) -∇ • (k∇p) + k c (p -p c ) = 0 ( 2
)
where u is the displacement vector, p is the interstitial pressure, ν is the Poisson's ratio, G is the shear modulus defined by E 2(1+ν) with E as Young's modulus, ρ is a density measure defined by (ρ f -ρ s ) with ρ f and ρ s as the density of the fluid and solid phases respectively, g is the gravitational unit vector, k c is the bronchi permeability, p c is the internal bronchi pressure and k the hydraulic conductivity. In addition, we should note that the interchange of air between the lung parenchyma and small bronchi has been approximated by lumping these exchange-effects using the organ-wide distributed term k c (p -p c ). It allows the simulation of air evacuation from small-scale lower bronchi structures distributed throughout the lung parenchyma, therefore producing local volume changes.

Geometry reconstruction

Using CGAL library (www.cgal.org), two finite element meshes of the lung with four-node-tetrahedral elements were generated (Fig. 2). The first mesh (L ct ) was generated from the lung parenchyma segmented in the CT image. The second mesh (L ldct ) was generated from the the thoracic cage segmented in the LDCT image. For the latter, the underlying assumption is that the original shape of the lung (i.e. before pneumothorax) matches the shape of the thoracic cage.

Material properties

As shown in previous studies, the lung parenchyma is highly heterogeneous [START_REF] Seyfi | Effect of gravity on subject-specific human lung deformation[END_REF]. The Young's modulus E, for instance, is considered to be lowest near the diaphragm and highest near the hilum. Similarly, we have made the approximation that the hydraulic conductivity k decreases significantly with the distance from the main airways. The natural anatomic structure of the lung with its decreasing bronchi diameter scale would result in a more restricted exchange, or a lower hydraulic conductivity. Similar to the hydraulic conductivity, we have adopted a distance-to-airway dependence with Young's modulus. As the arrangement of bronchi cartilage rings and subsequent structural elements is present, an approximation of tissue heterogeneity is proposed to decrease Young's modulus with increasing distance from these structures. We have used distance ranges of approximately 15 mm (0 -15 mm, 15 -30 mm, etc.) to vary the Young's modulus as an exponential function of this distance (Fig. 2 right). This relationship is represented by:

E(d) = αe -γd + β ( 3
)
where d is the distance to the airway tree and γ is the decay constant. Constants α and β are computed based on fit values of E min and E max such that E(d min ) = E max and E(d max ) = E min . An analogous exponential function was used for the hydraulic conductivity k values, based on fit values of k min and k max .

Initial alignment

Before nonrigid deformation, registration to a common reference frame was performed. We aligned the spine from the CT and LDCT images using rigid body registration guided by Mutual Information. This procedure was performed with the help of the multi-resolution, image-based registration techniques implemented in Elastix toolbox [START_REF] Klein | elastix: A Toolbox for Intensity-Based Medical Image Registration[END_REF].

Deformation compensation strategy

The deformation experienced by the lung from the CT configuration to the LDCT configuration might be understood as a combination of three interdependent events: (1) a change in the direction of gravity and its effects on internal structures; (2) a change in patient pose that affects the position and shape of the lung; and (3) the deflation of the lung as a result of pneumothorax. The true deformation is an intricate combination of the individual deformations resulting from each one of these events. In this preliminary study, we explored the feasibility of a linear biomechanical model to estimate these lung deformations. The three deformations were estimated independently before being superposed to produce the final post-pneumothorax registered CT image.

Deformations (1) and ( 2) dealt with the estimation of the pre-pneumothorax anatomy of the lung in the LDCT image configuration. Deformation (3) used the previously estimated anatomy to calculate the deformation associated with only the pneumothorax. It is worth mentioning that our approach does not employ a moving grid from stage to stage so geometric nonlinearities are neglected. We should also note that no hereditary approach accounting for strain history has been implemented. Rather, in this preliminary work, we have estimated each condition from a single mesh pose and deformation was a simple linear combination of deformation modes.

Gravity compensation

The direction of gravity with respect to the patient is different in the CT and LDCT configurations. We compensated for this difference by first estimating a gravity-free state of the lung in the CT configuration, and then recomputing the effects of gravity in the LDCT configuration. For the gravity-free state, the deformation was computed on L ct by applying a body force with the magnitude of gravity, but in the opposite direction, i.e. the right hand side of Eq. (1) becomes positive.

Tissue parameters were assigned as described in section 2.4, using empirically-found values within the range suggested in [START_REF] Seyfi | Effect of gravity on subject-specific human lung deformation[END_REF]. The fluid medium was disregarded during the simulation to compute this deformation. In addition, nodes on the boundaries of the parenchyma were allowed to slide along the surface, while assuming no friction between the visceral and parietal pleurae.

Change of pose compensation

Notable deformation is introduced in the thorax after the change of patient pose [START_REF] Alvarez | Lung deformation between preoperative CT and intraoperative CBCT for thoracoscopic surgery: a case study[END_REF]. We used a surface matching strategy to compensate for this deformation, where displacements were imposed on surface nodes of L ct to match the surface of L ldct .

First, the borders of the thoracic cage of the CT and LDCT images were matched using a BSpline-based Free-Form Deformation (FFD) registration algorithm guided by Mutual Information, as implemented in Elastix toolbox [START_REF] Klein | elastix: A Toolbox for Intensity-Based Medical Image Registration[END_REF]. From this, we obtained a displacement map describing the deformation necessary to align the surface of L ct to that of L ldct . Then, we imposed displacements on each node of the surface of L ct , according to the displacement map computed previously. Material properties were designated as in the previous step.

Pneumothorax compensation

We modeled the pneumothorax as being caused by differences in pressure between the surface and the interior of the parenchyma. The fluid mass conservation law associated with equation ( 2) allows for the evacuation of the fluid phase (air) via the bronchi and the subsequent volumetric contraction of tissue.A positive pressure (5 cm H 2 O) was assigned to the surface nodes of L ldct , while remaining nodes had a noflux boundary condition. Bronchi permeability was addressed with k c and p c , with pressure p c equals to zero inside the airways. Additionally, the deflation of the lung is constrained by the boundaries of the thoracic cage. Therefore, we assigned slip conditions to the nodes in contact with the thoracic cage, as visualized in the LDCT image.

Since, to the best of our knowledge, the material properties used in the governing equations have never been established for the pneumothorax, we approximated these values using the trust-region nonlinear optimization method implemented in MATLAB. The 7 parameters being fit for in this inverse model were: E max , E min , k max , k min , γ, ρ and k c . We iteratively minimized the difference between our modelestimated deformation and the true deformation measured in the LDCT image. More specifically, we used a measure of surface-to-surface proximity between the current (model) and target (LDCT) deformations as the optimality criteria. That is, we minimized the objective function Ω defined as:

Ω = 1 N N i=1 (x d,i -x t,i , n t,i ) 2 ( 4
)
where N is the number of deformed surface points, x d,i is an indexed deformed-surface point, x t,i is the corresponding closest target-surface point, n t,i is the normal of the target surface on x t,i , and (•,•) is the inner product operator. The objective function Ω is a variant of the point-to-point distance that takes into account the local shape of the target surface. We clarify that x d comes from the surface of L ldct after simulation, and x t comes from the surface of the deflated lung, in the LDCT image.

Deformation superposition

The computed deformations were superposed to obtain the post-pneumothorax registered CT image. First, deformations (1) and ( 2) were added on every node of L ct . The resulting deformation field was used to warp the CT image into the pre-pneumothorax, LDCT configuration. Then, this synthetic image was warped again using deformation (3) on L ldct to produce the final post-pneumothorax registered CT image.

EXPERIMENTAL RESULTS

Our approach was validated using Target Registration Errors (TRE) on a set of 22 paired anatomical landmarks. Twentyone of those landmarks were localized in clearly identifiable airway and vessel bifurcations, covering a large portion of the lung volume. Finally, an additional landmark was defined in the center of the lung nodule.

Figure 3 collects the resulting TREs after three different levels of deformation compensation. First, we applied the rigid initial alignment only, which provides an order of magnitude for the deformation. Then, we applied the pneumothorax compensation step only, to evaluate the relevance of accounting for the gravity and the change of pose. Finally, we evaluated our whole deformation compensation strategy. A graphical representation of the spatial distribution of post-deformation anatomical landmarks is depicted in Fig. 4. Each landmark is displayed with a color indicating final TRE. In particular, the landmark positioned in the nodule had a TRE of 28.39 mm after initial alignment, 8.06 mm after correcting for the pneumothorax only, and 3.64 mm after applying the whole deformation compensation strategy.

Finally, the parameters found in the optimization procedure were as follows: E max = 0.94 kPa, E min = 0.14 kPa, k max = 1.7 × 10 -12 m 3 s/kg, k min = 1.5 × 10 -13 m 3 s/kg, γ = 0.28, ρ = 165.96 kg/m 3 and k c = 2.51 × 10 -9 Pa/s.

DISCUSSION AND CONCLUSION

Preliminary results suggest that our proposed nonrigid registration approach could compensate for lung deformations occurring during VATS. The reported TREs are promising, considering that the necessary negative margins for nodule resection in VATS are approximately 15 mm [START_REF] Mohiuddin | Relationship between margin distance and local recurrence among patients undergoing wedge resection for small (≤2 cm) non-small cell lung cancer[END_REF]. In addition, we have also shown the importance of accounting for the change of pose of the patient, and gravity. To that end, it is interesting to note that our preliminary linear superposition of the three deformations resulted in improvements in terms of TRE.

As tissue parameters of the lung are patient-dependent, we did not know the property values required to induce the pneumothorax deformation. For this reason, we used optimization as the tool to estimate these tissue parameters, which ultimately dictated the required deformation. The found values were consistent with previous studies [START_REF] Seyfi | Effect of gravity on subject-specific human lung deformation[END_REF][START_REF] Dumpuri | An atlas-based method to compensate for brain shift: Preliminary results[END_REF]. We could then use these values to test the ability of the proposed model to reproduce the lung deformation occurring during VATS.

Future work will include validating our method with more clinical data, to reflect varying degrees of pneumothorax and different nodule locations. Also, the registration approach will have to be adapted for intraoperative CBCT images, by dealing with partial lung anatomy and increased image noise.

Fig. 1 .

 1 Fig. 1. Left lung containing a solid, solitary nodule, indicated with a circle. Left: Axial cut of the CT image with the patient in supine position. Right: Axial cut of the LDCT image after pneumothorax with the patient in lateral decubitus position. The pneumothorax is indicated by arrows.

Fig. 2 .

 2 Fig. 2. Left: Tetrahedral mesh of the lung parenchyma from the CT image (L ct ). Middle: Tetrahedral mesh of the thoracic cage from the LDCT image (L ldct ). Right: Example of mesh stratification on L ldct and the corresponding values of Young's modulus E (E min = 0.1 kPa and E max = 1 kPa).

Fig. 3 .

 3 Fig. 3. TREs after three levels of deformation compensation.

Fig. 4 .

 4 Fig. 4. Left: Spatial distribution of post-deformation landmarks represented with a color code indicating TREs. Right: 3D reconstructions of the ground truth tumor (in blue), deformation compensated tumor (in green) and rigidly transformed tumor (in black). The surfaces represent L ldct before and after simulation of pneumothorax.

Acknowledgments

This work was supported by the Région Bretagne through its Allocations de Recherche Doctorale (ARED) framework; the French National Research Agency (ANR) through the frameworks Investissements d'Avenir Labex CAMI (ANR-11-LABX-0004) and Infrastructure d'Avenir en Biologie et Santé (ANR-11-INBS-0006); and the National Institutes of Health -NINDS grant R01NS049251.