Laura Vinckenbosch

Céline Lacaux

Samy Tindel

Magalie Thomassin

Tiphaine Obara

Monte Carlo methods for light propagation in biological tissues

Keywords: Light propagation, equation of radiative transfer, Markov chain Monte Carlo methods, parameter estimation. 2010 Mathematics Subject Classification. Primary : 78M31; Secondary : 65C40, 60J20

Light propagation in turbid media is driven by the equation of radiative transfer. We give a formal probabilistic representation of its solution in the framework of biological tissues and we implement algorithms based on Monte Carlo methods in order to estimate the quantity of light that is received by an homogeneous tissue when emitted by an optic fiber. A variance reduction method is studied and implemented, as well as a Markov chain Monte Carlo method based on the Metropolis-Hastings algorithm. The resulting estimating methods are then compared to the so-called Wang-Prahl (or Wang) method. Finally, the formal representation allows to derive a non-linear optimization algorithm close to Levenberg-Marquardt that is used for the estimation of the scattering and absorption coefficients of the tissue from measurements.

Introduction

The results presented in this article have initially been motivated by several research projects grounded on Photodynamic therapy (PDT), which is a type of phototherapy used for treating several diseases such as acne, bacterial infection, viruses and some cancers. The aim of this treatment is to kill pathological cells with a photosensitive drug that is absorbed by the target cells and that is then activated by light. For appropriate wavelength and power, the light beam makes the photosensitizer produce singlet oxygen at high doses and induces the apoptosis and necrosis of the malignant cells. See [START_REF] Wilson | The physics of photodynamic therapy[END_REF][START_REF] Wilson | The physics, biophysics and technology of photodynamic therapy[END_REF] for a review on PDT.

The project that initiated this work focuses on an innovative application : the interstitial PDT for the treatment of high-grade brain tumors [START_REF] Benachour | Multifunctional peptide-conjugated hybrid silica nanoparticles for photodynamic therapy and mri[END_REF][START_REF] Bechet | Photodynamic therapy of malignant brain tumours: A complementary approach to conventional therapies[END_REF]. This strategy requires the installation of optical fibers to deliver light directly into the tumor tissue to be treated, while nanoparticles are used to carry the photosensitizer into the cancer cells.

Due to the complexity of interactions between physical, chemical and biological aspects and due to the high cost and the poor reproducibility of the experiments, mathematical and physical models must be developed to better control and understand PDT responses. In this new challenge, the two main questions to which these models should answer are:

1. What is the optimal shape, position and number of light sources in order to optimize the damage on malignant cells?

2. Is there a way to identify the physical parameters of the tissue which drive the light propagation?

The light propagation phenomenon involves three processes: absorption, emission and scattering that are described by the so-called equation of radiative transfer (ERT), see [START_REF] Chandrasekhar | Radiative transfer[END_REF]. In general, this equation does not admit any explicit solution, and its study relies on methods of approximation. One of them is its approximation by the diffusion equation and the use of finite elements methods to solve it numerically (see for example [START_REF] Arridge | A finite element approach for modeling photon transport in tissue[END_REF]). An other approach, which appeared in the 1970s, is the simulation of particle-transport with Monte Carlo (MC) method (see [START_REF] Bhan | Condensed history Monte Carlo methods for photon transport problems[END_REF][START_REF] Carter | Particle-transport simulation with the Monte Carlo method[END_REF][START_REF] Zhu | Review of Monte Carlo modeling of light transport in tissues[END_REF] and references therein). This method has been extended by several authors in order to deal with the special case of biological tissues and there is now a consensus in favor of the algorithm proposed by L. Wang and S. L. Jacques in [START_REF] Wang | Monte Carlo modeling of light transport in multi-layered tissues in standard C[END_REF], firstly described by S. A. Prahl in [START_REF] Prahl | Light transport in tissue[END_REF] and S. A. Prahl et al. in [22]. This method is based on a probabilistic interpretation of the trajectory of a photon. It is widely used and there exist now turnkey softwares based on this method. However, this method is time consuming in 3D and the associated softwares lie inside some kind of black boxes. Due to a slight lack of formalism, it is difficult to speed it up while controlling the estimation error, or to adapt it to inhomogenous tissues such as infiltrating gliomas. Finally, even though there exist several methods in order to estimate the optical parameters of the tissue (see for example [START_REF] Karlsson | Inverse Monte Carlo for estimation of scattering and absorption in liquid optical phantoms[END_REF][START_REF] Fredriksson | Inverse Monte Carlo method in a multilayered tissue model for diffuse reflectance spectroscopy[END_REF][START_REF] Bargo | In vivo determination of optical properties of normal and tumor tissue with white light reflectance and an empirical light transport model during endoscopy[END_REF][START_REF] Palmer | Monte Carlo-based inverse model for calculating tissue optical properties. Part I: Theory and validation on synthetic phantoms[END_REF]), one still misses formal representations that answer to the questions of identifiability.

In the current work, we wish to give a new point of view on simulation issues for ERT, starting from the very beginning. We first derive a rigorous probabilistic representation of the solution to ERT in homogeneous tissues, which will help us to propose an alternative MC method to Wang's algorithm [START_REF] Wang | Monte Carlo modeling of light transport in multi-layered tissues in standard C[END_REF]. Then we also propose a variance reduction method.

Interestingly enough, our formulation of the problem also allows us to design quite easily a Markov chain Monte Carlo (MCMC) method based on Metropolis-Hastings algorithm. We have compared both MC and MCMC algorithms, and our simulation results show that the plain MC method is still superior in case of an homogeneous tissue. However, MCMC methods induce quick mutations, which paves the way to very promising algorithms in the inhomogenous case. Finally we handle the inverse problem (of crucial importance for practitioners), consisting in estimating the optical coefficients of the tissue according to a series of measurements. Towards this aim, we derive a probabilistic representation of the variation of the fluence rate with respect to the absorption and scattering coefficients. This leads us to the implementation of a Levenberg-Marquardt type algorithm that gives an approximate solution to the inverse problem.

Our work should thus be seen as a complement to the standard algorithm described in [START_REF] Wang | Monte Carlo modeling of light transport in multi-layered tissues in standard C[END_REF]. Focusing on a rigorous formulation, it opens the way to a thorough analysis of convergence, generalizations to MCMC type methods and a mathematical formulation of the inverse problem.

The paper is organized as follows. We derive the probabilistic representation of the solution to ERT in Section 2. In Sections 3 and 4, we describe the MC and MCMC algorithms which are compared to Wang's algorithm in Section 5. Finally, the sensitivity of the measures with respect to the optical parameters of the medium, as well as their estimation are treated in Section 6.

2 Probabilistic representation of the fluence rate

The radiative transfer equation

Let D = R 3 be the set of positions in the biological homogeneous tissue and S 2 be the unit sphere in R 3 . Let us denote the optical parameters of the tissue by µ s > 0 for the scattering coefficient, µ a > 0 for the simple absorption coefficient and µ = µ s + µ a for the total absorption coefficient (or attenuation coefficient ). Moreover, let us denote by L e (x, ω) the emitted light from x in direction ω and by L(x, ω) the quantity of light at x in the direction ω. Then the equation of radiative transfer takes the following form (see e.g. [START_REF] Prahl | Light transport in tissue[END_REF][START_REF] Arvo | Transfer functions in global illumination[END_REF]):

L(x, ω) = L i (x, ω) + T L(x, ω), x ∈ D, ω ∈ S 2 , (2.1)
where L i (x, ω) is the incident volume emittance and T : L ∞ (R 3 × S 2 ; R) → L ∞ (R 3 × S 2 ; R) is the linear operator defined on the Banach space of essentially bounded real-valued functions ℓ : R 3 × S 2 → R, given by

T ℓ(x, ω) = µ s R+ dr exp(-µr) S 2 dσ(ω) f (ω, x -ωr, ω) ℓ(x -ωr, ω), (2.2)
with f the so-called bidirectional scattering distribution function and σ the uniform probability measure on the unit sphere S 2 . The incident volume emittance L i is also defined by applying a linear operator T i to L e :

L i = T i L e , with T i ℓ(x, ω) = +∞ 0 ℓ(x -rω, ω) exp(-µr)dr. (2.3)
In the following, we will denote the albedo coefficient by ρ := µs µ < 1. Moreover, since we consider an homogeneous biological tissue, the scattering function is given by the so-called Henyey-Greenstein function, see [START_REF] Henyey | Diffuse radiation in the galaxy[END_REF], that is

f (ω, x, ω) = f HG (ω, ω) = 1 -g 2 (1 + g 2 -2g ω, ω ) 3/2 , ω, ω ∈ S 2 , ∀x ∈ D, (2.4)
where the constant g ∈ [0, 1) is the anisotropy factor of the medium. The function ω → f HG (ω, ω) is a bounded and infinitely differentiable probability density function on S 2 with respect to the uniform probability σ. It only depends on the angle θ between ω and ω since ω, ω = cos(θ). The greater the anisotropic parameter g, the more likely the large values of cos(θ) and the less the scattering of the ray (see Fig. 1).

Neumann series expansion of the solution

In general, (2.1) admits no analytical solution and a classical way to express its solution is to expand it in Neumann series. This method is based on the next classical and general result.

Theorem 1 ([31] p.69). Let B be a Banach space equipped with a norm • and A a linear operator on B. If A < 1, then the Neumann series ∞ n=0 A n converges, the operator Id -A is invertible and for any x 0 ∈ B, the equation x = Ax + x 0 admits a unique solution given by

x = (Id -A) -1 x 0 = ∞ n=0 A n x 0 .
In order to apply Theorem 1 in our context, let us now bound the norm of the operator T defined above by (2.2).

Lemma 2. The operator T defined in (2.2), with f given by (2.4), satisfies T = ρ < 1, where we recall that we have set ρ := µs µ < 1.

Proof. Let ℓ ∈ L ∞ (R 3 × S 2 ; R). We have |T ℓ(x, ω)| µ s R+ dr exp(-µr) S 2 dσ(ω) |f HG (ω, ω) ℓ(x -ωr, ω)| µ s ℓ ∞ R+ dr exp(-µr) = µ s µ ℓ ∞ ,
since f HG is a density function on S 2 . Thus, T µs µ and since T 1 ≡ µs µ , we obtain T = µs µ and the proof is complete.

As a corollary of the previous considerations, we are able to derive an analytic expansion for the solution to equation (2.1): × S 2 ; R), then the radiative transfer equation (2.1) with a phase function given by (2.4) admits a unique solution L in L ∞ (R 3 × S 2 ; R). Moreover, L can be decomposed as

Corollary 3. If L e ∈ L ∞ (R
L = ∞ n=0 T n L i = ∞ n=0 [T n • T i ] L e
where T 0 ≡ Id and where for n ≥ 1, the linear operator T n • T i on L ∞ (R 3 × S 2 ; R) is given by:

[T n • T i ] ℓ(x, ω 0 ) = µ n s R n+1 + dr 0 • • • dr n exp   -µ n j=0 r j   (S 2 ) n dσ ⊗n (ω 1 , . . . , ω n ) n-1 j=0 f HG (ω j , ω j+1 ) ℓ x - n k=0 ω k r k , ω n . (2.5) Proof. Assume that L e ∈ L ∞ (R 3 × S 2 ; R). It is readily checked from the definition (2.3) of L i that we also have L i ∈ L ∞ (R 3 × S 2 ; R). Indeed, L i ∞ ≤ L e ∞ µ < +∞.
Hence, Theorem 1 and Lemma 2 provide the existence and uniqueness of the solution, as well as its expansion in Neumann series. Formula (2.5) is then found by induction.

Our next step is now to recast representation (2.5) into a probabilistic formula.

Probabilistic representation

The Neumann expansion of T enables us to express L = ∞ n=0 T n L i as an expectation. To this aim, we now introduce some notation. Let us define

A = ∞ n=0 M n , with M n = R n+1 + × S 2 n+1 . (2.6)
We denote by (r, ω) a generic element of A and by (r n , ω n ) a generic element of M n for n ∈ N with r n = (r 0 , . . . , r n ) and ω n = (ω 0 , . . . , ω n ). If (r, ω) ∈ A, we set

|r| = ∞ n=1 n 1 Mn (r, ω) (2.7)
and call it size or length of the path. For n ∈ N, let

G (n) x (r n , ω n ) = L e x - n k=0 ω k r k , ω n (2.8)
be defined on M n , and let

G x (r, ω) = ∞ n=0 G (n) x (r n , ω n ) 1 Mn (r, ω) (2.9)
be a function on A. Let Y = (R, W ) be a A-valued random variable defined on a probability space (Ω, F , P), whose law ν is given by

ν(F ) = ∞ n=0 (1 -ρ)ρ n ν n (F ∩ M n ) , (2.10)
where we recall that ρ = µs µ and where ν n is the probability measure on M n defined by

ν n (dr n , dω n ) = µ n+1 e -µ n j=0 rj n-1 j=0 f HG (ω j , ω j+1 ) dr n σ ⊗(n+1) (dω n ) (2.11)
with -1 j=0 a j = 1 by convention. Before we express L as an expectation involving G x and Y , let us state some properties of Y = (R, W ). Proposition 4. Let n ∈ N and let π be a permutation of {0, . . . , n}. Let us recall that ν defined by (2.10) is the distribution of Y = (R, W ). Then, conditionally on the event {Y ∈ M n }, the distribution of the variable (R 0 , . . . , R n , W 0 , . . . , W n ) is the probability measure ν n defined by (2.11), which satisfies ν n (dr 0 , . . . , dr n , dω 0 , . . . , dω n ) = ν n (dr π(0) , . . . , dr π(n) , dω n , dω n-1 , . . . , dω 0 ).

(2.12)

In other words, on the event {Y ∈ M n }, the random variables (R π(0) , . . . , R π(n) , W n , W n-1 , . . . , W 0 ) and (R 0 , . . . , R n , W 0 , . . . , W n ) have the same distribution ν n .

Proof. By definition of ν = L(Y ) and since (M p ) p∈N is a collection of pairwise disjoint sets, it is straightforward that ν n is the distribution of the variable (R 0 , . . . , R n , W 0 , . . . , W n ) on the set {Y ∈ M n }. Let us now emphasize that the phase function is symmetric, i.e. f HG (ω, •) = f HG (•, ω).

Then replacing f HG (ω j , ω j+1 ) by f HG (ω j+1 , ω j ) and using the invariance of the definition of ν n by permutations of the variables (r 0 , . . . , r n ), we obtain Equation (2.12).

Corollary 5. For all n ∈ N, conditionally on the events {Y ∈ M n } and {Y ∈ p≥n M p }, for any j = 0, . . . , n, the marginal distribution γ Wj of the direction W j is the uniform probability σ on S 2 . In particular, W 0 is uniformly distributed on S 2 .

Proof. Let n, j ∈ N such that j ≤ n. By proposition 4, on the event {Y ∈ M n }, the probability measure ν n defined by (2.11) is the distribution of (R n , W n ). Then, integrating the law ν n with respect to all variables except w j and using that f HG (ω, •) = f HG (•, ω) is a density function on S 2 , one obtains that on the event {Y ∈ M n } , W j is uniformly distributed on S 2 . Since this also holds replacing n by any p ≥ n and since M p , p ≥ n, are pairwise disjoint sets, this implies that on the event {Y ∈ p≥n M n }, W j is also uniformly distributed on S 2 . Proposition 6. The series ∞ n=0 T n L i can be expressed as

L(x, ω) = ∞ n=0 T n L i (x, ω) = 1 µ a E [G x (Y ) | W 0 = ω] . (2.13)
Proof. We shall relate [T n • T i ]L e with the measure ν n defined above, which will be sufficient for our purposes. To this aim, consider ω0 ∈ S 2 , and write a somehow more cumbersome version of formula (2.5) with ℓ = L e :

[T n • T i ] L e (x, ω0 ) = µ n s R n+1 + dr 0 • • • dr n exp   -µ n j=0 r j   (S 2 ) n dσ ⊗n (ω 1 , . . . , ω n ) f HG (ω 0 , ω 1 ) n-1 j=1 f HG (ω j , ω j+1 ) L e x -ω0 r 0 - n k=0 ω k r k , ω n .
Noting that µ n s = (1-ρ) ρ n µ n+1 µa , we thus get the following identity:

[T n • T i ] L e (x, ω0 ) = (1 -ρ) ρ n µ n+1 µ a R n+1 + dr 0 • • • dr n exp   -µ n j=0 r j   (S 2 ) n+1 dσ n,ω0 (ω 0 , . . . , ω n ) n-1 j=0 f HG (ω j , ω j+1 ) L e x - n k=0 ω k r k , ω n ,
where the measure σn,ω0 on S n+1 is given by σn,ω0 (dω 0 , . . . , dω n ) = δ ω0 (dω 0 ) ⊗ σ ⊗n (dω 1 , . . . , dω n ). Finally set ϕ n (r n , ω n ) = L e (x-n k=0 ω k r k , ω n ). Taking into account the identity above and (2.11), we easily get

[T n • T i ] L e (x, ω0 ) = (1 -ρ) ρ n µ a ν n (ϕ n | ω 0 = ω0 ) ,
from which our claim is straightforward.

Remark 7. In the following, we will call the random variable Y = (R, W ) a ray. Notice that it does not correspond exactly to a ray of light in the physics sense, since Y has a finite length (though random) and since a given realization of Y does not carry the information due to light absorption. Also notice that Y owns a complete probabilistic description which allows to exactly simulate it (see Proposition 10 for simulation considerations).

Model for light propagation

Observe that our formula (2.13) induces a Monte Carlo procedure to estimate L(x, ω) for each (x, ω) ∈ R 3 × S 2 based on the simulation of independent copies of Y . Nevertheless this procedure is time consuming. Indeed, assuming that the light is only emitted by an optical fiber, many realizations y of Y lead to a null contribution in the estimation of L(x, ω). Our aim is now to accelerate our simulation by means of a coupling between random variables corresponding to different (x, ω). Towards this aim, we now focus on an averaged model for light propagation.

Let thus V ⊂ R 3 be a cube whose center coincides with the origin. We discretize it into a partition of K smaller cubes (voxels in the image processing terminology) {V k , k = 0, . . . , K -1}, whose volume equals h 3 , h ∈ R + and such that the origin is the center of V 0 . Let us denote by x k the center of the voxel V k . We work under the following simplified assumption for the form of the light source: Hypothesis 8. We assume that the only emission of light in the domain V comes from the optical fiber. Let C 2α ⊂ S 2 denote the cone with opening angle 2α, whose summit is placed at the origin and whose axis follows -e 3 . The light source is defined by S = {(x, ω) : x ∈ V 0 , ω ∈ C 2α }. We assume that the emission of light satisfies

L e (x, ω) = c 1 V0×C 2α (x, ω) := c, if (x, ω) ∈ V 0 × C 2α , 0, otherwise, (2.14)
where c > 0 is a given constant.

This model remains close to reality and it is possible to refine it by weighting the light directions of the source in order to stick better to the shape of the fiber. With Hypothesis 8 in mind, we are interested in estimating the fluence rate at the center of the voxels V k , k = 0, that is the mean light intensity averaged in all directions

L(x k ) := S 2 L(x k , ω 0 ) σ(dω 0 ). (2.15)
This quantity admits a nice probabilistic representation.

Proposition 9. Let k ∈ {0, . . . , K -1} and let Y = (R, W ) be a random variable with distribution ν defined by (2.10). Then, the fluence rate L(x k ) at the center x k of the voxel V k , which is defined by (2.15), can also be expressed as

L(x k ) = c µ a P   x k - |R| j=0 R j W j ∈ V 0 , W |R| ∈ C 2α   (2.16)
where we recall that the length |R| of the ray Y is defined by (2.7).

Proof. Invoking Proposition 6 and by definition of L, we get

L(x k ) = 1 µ a S 2 E [G x k (Y )|W 0 = ω 0 ] σ(dω 0 )
where G x k is defined by (2.9). Since by Corollary 5, σ is the distribution of W 0 , the previous equation can be written as

L(x k ) = 1 µ a E [G x k (Y )] .
Then using equations (2.8) and (2.7), we get

G x k (Y ) = L e   x k - |R| j=0 R j W j , W |R|   .
Now applying Hypothesis 8, the random variable G x k (Y ) can be expressed as

G x k (Y ) = c 1 x k -|R| j=0 Rj Wj ∈V0, W |R| ∈C 2α ,
which finishes our proof. Now, instead of seeing Y as a ray starting at x k which possibly hits the light source V 0 × C 2α , we can imagine that it starts at the center of the light source and it possibly hits the voxel V k in any direction. This possibility stems from the invariance of ν stated in Proposition 4, and is exploited in the next result.

Proposition 10. For any 0 k K -1 we have :

L(x k ) = c(1 -cos α) 2µ a P (S N ∈ V k ) , (2.17)
where N ∼ NB(1, ρ) is a negative binomial random variable with parameter (1, ρ), that is a random variable whose law is given by P(N = n) = (1 -ρ)ρ n for all n ∈ N, and where for n ∈ N

S n = n i=0 R i W i with (R i ) i≥0
and (W i ) i≥0 satisfying the following assertions:

• (R i ) i≥0 is a sequence of independent identically distributed (i.i.d.) exponentially random variables of parameter µ (i.e such that E(R i ) = µ -1 ).

• W 0 is uniformly distributed on the cone C 2α .

• for any i ≥ 1, the conditional distribution of

W i given (W 0 , . . . , W i-1 ) = (ω 0 , . . . , ω i-1 ) is f HG (ω i-1 , ω i )σ(dω i ).
• N , (R i ) i 0 and (W i ) i 0 are independent.

Proof. Let k ∈ {0, . . . , K -1}. Notice that, if V k + x denotes the translation of the voxel V k by the vector x ∈ R 3 , then it is clear that V 0 + x k = V k . Therefore, we can rewrite (2.16) as

L(x k ) = c µ a P   |R| j=0 R j W j ∈ V k , W |R| ∈ C 2α   , (2.18)
where the distribution of Y = (R, W ) is the probability measure ν defined by (2.10). Then,

L(x k ) = c µ a +∞ n=0 (1 -ρ)ρ n P   n j=0 R j W j ∈ V k , W n ∈ C 2α  
where the distribution of (R 0 , . . . , R n , W 0 , . . . , W n ) is the probability ν n defined by (2.11). Therefore, applying Proposition 4, we get

L(x k ) = c µ a +∞ n=0 (1 -ρ)ρ n P   n j=0 R j W j ∈ V k , W 0 ∈ C 2α   .
By definition of ν n , we easily see that

P n i=0 R i W i ∈ V k , W 0 ∈ C 2α = σ(C 2α )P n i=0 R ′ i W ′ i ∈ V k ,
where (R ′ 0 , . . . , R ′ n , W ′ 0 , . . . , W ′ n ) is a random variable with distribution ν ′ n , defined by replacing in (2.11) the uniform probability σ(dω 0 ) on the sphere S 2 by the uniform probability

1 {ω ′ 0 ∈C 2α } σ(dω ′ 0 )/σ(C 2α ) on the cone C 2α . By definition of ν ′ n , this means that P n i=0 R i W i ∈ V k , W 0 ∈ C 2α = σ(C 2α )P(S n ∈ V k )
where S = (S p ) p≥0 is the random walk defined in the statement of the proposition. Thanks to the fact that N is a NB(1, ρ) random variable independent of the random walk S, we easily get relation (2.17).

Remark 11. The random variables N, R i and W 0 in Proposition 10 are simulated in a straightforward way. The simulation of the sequence (W i ) i 1 is obtained as follows. The direction W i of the i-th step of the random walk is sampled relatively to the direction W i-1 . We sample the spherical angles (Θ i , Φ i ) between the two directions according to the Henyey-Greenstein phase function. Namely, cos (Θ i ) = W i-1 , W i owns the following invertible cumulative distribution function:

F -1 (y) = 1 2g 1 + g 2 - 1 -g 2 1 -g + 2gy 2 , y ∈ [0, 1]
and Φ i , the azimuth angle of W i in the frame linked to W i-1 , is uniformly distributed on [0, 2π], see (2.4). To recover the cartesian coordinates of the directions, we inductively apply appropriate changes of frame. The corresponding formulas can be found in [21, p. 37].

Remark 12. Notice that S N does not depend on the voxel x k under consideration. This permits to use a single sample of realizations of this random variable in order to estimate the right hand side of (2.17) for all k ∈ {0, . . . , K -1} simultaneously. We call this improvement coupling, meaning that the random variables related to the Monte Carlo evaluations at different voxels are completely correlated.

Monte Carlo approach with variance reduction

In the last section, we have derived a probabilistic representation of L(x k ) for every voxel V k by means of the arrival position of a random walk (S n ) n 1 stopped at a negative binomial time. This classically means that L(x k ) can be approximated by MC methods. We first derive the expression of the approximate fluence rate by means of the basic MC method and then describe the variance reduction method that we implemented.

Proposition 13. Let us consider a random walk S = (S n ) n≥0 and a negative binomial random variable N as defined in Proposition 10. Let (S i , N i ) 1≤i≤M be M independent copies of (S, N ). Then, for k = 0, . . . , K -1,

L MC (x k ) := c (1 -cos α) 2µ a M M i=1 1 S i N i ∈V k (3.1)
is an unbiased and strongly consistent estimator of L(x k ).

Proof. This statement follows simply from the discussion of the previous section and the law of large numbers (LNN).

In addition to this proposition, let us highlight the fact that central limit theorem provides confidence intervals for the estimators L MC (x k ). Furthermore, owing to Remark 12, the family (S i , N i ) 1≤i≤M enables to estimate the fluence rate L(x k ) for all k = 0, . . . , K -1 at once. The reader should be aware of the fact that the quantity ρ is large in general in biological tissues, which means that the size of the ray, given by N ∼ NB(1, ρ), will often be large. Typical values of the parameters are provided in Table 1. Therefore, sampling a ray is relatively time consuming and it is necessary to improve the basic Monte Carlo algorithm in order to reduce the variance of the estimates.

µ s µ a µ ρ g
healthy tissue 280 cm -1 0.57 cm -1 280.57 cm -1 0.998 0.9 tumor 73 cm -1 1.39 cm -1 74.39 cm -1 0.981 0.9

Table 1: Values given by [START_REF] Angell-Petersen | Determination of fluence rate and temperature distributions in the rat brain; implications for photodynamic therapy[END_REF] for the optical parameters in the rat brain for a wavelength λ = 632 nm (red light).

Furthermore, because of the formulation (2.18), only the last point of each whole ray is used in the estimation. It is however possible to take into account more points of the rays and still have unbiased estimators. Finally, the angular symmetry of the problem, allows us to replicate observed rays by applying rotation. We took these two considerations into account and named the resulting method Monte Carlo with some points (MC-SOME). The idea is to firstly draw some random walks which share the same initial direction and to pick a given number of points of each walk. Then, we apply rotations to that set of points with respect to different initial directions. We finally count the number of points in each voxel. This artificially increase the size of the samples and thus reduce the variance of our estimation of L(x k ). Specifically, the resulting estimator is given by : Definition 14 (MC-SOME). Let M, M points , M rot ∈ N * be the parameters of the method. Let us assume that the following assertions hold:

• (W j 0 ) 1≤j≤Mrot are i.i.d. copies of W 0 ∼ U(C 2α ).
• N ℓ i , 1 i M , 1 ℓ M points , are i.i.d. copies of a negative binomial random variable with parameter (1, ρ).

• S i n n≥0 , 1 ≤ i ≤ M , are i.i.d. copies of the random walk S defined in Proposition 10, all sharing the same initial direction W 1 0 .

• The sequences N ℓ i i,ℓ and W j 0 , S i n i,j are independent.

Let S i,j N ℓ i denotes the N ℓ i -th point of the i-th random walk S i n n≥0 after a rotation corresponding to the j-th initial direction W 0 = W j 0 . Then, for k ∈ {0, . . . , K -1}, the MC-SOME estimate of L(x k ) is defined by

L MC-SOME (x k ) = c (1 -cos α) 2µ a M rot M points M M i=1 Mpoints ℓ=1 Mrot j=1 1 S i,j N ℓ i ∈V k . (3.2)
This estimator is unbiased and strongly consistent. Its construction is illustrated in Fig. 2 and it will be compared to Wang's algorithm estimator in Section 5.

Figure 2: Description of MC-SOME method. In this example, the grey path is a rotation of the black one with respect to its initial direction ω j+1 0 and M = 4.

Remark 15. Choosing only some points (S N ℓ , 1 ℓ M points ) of the path instead of all (S i , 1 i n) provides a more efficient estimation. Indeed, it would take a lot of time to run over all voxels V k in order to evaluate the indicator functions 1 {Si∈V k } . Moreover, the information brought by close points is in a sense redundant. Choosing the points according to a negative binomial law maintain the estimator unbiased, while speeding up the estimation.

A Metropolis-Hastings algorithm for light propagation

Inspired by results in computer graphics, see [START_REF] Jensen | Monte Carlo ray tracing[END_REF][START_REF] Shirley | Monte Carlo and quasi-Monte Carlo methods for computer graphics[END_REF][START_REF] Veach | Robust Monte Carlo methods for light transport simulation[END_REF], we implemented a Metropolis Hastings algorithm which is a Markov chain Monte Carlo method (MCMC) by random walk [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF][START_REF] Roberts | General state space Markov chains and MCMC algorithms[END_REF]. We shall first discuss general principles and then practical implementation issues.

For simplicity reasons, by slightly abusing the notations, we identify the stopped random walk S = (S n ) 0 n N of Proposition 10 with the ray (R 0 , . . . , R N , W 0 , . . . , W N ) which defines it. The law of this ray, which will still be denoted by ν, is given by replacing in (2.10) the uniform measure σ(dω 0 ) on the sphere by the uniform measure on the cone C 2α .

A realization of the walk S stopped at N = n will be indifferently referred to as (S 0 , . . . , S n ), as (r 0 , . . . , r n , ω 0 , . . . , ω n ) or as (r 0 , θ 0 , ϕ 0 , . . . , r n , θ n , ϕ n ) where (θ 0 , ϕ 0 ) are the spherical coordinates of ω 0 and for 1 i n, cos(θ i ) = ω i-1 , ω i and ϕ i is the azimuth angle of ω i in the frame linked to ω i-1 .

General principle

For a given ω 0 ∈ C 2α and for 0 k K -1, we are willing to estimate the conditional probability P (S N ∈ V k | W 0 = ω 0 ) by generating a Markov chain whose steady-state measure is the conditional distribution ν| W0=ω0 and by applying LNN for ergodic Markov chains. We then combined this estimation with the classical LNN sampling the initial direction W 0 in C 2α to obtain an estimate of L(x k ) viewed as in (2.17).

An overview of the MCMC dynamics in this context is the following: Let ω 0 ∈ C 2α be a fixed initial direction. The Markov chain starts at time t = 1 in the state S(1) ∈ A with W 0 = ω 0 . At each time t ∈ N * , a move (mutation) is proposed from the current state S(t) to the state S ′ (t) according to a proposal density q(S(t), •) and such that the initial direction of S ′ (t) is still ω 0 . The chain then jumps to S ′ (t) with acceptance probability α(S(t), S ′ (t)) or stays in S(t) with probability 1 -α(S(t), S ′ (t)). This is described in pseudo-code in Algorithm 1. The MCMC simulation generates a Markov chain {S(t); t 1} on the space of rays A whose steady-state measure is the desired distribution ν| W0=ω0 (see [26, §2.3.1]).

Algorithm 1 Metropolis-Hastings algorithm for light propagation

Initialization: draw ω 0 uniformly on C 2α , draw S(1) according to ν| W0=ω0 for t = 1 to T -1 do S ′ (t) ∼ q(S(t), •) α(S(t), S ′ (t)) ← min 1, ν| W0=ω0 (S ′ (t)) q(S ′ (t), S(t)) ν| W0=ω0 (S(t)) q(S(t), S ′ (t)) if Rand()< α(S(t), S ′ (t)) then S(t + 1) ← S ′ (t) else S(t + 1) ← S(t) end if end for If the ray S(t) = (S 0 (t), . . . , S Nt (t)) denotes the position of the chain at a given time t, then, for We can then sample the law of W 0 to recover an estimate of P (S N ∈ V k ). Let M rot ∈ N * and let (ω 1 0 , . . . , ω Mrot 0 ) be a sample of i.i.d. initial directions drawn according to U C 2α . For i = 1, . . . , M rot and t = 1, . . . , T , let S (i) (t) denote the rotation of the random walk S(t) with respect to the initial direction ω i 0 , see Fig. 2. For k ∈ {0, . . . , K -1}, our Metropolis-Hastings estimator of L(x k ) is defined by

0 k K -1, lim T →∞ 1 T T t=1 1 {SN t (t)∈V k } = P (S N ∈ V k | W 0 = ω 0 ) almost surely, ( 4
L MH (x k ) = c (1 -cos α) 2µ a M rot T Mrot i=1 T t=1 1 S (i) N t (t)∈V k . (4.2)
This estimator is strongly consistent on condition that the proposal density q(•, •) of Algorithm 1 provides a Harris recurrent chain.

Mutation strategy

One of the delicate issues in the implementation of MCMC methods is the choice of a convenient proposal density q. In our case, we have tested a mixture of mutations of two types for the Metropolis-Hastings algorithm: rotation-translation and deletion-addition. In order to describe them, let us introduce some notations. The method uses a perturbed phase function f εg HG , ε ∈ [-1, 1] which has an anisotropic coefficient εg instead of g. It also uses two coprime integers 1 j < J, which denote respectively the small and the big size of a length change.

Remark 16. The need for these two numbers to be coprime will be made clearer in the proof of Proposition 19, as it will ensure that the mutations produce paths of arbitrary length.

Definition 17 (Mutation rule).

Let us assume that, at time t ∈ N * , the current ray is given by S(t) = (r 0 , θ 0 , ϕ 0 , . . . , r nt , θ nt , ϕ nt ). Our proposition for the next move from S(t) to S(t + 1) is the following. (i) With probability 1 2 , the mutation is of type deletion-addition. Draw ∆(n t ) according to the following law that depends on the size of the current ray

∆(n) ∼      U ({-J, -j, j, J}) , if n J, U ({-j, j, J}) , if j n < J, U ({j, J}) , if 0 n < j.
• If ∆(n t ) < 0, then delete the last |∆(n t )| edges of S(t). The proposed path is S ′ (t) = r 0 , θ 0 , ϕ 0 , . . . , r nt-∆(nt) , θ nt-∆(nt) , ϕ nt-∆(nt) .

• If ∆(n t ) > 0, then add ∆(n t ) new edges at the end of S(t):

-Draw r new nt+1 , . . . , r new nt+∆(nt) i.i.d. according to E(µ).

-Draw θ new nt+1 , . . . , θ new nt+∆(nt) i.i.d. according to f εg HG .

-Draw ϕ new nt+1 , . . . , ϕ new nt+∆(nt) i.i.d. uniformly on [0, 2π].

The proposed path is (ii) With probability 1 2 , the mutation is of type rotation-translation. Choose an index i uniformly over {0, . . . , n t }.

S ′ (t) = r 0 , θ 0 , ϕ 0 , . . . ,
• If i = 0, then make a rotation of the path at the i-th edge:

-Draw a new angle θ new The proposed path is (see Fig. 3)

S ′ (t) = (r 0 , θ 0 , ϕ 0 , . . . , r i , θ new i , ϕ new i , . . . , r nt , θ nt , ϕ nt ) .
• If i = 0, then translate from the initial edge:

-Draw a new edge length r new 0 according to E(µ).

The proposed path is S ′ (t) = (r new 0 , θ 0 , ϕ 0 , r 1 , θ 1 , ϕ 1 , . . . , r nt , θ nt , ϕ nt ) .

Remark 18. The initial direction is fixed, thus mutations of (θ 0 , ϕ 0 ) are forbidden. However, the length r 0 has to change from time to time. This is ensured by the translations and this is why the initial edge need to be considered separately in the mutation rule. Let us now compute the proposal density q(S, •) of this mutation rule. For m ∈ N * , let

ζ(m) =      1 4 , if m J, 1 3 , if j m < J, 1
2 , if 0 m < j. Assume that S ′ is a mutation of S, denote by i the first index where there is a difference between them, denote by n ′ and n their respective length and set ∆ = n ′ -n. We have

• if ∆ = 0 and i 1, then q(S, S ′ ) = 1 2 1 n+1 1 2π f εg HG (θ ′ i );
• if ∆ = 0 and i = 0, then q(S, S ′ ) = 1 2 1 n+1 µe -µr ′ 0 ;

• if ∆ < 0, then q(S, S ′ ) = 1 2 ζ(n);

• if ∆ > 0, then q(S, S ′ ) = 1 2 ζ(n)e -µ ∆ k=1 r ′ n+k µ 2π ∆ ∆ k=1 f εg HG θ ′ n+k .
From these formulas, it is straightforward to recover the acceptance probability.

The idea behind this mixture of mutations is to find a compromise between large jumping size of the Markov chain which implies a lot of "burnt" samples, and smaller jumps which provide more correlated samples, hence a worse convergence. The rotations lead to a good exploration of the domain at low cost, whereas the addition-deletion mutations ensure the visit of the whole state space A with W 0 = ω 0 . The use of a perturbed phase function decreases the acceptance probability of the mutations and thus, increases the number of samples needed in order to converge to the invariant measure. But, it allows a better exploration of the domain and this why the parameter ε, as well as the sizes j, J, need to be adapted on a case by case basis. Finally, we can prove that, with this rule of mutations, Algorithm 1 produces a Markov chain that satisfies the LLN. This guarantees the convergence of the estimator defined in (4.2).

Proposition 19. If the chain (S(t)) t∈N * is obtained by Algorithm 1 with the mutation rule given in Definition 17, then it is Harris positive with respect to the measure ν| W0=ω0 and the estimator L MH (x k ) defined in (4.2) is strongly consistent for all 0 k K -1.

Proof. The fact that ν| W0=ω0 is an invariant measure of (S(t)) t∈N * is an inherent property of Metropolis-Hastings algorithm ( [START_REF] Roberts | Harris recurrence of Metropolis-within-Gibbs and trans-dimensional Markov chains[END_REF][START_REF] Tierney | Markov chains for exploring posterior distributions[END_REF]). The Harris recurrence is then obtained by checking that the chain is irreducible with respect to ν| W0=ω0 , see [START_REF] Tierney | Markov chains for exploring posterior distributions[END_REF]Corollary 2]. Let τ A = inf {t ∈ N * : S(t) ∈ A} denote the hitting time of any A ⊂ A such that ν| W0=ω0 (A) > 0. We must demonstrate that (S(t)) t∈N * is irreducible with respect to ν| W0=ω0 , that is,

P s (τ A < +∞) > 0, for all s ∈ A, (4.3)
where P s (S(1) = s) = 1. Furthermore, notice that it is sufficient to check this property for subsets A of the type

A = {ω 0 } × I 0 × n i=1 (I i × E i ) , (4.4)
where n ∈ N * and where for all 1 i n, the sets I i ⊂ R + and E i ⊂ S 2 are all sets of positive Lebesgue measure.

In order to prove relation (4.3) for sets of the form (4.4), consider the conditional measure ν n | W0=ω0 using (2.11). By assumptions on A, we see easily that ν| W0=ω0 (A) > 0. Now, notice that by the Markov property, if τ A and τ {ω0}×I0 denote respectively the time for the chain to be in A, resp. in {ω 0 } × I 0 , then we have

P s (τ A < +∞) P s τ {ω0}×I0 < +∞ P {ω0}×I0 (τ A < +∞) ,
for all s ∈ A.

Now we can lower bound the right hand side of this relation in the following way:

(i) We have that P s τ {ω0}×I0 < +∞ is greater than the probability of deleting all the edges of s except (r 0 , ω 0 ) and of modifying its length so that r ′ 0 ∈ I 0 . This probability is strictly positive, as well as its acceptance. Indeed, we use here the fact that j and J are coprime (through Bezout's lemma) plus elementary relations for uniform distributions to assert that the probability of deleting all the edges is strictly positive. The positivity of acceptance is due to absolute continuity properties of q(s, •). (ii) The same kind of argument works in order to lower bound P {ω0}×I0 (τ A < +∞). Namely, this quantity is greater than the probability to construct directly a ray s ∈ A, which is itself strictly positive. Indeed, since j and J are mutually prime, it is possible to construct a ray of any desired length. Moreover, at each step, the probability of adding an edge (r ′ i , ω ′ i ) ∈ I i × E i , as well as its acceptance are always strictly positive. We have thus obtained that P s τ {ω0}×I0 < +∞ P {ω0}×I0 (τ A < +∞) > 0, which concludes the proof.

Remark 20. The process (N t ) t 1 that gives the length of the ray S(t) at time t behaves like a birthdeath process with time inhomogeneous rates. If there exists an invariant measure of the process (N t ) t 1 , then it must coincide with the negative binomial distribution NB(1, ρ) that drives the length of a path S ∼ ν. This provides an easy criterium in order to check that the chain has already mixed, for example with a chi-squared test on the empirical distribution of (N t ) t 1 .

Simulation and comparison of the methods

In this section, we compare the estimates of the fluence rate L(x k ) provided by three methods: Monte Carlo with Wang-Prahl algorithm (denoted by WANG, see [START_REF] Prahl | Light transport in tissue[END_REF][START_REF] Wang | Monte Carlo modeling of light transport in multi-layered tissues in standard C[END_REF]), MC-SOME (see (3.2)) and the Metropolis-Hastings (MH) (see (4.2)) with the mutation rule given in Definition 17. We tested the methods in different settings. Here, we present results in a framework corresponding to a healthy homogeneous rat brain tissue. We chose to follow [START_REF] Angell-Petersen | Determination of fluence rate and temperature distributions in the rat brain; implications for photodynamic therapy[END_REF] for the values of the optical parameters (see Table 1). Other values for rat or human brain can be found in [START_REF] Bevilacqua | In vivo local determination of tissue optical properties: applications to human brain[END_REF][START_REF] Cheong | A review of the optical properties of biological tissues[END_REF][START_REF] Johns | Determination of reduced scattering coefficient of biological tissue from a needle-like probe[END_REF]. The volume of the cube V equals 8 cm 3 , that is V = [-1, 1] 3 . It is discretized into voxels whose volume is (0.04) 3 cm 3 . The half-opening angle of the optical fiber was set to α = π 10 and the constant c in (2.16) was set c = 1. We chose the following simulation parameters for the three methods so that they need the same amount of computational time. Those are WANG: M = 6000 photons trajectory. MC-SOME: M = 30000 rays, M points = 40 points chosen in each ray and M rot = 30 rotations with respect to the initial direction.

MH j = 10, J = 21, ε = 0.9, T = 250000 steps of the chain and M rot = 30 rotations with respect to the initial direction.

In Fig. 4, we picture for each methods, a zoom of the contour plot of the estimates in the plane x = -0.04. In Fig. 6, we compare the estimates along several lines of voxels (ℓ i ) i=1,...,6 which are parallel to the y-axis and pass through the points (0, 0, -0.08), (0, 0, -0.12), (0, 0, -0.16), (0, 0, -0.4), (0, 0, -0.48) and (0, 0, -0.6) respectively (see Fig. 5). We notice that MC-SOME gives more consistant estimates than the two other algorithms whose estimates are more noisy. Moreover, it seems that the MH-estimates are not very symmetric. Perhaps because the algorithm had not converged yet. Let us conclude this section by studying the accuracy of the methods by mean of 50 independent replicates of these estimates. In Fig. 7, the boxplots compare the dispersion of the 50 estimates of each method in six voxels (v i ) i=1,...,6 such that (see Fig. 5) (0, 0.2, 0) ∈ v 1 , (0, 0.6, 0) ∈ v 2 , (0, 0, -0.2) ∈ v 3 , (0, 0, -0.6) ∈ v 4 , (0, 0.2, -0.2) ∈ v 5 , (0, 0.6, -0.6) ∈ v 6 .

(

On one hand, we see that MC-SOME is much more consistent than WANG, because of the variance reduction that we use in MC-SOME. On the other hand, MH gives spread estimates. This is due to the simulations for which the Markov chain has not yet converged. In Table 2, we have the mean of the estimates and their mean square error in each of the 6 voxels.

Mean

Mean Square Error

Wang MC-SOME M-H Wang MC-SOME M-H v 1 1.203 • 10 -5

Inverse problem and sensitivity

For biologists, it is of considerable practical importance to have good estimates of the optical coefficients of the tissue they consider. One way to do this estimate is to compare simulated data with measurements of the fluence rate in the tissue and adjust the optical parameters of the simulation until obtaining values close to the measurements. Thanks to the probabilistic representation (2.16), this problem can be numerically solved as we shall see.

Sensitivity of the measurements

As a preliminary step towards a good resolution of the inverse problem, we first observe how fluence rate measurements vary with respect to the optical parameters g, µ s and µ a . To this aim, we built a small database of simulations for different values of the parameters and then compared the estimated fluence rate. The estimates are computed by resorting to MC-SOME, which is the best performing method among the three we have implemented according to Section 5. First, we choose a reference simulation obtained with reference parameters (g * , µ * a , µ * s ). Then, we pick n 1 voxels and consider their respective fluence rate estimates as measurements. That is, we choose n voxel centers (x ki ) i=1,...,n and define

m i = L(x ki ; g * , µ * a , µ * s ), i = 1, . . . , n,
where we recall that L(x ki ) is defined by (2.15) with γ W0 = 1 {ω ′ 0 ∈C 2α } σ(dω ′ 0 )/σ(C 2α ) and where we stress the dependence on the optical coefficients by writing L(x ki ; g * , µ * a , µ * s ) ≡ L(x ki ). Now for each Fluence rate

Estimates in voxel v possible triplet of parameters (g, µ a , µ s ), we compute the normalized quadratic error (or evaluation error )

J(µ a , µ s , g) = 1 2 n i=1 L(x ki ; g, µ a , µ s ) -m i m i 2 . (6.1)
For the dataset of simulations, we use the same settings as in Section 5 (|V | = 8 cm 3 , the volume of a voxel is (0.04) 3 cm 3 and α = π 10 ). The variable parameters are: g, µ a and µ s . Their values are given in Table 3. This choice is motivated by [START_REF] Angell-Petersen | Determination of fluence rate and temperature distributions in the rat brain; implications for photodynamic therapy[END_REF][START_REF] Bevilacqua | In vivo local determination of tissue optical properties: applications to human brain[END_REF][START_REF] Cheong | A review of the optical properties of biological tissues[END_REF][START_REF] Johns | Determination of reduced scattering coefficient of biological tissue from a needle-like probe[END_REF]. The anisotropy parameter g does not vary a lot between tissue type (healthy or tumorous) and it is often even hidden in a reduction of the scattering coefficient µ ′ s = µ s (1 -g). For this reason, we chose only three values in a small range of common values. Concerning the other parameters, we chose five values in intervals covering values corresponding to healthy and tumorous brain tissues according to [START_REF] Angell-Petersen | Determination of fluence rate and temperature distributions in the rat brain; implications for photodynamic therapy[END_REF][START_REF] Johns | Determination of reduced scattering coefficient of biological tissue from a needle-like probe[END_REF]. Table 3: Values of the optical parameters for the study of sensitivity.

Figures 8 and9 give different representation of the variation of the error J(µ a , µ s , g) with respect to the optical parameters. The real values are (µ * a , µ * s , g * ) = (0.75, 105, 0.9) and we set n = 3, x k1 ∈ v 2 , x k2 ∈ v 4 , x k3 ∈ v 6 respectively (see Fig. 5). We see that the sensitivity in the parameters µ s and g is very low compared to the sensitivity in µ a . In Fig. 8, we see that a wrong value of µ a has strong effects on the error function and that it becomes then almost impossible to see any tendency for the anisotropy parameter g. Notice also that an undervaluation of µ a is worse than an overvaluation in terms of the error. Figure 8: Colormap of the quadratic error (µ a , µ s ) → J(µ a , µ s , g) for three values of g, where µ s is displayed on the x axis and µ a on the y axis.

Parameters estimation

This section is devoted to the estimation of the parameters µ a and µ s only. Indeed, we have seen in the last section that the sensitivity of the fluence rate with respect to the anisotropic parameter g is low. Moreover, simulations do not show any monotonicity or tendency in the error for this parameter because, in our settings, the Monte Carlo error prevails over the evaluation error. In addition, for our purpose, the uncertainty about g is small in front of the uncertainty of the two other parameters (see [START_REF] Angell-Petersen | Determination of fluence rate and temperature distributions in the rat brain; implications for photodynamic therapy[END_REF]). We shall thus suppose in the sequel that g is known. With these preliminary considerations in mind, our goal is to solve the following nonlinear least square minimization problem: Find (µ s , µ a ) in order to minimize

J(µ s , µ a ) = 1 2 n i=1 L(x ki ; µ s , µ a ) -m i m i 2 , (6.2)
where (m i ) i=1,...,n are measurements in n different voxels centered at (x ki ) i=1,...,n .

The optimization method that we implemented to solve this problem is based on the Levenberg-Marquardt algorithm (see [START_REF] Fletcher | Practical methods of optimization[END_REF]). This gradient descent algorithm involves the computation of the gradient, as well as the Hessian matrix of the score function J. It is described in pseudo-code in Algorithm 2. In this description, we have set H k = Hess(J)(µ k s , µ k a ) for k 0. The term Diag(H k ) is the diagonal matrix of H k , λ k is the so-called damping factor which may be either constant or corrected at each step, and τ k ∈ R + controls the step size of each iteration.

Algorithm 2 Gradient descent algorithm for the estimation of µ a and µ s Input: measurements (m i ) i=1,...,n , initial couple (µ 0 s , µ 0 a ), precision ε > 0.

1: k ← 0 2: while J(µ k s , µ k a ) > ε do 3: (µ k+1 s , µ k+1 a ) ← (µ k s , µ k a ) -τ k [H k + λ k Diag(H k )] -1 ∇J(µ k s , µ k a ) 4: k ← k + 1 5: end while Output: an approximation (µ k s , µ k a ) of the real parameters (µ * s , µ * a )
The simple form of the objective function in (6.2) allows to express the term on the right hand side of line 3 in Algorithm 2 explicitly as a function of the partial derivatives of L. Indeed, the gradient of J is given by

∇J(•) = n i=1 L(x ki ; •) -m i m 2 i ∇L(x ki ; •), (6.3)
and its Hessian matrix is given by Hess

(J)(•) = n i=1 L(x ki ; •) -m i m 2 i Hess(L)(x i ; •) + 1 m 2 i ∇L(x ki ; •)∇ t L(x ki ; •) . (6.4)
Moreover, as stated in the following proposition, the formal representation in Proposition 6 allows to also use the Monte Carlo method MC-SOME in order to estimate the first order and the second order partial derivatives of L which can be expressed similarly to (2.18). Proposition 21. The partial derivatives of L(x ki ; µ s , µ a ) can be expressed as the expectation of fully simulable random variables. Using the same notations as in (2.17), they are given by

∂L ∂µ a (x ki ; µ s , µ a ) = - c(1 -cos(α)) 2µ a E   1 {SN ∈V k i } N j=0 R j   (6.5) and ∂L ∂µ s (x ki ; µ s , µ a ) = c(1 -cos(α)) 2µ a E   1 {SN ∈V k i }   N µ s - N j=0 R j     . (6.6)
Proof. We start by differentiating term-by-term the Neumann series of Corollary 3. Let us first note that by definition of L i we have:

∂L i ∂µ a (x, ω 0 ; µ s , µ a ) = - +∞ 0 r exp(-(µ s + µ a )r)L e (x -rw, w)dr.
Moreover, for n 1, by definition of T n • T i (see (2.5)), we also have:

∂ [T n • T i ] L e ∂µ a (x, ω 0 ; µ s , µ a ) = µ n s R n + dr 0 • • • dr n   - n j=0 r j   exp   -(µ s + µ a ) n j=0 r j   (S 2 ) (n+1) dσ ⊗(n+1) (ω 0 , . . . , ω n ) n-1 j=0 f HG (ω j , ω j+1 ) L e x - n k=0 r k ω k , ω n .
Looking back at Section 2.3 and using the same notations, we deduce that

∞ n=0 ∂ [T n • T i ] L e ∂µ a (x, ω 0 ; µ s , µ a ) = - A ν(dr, dω) G x (r, ω) |r| i=0 r i ,
where we recall that |r| stands for the size of r. Assuming that the left-hand side coincides with the partial derivative ∂L ∂µa , then (6.5) is found just like (2.17) and the same arguments provide (6.6), considering that for all n ≥ 0

∂ [T n • T i ] L e ∂µ s = n µ s [T n • T i ] L e + ∂ [T n • T i ] L e ∂µ a .
To conclude the proof, notice that the match between the partial derivatives and the term-by-term differentiation of the Neumann series is ensured by the fact that the operator T n • T i is infinitely continuously differentiable for all n and by the uniform convergence of the corresponding sequences of truncated sums

s m = m n=0 ∂[T n • T i ]L e
∂µ s (x, ω 0 ; µ s , µ a ), m 0.

Remark 22. Similar formula to (6.5) and (6.6) can be easily found for the second order derivatives Implementation and discussion. The randomness coming from Monte Carlo estimation of the score J, of its gradient and of its Hessian matrix during the run of the algorithm, makes a precise estimation of the real values of (µ * a , µ * s ) difficult. We observed that far from the real value µ * a , the eigenvalues of the Hessian matrix Ĥk are very small and that their sign can vary a lot because of the volatility of the estimates. Conversely, near the real value µ * a , the estimate Ĥk is more robust and its eigenvalues are almost always both positive, which legitimates the quadratic approximation of Levenberg-Marquardt algorithm. For these reasons, we implemented a hybrid algorithm which chooses between the Levenberg-Marquardt descent and the classic steepest gradient descent depending on the sign of the eigenvalues of H k . If they are both positive, one moves to the next point following (6.7), else one makes a move in the opposite direction of the gradient, -∇J.

In Fig. 10 and 11, we can see two examples of descent of our algorithm. The settings are the following: (1) we choose n = 3 positions for the measurements: x k1 ∈ v 2 , x k2 ∈ v 4 , x k3 ∈ v 6 (see Fig. 5) and the values of the measurements m 1 , m 2 and m 3 are taken from the database of simulations described in Sec. 6.1 with the desired parameters, (2) the anisotropy factor is set to g = 0.9, (3) the damping parameter is constant λ = 0.01, (4) the precision parameter (see Algo. 2) is set to ε = 0.005, (5) the sequence (τ k ) k 1 controlling the step size of each iteration is decreasing in k and depends on the score of the iteration, as well as on the sign of the eigenvalues of H k .

In Fig. 10, the reference parameter are (µ * a = 1, µ * s = 75). Notice the oscillations around µ * a . Those descent zigzags near the real value of µ a are also apparent in the other descent in Fig. 11 for which (µ * a = 1, µ * s = 105). They correspond to iteration where the descent is done according to the classic steepest descent. The iterations for which one moves more vertically correspond, as for them, to the case where the descent is done according to (6.7). As we can see, a satisfying estimate of µ * a comes up rapidly, whereas µ * s is more difficult to approach. This is due to the low sensitivity of J with respect to µ s discussed in the previous section. ). The starting point is [START_REF] Arridge | A finite element approach for modeling photon transport in tissue[END_REF]90). The final estimate is (1.17, 107.04) with a score equal to J = 0.0048.

Figure 1 :

 1 Figure 1: Density function of cos(θ) = ω, ω when ω is scattered according to the Henyey-Greenstein function f HG (ω, •) for several values of the anisotropy parameter g.

i

  according to f εg HG . -Draw a new angle ϕ new i uniformly on [0, 2π].

Figure 3 :

 3 Figure 3: Metropolis-Hastings algorithm. Example of a mutation by rotation.

Figure 4 :

 4 Figure 4: Contour plots of the fluence rate estimates in the plane x = -0.04 for WANG, MC-SOME and MH.

Figure 5 : 6 Figure 6 :

 566 Figure 5: Choice of six particular voxels and position of the lines (ℓ i ) i=1,...,6 in the cube V .

Figure 7 :

 7 Figure 7: Boxplots of 50 fluence rate estimates in the voxels (v i ) i=1,...,6 with WANG, MC-SOME and MH.

µaµaµa 135 Figure 9 :

 1359 Figure 9: Quadratic error µ a → J(µ a , µ s , g) for three values of g.

  of L(x ki ; µ s , µ a ) in(2.18) and its partial derivatives allows us to estimate the score J(µ s , µ a ), its gradient and its Hessian matrix by Monte Carlo methods. A sole sample (y 1 , . . . , y n ) ∈ A n of n observations of the random ray Y can be used to estimate the expectations in L, ∇L(x ki ; •) and Hess(L)(x i ; •) at the same time. We denote these estimates by L, ∇L(x ki ; •) and Hess(L)(x i ; •) and the corresponding score by J. The updating rule at line 3 in Algorithm 2 becomes then(µ k+1 s , µ k+1 a ) = (µ k s , µ k a ) -τ k H k + λ k Diag( H k ) -1∇J(µ k s , µ k a ). (6.7)

Figure 10 :Figure 11 :

 1011 Figure 10: Parameters estimation with an adaptation of Levenberg-Marquardt descent algorithm. The real value of (µ * a , µ * s ) is (1, 75). The starting point is (2, 90). The final estimate is (1.07, 76.02) with a score equal to J = 0.0035.

  3

  r nt , θ nt , ϕ nt , r new

	nt+1 , θ new nt+1 , ϕ new nt+1 , . . . , r new nt+∆(nt) , θ new nt+∆(nt) , ϕ new nt+∆(nt) .

  1.2366 • 10 -5 1.3083 • 10 -5 1.4484 • 10 -12 3.9108 • 10 -13 4.5158 • 10 -12 v 2 2.4309 • 10 -7 2.7177 • 10 -7 1.6545 • 10 -7 6.3788 • 10 -15 2.5597 • 10 -15 1.4411 • 10 -14 v 3 2.0014 • 10 -5 2.0033 • 10 -5 1.8337 • 10 -5 2.503 • 10 -12 3.0408 • 10 -13 1.6200 • 10 -11 v 4 3.4947 • 10 -7 3.5713 • 10 -7 2.1497 • 10 -7 1.4136 • 10 -14 8.1737 • 10 -16 1.5819 • 10 -14 v 5 6.7519 • 10 -6 6.6047 • 10 -6 5.4786 • 10 -6 8.0457 • 10 -13 4.4781 • 10 -14 1.1378 • 10 -12 v 6 3.9952 • 10 -8 4.217 • 10 -8 1.8755 • 10 -8 8.4058 • 10 -16 6.6977 • 10 -17 6.5413 • 10 -16

Table 2 :

 2 Mean and mean square error of the 50 fluence rate estimates in 6 voxels with WANG, MC-SOME and MH.